program piform1

!   This evaluates the pi summation formulas given in DHB's paper "A catalogue of
!   mathematical formulas involving Pi, with analysis," available at
!     https://www.davidhbailey.com/dhbpapers/pi-formulas.pdf

!   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 20000 (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 piform1" if using GNU fortran or similar script for
!      another compiler.
!   6. Type "./piform1" > piform1.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 = 20000, neps1 = -ndp1, nwds1 = int (ndp1 / mpdpw + 2))
real (8) second, tmx, tm0, tm1
type (mp_real) eps1, mplg1, mppi1
external second

eps1 = mpreal (10.d0, nwds1) ** neps1
mplg1 = mplog2 (nwds1)
mppi1 = mppi (nwds1)

!  Evaluate infinite series formulas.

!  Formula 02:

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

!  Formula 03:

write (6, 1) 3
tm0 = second ()
call form03 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 04:

write (6, 1) 4
tm0 = second ()
call form04 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 05:

write (6, 1) 5
tm0 = second ()
call form05 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 06:

write (6, 1) 6
tm0 = second ()
call form06 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 07:

write (6, 1) 7
tm0 = second ()
call form07 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 08:

write (6, 1) 8
tm0 = second ()
call form08 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 09:

write (6, 1) 9
tm0 = second ()
call form09 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 11:

write (6, 1) 11
tm0 = second ()
call form11 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 12:

write (6, 1) 12
tm0 = second ()
call form12 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 13:

write (6, 1) 13
tm0 = second ()
call form13 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 14:

write (6, 1) 14
tm0 = second ()
call form14 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 15:

write (6, 1) 15
tm0 = second ()
call form15 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 16:

write (6, 1) 16
tm0 = second ()
call form16 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 17:

write (6, 1) 17
tm0 = second ()
call form17 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 18:

write (6, 1) 18
tm0 = second ()
call form18 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 19:

write (6, 1) 19
tm0 = second ()
call form19 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 21:

write (6, 1) 21
tm0 = second ()
call form21 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 22:

write (6, 1) 22
tm0 = second ()
call form22 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 23:

write (6, 1) 23
tm0 = second ()
call form23 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 24:

write (6, 1) 24
tm0 = second ()
call form24 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 25:

write (6, 1) 25
tm0 = second ()
call form25 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 26:

write (6, 1) 26
tm0 = second ()
call form26 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 27:

write (6, 1) 27
tm0 = second ()
call form27 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 28:

write (6, 1) 28
tm0 = second ()
call form28 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 29:

write (6, 1) 29
tm0 = second ()
call form29 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 30:

write (6, 1) 30
tm0 = second ()
call form30 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 31:

write (6, 1) 31
tm0 = second ()
call form31 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 32:

write (6, 1) 32
tm0 = second ()
call form32 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 33:

write (6, 1) 33
tm0 = second ()
call form33 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 34:

write (6, 1) 34
tm0 = second ()
call form34 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 35:

write (6, 1) 35
tm0 = second ()
call form35 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 36:

write (6, 1) 36
tm0 = second ()
call form36 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 37:

write (6, 1) 37
tm0 = second ()
call form37 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 38:

write (6, 1) 38
tm0 = second ()
call form38 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 39:

write (6, 1) 39
tm0 = second ()
call form39 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 40:

write (6, 1) 40
tm0 = second ()
call form40 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 41:

write (6, 1) 41
tm0 = second ()
call form41 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 42:

write (6, 1) 42
tm0 = second ()
call form42 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 43:

write (6, 1) 43
tm0 = second ()
call form43 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 44:

write (6, 1) 44
tm0 = second ()
call form44 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 45:

write (6, 1) 45
tm0 = second ()
call form45 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 46:

write (6, 1) 46
tm0 = second ()
call form46 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 48:

write (6, 1) 48
tm0 = second ()
call form48 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

!  Formula 49:

write (6, 1) 49
tm0 = second ()
call form49 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

111 continue

!  Formula 50:

write (6, 1) 50
tm0 = second ()
call form50 (nwds1, eps1, mplg1, mppi1)
tm1 = second ()
write (6, 2) tm1 - tm0

stop
end

subroutine form02 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
is = 1
s1 = one1 / 2.d0
t2 = mpreal (2.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 4.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
is = 1
s2 = one1 / 3.d0
t2 = mpreal (3.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 9.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s2 = s2 + t3
  if (abs (t3) < eps1) goto 102
enddo

102 continue

k2 = k
write (6, 1) k1, k2
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 + s2 - 0.25d0 * mppi1)
return
end

subroutine form03 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
is = 1
s1 = one1 / 5.d0
t2 = mpreal (5.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 25.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
is = 1
s2 = one1 / 239.d0
t2 = mpreal (239.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 57121.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s2 = s2 + t3
  if (abs (t3) < eps1) goto 102
enddo

102 continue

k2 = k
write (6, 1) k1, k2
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 4.d0 * s1 - s2 - 0.25d0 * mppi1)
return
end

subroutine form04 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
is = 1
s1 = one1 / 2.d0
t2 = mpreal (2.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 4.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
is = 1
s2 = one1 / 5.d0
t2 = mpreal (5.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 25.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s2 = s2 + t3
  if (abs (t3) < eps1) goto 102
enddo

102 continue

k2 = k
is = 1
s3 = one1 / 8.d0
t2 = mpreal (8.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 64.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s3 = s3 + t3
  if (abs (t3) < eps1) goto 103
enddo

103 continue

k3 = k
write (6, 1) k1, k2, k3
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 + s2 + s3 - 0.25d0 * mppi1)
return
end

subroutine form05 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
is = 1
s1 = one1 / 4.d0
t2 = mpreal (4.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 16.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
is = 1
s2 = one1 / 20.d0
t2 = mpreal (20.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 400.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s2 = s2 + t3
  if (abs (t3) < eps1) goto 102
enddo

102 continue

k2 = k
is = 1
s3 = one1 / 1985.d0
t2 = mpreal (1985.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 3940225.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s3 = s3 + t3
  if (abs (t3) < eps1) goto 103
enddo

103 continue

k3 = k
write (6, 1) k1, k2, k3
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 3.d0 * s1 + s2 + s3 - 0.25d0 * mppi1)
return
end

subroutine form06 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
is = 1
s1 = one1 / 49.d0
t2 = mpreal (49.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 2401.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
is = 1
s2 = one1 / 57.d0
t2 = mpreal (57.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 3249.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s2 = s2 + t3
  if (abs (t3) < eps1) goto 102
enddo

102 continue

k2 = k
is = 1
s3 = one1 / 239.d0
t2 = mpreal (239.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 57121.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s3 = s3 + t3
  if (abs (t3) < eps1) goto 103
enddo

103 continue

k3 = k
is = 1
s4 = one1 / 110443.d0
t2 = mpreal (110443.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 12197656249.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s4 = s4 + t3
  if (abs (t3) < eps1) goto 104
enddo

104 continue

k4 = k
write (6, 1) k1, k2, k3, k4
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 12.d0*s1 + 32.d0*s2 - 5.d0*s3 + 12.d0*s4 - 0.25d0 * mppi1)
return
end

subroutine form07 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
is = 1
s1 = one1 / 57.d0
t2 = mpreal (57.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 3249.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
is = 1
s2 = one1 / 239.d0
t2 = mpreal (239.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 57121.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s2 = s2 + t3
  if (abs (t3) < eps1) goto 102
enddo

102 continue

k2 = k
is = 1
s3 = one1 / 682.d0
t2 = mpreal (682.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 465124.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s3 = s3 + t3
  if (abs (t3) < eps1) goto 103
enddo

103 continue

k3 = k
is = 1
s4 = one1 / 12943.d0
t2 = mpreal (12943.d0, nwds1)

do k = 1, itr
  is = - is
  t2 = 167521249.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s4 = s4 + t3
  if (abs (t3) < eps1) goto 104
enddo

104 continue

k4 = k
write (6, 1) k1, k2, k3, k4
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 44.d0*s1 + 7.d0*s2 - 12.d0*s3 + 24.d0*s4 - 0.25d0 * mppi1)
return
end

subroutine form08 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
is = 1
s1 = one1
t2 = one1

do k = 1, itr
  is = - is
  t2 = 3.d0 * t2
  t3 = is / ((2*k + 1)* t2)
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, sqrt (12.d0 * one1) * s1 - mppi1)
return
end

subroutine form09 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = -one1 / 12.d0
t1 = one1
t2 = 4.d0 * one1

do k = 1, itr
  t1 = t1 * 2 * (2*k - 1) / k
  t2 = 16.d0 * t2
  t3 = t1 / ((2*k + 3) * (2*k - 1) * t2)
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 3.d0 * sqrt (3.d0 * one1) / 4.d0 - 24.d0 * s1 - mppi1)
return
end

subroutine form11 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = 1103.d0 * one1
t1 = one1
t2 = one1

do k = 1, itr
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3)
  t2 = t2 * 24591257856.d0 * k**2 * k**2
  t3 = t1 * (1103.d0 + 26390.d0 * k) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 2.d0 * sqrt (2.d0 * one1) / 9801.d0 * s1 - 1.d0 / mppi1)
return
end

subroutine form12 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = 13591409.d0 * one1
t1 = one1
t2 = one1
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * (6*k) * (6*k - 1) * (6*k - 2) * (6*k - 3) * (6*k - 4) * (6*k - 5)
  t2 = t2 * (3*k) * (3*k - 1) * (3*k - 2) * k**3 * 640320 * 640320 * 640320
  t3 = is * t1 * (13591409.d0 + 545140134.d0 * k) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 12.d0 * s1 / (640320.d0 * sqrt (640320.d0*one1)) - 1.d0 / mppi1)
return
end

subroutine form13 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = 4.d0 - 0.5d0 - one1 / 5.d0 - one1 / 6.d0
t2 = one1

do k = 1, itr
  t2 = 16.d0 * t2
  t3 = (4.d0 * one1 / (8*k + 1) - 2.d0 * one1 / (8*k + 4) - one1 / (8*k + 5) - one1 / (8*k + 6)) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form14 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = one1
t2 = one1
is = 1

do k = 1, itr
  is = - is
  t2 = 4.d0 * t2
  t3 = is / ((2*k + 1) * t2)
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
is = 1
s2 = 32.d0 + 4.d0 + one1 / 3.d0
t2 = one1

do k = 1, itr
  is = - is
  t2 = 1024.d0 * t2
  t3 = is * (32.d0 * one1 / (4*k + 1) + 8.d0 * one1 / (4*k + 2) + one1 / (4*k + 3)) / t2
  s2 = s2 + t3
  if (abs (t3) < eps1) goto 102
enddo

102 continue

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

subroutine form15 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = 3.d0 + one1 / 3.d0
t2 = one1
is = 1

do k = 1, itr
  is = - is
  t2 = 4.d0 * t2
  t3 = is * (2.d0 * one1 / (4*k + 1) + 2.d0 * one1 / (4*k + 2) + one1 / (4*k + 3)) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form16 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = 16.d0 - 6.d0 - 8.d0 * one1 / 9.d0 - 6.d0 * one1 / 16.d0 + one1 / 25.d0
t2 = one1

do k = 1, itr
  t2 = 64.d0 * t2
  t3 = (16.d0 * one1 / (6*k + 1)**2 - 24.d0 * one1 / (6*k + 2)**2 &
    - 8.d0 * one1 / (6*k + 3)**2 - 6.d0 * one1 / (6*k + 4)**2 &
    + one1 / (6*k + 5)**2) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form17 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = 81.d0 - 27.d0 - 9.d0 * one1 / 4.d0 - 2.d0 - 3.d0 * one1 / 7.d0 - one1 / 4.d0 - one1 / 10.d0
t2 = one1

do k = 1, itr
  t2 = 729.d0 * t2
  t3 = (81.d0 * one1 / (12*k + 1) - 54.d0 * one1 / (12*k + 2) &
    - 9.d0 * one1 / (12*k + 4) - 12.d0 * one1 / (12*k + 6) &
    - 3.d0 * one1 / (12*k + 7) - 2.d0 * one1 / (12*k + 8) &
    - one1 / (12*k+10)) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 / 9.d0 - mppi1 * sqrt (3.d0 * one1))
return
end

subroutine form18 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = 243.d0 - 405.d0 * one1 / 4.d0 - 81.d0 * one1 / 16.d0 - 27.d0 * one1 / 25.d0 &
  - 2.d0 - 9.d0 * one1 / 49.d0 - 9.d0 * one1 / 64.d0 - 5.d0 * one1 / 100.d0 &
  + one1 / 121.d0
t2 = one1

do k = 1, itr
  t2 = 729.d0 * t2
  t3 = (243.d0 * one1 / (12*k + 1)**2 - 405.d0 * one1 / (12*k + 2)**2 &
    - 81.d0 * one1 / (12*k + 4)**2 - 27.d0 * one1 / (12*k + 5)**2 &
    - 72.d0 * one1 / (12*k + 6)**2 - 9.d0 * one1 / (12*k + 7)**2 &
    - 9.d0 * one1 / (12*k + 8)**2 - 5.d0 * one1 / (12*k + 10)**2 &
    + one1 / (12*k + 11)**2) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form19 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 4.d0 * t1
  t2 = t2 * (4*k - 1) * (4*k - 2) * (4*k - 3) / k / (2*k - 1)**2
  t3 = k * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, 12.d0 * s1 - (3.d0 * mppi1 + 8.d0))
return
end

subroutine form21 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (8.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 8.d0 * t1
  t2 = t2 * (6*k) * (6*k - 1) * (6*k - 2) * (6*k - 3) * (6*k - 4) * (6*k - 5) &
    / (3*k)**2 / (3*k - 1)**2 / (3*k - 2)**2
  t3 = (126.d0 * k**2 - 24.d0 * k + 8.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (15.d0 * mppi1 + 52.d0))
return
end

subroutine form22 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (-16.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 16.d0 * t1
  t2 = t2 * (8*k) * (8*k - 1) * (8*k - 2) * (8*k - 3) * (8*k - 4) * (8*k - 5) &
    * (8*k - 6) * (8*k - 7) / (4*k)**2 / (4*k - 1)**2 / (4*k - 2)**2 / (4*k - 3)**2
  t3 = mpreald (1920.d0 * k**2 * k - 928.d0 * k**2 + 424.d0 * k - 16.d0, nwds1) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (105.d0 * mppi1 + 304.d0))
return
end

subroutine form23 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = one1
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 8.d0 * t1
  t2 = t2 * 3 * (3*k) * (3*k - 1) * (3*k - 2) / (2*k**2) / (2*k - 1)
  t3 = (49.d0 * k + 1.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (16.d0 * mppi1 * sqrt (3.d0 * one1) + 81.d0))
return
end

subroutine form24 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (338.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 8.d0 * t1
  t2 = t2 * 3 * (3*k) * (3*k - 1) * (3*k - 2) / (2*k**2) / (2*k - 1)
  t3 = (-245.d0 * k + 338.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (162.d0 - 6.d0 * mppi1 * sqrt (3.d0 * one1) &
  - 18.d0 * log (3.d0 * one1)))
return
end

subroutine form25 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (-6.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t2 = t2 * 2 * (3*k) * (3*k - 1) * (3*k - 2) / (2*k**2) / (2*k - 1)
  t3 = (50.d0 * k - 6.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form26 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * 4 * (2*k)**2 * (2*k - 1)**2 * (3*k) * (3*k - 1) * (3*k - 2)
  t2 = t2 * (6*k) * (6*k - 1) * (6*k - 2) * (6*k - 3) * (6*k - 4) * (6*k - 5) * k
  t3 = is * (201.d0 - 952.d0 * k) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (15.d0 * mppi1 + 42.d0))
return
end

subroutine form27 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (-17.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = t1 * 8 * (2*k)**2 * (2*k - 1)**2 * (3*k) * (3*k - 1) * (3*k - 2)
  t2 = t2 * (6*k) * (6*k - 1) * (6*k - 2) * (6*k - 3) * (6*k - 4) * (6*k - 5) * k
  t3 = (350.d0 * k - 17.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (15.d0 * mppi1 * sqrt (2.d0 * one1) + 27.d0))
return
end

subroutine form28 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * 27 * (2*k)**2 * (2*k - 1)**2 * (3*k) * (3*k - 1) * (3*k - 2)
  t2 = t2 * (6*k) * (6*k - 1) * (6*k - 2) * (6*k - 3) * (6*k - 4) * (6*k - 5) * k
  t3 = is * (81.d0 - 1080.d0 * k) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (40.d0 * mppi1 * sqrt (3.d0 * one1) + 243.d0))
return
end

subroutine form29 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * (2*k)**2 * (2*k - 1)**2 * (3*k) * (3*k - 1) * (3*k - 2)
  t2 = t2 * 3 * (6*k) * (6*k - 1) * (6*k - 2) * (6*k - 3) * (6*k - 4) * (6*k - 5) * k
  t3 = is * (4123.d0 - 22100.d0 * k) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (20.d0 * mppi1 * sqrt (3.d0 * one1) + 89.d0))
return
end

subroutine form30 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1
is = 1

do k = 1, itr
  is = - is
  t2 = t2 * 2 * (5*k) * (5*k - 1) * (5*k - 2) * (5*k - 3) * (5*k - 4) &
    / (2*k) / (2*k - 1) / (3*k) / (3*k - 1) / (3*k - 2)
  t3 = is * mpreald (89012.d0 * k**2 * k - 77362.d0 * k**2 + 482.d0 * k + 3028.d0, nwds1) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (15.d0 * mppi1 + 240.d0 * mplg1 - 528.d0))
return
end

subroutine form31 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 9.d0 * t1
  t2 = t2 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3) &
    / k / (3*k) / (3*k - 1) / (3*k - 2)
  t3 = t1 * mpreald (2743.d0 * k * k - 130971.d0 * k - 12724.d0, nwds1) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (24516.d0 - 360.d0 * mppi1 * sqrt (3.d0 * one1)))
return
end

subroutine form32 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 8.d0 * t1
  t2 = t2 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3) &
    / k / (3*k) / (3*k - 1) / (3*k - 2)
  t3 = t1 * mpreald (430.d0 * k * k - 6240.d0 * k - 520.d0, nwds1) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (45.d0 * mppi1 + 1164.d0))
return
end

subroutine form33 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 3.d0 * t1
  t2 = t2 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3) &
    / k / (3*k) / (3*k - 1) / (3*k - 2)
  t3 = t1 * mpreald (7175.d0 * k * k - 15215.d0 * k + 480.d0, nwds1) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (40.d0 * mppi1 * sqrt (3.d0 * one1) + 1872.d0))
return
end

subroutine form34 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 9.d0 * t1
  t2 = t2 * 8 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3) &
    / k / (3*k) / (3*k - 1) / (3*k - 2)
  t3 = t1 * mpreald (5692.d0 + 6335.d0 * k - 5415.d0 * k * k, nwds1) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (288.d0 * mppi1 * sqrt (3.d0 * one1) &
  - 576.d0 * mplg1 + 324.d0))
return
end

subroutine form35 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = 9.d0 * t1
  t2 = t2 * 8 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3) &
    / k / (3*k) / (3*k - 1) / (3*k - 2)
  t3 = t1 * mpreald (7517.d0 + 1145.d0 * k + 18050.d0 * k * k, nwds1) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - (1008.d0 * mppi1 * sqrt (3.d0 * one1) &
  - 576.d0 * mplg1 + 7587.d0))
return
end

subroutine form36 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (5.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = t1 * 2 * (2*k - 1) / k
  t2 = 4096.d0 * t2
  t3 = t1**3 * (42.d0 * k + 5.d0) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form37 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (1.5d0, nwds1)
t1 = one1
t2 = mpreal (2.d0, nwds1)
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3)
  t2 = t2 * 1024.d0 * k**2 * k**2
  t3 = is * (20.d0 * k + 3.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form38 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (23.d0, nwds1) / 18.d0
t1 = one1
t2 = mpreal (18.d0, nwds1)
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3)
  t2 = t2 * 82944.d0 * k**2 * k**2
  t3 = is * (260.d0 * k + 23.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form39 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (1123.d0, nwds1) / 882.d0
t1 = one1
t2 = mpreal (882.d0, nwds1)
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3)
  t2 = t2 * 199148544.d0 * k**2 * k**2
  t3 = is * (21460.d0 * k + 1123.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form40 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = one1 / 3.d0
t1 = one1
t2 = mpreal (3.d0, nwds1)

do k = 1, itr
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3)
  t2 = t2 * 2304.d0 * k**2 * k**2
  t3 = (8.d0 * k + 1.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

k1 = k
write (6, 1) k1
1 format ('Iterations, error =',4i10)
call mpwrite (6, 80, 60, s1 - 2.d0 / (mppi1 * sqrt (3.d0 * one1)))
return
end

subroutine form41 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = one1 / 9.d0
t1 = one1
t2 = mpreal (9.d0, nwds1)

do k = 1, itr
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3)
  t2 = t2 * 20736.d0 * k**2 * k**2
  t3 = (10.d0 * k + 1.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form42 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (0.75d0, nwds1)
t1 = one1
t2 = mpreal (4.d0, nwds1)
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3)
  t2 = t2 * 12288.d0 * k**2 * k**2
  t3 = is * (28.d0 * k + 3.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form43 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (41.d0, nwds1) / 72.d0
t1 = one1
t2 = mpreal (72.d0, nwds1)
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3)
  t2 = t2 * 6635520.d0 * k**2 * k**2
  t3 = is * (644.d0 * k + 41.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form44 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (3.d0, nwds1) / 49.d0
t1 = one1
t2 = mpreal (49.d0, nwds1)

do k = 1, itr
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3)
  t2 = t2 * 614656.d0 * k**2 * k**2
  t3 = (40.d0 * k + 3.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form45 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (3.d0, nwds1)
t1 = one1
t2 = one1

do k = 1, itr
  t1 = t1 * (4*k) * (4*k - 1) * (4*k - 2) * (4*k - 3) &
    / (2*k)**2 / (2*k - 1)**2 * 16 * (2*k - 1)**2 * (2*k - 1)**2 / k**2 / k**2
  t2 = t2 * 65536.d0
  t3 = (120.d0 * k * k + 34.d0 * k + 3.d0) * t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form46 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (13.d0, nwds1)
t1 = one1
t2 = one1
is = 1

do k = 1, itr
  is = - is
  t1 = t1 * 2 * (2*k - 1) / k
  t2 = t2 * 1048576.d0
  t3 = is * (820.d0 * k * k + 180.d0 * k + 13.d0) * t1**5 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form48 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = one1
t1 = one1
t2 = one1

do k = 1, itr
  t1 = t1 * 2 * (2*k - 1) / k
  t2 = t2 * 256.d0
  t3 = (6.d0 * k + 1.d0) * t1**3 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form49 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (2.d0, nwds1)
t1 = mpreal (2.d0, nwds1)
t2 = one1

do k = 1, itr
  t1 = 2.d0 * t1
  t2 = t2 * 2 * (2*k - 1) / k
  t3 = t1 / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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

subroutine form50 (nwds1, eps1, mplg1, mppi1)
use mpmodule
implicit none
integer is, itr, k, k1, k2, k3, k4, nwds1
parameter (itr = 1000000)
type (mp_real) eps1, mplg1, mppi1, one1, s1, s2, s3, s4, t1, t2, t3

one1 = mpreal (1.d0, nwds1)
s1 = mpreal (9.d0, nwds1) / 1000.d0
t1 = one1
t2 = mpreal (1000.d0, nwds1)

do k = 1, itr
  t1 = t1 * (6*k) * (6*k - 1) * (6*k - 2) * (6*k - 3) * (6*k - 4) * (6*k - 5)
  t2 = t2 * 1000000.d0 * k**2 * k**2 * k**2
  t3 = t1 * (532.d0 * k**2 + 126.d0 * k + 9.d0) / t2
  s1 = s1 + t3
  if (abs (t3) < eps1) goto 101
enddo

101 continue

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