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