! { dg-do compile }
! { dg-options "-O2 -fdump-tree-optimized" }
!
! PR fortran/54556
!
! Contributed by Joost VandeVondele
!
MODULE parallel_rng_types

  IMPLICIT NONE

  ! Global parameters in this module
  INTEGER, PARAMETER :: dp=8

  TYPE rng_stream_type
    PRIVATE
    CHARACTER(LEN=40)             :: name
    INTEGER                       :: distribution_type
    REAL(KIND=dp), DIMENSION(3,2) :: bg,cg,ig
    LOGICAL                       :: antithetic,extended_precision
    REAL(KIND=dp)                 :: buffer
    LOGICAL                       :: buffer_filled
  END TYPE rng_stream_type

  REAL(KIND=dp), DIMENSION(3,3) :: a1p0,a1p76,a1p127,&
                                   a2p0,a2p76,a2p127,&
                                   inv_a1,inv_a2

  INTEGER, PARAMETER          :: GAUSSIAN = 1,&
                                 UNIFORM  = 2

  REAL(KIND=dp), PARAMETER :: norm  = 2.328306549295727688e-10_dp,&
                              m1    = 4294967087.0_dp,&
                              m2    = 4294944443.0_dp,&
                              a12   = 1403580.0_dp,&
                              a13n  = 810728.0_dp,&
                              a21   = 527612.0_dp,&
                              a23n  = 1370589.0_dp,&
                              two17 = 131072.0_dp,&            ! 2**17
                              two53 = 9007199254740992.0_dp,&  ! 2**53
                              fact  = 5.9604644775390625e-8_dp ! 1/2**24


CONTAINS

  FUNCTION rn32(rng_stream) RESULT(u)

    TYPE(rng_stream_type), POINTER           :: rng_stream
    REAL(KIND=dp)                            :: u

    INTEGER                                  :: k
    REAL(KIND=dp)                            :: p1, p2

! -------------------------------------------------------------------------
! Component 1

    p1 = a12*rng_stream%cg(2,1) - a13n*rng_stream%cg(1,1)
    k = INT(p1/m1)
    p1 = p1 - k*m1
    IF (p1 < 0.0_dp) p1 = p1 + m1
    rng_stream%cg(1,1) = rng_stream%cg(2,1)
    rng_stream%cg(2,1) = rng_stream%cg(3,1)
    rng_stream%cg(3,1) = p1

    ! Component 2

    p2 = a21*rng_stream%cg(3,2) - a23n*rng_stream%cg(1,2)
    k = INT(p2/m2)
    p2 = p2 - k*m2
    IF (p2 < 0.0_dp) p2 = p2 + m2
    rng_stream%cg(1,2) = rng_stream%cg(2,2)
    rng_stream%cg(2,2) = rng_stream%cg(3,2)
    rng_stream%cg(3,2) = p2

    ! Combination

    IF (p1 > p2) THEN
      u = (p1 - p2)*norm
    ELSE
      u = (p1 - p2 + m1)*norm
    END IF

    IF (rng_stream%antithetic) u = 1.0_dp - u

  END FUNCTION rn32

! *****************************************************************************
  FUNCTION rn53(rng_stream) RESULT(u)

    TYPE(rng_stream_type), POINTER           :: rng_stream
    REAL(KIND=dp)                            :: u

    u = rn32(rng_stream)

    IF (rng_stream%antithetic) THEN
      u = u + (rn32(rng_stream) - 1.0_dp)*fact
      IF (u < 0.0_dp) u = u + 1.0_dp
    ELSE
      u = u + rn32(rng_stream)*fact
      IF (u >= 1.0_dp) u = u - 1.0_dp
    END IF

  END FUNCTION rn53

END MODULE

! { dg-final { scan-module-absence "parallel_rng_types" "IMPLICIT_PURE" } }
! { dg-final { scan-tree-dump-times "rn32 \\(rng_stream" 3 "optimized" } }
! { dg-final { cleanup-tree-dump "optimized" } }