program piqp c This program employs the recently discovered digit extraction scheme c to produce hex digits of pi. This code is valid up to ic = 2^24 on c systems with IEEE arithmetic. real*8 pid, s1, s2, s3, s4, series parameter (ic = 100000, nhx = 12) character*1 chx(nhx) c ic is the hex digit position -- output begins at position ic + 1. s1 = series (1, ic) s2 = series (4, ic) s3 = series (5, ic) s4 = series (6, ic) pid = mod (4.d0 * s1 - 2.d0 * s2 - s3 - s4, 1.d0) + 1.d0 call ihex (pid, nhx, chx) write (6, 1) ic, pid, chx 1 format ('Pi hex digit computation'/'position =',i12,' + 1'/ $ f20.15/12a1) stop end subroutine ihex (x, nhx, chx) c This returns, in chx, the first nhx hex digits of the fraction of x. real*8 x, y character*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.d0 * (y - aint (y)) chx(i) = hx(int(y)) enddo return end function series (m, ic) c This routine evaluates the series sum_k 16^(ic-k)/(8*k+m) c using the modular exponentiation technique. real*8 ak, eps, expm, p, s, series, t parameter (eps = 1d-17) s = 0.d0 c Sum the series up to ic. do k = 0, ic - 1 ak = 8 * k + m p = ic - k t = expm (p, ak) s = mod (s + t / ak, 1.d0) enddo c Compute a few terms where k >= ic. do k = ic, ic + 100 ak = 8 * k + m t = 16.d0 ** (ic - k) / ak if (t .lt. eps) goto 100 s = mod (s + t, 1.d0) enddo 100 continue series = s return end function expm (p, ak) c expm = 16^p mod ak. This routine uses the left-to-right binary c exponentiation scheme. It is valid for ak <= 2^24. parameter (ntp = 25) real*8 ak, expm, p, p1, pt, r, tp(ntp) save tp data tp/ntp*0.d0/ c If this is the first call to expm, 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 if (ak .eq. 1.d0) then expm = 0.d0 return endif c Find the greatest 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) p1 = p r = 1.d0 c Perform binary exponentiation algorithm modulo ak. 110 continue if (p1 .ge. pt) then r = mod (16.d0 * r, ak) p1 = p1 - pt endif pt = 0.5d0 * pt if (pt .ge. 1.d0) then r = mod (r * r, ak) goto 110 endif expm = mod (r, ak) return end