c   This function 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.

      function pidf (ic)

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

      real*8 s1, s2, s3, s4, series, pidf

      s1 = series (1, ic)
      s2 = series (4, ic)
      s3 = series (5, ic)
      s4 = series (6, ic)
c      pidf = mod (4.d0 * s1 - 2.d0 * s2 - s3 - s4, 1.d0) + 1.d0
      pidf = mod (4.d0 * s1 - 2.d0 * s2 - s3 - s4, 1.d0)
      if (pidf .lt. 0) then
        pidf = pidf + 1.d0
      endif

      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
