u64_mul Function

public pure function u64_mul(a, b) result(r)

(a * b) reduced modulo 2**64, without ever overflowing a signed 64-bit integer.

Schoolbook multiplication on 16-bit limbs: each partial product is below 232 and at most four of them plus a carry are accumulated, so every intermediate stays below 235 and signed overflow is impossible regardless of the sign of the inputs. Limbs are extracted with ibits (bit-model, sign-safe) and the result is reassembled with ior/ishft.

Arguments

Type IntentOptional Attributes Name
integer(kind=int64), intent(in) :: a
integer(kind=int64), intent(in) :: b

Return Value integer(kind=int64)


Called by

proc~~u64_mul~~CalledByGraph proc~u64_mul u64_mul proc~mix64 mix64 proc~mix64->proc~u64_mul proc~pcg32_step pcg32_step proc~pcg32_step->proc~u64_mul proc~stream_seed stream_seed proc~stream_seed->proc~u64_mul proc~stream_seed->proc~mix64 proc~pcg32_next pcg32_t%pcg32_next proc~pcg32_next->proc~pcg32_step proc~pcg32_next_u64 pcg32_t%pcg32_next_u64 proc~pcg32_next_u64->proc~pcg32_step proc~pcg32_seed pcg32_t%pcg32_seed proc~pcg32_seed->proc~pcg32_step proc~splitmix64_next splitmix64_t%splitmix64_next proc~splitmix64_next->proc~mix64 proc~stream_for_pcg32 stream_for_pcg32 proc~stream_for_pcg32->proc~stream_seed proc~stream_for_pcg32->proc~splitmix64_next proc~stream_for_splitmix64 stream_for_splitmix64 proc~stream_for_splitmix64->proc~stream_seed interface~stream_for stream_for interface~stream_for->proc~stream_for_pcg32 interface~stream_for->proc~stream_for_splitmix64 proc~pcg32_next_below pcg32_next_below proc~pcg32_next_below->proc~pcg32_next_u64 proc~pcg32_next_exponential_int pcg32_next_exponential_int proc~pcg32_next_exponential_int->proc~pcg32_next_u64 proc~pcg32_next_real_dp pcg32_next_real_dp proc~pcg32_next_real_dp->proc~pcg32_next_u64 proc~pcg32_uniform_below pcg32_uniform_below proc~pcg32_uniform_below->proc~pcg32_next_u64 proc~splitmix64_next_u64 splitmix64_t%splitmix64_next_u64 proc~splitmix64_next_u64->proc~splitmix64_next interface~next_below next_below interface~next_below->proc~pcg32_next_below proc~splitmix64_next_below splitmix64_next_below interface~next_below->proc~splitmix64_next_below interface~next_exponential_int next_exponential_int interface~next_exponential_int->proc~pcg32_next_exponential_int proc~splitmix64_next_exponential_int splitmix64_next_exponential_int interface~next_exponential_int->proc~splitmix64_next_exponential_int interface~next_real_dp next_real_dp interface~next_real_dp->proc~pcg32_next_real_dp proc~splitmix64_next_real_dp splitmix64_next_real_dp interface~next_real_dp->proc~splitmix64_next_real_dp proc~pcg32_next_bernoulli_ppm pcg32_next_bernoulli_ppm proc~pcg32_next_bernoulli_ppm->proc~pcg32_uniform_below proc~pcg32_next_poisson_int pcg32_next_poisson_int proc~pcg32_next_poisson_int->proc~pcg32_next_exponential_int proc~pcg32_next_range pcg32_next_range proc~pcg32_next_range->proc~pcg32_uniform_below proc~splitmix64_next_below->proc~splitmix64_next_u64 proc~splitmix64_next_exponential_int->proc~splitmix64_next_u64 proc~splitmix64_next_real_dp->proc~splitmix64_next_u64 proc~splitmix64_uniform_below splitmix64_uniform_below proc~splitmix64_uniform_below->proc~splitmix64_next_u64 interface~next_bernoulli_ppm next_bernoulli_ppm interface~next_bernoulli_ppm->proc~pcg32_next_bernoulli_ppm interface~next_poisson_int next_poisson_int interface~next_poisson_int->proc~pcg32_next_poisson_int interface~next_range next_range interface~next_range->proc~pcg32_next_range proc~pcg32_next_exponential_dp pcg32_next_exponential_dp proc~pcg32_next_exponential_dp->interface~next_real_dp proc~pcg32_next_normal_dp pcg32_next_normal_dp proc~pcg32_next_normal_dp->interface~next_real_dp proc~splitmix64_next_bernoulli_ppm splitmix64_next_bernoulli_ppm proc~splitmix64_next_bernoulli_ppm->proc~splitmix64_uniform_below proc~splitmix64_next_exponential_dp splitmix64_next_exponential_dp proc~splitmix64_next_exponential_dp->interface~next_real_dp proc~splitmix64_next_normal_dp splitmix64_next_normal_dp proc~splitmix64_next_normal_dp->interface~next_real_dp proc~splitmix64_next_poisson_int splitmix64_next_poisson_int proc~splitmix64_next_poisson_int->proc~splitmix64_next_exponential_int proc~splitmix64_next_range splitmix64_next_range proc~splitmix64_next_range->proc~splitmix64_uniform_below

Variables

Type Visibility Attributes Name Initial
integer(kind=int64), private :: a0
integer(kind=int64), private :: a1
integer(kind=int64), private :: a2
integer(kind=int64), private :: a3
integer(kind=int64), private :: b0
integer(kind=int64), private :: b1
integer(kind=int64), private :: b2
integer(kind=int64), private :: b3
integer(kind=int64), private :: carry
integer(kind=int64), private :: r0
integer(kind=int64), private :: r1
integer(kind=int64), private :: r2
integer(kind=int64), private :: r3
integer(kind=int64), private :: t

Source Code

   pure function u64_mul(a, b) result(r)
      !! `(a * b)` reduced modulo 2**64, without ever overflowing a signed
      !! 64-bit integer.
      !!
      !! Schoolbook multiplication on 16-bit limbs: each partial product is
      !! below 2**32 and at most four of them plus a carry are accumulated, so
      !! every intermediate stays below 2**35 and signed overflow is impossible
      !! regardless of the sign of the inputs. Limbs are extracted with `ibits`
      !! (bit-model, sign-safe) and the result is reassembled with `ior`/`ishft`.
      integer(int64), intent(in) :: a
      integer(int64), intent(in) :: b
      integer(int64) :: r

      integer(int64) :: a0, a1, a2, a3
      integer(int64) :: b0, b1, b2, b3
      integer(int64) :: t, carry
      integer(int64) :: r0, r1, r2, r3

      a0 = ibits(a, 0, 16)
      a1 = ibits(a, 16, 16)
      a2 = ibits(a, 32, 16)
      a3 = ibits(a, 48, 16)
      b0 = ibits(b, 0, 16)
      b1 = ibits(b, 16, 16)
      b2 = ibits(b, 32, 16)
      b3 = ibits(b, 48, 16)

      t = a0*b0
      r0 = iand(t, MASK16)
      carry = ibits(t, 16, 48)

      t = a0*b1 + a1*b0 + carry
      r1 = iand(t, MASK16)
      carry = ibits(t, 16, 48)

      t = a0*b2 + a1*b1 + a2*b0 + carry
      r2 = iand(t, MASK16)
      carry = ibits(t, 16, 48)

      t = a0*b3 + a1*b2 + a2*b1 + a3*b0 + carry
      r3 = iand(t, MASK16)

      r = ior(ior(ishft(r3, 48), ishft(r2, 32)), ior(ishft(r1, 16), r0))
   end function u64_mul