program piform3

!   This evaluates the pi iteration formulas given in DHB's paper.

!   To run this program:

!   1. Download the latest version of MPFUN-MPFR (faster) or MPFUN-Fort (slower but
!      easier to compile) library from https://www.davidhbailey.com/dhbsoftware/.
!      Unpack the file if one's computer did not do this automatically.
!   2. If using the MPFR-MPFR package, you must first install both the GMP package
!      and the MPFR package. See the README.txt file in the main directory of the
!      MPFUN-MPFR package for details on installing GMP and MPFR.
!   3. Go to the fortran subdirectory of either MPFUN-MPFR or MPFUN-Fort. Edit the
!      file mpmodf.f90, setting the value of the parameter mpipl to 100000 (or more).
!   4. Type "./gnu-complib2.scr" if using the GNU gfortran compiler or one of the other
!      similar scripts if using another compiler.
!   5. Type "./gnu-complink2.scr piform3" if using GNU fortran or similar script for
!      another compiler.
!   6. Type "./piform3" > piform3.txt to execute this program and save output.

!   For other details on using MPFUN-MPFR and MPFUN-Fort, see the README.txt file
!   in the main MPFUN-MPFR or MPFUN-Fort directory, or the technical paper at
!   https://www.davidhbailey.com/dhbpapers/mpfun2015.pdf.

!   David H Bailey   22 Apr 2020

use mpmodule
implicit none
integer ndp1, neps1, nwds1
parameter (ndp1 = 100000, neps1 = -ndp1, nwds1 = int (ndp1 / mpdpw + 2))
real (8) second, tmx, tm0, tm1
type (mp_real) eps1, mppi1
external second

mppi1 = mppi (nwds1)

!  Perform iterations.

!  Formula 68:

write (6, 1) 68
1 format ('Formula',i3)
tm0 = second ()
call form68 (nwds1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0
2 format ('CPU time =',f12.6)

write (6, 1) 69
tm0 = second ()
call form69 (nwds1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

write (6, 1) 70
tm0 = second ()
call form70 (nwds1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

write (6, 1) 71
tm0 = second ()
call form71 (nwds1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

write (6, 1) 72
tm0 = second ()
call form72 (nwds1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

stop
end

subroutine form68 (nwds1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mppi1, one1, a0, a1, b0, b1

one1 = mpreal (1.d0, nwds1)
eps1 = mpreal (0.5d0, nwds1) ** (nwds1 * mpnbt - 20)
a0 = 2.d0 * sqrt (3.d0 * one1)
b0 = mpreal (3.d0, nwds1)

do k = 0, itr

  if (mod (k, 1000) == 0) write (6, *) k
  
  a1 = 2.d0 * a0 * b0 / (a0 + b0)
  b1 = sqrt (a1 * b0)
  if (abs (a1 - a0) < eps1) goto 101
  a0 = a1
  b0 = b1
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, a1 - mppi1)
return
end

subroutine form69 (nwds1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mppi1, one1, a0, a1, b0, b1, c1, p1, s0, s1

one1 = mpreal (1.d0, nwds1)
eps1 = mpreal (0.5d0, nwds1) ** (nwds1 * mpnbt / 2)
a0 = one1
b0 = 1.d0 / sqrt (2.d0 * one1)
s0 = mpreal (0.5d0, nwds1)

do k = 0, itr
  a1 = 0.5d0 * (a0 + b0)
  b1 = sqrt (a0 * b0)
  c1 = a1**2 - b1**2
  s1 = s0 - 2.d0**(k+1) * c1
  p1 = 2.d0 * a1**2 / s1
  if (abs (a1 - a0) < eps1) goto 101
  a0 = a1
  b0 = b1
  s0 = s1
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, p1 - mppi1)
return
end

subroutine form70 (nwds1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mppi1, one1, a0, a1, r0, r1, s0, s1

one1 = mpreal (1.d0, nwds1)
eps1 = mpreal (0.5d0, nwds1) ** (nwds1 * mpnbt / 3)
a0 = one1 / 3.d0
s0 = 0.5d0 * (sqrt (3.d0 * one1) - 1.d0)

do k = 0, itr
  r1 = 3.d0 / (1.d0 + 2.d0 * mpnrt (1.d0 - s0**3, 3))
  s1 = 0.5d0 * (r1 - 1.d0)
  a1 = r1**2 * a0 - 3**k * (r1**2 - 1.d0)
  if (abs (a1 - a0) < eps1) goto 101
  a0 = a1
  s0 = s1
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 1.d0 / a1 - mppi1)
return
end

subroutine form71 (nwds1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mppi1, one1, a0, a1, t1, y0, y1

one1 = mpreal (1.d0, nwds1)
eps1 = mpreal (0.5d0, nwds1) ** (nwds1 * mpnbt / 4)
a0 = 6.d0 - 4.d0 * sqrt (2.d0 * one1)
y0 = sqrt (2.d0 * one1) - 1.d0

do k = 0, itr
  t1 = mpnrt (1.d0 - y0**4, 4)
  y1 = (1.d0 - t1) / (1.d0 + t1)
  a1 = a0 * (1.d0 + y1)**4 - 2.d0**(2*k+3) * y1 * (1.d0 + y1 + y1**2)
  if (abs (a1 - a0) < eps1) goto 101
  a0 = a1
  y0 = y1
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 1.d0 / a1 - mppi1)
return
end

subroutine form72 (nwds1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mppi1, one1, a0, a1, r0, r1, s0, s1, t1, u1, v1, w1

one1 = mpreal (1.d0, nwds1)
eps1 = mpreal (0.5d0, nwds1) ** (nwds1 * mpnbt / 9)
a0 = one1 / 3.d0
r0 = 0.5d0 * (sqrt (3.d0 * one1) - 1.d0)
s0 = mpnrt (1.d0 - r0**3, 3)

do k = 0, itr
  t1 = 1.d0 + 2.d0 * r0
  u1 = mpnrt (9.d0 * r0 * (1.d0 + r0 + r0**2), 3)
  v1 = t1**2 + t1 * u1 + u1**2
  w1 = 27.d0 * (1.d0 + s0 + s0**2) / v1
  a1 = w1 * a0 + 3.d0**(2*k) * (1.d0 - w1) / 3.d0
  s1 = (1.d0 - r0)**3 / ((t1 + 2.d0 * u1) * v1)
  r1 = mpnrt (1.d0 - s1**3, 3)
  if (k > 0 .and. abs (a1 - a0) < eps1) goto 101
  a0 = a1
  r0 = r1
  s0 = s1
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 1.d0 / a1 - mppi1)
return
end

