pf_random_normal_portable_at Interface

public interface pf_random_normal_portable_at

A standard normal that is identical on every platform, compiler and flag set, not merely for a given libm.

Same arguments and same distribution as pf_random_normal_at; a different algorithm, and so a different value at the same coordinate. Marsaglia's polar method: draw a point in the square, reject it unless it lands in the unit disc, and map the survivor through sqrt(-2*log(s)/s). Every step is IEEE + - * / and sqrt -- both correctly rounded by the standard -- over parquet_expkey's frozen logarithm, with rounding barriers on the two squares and on the argument of the sqrt so that no compiler may fuse or regroup them.

It costs more than the Ziggurat, in two ways. It rejects 21.5 % of candidate pairs rather than 1.5 % of single draws, it discards the second variate the method produces (a stream that kept it would have hidden state, so %rewind would stop being exact), and its logarithm is about 3x libm's. Reach for it when a stored result must reproduce across machines; take pf_random_normal_at otherwise.

Not pure, because the frozen transform's rounding barrier is a volatile local and a pure procedure may not have one. Elemental use still works, but a caller's own pure procedure and a do concurrent body cannot reach it -- pf_random_normal_at can.


Module Procedures

private impure elemental function pf_random_normal_portable_at_i32(seed, i, draw) result(r)

pf_random_normal_portable_at for an integer(int32) stream index.

impure because exp_key is: the frozen transform's rounding barrier is a volatile local, and a pure procedure may not have one. Elemental use is unaffected.

Arguments

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

the stream family's seed
integer(kind=int32), intent(in) :: i

stream index; sign-extends, so any value is valid
integer(kind=int64), intent(in), optional :: draw

1-based value index; absent means 1

Return Value real(kind=real64)

a standard normal draw

private impure elemental function pf_random_normal_portable_at_i64(seed, i, draw) result(r)

pf_random_normal_portable_at for an integer(int64) stream index.

Arguments

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

the stream family's seed
integer(kind=int64), intent(in) :: i

stream index; every value is valid
integer(kind=int64), intent(in), optional :: draw

1-based value index; absent means 1

Return Value real(kind=real64)

a standard normal draw