16 Jan 2024 12:33
(Edited: 3 Feb 2024 7:53)
#30957
The late George Marsaglia posted his SuperKiss-64 RNG code in the 'General' forum in 2009 ( https://forums.silverfrost.com/Forum/Topic/1217 ). At the time (and, perhaps, even now), 32-bit code generated by FTN95 used the FPU for 8-byte integer operations, which leads to lots of FPU exceptions, so we use only /64 for this code.
Marsaglia's test program generates a billion 8-byte random integers and compares the last of those to a known value. With FTN95 8.972 and 9.0, the test fails. With 8.9, 8.92, 8.95, the test passes. Here is a very short program that contains a small portion of Marsaglia's code, just enough to display the bug.
! Suitable only for 64-bit compilation
! Some of the last digits are wrong with Versions 8.972, 9.0
! Correct output with Versions 8.9, 8.92, 8.95, as follows
! i h z carry q(i)
! 1 1 -2945433577893571930 863862222379 5890867155787143858
! 2 1 6661552707265476373 1293535672257 5123638659178598868
! 3 1 -4642211112981496864 93621267173 -9162321847746557890
! 4 1 -8101377823381285518 1808705431682 -2243988426946980582
!
program S64BITOPS
implicit None
Integer, Parameter :: QSIZ = 4
Integer, Parameter :: i8 = selected_int_kind(18)
Integer (Kind=i8) :: i, z, h, q(QSIZ)
Integer (Kind=i8), Parameter :: CMUL = 6906969069_i8, COFFS = 123_i8
Integer (Kind=i8) :: carry = 36243678541_i8, &
xcng = 12367890123456_i8, xs = 521288629546311_i8
Do i = 1, qsiz ! fill Q with Congruential + Xorshift
xcng = xcng*CMUL + COFFS ! ignore integer overflow here
xs = ieor(xs, ishft(xs,13))
xs = ieor(xs, ishft(xs,-17))
xs = ieor(xs, ishft(xs,43))
q(i) = xcng + xs
End Do
print *,' i h z carry q(i)'
Do i = 1, QSIZ
h = iand(carry, 1_i8)
z = ishft(ishft(q(i),41), -1) + &
ishft(ishft(q(i),39), -1) + &
ishft(carry, -1)
carry = ishft(q(i), -23) + ishft(q(i), -25) + ishft(z, -63)
q(i) = not(ishft(z,1)+h)
print '(2i3,i21,i15,i21)',i,h,z,carry,q(i)
End Do
end program