module quadglobal
use mpmodule
implicit none
type (mp_real) mplg1, mppi1, sqrtpi1
end module

program piform2

!   This evaluates the 1-D integral pi 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 4000 (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 piform2" if using GNU fortran or similar script for
!      another compiler.
!   6. Type "./piform2" > piform2.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
use quadglobal
implicit none
integer i, is, j, k, ndp1, ndp2, neps1, neps2, nq1, nq2, nwds1, nwds2, n1
parameter (ndp1 = 4000, ndp2 = 8000, neps1 = -ndp1, neps2 = -ndp2, &
  nq1 = 13, nq2 = 12 * 2 ** nq1, nwds1 = int (ndp1 / mpdpw + 2), &
  nwds2 = int (ndp2 / mpdpw + 2))
real (8) dplog10q, d1, d2, second, tm0, tm1, tm2
type (mp_real) dfun2, err, fun51, fun52, fun53, fun54, fun55, fun56, &
  fun57, fun58, fun59, fun60, fun61, fun62, fun63, fun64, fun65, fun66, &
  fun67, fun68, fun69, fun70, quades, quadss, quadts, &
  one1, one2, mppi2, t1, t2, t3, &
  wkes(-1:nq2), xkes(-1:nq2), wkss(-1:nq2), xkss(-1:nq2), &
  wkts(-1:nq2), xkts(-1:nq2), x1, x2, zero1, zero2
external dfun2, fun51, fun52, fun53, fun54, fun55, fun56, fun57, fun58, &
  fun59, fun60, fun61, fun62, fun63, fun64, fun65, fun66, fun67, fun68, &
  fun69, fun70, quades, quadss, quadts, second

!   Check to see if default precision is high enough.

if (ndp2 > mpipl) then
  write (6, '("Increase default precision in module MPFUNF.")')
  stop
endif

!   Compute 0, 1, log(2) and pi to low (nwds1) and high precision (nwds2) words.

one1 = mpreal (1.d0, nwds1)
one2 = mpreal (1.d0, nwds2)
mplg1 = mplog2 (nwds1)
mppi1 = mppi (nwds1)
mppi2 = mppi (nwds2)
sqrtpi1 = sqrt (mppi1)
zero1 = mpreal (0.d0, nwds1)
zero2 = mpreal (0.d0, nwds2)
write (6, 1) nq1, ndp1, ndp2, neps1, neps2
1 format ('Quadts test:  Quadlevel =',i6/'Digits1 =',i6,'  Digits2 =',i6, &
  '  Epsilon1 =',i6,'  Epsilon2 =',i6)

!   Initialize quadrature tables wk and xk (weights and abscissas).

tm0 = second ()
call initqes (nq1, nq2, nwds1, neps2, wkes, xkes)
call initqss (nq1, nq2, nwds1, neps2, wkss, xkss)
call initqts (nq1, nq2, nwds1, neps2, wkts, xkts)
tm1 = second ()
tm2 = tm1 - tm0
write (6, 2) tm1 - tm0
2 format ('Quadrature initialization completed: cpu time =',f12.6)

!   Begin quadrature calculations.

write (6, 11) 51
11 format (/'Formula',i6)
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun51, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
3 format ('Quadrature completed: CPU time =',f12.6/'Result =')
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = 0.25d0 * mppi1
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1
4 format ('Actual error =',f10.6,'x10^',i6)

write (6, 11) 52
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun52, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = 0.25d0 * mppi1
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 53
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun53, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = (mppi1 - 2.d0) / 4.d0
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 54
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun54, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = mppi1 * (mppi1 - 12.d0) / 48.d0 + mplg1 / 2.d0
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 55
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun55, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = mpreal (22.d0, nwds1) / 7.d0 - mppi1
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 56
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun56, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = 0.125d0 * mppi1
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 57
x1 = zero2
tm0 = second ()
t1 = quades (fun57, x1, nq1, nq2, nwds1, nwds2, neps1, wkes, xkes)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = mppi1 * (1.d0 + 2.d0 * mplg1) / 8.d0
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 58
x1 = zero2
tm0 = second ()
t1 = quades (fun58, x1, nq1, nq2, nwds1, nwds2, neps1, wkes, xkes)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = 4.d0 * mppi1 * mplg1**2 + mppi1**3 / 3.d0
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 59
x1 = zero2
x2 = 0.5d0 * mppi2
tm0 = second ()
t1 = quadts (fun59, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = mppi1 * mplg1
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 60
x1 = zero2
x2 = 0.5d0 * mppi2
tm0 = second ()
t1 = quadts (fun60, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = mppi1**3 / 24.d0 + mppi1 * mplg1**2 / 2.d0
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 61
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun61, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = 8.d0 * mppi1**3 / (81.d0 * sqrt (3.d0 * one1))
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 62
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun62, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = 0.5d0 * mppi1 - mplg1
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 63
x1 = zero2
tm0 = second ()
t1 = quades (fun63, x1, nq1, nq2, nwds1, nwds2, neps1, wkes, xkes)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = sqrt (mppi1)
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 64
x1 = zero2
tm0 = second ()
x2 = one2
! t1 = quades (fun64, x1, nq1, nq2, nwds1, nwds2, neps1, wkes, xkes)
! t1 = 2.d0 * t1 / sqrtpi
t1 = quadts (fun64, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
t1 = 3.d0 * (t1 - mppi1 / 4.d0)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = - mppi1 / 4.d0 + 3.d0 * log (2.d0 + sqrt (3.d0 * one1)) / 2.d0
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 65
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun65, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
t1 = 0.5d0 * (t1 - mppi1 / 4.d0)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = -mppi1 / 24.d0 + sqrt (3.d0 * one1) / 4.d0 + log (2.d0 + sqrt (3.d0 * one1)) / 2.d0
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 66
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun66, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
t1 = one1 / 5.d0 * (t1 - mppi1 / 4.d0)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = -mppi1 / 60.d0 + 2.d0 * sqrt (3.d0 * one1) / 5.d0 &
  + 7.d0 * log (2.d0 + sqrt (3.d0 * one1)) / 20.d0
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 11) 67
x1 = zero2
x2 = one2
tm0 = second ()
t1 = quadts (fun67, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)
tm1 = second ()
tm2 = tm2 + (tm1 - tm0)
write (6, 3) tm1 - tm0
! call mpwrite (6, ndp1 + 20, ndp1, t1)
t2 = 5.d0 - mppi1**2 - 4.d0 * mplg1 + 16.d0 * mplg1**2
call decmdq (t2 - t1, d1, n1)
write (6, 4) d1, n1

write (6, 99) tm2
99 format ('Total CPU time =',f12.6)

stop
end

function fun51 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun51, t1, t2
type (mp_real) t

t1 = mpreal (t, nwds1)
fun51 = 1.d0 / (1.d0 + t1**2)
return
end

function fun52 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun52, t1
type (mp_real) t

t1 = mpreal (t, nwds1)
fun52 = sqrt (1.d0 - t1**2)
return
end

function fun53 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun53, t1
type (mp_real) t

t1 = mpreal (t, nwds1)
fun53 = t1 * atan (t1)
return
end

function fun54 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun54, t1, t2
type (mp_real) t

t1 = mpreal (t, nwds1)
fun54 = log (t1) * atan (t1)
return
end

function fun55 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun55, t1, t2
type (mp_real) t

t1 = mpreal (t, nwds1)
t2 = t1**2
fun55 = t2**2 * (1.d0 - t1)**4 / (1.d0 + t2)
return
end

function fun56 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun56, t1, t2, t3
type (mp_real) t

t1 = mpreal (t, nwds1)
t2 = t1**2
t3 = mpreal (1.d0 - t**4, nwds1)
if (t3 > 0.d0) then
  fun56 = t2 / ((1.d0 + t2**2) * sqrt (t3))
else
  fun56 = mpreal (0.d0, nwds1)
endif
return
end

function fun57 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun57, t1, t2
type (mp_real) t

t1 = mpreal (t, nwds1)
t2 = exp (-t1)
fun57 = t1 * t2 * sqrt (1.d0 - t2**2)
return
end

function fun58 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun58, t1, t2
type (mp_real) t

t1 = mpreal (t, nwds1)
t2 = mpreal (exp (t1) - 1.d0, nwds1)
if (t2 > 0.d0) then
  fun58 = t1**2 / sqrt (t2)
else
  fun58 = mpreal (0.d0, nwds1)
endif
return
end

function fun59 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun59, t1
type (mp_real) t

t1 = mpreal (t, nwds1)
fun59 = t1**2 / sin (t1)**2
return
end

function fun60 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun60, t1
type (mp_real) t

t1 = mpreal (t, nwds1)
fun60 = log (cos (t1)) ** 2
return
end

function fun61 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun61, t1, t2
type (mp_real) t

t1 = mpreal (t, nwds1)
fun61 = log (t1)**2 / (t1**2 + t1 + 1.d0)
return
end

function fun62 (t, nwds1, nwds2)
use mpmodule
implicit none
integer is, k, nwds1, nwds2
type (mp_real) eps1, fun62, t1, t2, t3, t4, t5
type (mp_real) t

eps1 = mpreal (0.5d0, nwds1) ** (nwds1 * mpnbt)
t1 = mpreal (t, nwds1)
t2 = t1**2
if (t2 > 0.5d0**100) then
  fun62 = log (1.d0 + t2) / t2
else
  t3 = mpreal (1.d0, nwds1)
  t4 = mpreal (1.d0, nwds1)
  is = 1
  
  do k = 2, 1000000
    is = - is
    t4 = t2 * t4
    t5 = is * t4 / k
    t3 = t3 + t5
    if (abs (t5) < eps1) goto 101
  enddo

101 continue

  fun62 = t3
endif
return
end

function fun63 (t, nwds1, nwds2)
use mpmodule
implicit none
integer nwds1, nwds2
type (mp_real) fun63, t, t1, t2, t3

t1 = mpreal (t, nwds1)
fun63 = 1.d0 / (sqrt (t1) * exp (t1))
return
end

function fun64 (t, nwds1, nwds2)
use mpmodule
use quadglobal
implicit none
integer nwds1, nwds2
type (mp_real) fun64, t, t1, t2, t3

t1 = mpreal (t, nwds1)
t2 = t1**2 + 1.d0
fun64 = sqrt (t2 + 1.d0) / t2
return
end

function fun65 (t, nwds1, nwds2)
use mpmodule
use quadglobal
implicit none
integer is, j, k, nwds1, nwds2
type (mp_real) eps, fun65, t, t0, t1, t2, t3, t4, t5

t1 = t
t2 = t1**2 + 1.d0
fun65 = (t2 + 1.d0) * sqrt (t2 + 1.d0) / t2
return
end

function fun66 (t, nwds1, nwds2)
use mpmodule
use quadglobal
implicit none
integer is, j, k, nwds1, nwds2
type (mp_real) eps, fun66, t, t0, t1, t2, t3, t4, t5

t1 = t
t2 = t1**2 + 1.d0
fun66 = (t2 + 1.d0)**2 * sqrt (t2 + 1.d0) / t2
return
end

function fun67 (t, nwds1, nwds2)
use mpmodule
use quadglobal
implicit none
integer nwds1, nwds2
type (mp_real) fun67, t, t1, t2, t3, t4

t1 = mpreal (t, nwds1)
t2 = t1**2
t3 = t1 * t2
t4 = t1 * t3
if (t1 > 0.d0 .and. t1 < 1.d0) then
  fun67 = (t1*(-7 + 4*mplg1) + t2*(1 + 20*mplg1) + t3*(3 + 12*mplg1) + &
    t4*(3 - 4*mplg1) - 4*(1 + t1)*(-1 + 4*t1 + t2)*log(1 + t1))&
    /(t1*(t1 - 1)*(t1 + 1)**2)
else
  fun67 = mpreal (0.d0, nwds1)
endif
return
end

subroutine initqes (nq1, nq2, nwds1, neps2, wkes, xkes)

!   This subroutine initializes the quadrature arrays xkes and wkes for quades.
!   The argument nq2 is the space allocated for wkes and xkes in the calling
!   program.  By default it is set to 12 * 2^nq1.  If initqes outputs the message
!   "Table space parameter is too small", adjust nq2.  Also, if quades outputs
!   the message "Terms too large", adjust nq1 and neps2 as necessary in the call
!   to initqes. The argument neps2 controls termination of the loop below, which
!   ends when wkes(k) * 10^(neps2) > 1.

!   The wkes and xkes arrays are computed based on the transformation
!   t = exp (pi/2 * sinh (x)).  See comments below.

!   Both initqes and quades are 100% THREAD SAFE -- all requisite parameters
!   and arrays are passed through subroutine arguments. 

!   David H Bailey   31 Jan 2020

use mpmodule
implicit none
integer i, ierror, iprint, j, k, k1, ndebug, neps2, nq1, nq2, nwds1
parameter (iprint = 1024, ndebug = 2)
type (mp_real) eps2, h, p2, t1, t2, t3, t4, u1, u2, &
  wkes(-1:nq2), xkes(-1:nq2)

write (6, 1)
1 format ('initqes: Exp-sinh quadrature initialization')

eps2 = mpreal (10.d0, nwds1) ** neps2
p2 = 0.5d0 * mppi (nwds1)
h = mpreal (0.5d0 ** nq1, nwds1)
wkes(-1) = mpreal (dble (nq1), nwds1)

do k = 0, nq2
  if (ndebug >= 2 .and. mod (k, iprint) == 0) write (6, *) k, nq2
    t1 = mpreal (dble (k) * h, nwds1)

!   xkes(k) = exp (u1)
!   wkes(k) = exp (u1) * u2
!   where u1 = pi/2 * sinh (t1) and u2 = pi/2 * cosh (t1)

  t2 = exp (t1)
  u1 = 0.5d0 * p2 * (t2 - 1.d0 / t2)
  u2 = 0.5d0 * p2 * (t2 + 1.d0 / t2)
  xkes(k) = exp (u1)
  wkes(k) = xkes(k) * u2

  if (wkes(k) * eps2 > 1.d0) goto 100
enddo

write (6, 2) nq2
2 format ('initqes: Table space parameter is too small; value =',i8)
stop

100 continue

xkes(-1) = mpreal (dble (k), nwds1)
if (ndebug >= 2) then
  write (6, 3) k
3 format ('initqes: Table spaced used =',i8)
endif

return
end

function quades (fun, x1, nq1, nq2, nwds1, nwds2, neps1, wkes, xkes)

!   This routine computes the integral of the function fun on the interval
!   (x1, inf) with a target tolerance of 10^neps1. The quadrature level is
!   progressively increased (approximately doubling the work with each level)
!   until level nq1 has been performed or the target tolerance has been met.
!   nq2 is the size of the wkes and xkes arrays, which must first be 
!   initialized by calling initqes. If quades outputs the message "Terms too
!   large", adjust nq1 and neps2 as necessary in the call to initqes.

!   Both initqes and quades are 100% THREAD SAFE -- all requisite parameters
!   and arrays are passed through subroutine arguments. 

!   For some functions, it is important that the endpoint x1 be
!   computed to high precision (nwds2 words) in the calling program, that
!   these high-precision values be passed to quades (where scaled abscissas
!   are calculated to high precision), and that the function definition
!   itself uses these high-precision scaled abscissas in any initial
!   subtractions or other sensitive operations involving the input argument.
!   Otherwise the accuracy of the quadrature result might only be half as
!   high as it otherwise could be. See the function definitions of fun06,
!   fun07, fun09 and fun10 for examples on how this is done.  Otherwise the
!   function evaluation can and should be performed with low precision
!   (nwds1 words) for faster run times. The two precision levels (nwds1
!   and nwds2) are specified by the user in the calling program.

!   David H Bailey  31 Jan 2020

use mpmodule
implicit none
integer i, ierror, ip(0:100), iz1, iz2, izx, j, k, k1, k2, n, ndebug, &
  nds, neps1, nq1, nq2, nqq1, nwds1, nwds2
parameter (izx = 5, ndebug = 2)
logical log1
real (8) d1, d2, d3, d4, dplog10q
type (mp_real) c10, eps1, eps2, epsilon1, err, fun, h, &
  quades, tsum, s1, s2, s3, t1, t2, t3, tw1, tw2, twi1, twi2, twmx, &
  wkes(-1:nq2), xkes(-1:nq2)
type (mp_real) x1, xx1, xx2
external fun, dplog10q

epsilon1 = mpreal (10.d0, nwds1) ** neps1
tsum = mpreal (0.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
s2 = mpreal (0.d0, nwds1)
h = mpreal (1.d0, nwds1)
c10 = mpreal (10.d0, nwds1)

if (wkes(-1) < dble (nq1)) then
  write (6, 1) nq1
1 format ('quades: quadrature arrays have not been initialized; nq1 =',i6)
  goto 140
endif
nqq1 = dble (wkes(-1))
n = dble (xkes(-1))

do k = 0, nqq1
  ip(k) = 2 ** k
enddo

do k = 1, nq1
  h = 0.5d0 * h
  s3 = s2
  s2 = s1
  k1 = ip(nqq1-k)
  k2 = ip(nqq1-k+1)
  iz1 = 0
  iz2 = 0
  twmx = mpreal (0.d0, nwds1)

!   Evaluate function at level k in x, avoiding unnecessary computation.

  do i = 0, n, k1
    if (mod (i, k2) /= 0 .or. k == 1) then

!   These next few lines, which scale the abscissas, must be performed in
!   high precision (nwds2 words) to ensure full accuracy in the quadrature
!   results, even though the abscissas xkes(i) were computed in low precision.

      xx1 = x1 + mpreal (xkes(i), nwds2)
      xx2 = x1 + 1.d0 / mpreal (xkes(i), nwds2)
      log1 = xx1 > x1
  
!   The remaining computations are performed in low precision (nwds1 words).

      if (iz1 < izx) then
        t1 = fun (xx1, nwds1)
        tw1 = t1 * wkes(i)
        twi1 = abs (tw1)
        if (twi1 < epsilon1) then
          iz1 = iz1 + 1
        else
          iz1 = 0
        endif
      else
        t1 = mpreal (0.d0, nwds1)
        tw1 = mpreal (0.d0, nwds1)
      endif

      if (i > 0 .and. log1 .and. iz2 < izx) then
        t2 = fun (xx2, nwds1)
        tw2 = t2 * wkes(i) / xkes(i)**2
        twi2 = abs (tw2)
        if (twi2 < epsilon1) then
          iz2 = iz2 + 1
        else
          iz2 = 0
        endif
      else
        t2 = mpreal (0.d0, nwds1)
        tw2 = mpreal (0.d0, nwds1)
      endif

      tsum = tsum + tw1 + tw2
      twmx = max (twmx, abs (tw1), abs (tw2))
    endif
  enddo

!   Compute s1 = current integral approximation and err = error estimate.
!   Tsum is the sum of all tw1 and tw2 from the loop above.
!   Twmx is the largest absolute value of tw1 and tw2 from the loop above.
!   Twi1 and twi2 are the final nonzero values of abs(tw1) and abs(tw2).

  s1 =  h * tsum
  eps1 = twmx * epsilon1
  eps2 = abs (max (twi1, twi2) / s1)
  d1 = dplog10q (abs ((s1 - s2) / s1))
  d2 = dplog10q (abs ((s1 - s3) / s1))
  d3 = dplog10q (eps1) - 1.d0
  d4 = dplog10q (eps2) - 1.d0

  if (k <= 2) then
    err = mpreal (1.d0, nwds1)
  elseif (d1 .eq. -999999.d0) then
    err = mpreal (0.d0, nwds1)
  else
    err = c10 ** nint (min (0.d0, max (d1 ** 2 / d2, 2.d0 * d1, d3, d4)))
  endif

!   Output current integral approximation and error estimate, to 60 digits.

  if (ndebug >= 2) then
    write (6, 2) k, nq1, nint (dplog10q (abs (err)))
2   format ('quades: Iteration',i3,' of',i3,'; est error = 10^',i7, &
      '; approx value =')
    call mpwrite (6, 80, 60, s1)
  endif

  if (k >= 3 .and. iz1 == 0 .and. iz2 == 0) then
    write (6, 3)
3   format ('quades: Terms too large -- adjust neps2 in call to initqes.')
    goto 140
  endif

  if (k >= 3 .and. err < eps1) then
    write (6, 4) nint (dplog10q (abs (err)))
4   format ('quades: Estimated error = 10^',i7)
    goto 140
  endif

  if (k >= 3 .and. err < eps2) then
    write (6, 5) nint (dplog10q (abs (err)))
5   format ('quades: Estimated error = 10^',i7/&
    'Adjust nq1 and neps2 in initqes for greater accuracy.')
    goto 140
  endif
enddo

140 continue

quades = s1
return
end

subroutine initqss (nq1, nq2, nwds1, neps2, wkss, xkss)

!   This subroutine initializes the quadrature arrays xkss and wkss for quadss.
!   The argument nq2 is the space allocated for wkss and xkss in the calling
!   program. By default it is set to 12 * 2^nq1. If initqss outputs the message
!   "Table space parameter is too small", adjust nq2. Also, if quadss outputs
!   the message "Terms too large", adjust nq1 and neps2 as necessary in the call
!   to initqss. The argument neps2 controls termination of the loop below, which
!   ends when wkss(k) * 10^(neps2) > 1.

!   The wkss and xkss arrays are computed based on the transformation
!   t = sinh (pi/2 * sinh (x)).  See comments below.

!   Both initqss and quadss are 100% THREAD SAFE -- all requisite parameters
!   and arrays are passed through subroutine arguments. 

!   David H Bailey   31 Jan 2020

use mpmodule
implicit none
integer i, ierror, iprint, j, k, k1, ndebug, neps2, nq1, nq2, nwds1
parameter (iprint = 1024, ndebug = 2)
type (mp_real) eps2, h, p2, t1, t2, t3, t4, u1, u2, &
  wkss(-1:nq2), xkss(-1:nq2)

write (6, 1)
1 format ('initqss: Sinh-sinh quadrature initialization')

eps2 = mpreal (10.d0, nwds1) ** neps2
p2 = 0.5d0 * mppi (nwds1)
h = mpreal (0.5d0 ** nq1, nwds1)
wkss(-1) = mpreal (dble (nq1), nwds1)

do k = 0, nq2
  if (ndebug >= 2 .and. mod (k, iprint) == 0) write (6, *) k, nq2
    t1 = mpreal (dble (k) * h, nwds1)

!   xkss(k) = sinh (u1)
!   wkss(k) = cosh (u1) * u2
!   where u1 = pi/2 * sinh (t1) and u2 = pi/2 * cosh (t1)

  t2 = exp (t1)
  u1 = 0.5d0 * p2 * (t2 - 1.d0 / t2)
  u2 = 0.5d0 * p2 * (t2 + 1.d0 / t2)
  t3 = exp (u1)
  xkss(k) = 0.5d0 * (t3 - 1.d0 / t3)
  wkss(k) = 0.5d0 * (t3 + 1.d0 / t3) * u2

  if (wkss(k) * eps2 > 1.d0) goto 100
enddo

write (6, 2) nq2
2 format ('initqss: Table space parameter is too small; value =',i8)
stop

100 continue

xkss(-1) = mpreal (dble (k), nwds1)
if (ndebug >= 2) then
  write (6, 3) k
3 format ('initqss: Table spaced used =',i8)
endif

return
end

function quadss (fun, nq1, nq2, nwds1, neps1, wkss, xkss)

!   This routine computes the integral of the function fun on the interval
!   (-inf, inf) with a target tolerance of 10^neps1. The quadrature level is
!   progressively increased (approximately doubling the work with each level)
!   until level nq1 has been performed or the target tolerance has been met.
!   nq2 is the size of the wkss and xkss arrays, which must first be 
!   initialized by calling initqss. If quadss outputs the message "Terms too
!   large", adjust nq1 and neps2 as necessary in the call to initqss.

!   Both initqss and quadss are 100% THREAD SAFE -- all requisite parameters
!   and arrays are passed through subroutine arguments. 

!   David H Bailey  31 Jan 2020

use mpmodule
implicit none
integer i, ierror, ip(0:100), iz1, iz2, izx, j, k, k1, k2, n, ndebug, &
  nds, neps1, nq1, nq2, nqq1, nwds1
parameter (izx = 5, ndebug = 2)
real (8) d1, d2, d3, d4, dplog10q
type (mp_real) c10, eps1, eps2, epsilon1, err, fun, h, &
  quadss, tsum, s1, s2, s3, t1, t2, t3, tw1, tw2, twi1, twi2, twmx, &
  wkss(-1:nq2), xkss(-1:nq2)
external fun, dplog10q

epsilon1 = mpreal (10.d0, nwds1) ** neps1
tsum = mpreal (0.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
s2 = mpreal (0.d0, nwds1)
h = mpreal (1.d0, nwds1)
c10 = mpreal (10.d0, nwds1)

if (wkss(-1) < dble (nq1)) then
  write (6, 1) nq1
1 format ('quadss: quadrature arrays have not been initialized; nq1 =',i6)
  goto 140
endif
nqq1 = dble (wkss(-1))
n = dble (xkss(-1))

do k = 0, nqq1
  ip(k) = 2 ** k
enddo

do k = 1, nq1
  h = 0.5d0 * h
  s3 = s2
  s2 = s1
  k1 = ip(nqq1-k)
  k2 = ip(nqq1-k+1)
  iz1 = 0
  iz2 = 0
  twmx = mpreal (0.d0, nwds1)

!   Evaluate function at level k in x, avoiding unnecessary computation.

  do i = 0, n, k1
    if (mod (i, k2) /= 0 .or. k == 1) then
      if (iz1 < izx) then
        t1 = fun (xkss(i), nwds1)
        tw1 = t1 * wkss(i)
        twi1 = abs (tw1)
        if (twi1 < epsilon1) then
          iz1 = iz1 + 1
        else
          iz1 = 0
        endif
      else
        t1 = mpreal (0.d0, nwds1)
        tw1 = mpreal (0.d0, nwds1)
      endif

      if (i > 0 .and. iz2 < izx) then
        t2 = fun (-xkss(i), nwds1)
        tw2 = t2 * wkss(i)
        twi2 = abs (tw2)
        if (twi2 < epsilon1) then
          iz2 = iz2 + 1
        else
          iz2 = 0
        endif
      else
        t2 = mpreal (0.d0, nwds1)
        tw2 = mpreal (0.d0, nwds1)
      endif

      tsum = tsum + tw1 + tw2
      twmx = max (twmx, abs (tw1), abs (tw2))
    endif
  enddo

!   Compute s1 = current integral approximation and err = error estimate.
!   Tsum is the sum of all tw1 and tw2 from the loop above.
!   Twmx is the largest absolute value of tw1 and tw2 from the loop above.
!   Twi1 and twi2 are the final nonzero values of abs(tw1) and abs(tw2).

  s1 =  h * tsum
  eps1 = twmx * epsilon1
  eps2 = abs (max (twi1, twi2) / s1)
  d1 = dplog10q (abs ((s1 - s2) / s1))
  d2 = dplog10q (abs ((s1 - s3) / s1))
  d3 = dplog10q (eps1) - 1.d0
  d4 = dplog10q (eps2) - 1.d0

  if (k <= 2) then
    err = mpreal (1.d0, nwds1)
  elseif (d1 .eq. -999999.d0) then
    err = mpreal (0.d0, nwds1)
  else
    err = c10 ** nint (min (0.d0, max (d1 ** 2 / d2, 2.d0 * d1, d3, d4)))
  endif

!   Output current integral approximation and error estimate, to 60 digits.

  if (ndebug >= 2) then
    write (6, 2) k, nq1, nint (dplog10q (abs (err)))
2   format ('quadss: Iteration',i3,' of',i3,'; est error = 10^',i7, &
      '; approx value =')
    call mpwrite (6, 80, 60, s1)
  endif

  if (k >= 3 .and. iz1 == 0 .and. iz2 == 0) then
    write (6, 3)
3   format ('quadss: Terms too large -- adjust neps2 in call to initqss.')
    goto 140
  endif

  if (k >= 3 .and. err < eps1) then
    write (6, 4) nint (dplog10q (abs (err)))
4   format ('quadss: Estimated error = 10^',i7)
    goto 140
  endif

  if (k >= 3 .and. err < eps2) then
    write (6, 5) nint (dplog10q (abs (err)))
5   format ('quadss: Estimated error = 10^',i7/&
    'Adjust nq1 and neps2 in initqss for greater accuracy.')
    goto 140
  endif
enddo

140 continue

quadss = s1
return
end

subroutine initqts (nq1, nq2, nwds1, neps2, wkts, xkts)

!   This subroutine initializes the quadrature arrays xkts and wkts for quadts.
!   The argument nq2 is the space allocated for wkts and xkts in the calling
!   program.  By default it is set to 12 * 2^nq1.  If initqts outputs the message
!   "Table space parameter is too small", adjust nq2.  Also, if quadts outputs
!   the message "Terms too large", adjust nq1 and neps2 as necessary in the call
!   to initqts. The argument neps2 controls termination of the loop below, which
!   ends when wkts(k) < 10^(neps2).

!   The wkts and xkts arrays are computed based on the transformation
!   t = tanh (pi/2 * sinh (x)).  Note however that xkts contains one minus the
!   conventional abscissas, in order to conserve precision. See comments below.

!   Both initqts and quadts are 100% THREAD SAFE -- all requisite parameters
!   and arrays are passed through subroutine arguments. 

!   David H Bailey   31 Jan 2020

use mpmodule
implicit none
integer i, ierror, iprint, j, k, k1, ndebug, neps2, nq1, nq2, nwds1
parameter (iprint = 1024, ndebug = 2)
type (mp_real) eps2, h, p2, t1, t2, t3, t4, u1, u2, wkts(-1:nq2), xkts(-1:nq2)

write (6, 1)
1 format ('initqts: Tanh-sinh quadrature initialization')

eps2 = mpreal (10.d0, nwds1) ** neps2
p2 = 0.5d0 * mppi (nwds1)
h = mpreal (0.5d0 ** nq1, nwds1)
wkts(-1) = mpreal (dble (nq1), nwds1)

do k = 0, nq2
  if (ndebug >= 2 .and. mod (k, iprint) == 0) write (6, *) k, nq2
  t1 = mpreal (dble (k) * h, nwds1)

!   xkts(k) = 1 - tanh (u1) = 1 /(e^u1 * cosh (u1))
!   wkts(k) = u2 / cosh (u1)^2
!   where u1 = pi/2 * cosh (t1), u2 = pi/2 * sinh (t1)

  t2 = exp (t1)
  u1 = 0.5d0 * p2 * (t2 + 1.d0 / t2)
  u2 = 0.5d0 * p2 * (t2 - 1.d0 / t2)
  t3 = exp (u2)
  t4 = 0.5d0 * (t3 + 1.d0 / t3)
  xkts(k) = 1.d0 / (t3 * t4)
  wkts(k) = u1 / t4 ** 2

  if (wkts(k) < eps2) goto 100
enddo

write (6, 2) nq2
2 format ('initqts: Table space parameter is too small; value =',i8)
stop

100 continue

xkts(-1) = mpreal (dble (k), nwds1)
if (ndebug >= 2) then
  write (6, 3) k
3 format ('initqts: Table spaced used =',i8)
endif

return
end

function quadts (fun, x1, x2, nq1, nq2, nwds1, nwds2, neps1, wkts, xkts)

!   This routine computes the integral of the function fun on the interval
!   (x1, x2) with a target tolerance of 10^neps1. The quadrature level is
!   progressively increased (approximately doubling the work with each level)
!   until level nq1 has been performed or the target tolerance has been met.
!   nq2 is the size of the wkts and xkts arrays, which must first be 
!   initialized by calling initqts. If quadts outputs the message "Terms too
!   large", adjust nq1 and neps2 as necessary in the call to initqts.

!   Both initqts and quadts are 100% THREAD SAFE -- all requisite parameters
!   and arrays are passed through subroutine arguments. 

!   For some functions, it is important that the endpoint x1 be
!   computed to high precision (nwds2 words) in the calling program, that
!   these high-precision values be passed to quadts (where scaled abscissas
!   are calculated to high precision), and that the function definition
!   itself uses these high-precision scaled abscissas in any initial
!   subtractions or other sensitive operations involving the input argument.
!   Otherwise the accuracy of the quadrature result might only be half as
!   high as it otherwise could be. See the function definitions of fun06,
!   fun07, fun09 and fun10 for examples on how this is done.  Otherwise the
!   function evaluation can and should be performed with low precision
!   (nwds1 words) for faster run times. The two precision levels (nwds1
!   and nwds2) are specified by the user in the calling program.

!   David H Bailey  31 Jan 2020

use mpmodule
implicit none
integer i, ierror, ip(0:100), iz1, iz2, izx, j, k, k1, k2, n, ndebug, &
  nds, neps1, nq1, nq2, nqq1, nwds1, nwds2
parameter (izx = 5, ndebug = 2)
logical log1, log2
real (8) d1, d2, d3, d4, dplog10q
type (mp_real) c10, eps1, eps2, epsilon1, err, fun, h, &
  quadts, tsum, s1, s2, s3, t1, t2, t3, tw1, tw2, twi1, twi2, twmx, &
  wkts(-1:nq2), xkts(-1:nq2)
type (mp_real) ax, bx, x1, x2, xki, xt1, xx1, xx2
external fun, dplog10q

!  These two lines are performed in high precision (nwds2 words).

ax = 0.5d0 * (x2 - x1)
bx = 0.5d0 * (x2 + x1)

!  The remaining initialization is performed in low precision (nwds1 words).

epsilon1 = mpreal (10.d0, nwds1) ** neps1
tsum = mpreal (0.d0, nwds1)
s1 = mpreal (0.d0, nwds1)
s2 = mpreal (0.d0, nwds1)
h = mpreal (1.d0, nwds1)
c10 = mpreal (10.d0, nwds1)

if (wkts(-1) < dble (nq1)) then
  write (6, 1) nq1
1 format ('quadts: Quadrature arrays have not been initialized; nq1 =',i6)
  goto 140
endif
nqq1 = dble (wkts(-1))
n = dble (xkts(-1))

do k = 0, nqq1
  ip(k) = 2 ** k
enddo

do k = 1, nq1
  h = 0.5d0 * h
  s3 = s2
  s2 = s1
  k1 = ip(nqq1-k)
  k2 = ip(nqq1-k+1)
  iz1 = 0
  iz2 = 0
  twmx = mpreal (0.d0, nwds1)

!   Evaluate function at level k in x, avoiding unnecessary computation.

  do i = 0, n, k1  
    if (mod (i, k2) /= 0 .or. k == 1) then

!   These next few lines, which scale the abscissas, must be performed in
!   high precision (nwds2 words) to ensure full accuracy in the quadrature
!   results, even though the abscissas xkts(i) were computed in low precision.

      xki = xkts(i)
      xt1 = 1.d0 - mpreal (xki, nwds2)
      xx1 = - ax * xt1 + bx
      xx2 = ax * xt1 + bx
      log1 = xx1 > x1
      log2 = xx2 < x2      

!   The remaining computations are performed in low precision (nwds1 words).

      if (log1 .and. iz1 < izx) then
        t1 = fun (xx1, nwds1, nwds2)
        tw1 = t1 * wkts(i)
        twi1 = abs (tw1)
        if (twi1 < epsilon1) then
          iz1 = iz1 + 1
        else
          iz1 = 0
        endif
      else
        t1 = mpreal (0.d0, nwds1)
        tw1 = mpreal (0.d0, nwds1)
      endif

      if (i > 0 .and. log2 .and. iz2 < izx) then
        t2 = fun (xx2, nwds1, nwds2)
        tw2 = t2 * wkts(i)
        twi2 = abs (tw2)
        if (twi2 < epsilon1) then
          iz2 = iz2 + 1
        else
          iz2 = 0
        endif
      else
        t2 = mpreal (0.d0, nwds1)
        tw2 = mpreal (0.d0, nwds1)
      endif

      tsum = tsum + tw1 + tw2
      twmx = max (twmx, abs (tw1), abs (tw2))
    endif
  enddo

!   Compute s1 = current integral approximation and err = error estimate.
!   Tsum is the sum of all tw1 and tw2 from the loop above.
!   Twmx is the largest absolute value of tw1 and tw2 from the loop above.
!   Twi1 and twi2 are the final nonzero values of abs(tw1) and abs(tw2).

  s1 =  mpreal (ax, nwds1) * h * tsum
  eps1 = twmx * epsilon1
  eps2 = abs (max (twi1, twi2) / s1)
  d1 = dplog10q (abs ((s1 - s2) / s1))
  d2 = dplog10q (abs ((s1 - s3) / s1))
  d3 = dplog10q (eps1) - 1.d0
  d4 = dplog10q (eps2) - 1.d0

  if (k <= 2) then
    err = mpreal (1.d0, nwds1)
  elseif (d1 .eq. -999999.d0) then
    err = mpreal (0.d0, nwds1)
  else
    err = c10 ** nint (min (0.d0, max (d1 ** 2 / d2, 2.d0 * d1, d3, d4)))
  endif

!   Output current integral approximation and error estimate, to 60 digits.

  if (ndebug >= 2) then
    write (6, 2) k, nq1, nint (dplog10q (abs (err)))
2   format ('quadts: Iteration',i3,' of',i3,'; est error = 10^',i7, &
      '; approx value =')
    call mpwrite (6, 80, 60, s1)
  endif

  if (k >= 3 .and. iz1 == 0 .and. iz2 == 0) then
    write (6, 3)
3   format ('quadts: Terms too large -- adjust neps2 in call to initqts.')
    goto 140
  endif

  if (k >= 3 .and. err < eps1) then
    write (6, 4) nint (dplog10q (abs (err)))
4   format ('quadts: Estimated error = 10^',i7)
    goto 140
  endif

  if (k >= 3 .and. err < eps2) then
    write (6, 5) nint (dplog10q (abs (err)))
5   format ('quadts: Estimated error = 10^',i7/&
    'Adjust nq1 and neps2 in initqts for greater accuracy.')
    goto 140
  endif
enddo

140 continue

quadts = s1
return
end

function dplog10q (a)

!   For input MP value a, this routine returns a DP approximation to log10 (a).

use mpmodule
implicit none
integer ia
real (8) da, dplog10q, t1
type (mp_real) a

call mpmdi (a, da, ia)
if (da .eq. 0.d0) then
  dplog10q = -999999.d0
else
  dplog10q = log10 (abs (da)) + ia * log10 (2.d0)
endif

100 continue
return
end

subroutine decmdq (a, b, ib)

!   For input MP value a, this routine returns DP b and integer ib such that 
!   a = b * 10^ib, with 1 <= abs (b) < 10 for nonzero a.

use mpmodule
implicit none
integer ia, ib
real (8) da, b, t1, xlt
parameter (xlt = 0.3010299956639812d0)
type (mp_real) a

call mpmdi (a, da, ia)
if (da .ne. 0.d0) then
  t1 = xlt * ia + log10 (abs (da))
  ib = t1
  if (t1 .lt. 0.d0) ib = ib - 1
  b = sign (10.d0 ** (t1 - ib), da)
else
  b = 0.d0
  ib = 0
endif

return
end
