program piqp

!   This program employs the recently discovered digit extraction scheme due
!   to Plouffe and P. Borwein to produce hex digits of pi.
!   On IBM RS6000 systems, compile with -O3 -qstrict.

!   David H. Bailey     950919

integer*8 ic, ic4
real*8 second, tm0, tm1
real*16 pid, s1, s2, s3, s4, series
parameter (ic = 9999999, nbt = 96)
! parameter (ic = 10000000, nbt = 96)
integer ibt(nbt)
character*1 chx(nbt/4)

!   ic is the hex digit position -- output begins at position ic + 1.

write (6, 1) ic
1 format ('Pi hex digit computation'/'position =',i12,' + 1')
ic4 = 4 * ic

!   Evaluate the four series used in the formula.

tm0 = second ()
s1 = series (1, ic4)
call ibits (s1, nbt, ibt, chx)
write (6, 2) 1, s1, chx, ibt
2 format ('Sum of series',i2,' ='/f34.30/24a1/48i1.1/48i1.1)
s2 = series (4, ic4)
call ibits (s2, nbt, ibt, chx)
write (6, 2) 2, s2, chx, ibt
s3 = series (5, ic4)
call ibits (s3, nbt, ibt, chx)
write (6, 2) 3, s3, chx, ibt
s4 = series (6, ic4)
call ibits (s4, nbt, ibt, chx)
write (6, 2) 4, s4, chx, ibt

!   Compute hex digits of pi^2 and log(2)^2 using P-B's formulas.

pid = mod (4.q0 * s1 - 2.q0 * s2 - s3 - s4, 1.q0)
if (pid .lt. 0.q0) pid = pid + 1.q0
tm1 = second ()

!   Output results.

write (6, 3) tm1 - tm0
3 format ('Run complete.  CPU time =',f12.4)
call ibits (pid, nbt, ibt, chx)
write (6, 4) pid, chx, ibt
4 format ('Pi ='/f34.30/24a1/48i1.1/48i1.1)

stop
end

subroutine ibits (x, nbt, ibt, chx)

!   This returns in ibt the first nbt bits of the fraction of the quad
!   precision number x, and returns the first nbt/4 hex digits in chx.
!   nbt must be divisible by 4 and must not exceed 104.

real*16 x, y
integer ibt(nbt)
character*1 chx(nbt/4), 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, nbt
  y = 2.d0 * (y - aint (y))
  ibt(i) = y
enddo

do i = 1, nbt / 4
  ii = 8 * ibt(4*i-3) + 4 * ibt(4*i-2) + 2 * ibt(4*i-1) + ibt(4*i)
  chx(i) = hx(ii)
enddo

return
end

function series (m, ic)

!   This routine evaluates the series  sum_k 2^(ic-4*k)/(8*k+m)
!   using the Plouffe-Borwein trick and some advanced exponentiation schemes.

integer*8 ic, j, k
real*8 aj, ak, eps, expm1, expm2, p, pt, tp, t12, t24, t48
real*16 s, series, t
parameter (ntp = 50, eps = 1d-33, t12 = 4096.d0, t24 = t12 * t12, &
  t48 = t24 * t24)
dimension tp(ntp)
save tp
data tp/ntp*0.d0/

!   If this is the first call, fill the power of two table tp.

if (tp(1) .eq. 0.d0) then
  tp(1) = 1.d0

  do i = 2, ntp
    tp(i) = 2.d0 * tp(i-1)
  enddo
endif

s = 0.q0

!  Sum the series until k exceeds ic/4.

do k = 0, (ic - 1) / 4 
  ak = 8 * k + m
  if (ak .eq. 1.d0) goto 110
  p = ic - 4 * k

!   Find the largest power of two less than or equal to p.

  do i = 1, ntp
    if (tp(i) .gt. p) goto 100
  enddo

100 pt = tp(i-1)

!   Perform one of three modular exponentiation schemes, depending on ak.

  if (ak .le. t24) then
    t = expm1 (p, pt, ak)
  elseif (ak .le. t48) then
    t = expm2 (p, pt, ak)
  endif

  s = mod (s + t / ak, 1.q0)

110 continue
enddo

!   Complete the computation by manually computing a few terms where k > ic/2.

do j = k, k + 100
  aj = 8 * j + m
  t = 2.q0 ** (ic - 4 * j) / aj
  if (t .lt. eps) goto 120
  s = mod (s + t, 1.q0)
enddo

120 continue

series = s
return
end

function expm1 (p, pt, ak)

!   expm2 = 2^p mod ak.    pt is the largest power of two <= p.
!   This routine uses a straightforward DP left-to-right binary exponentiation
!   scheme.  This routine is valid for  ak <= 2^24.

real*8 ak, expm1, p, p1, pt, pt1, r

p1 = p
pt1 = pt
r = 1.d0

100 continue

if (p1 .ge. pt1) then
  r = 2.d0 * r
  p1 = p1 - pt1
endif
pt1 = 0.5d0 * pt1
if (pt1 .ge. 1.d0) then
  r = mod (r * r, ak)
  goto 100
endif

expm1 = mod (r, ak)

return
end

function expm2 (p, pt, ak)

!   expm2 = 2^p mod ak.    pt is the largest power of two <= p.
!   This routine uses a straightforward DP left-to-right binary exponentiation
!   scheme.  This routine is valid for  ak <= 2^48.

!   This routine employs an advanced trick for obtaining quad precision
!   products that only works on certain systems, including the IBM RS6000 and 
!   the SGI R8000.

real*8 ak, expm2, p, p1, pt, pt1, r, t1, t2, t3, t4, t5

p1 = p
pt1 = pt
r = 1.d0

100 continue

if (p1 .ge. pt1) then
  r = 2.d0 * r
  p1 = p1 - pt1
endif
pt1 = 0.5d0 * pt1
if (pt1 .ge. 1.d0) then
!   r = mod (r * r, ak)
  t1 = r * r
  t2 = r * r - t1
  t3 = anint (t1 / ak)
  t4 = ak * t3
  t5 = ak * t3 - t4
  r = (t1 - t4) + (t2 - t5)
  goto 100
endif

!   Make sure result is in the range [0, ak).

r = mod (r, ak)
if (r .lt. 0.d0) r = r + ak
expm2 = r

return
end