Compiler Explorer

Source code

Module prec
      Use iso_fortran_env
      Implicit None
      Integer wp
      Parameter (wp=real64)
    End Module

Program piqp_coarray

! This program employs the recently discovered digit extraction scheme
      ! to produce hex digits of pi.  This code is valid up to ic = 2^24 on
      ! systems with IEEE arithmetic.
      
      Use prec
      Implicit None
      Real (Kind=wp) :: pid, s1, s2, s3, s4, s[ * ]
      Integer :: i, length, status, ic, iextra, nhx
      Parameter (nhx=12)
      Character chx(nhx)
      Character (Len=256) arg
      ! ic is the hex digit position -- output begins at position ic + 1.

If (command_argument_count()/=2) Then
        If (this_image()==1) Then
          Write (error_unit, '(A)') '<prog> hex-digit-position safety-iters'
          Error Stop 1
        End If
      End If
      Call get_command_argument(1, arg, length, status)
      Read (arg, '(I256)') ic
      Call get_command_argument(2, arg, length, status)
      Read (arg, '(I256)') iextra

s1 = series(1, ic)
      s2 = series(4, ic)
      s3 = series(5, ic)
      s4 = series(6, ic)
      s = 4.0_wp*s1 - 2.0_wp*s2 - s3 - s4

Sync All
      If (this_image()==1) Then
        Do i = 2, num_images()
          s = s + s[ i ]
        End Do
        pid = mod(s, 1.0_wp) + 1.0_wp
        Call aihex(pid, nhx, chx)
        Write (output_unit, 100) ic, pid, chx
      End If
      Stop
100   Format ('Pi hex digit computation', /, 'position =', I12, ' + 1', /, &
          F20.15, /, 12A1)
      
    Contains
      Subroutine aihex(x, nhx, chx)
        ! This returns, in chx, the first nhx hex digits of the fraction of x.
        Real (Kind=wp) :: x, y
        Integer nhx, i
        Character (Len=1) chx(nhx), hx(0:15)
        Data hx/'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B', &
          'C', 'D', 'E', 'F'/

y = abs(x)
        Do i = 1, nhx
          y = 16.0_wp*(y-aint(y))
          chx(i) = hx(int(y))
        End Do
        Return
      End Subroutine

Function series(m, ic) Result(s)
        ! This routine evaluates the series  sum_k 16^(ic-k)/(8*k+m)
        ! using the modular exponentiation technique.
        Integer m, ic, k
        Real (Kind=wp) :: ak, p, s, t, corr, yy, tt, eps
        Parameter(eps=1.d-16)

s = 0.0_wp
        corr = 0.0_wp

! Sum the series up to ic.
        Do k = this_image() - 1, ic, num_images()
          ak = 8 * k + m
          p = ic - k
          if (p == 0.0_wp) then
            t = 1.0_wp
          else if (ak == 1.0_wp) then
            t = 0.0_wp
          else
            t = expm1(p, ak)
          end if
          s = mod(s + t / ak,1.0_wp)
          ! Kahan summation instead of above:
          ! yy = t/ak - corr
          ! tt = s + yy
          ! corr = (tt-s) - yy
          ! s = mod(tt, 1.0_wp)
          ! End Kahan
        End Do
        ! We don't do the extra (ca. 100) iterations
        ! If (this_image()==1) then
        !  Do k = ic + 1, ic+100
        !    ak = 8 * k + m
        !    t = 16.0_wp ** (ic - k)
        !    if (t < eps) Exit
        !    s = mod(s + t / ak, 1.0_wp)
        !  End Do
        ! End if
        Return
      End Function

Function expm1(p, ak)
        ! expm1 = 16^p mod ak.  This routine uses the left-to-right binary
        ! exponentiation scheme.  It is valid for  ak <= 2^24.
        Integer :: ntp, i
        Parameter (ntp=30)
        Real (Kind=wp) :: ak, expm1, p, p1, pt, r, tp(ntp), zero, one, two, &
          half, sixteen
        Parameter (zero=0.0_wp, one=1.0_wp, two=2.0_wp, half=0.5_wp, &
          sixteen=16.0_wp)
        Save tp
        Data tp/ntp*0.0_wp/

If (ak==one) Then
          expm1 = zero
          Return
        End If
        ! If this is the first call to expm1, fill the power of two table tp.
        If (tp(1)==zero) Then
          tp(1) = one
          Do i = 2, ntp
            tp(i) = two*tp(i-1)
          End Do
        End If
        ! Find the greatest power of two less than or equal to p.
        Do i = 1, ntp
          If (tp(i)>p) Exit
        End Do
        pt = tp(i-1)
        p1 = p
        r = one
        ! Perform binary exponentiation algorithm modulo ak.
        Do
          If (p1>=pt) Then
            r = mod(sixteen*r, ak)
            p1 = p1 - pt
          End If
          pt = half*pt
          If (pt>=one) Then
            r = mod(r*r, ak)
            Cycle
          Else
            Exit
          End If
        End Do
        expm1 = mod(r, ak)
        Return
      End Function
    End Program