Prepared for quad

src/globals.f90

@@ -4,11 +4,11 @@ MODULE globals
 
 
   integer, parameter :: prec = selected_real_kind(15,32)  
-  real(kind=prec), parameter :: zero = 1e-15_prec     ! values smaller than this count as zero
-  real(kind=prec), parameter :: pi = 3.14159265358979323846_prec
+  real(kind=prec), parameter :: zero = 10._prec**(-precision(1._prec))     ! values smaller than this count as zero
+  real(kind=prec), parameter :: pi = 3.1415926535897932384626433832795028841971693993751_prec
 
   ! The following parameters control the accuracy of the evaluation
-  real(kind=prec), protected :: MPLdel = 1e-15_prec      ! if the MPL sum changes less then del it is truncated.
+  real(kind=prec), protected :: MPLdel = zero            ! if the MPL sum changes less then del it is truncated.
   integer, protected :: PolylogInfinity = 1000           ! expansion order for Polylogs
   real(kind=prec), protected :: HoelderCircle = 1.1_prec ! when to apply Hoelder convolution?
   integer, parameter :: PolyLogCacheSize(2) = (/ 5, 100 /)

src/maths_functions.f90

@@ -8,10 +8,11 @@ MODULE maths_functions
     module procedure polylog1, polylog2
   end interface polylog
 
-  real(kind=prec), parameter :: zeta(2:10) = (/1.6449340668482262_prec, 1.2020569031595942_prec, 1.0823232337111381_prec, &
-               1.03692775514337_prec, 1.0173430619844488_prec, 1.008349277381923_prec, &
-               1.0040773561979441_prec, 1.0020083928260821_prec, 1.000994575127818_prec/)
-
+  real(kind=prec), parameter :: zeta(2:10) = (/ 1.6449340668482264364724151666460251892189499012068_prec,  &
+      1.2020569031595942853997381615114499907649862923405_prec, 1.0823232337111381915160036965411679027747509519187_prec,  &
+      1.0369277551433699263313654864570341680570809195019_prec, 1.0173430619844491397145179297909205279018174900329_prec,  &
+      1.0083492773819228268397975498497967595998635605652_prec, 1.0040773561979443393786852385086524652589607906499_prec,  &
+      1.0020083928260822144178527692324120604856058513949_prec, 1.0009945751278180853371459589003190170060195315645_prec /)
   type el
     type(inum) :: c
     complex(kind=prec) ans
@@ -38,7 +39,7 @@ CONTAINS
 
    !! Dilogarithm for arguments x < = 1.0
 
-   real (kind=prec):: X,Y,T,S,A,PI3,PI6,ZERO,ONE,HALF,MALF,MONE,MTWO
+   real (kind=prec):: X,Y,T,S,A,ZERO,ONE,HALF,MALF,MONE,MTWO
    real (kind=prec):: C(0:18),H,ALFA,B0,B1,B2,LI2_OLD
    real (kind=prec):: Li2
    integer :: i
@@ -46,7 +47,6 @@ CONTAINS
    DATA ZERO /0.0_prec/, ONE /1.0_prec/
    DATA HALF /0.5_prec/, MALF /-0.5_prec/ 
    DATA MONE /-1.0_prec/, MTWO /-2.0_prec/
-   DATA PI3 /3.289868133696453_prec/, PI6 /1.644934066848226_prec/
 
    DATA C( 0) / 0.4299669356081370_prec/
    DATA C( 1) / 0.4097598753307711_prec/
@@ -75,28 +75,28 @@ CONTAINS
    endif    
 
    IF(X > 0.999999_prec) THEN
-    LI2_OLD=PI6
+    LI2_OLD=zeta(2)
     Li2 = Real(LI2_OLD,prec)
     RETURN
    ELSE IF(abs(x-MONE) < zero) THEN
-    LI2_OLD=MALF*PI6
+    LI2_OLD=MALF*zeta(2)
     RETURN
    END IF
    T=-X
    IF(T .LE. MTWO) THEN
     Y=MONE/(ONE+T)
     S=ONE
-    A=-PI3+HALF*(LOG(-T)**2-LOG(ONE+ONE/T)**2)
+    A=-2*zeta(2)+HALF*(LOG(-T)**2-LOG(ONE+ONE/T)**2)
    ELSE IF(T .LT. MONE) THEN
     Y=MONE-T
     S=MONE
     A=LOG(-T)
-    A=-PI6+A*(A+LOG(ONE+ONE/T))
+    A=-zeta(2)+A*(A+LOG(ONE+ONE/T))
    ELSE IF(T .LE. MALF) THEN
     Y=(MONE-T)/T
     S=ONE
     A=LOG(-T)
-    A=-PI6+A*(MALF*A+LOG(ONE+T))
+    A=-zeta(2)+A*(MALF*A+LOG(ONE+T))
    ELSE IF(T .LT. ZERO) THEN
     Y=-T/(ONE+T)
     S=MONE
@@ -108,7 +108,7 @@ CONTAINS
    ELSE
     Y=ONE/T
     S=MONE
-    A=PI6+HALF*LOG(T)**2
+    A=zeta(2)+HALF*LOG(T)**2
    END IF
 
    H=Y+Y-ONE