doi:10.1016/0010-4655(95)00005-z

ranshirand.f90

program ranshirand
    common /ranbuf/ mbuff(130)

    real*4 a(1024)                                                                      ! size of array to fill with random numbers

    call ranvin(1337,100)                                                               ! seed and number of generations to warm up mbuff for
    call ranshi(size(a),a)
    print *, a

    contains
        subroutine ranshi(n, a)                                                         !1 put n random
                                                                                        !2 numbers into array a
            parameter(nb = 128, nbit = 32, is = 17, amin = 1.e-10)                      !3
            parameter(nc = nb + 2, mask =nb -1, scalin= 2.**(1-nbit))                   !4
            common /ranbuf/ mbuff(nc)                                                   !5 buffer
            real*4 a(*)                                                                 !6 output randoms

            iangle = mbuff(nb + 1)                                                      !7 load red angle
            ic = mbuff(nb + 2)                                                          !8 load red spin

            do i = 1, n                                                                 !9 vector-like loop
                ir = mbuff( iangle )                                                    !10 hit ball (iangle)
                a(i)= float(ishft(ir,-1)) * scalin + amin                               !11 avoid zero output
                iro = ishftc( ir, is, nbit)                                             !12 circ. left shift by
                                                                                        !   is bits, see below
                mbuff( iangle ) = ieor( iro, ic)                                        !13 replace black spin
                ic = ir                                                                 !14 replace red spin
                iangle = iand(mask , ir + i ) + 1                                       !15 boost spin and
                                                                                        !16 mask for new angle
            enddo                                                                       !17 end of loop
            mbuff(nb + 1) = iangle                                                      !18 store red angle
            mbuff(nb + 2) = ic                                                          !19 store red spin

            return                                                                      !20
            end                                                                         !21

        subroutine ranvin(iseed, iwarm)

            !! initialize buffer !!

            real*8 rtval, modu, modi, aa, ac, xint, scale
            parameter(modu=1771875.0d0, aa = 2416.d0, ac = 37444.d0)
            parameter(nb = 128, nc =nb + 2, nbit = 32)
            parameter( scale = 2.d0 ** nbit )
            common /ranbuf/ mbuff(nc)
            
            !! initialize buffer !!
            real*4 a(nb)
            modi = 1.d0 / modu
            rtval = float(4*iseed+1)
            do j = 1, nc
                xint = rtval * aa + ac
                ntval = xint * modi
                rtval = xint - ntval * modu
                rnanf = rtval * modi
                mbuff(j) = scale * (rnanf - 0.5)
            enddo
            mbuff(nb+1)= iand(mbuff(nb+1), nb- 1) + 1                                   ! restrict red angle
            do j = 1, iwarm                                                             ! warm that shit up
                call ranshi (nb, a)
            enddo
            return 
            end

end program