scale_by_mean Function

private pure function scale_by_mean(x_q32, mean) result(r)

Multiply a Q32.32 value by an integer mean and drop the fraction.

Done in two halves so that nothing overflows. The integer part of x_q32 is at most 9ln(4096), under 75, so the first product cannot exceed 75 * huge(int32). The fractional part is below 232 and mean is at most huge(int32), so the second product is at most (232 - 1) * (2*31 - 1); adding the rounding term below still leaves it under huge(int64), with about 4e9 to spare.

The fraction is rounded to nearest rather than truncated. Truncating would lose on average half a unit per draw, which is invisible at a mean of 4000 ticks and a 23% error at a mean of 2.

Rounding shrinks that but does not remove it. Discretising an exponential to integers leaves E[round(mean*X)] at exp(-1/(2*mean)) / (1 - exp(-1/mean)), which is 4.0% low at a mean of 1, 1.0% at 2, 0.46% at 3 and 0.04% at 10. Above about 100 it is gone. That is inherent to an integer result, not a defect in the scaling, but a caller working in units that small should know the mean it gets is not quite the mean it asked for.

Arguments

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

Deviate in Q32.32 fixed point.
integer(kind=int64), intent(in) :: mean

Mean, in [0, EXPONENTIAL_MAX_MEAN].

Return Value integer(kind=int64)


Calls

proc~~scale_by_mean~~CallsGraph proc~scale_by_mean scale_by_mean proc~u64_shr u64_shr proc~scale_by_mean->proc~u64_shr

Called by

proc~~scale_by_mean~~CalledByGraph proc~scale_by_mean scale_by_mean proc~pcg32_next_exponential_int pcg32_next_exponential_int proc~pcg32_next_exponential_int->proc~scale_by_mean proc~splitmix64_next_exponential_int splitmix64_next_exponential_int proc~splitmix64_next_exponential_int->proc~scale_by_mean interface~next_exponential_int next_exponential_int interface~next_exponential_int->proc~pcg32_next_exponential_int interface~next_exponential_int->proc~splitmix64_next_exponential_int proc~pcg32_next_poisson_int pcg32_next_poisson_int proc~pcg32_next_poisson_int->proc~pcg32_next_exponential_int proc~splitmix64_next_poisson_int splitmix64_next_poisson_int proc~splitmix64_next_poisson_int->proc~splitmix64_next_exponential_int interface~next_poisson_int next_poisson_int interface~next_poisson_int->proc~pcg32_next_poisson_int interface~next_poisson_int->proc~splitmix64_next_poisson_int

Source Code

   pure function scale_by_mean(x_q32, mean) result(r)
      !! Multiply a Q32.32 value by an integer mean and drop the fraction.
      !!
      !! Done in two halves so that nothing overflows. The integer part of
      !! `x_q32` is at most 9*ln(4096), under 75, so the first product cannot
      !! exceed 75 * huge(int32). The fractional part is below 2**32 and
      !! `mean` is at most huge(int32), so the second product is at most
      !! (2**32 - 1) * (2**31 - 1); adding the rounding term below still
      !! leaves it under huge(int64), with about 4e9 to spare.
      !!
      !! The fraction is rounded to nearest rather than truncated. Truncating
      !! would lose on average half a unit per draw, which is invisible at a
      !! mean of 4000 ticks and a 23% error at a mean of 2.
      !!
      !! Rounding shrinks that but does not remove it. Discretising an
      !! exponential to integers leaves `E[round(mean*X)]` at
      !! `exp(-1/(2*mean)) / (1 - exp(-1/mean))`, which is 4.0% low at a mean
      !! of 1, 1.0% at 2, 0.46% at 3 and 0.04% at 10. Above about 100 it is
      !! gone. That is inherent to an integer result, not a defect in the
      !! scaling, but a caller working in units that small should know the
      !! mean it gets is not quite the mean it asked for.
      integer(int64), intent(in) :: x_q32
         !! Deviate in Q32.32 fixed point.
      integer(int64), intent(in) :: mean
         !! Mean, in `[0, EXPONENTIAL_MAX_MEAN]`.
      integer(int64) :: r

      r = u64_shr(x_q32, 32_default_int)*mean &
          + u64_shr(iand(x_q32, Q32_MASK)*mean + Q32_HALF, 32_default_int)
   end function scale_by_mean