Attempt of porting the PCG random number generator by Melissa O'Neill to Fortran. The code is not tested. Signed integers are used for the state.

pcg32.f90

!
! PCG Random Number Generation for Fortran
! Ported from the minimal C version by Melissa O'Neill
!
! Copyright 2020 Ivan Pribec <ivan.pribec@gmail.com>
!

!
! Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
! For additional information about the PCG random number generation scheme,
! including its license and other licensing options, visit
!
!       http://www.pcg-random.org
!

!
! This code is derived from the full C implementation, which is in turn
! derived from the canonical C++ PCG implementation. The C++ version
! has many additional features and is preferable if you can use C++ in
! your project.
!

module pcg32

  use iso_fortran_env, only: int32, int64
  implicit none

  integer, parameter :: dp = kind(1.0d0)

  type pcg32_random_t
    private
    integer(int64) :: state
    integer(int64) :: inc
  end type

  ! If you *must* statically initialize it, here's one.
  type(pcg32_random_t), parameter :: pcg32_initializer = pcg32_random_t( &
      -8846114313915602277_int64, &
      -2720673578348880933_int64)

  ! z'853c49e6748fea9b'
  ! b'1000010100111100010010011110011001110100100011111110101010011011'
  ! decimal: 9600629759793949339
  ! 64-bit signed int: -8846114313915602277

  ! z'da3e39cb94b95bdb'
  ! b'1101101000111110001110011100101110010100101110010101101111011011'
  ! decimal: 15726070495360670683
  ! 64-bit signed int: -2720673578348880933

  ! State of the random number generator
  type(pcg32_random_t) :: pcg32_global = pcg32_initializer

contains

  !
  ! Seed the rng. Specified in two parts, state initializer and a
  ! sequence selection constant (a.k.a stream id)
  !
  subroutine pcg32_srandom(initstate, initseq)
    integer(int64), intent(in) :: initstate, initseq
    call pcg32_srandom_r(pcg32_global, initstate, initseq)
  end subroutine

  subroutine pcg32_srandom_r(rng,initstate,initseq)
    type(pcg32_random_t), intent(out) :: rng
    integer(int64), intent(in) :: initstate, initseq
    integer(int32) :: discard
    rng%state = 0
    rng%inc = ior(shiftl(initseq, 1_int64),1_int64)
    discard = pcg32_random_r(rng)
    rng%state = rng%state + initstate
    discard = pcg32_random_r(rng)
  end subroutine

  !
  ! Generate a uniformly distributed 32-bit random number
  !
  integer(int32) function pcg32_random_r(rng)
    type(pcg32_random_t), intent(inout) :: rng
    integer(int64) :: oldstate
    integer(int32) :: xorshifted, rot
    oldstate = rng%state
    rng%state = oldstate*6364136223846793005_int64 + rng%inc
    xorshifted = int(shiftr(ieor(shiftr(oldstate,18_int64),oldstate),27_int64),kind=int32)
    rot = int(shiftr(oldstate,59_int64),kind=int32)
    pcg32_random_r = ior(shiftr(xorshifted,rot),shiftl(xorshifted,iand(-rot,31_int32)))
  end function
  integer(int32) function pcg32_random()
    pcg32_random = pcg32_random_r(pcg32_global)
  end function

  !
  ! Generate a 32-bit floating point number between 0 and 1
  !
  ! WARNING: This procedure is not verified!
  !
  real(dp) function pcg32_random_double_r(rng)
    type(pcg32_random_t), intent(inout) :: rng
    integer(int32) :: ix
    real(dp) :: x
    ix = pcg32_random_r(rng)
    
    ! We add huge(ix) to recover the original signed integer value.
    x = real(ix,dp) + huge(ix) + 1

    ! See http://computer-programming-forum.com/49-fortran/42bffeb406269a1a.htm
    ! for a discussion of the scale function, supposed to be equivalent to ldexp in C.
    pcg32_random_double_r = scale(x,-32)
  end function
  real(dp) function pcg32_random_double()
    pcg32_random_double = pcg32_random_double_r(pcg32_global)
  end function

  subroutine random_double(x)
    real(dp), intent(out) :: x(:)
    integer :: i
    do i = 1, size(x)
      x(i) = pcg32_random_double()
    end do
  end subroutine

end module

program test_pcg

  use pcg32
  implicit none
  integer :: i

  real(dp), allocatable :: A(:)
  real(dp) :: t0, t1

  allocate(A(100000000))

  do i = 1, 1000
    print *, pcg32_random_double()
  end do

  call cpu_time(t0)
  call random_double(A)
  call cpu_time(t1)
  print *, "time = ", t1 - t0
  print *, all(A >= 0.0_dp), all(A < 1.0_dp), minval(A), maxval(A)

  call cpu_time(t0)
  call random_number(A)
  call cpu_time(t1)
  print *, "time = ", t1 - t0
  print *, all(A >= 0.0_dp), all(A < 1.0_dp), minval(A), maxval(A)

end program