Commit 9dbf60e8 authored by ulrich_y's avatar ulrich_y

Follow prec standard strictly

Test will not compile because they are still inprecise
parent d2070230
......@@ -186,7 +186,6 @@ addflag() {
}
gnuflags() {
eval addflag FFLAGS -fdefault-real-8 # default double is 8 byte
eval addflag FFLAGS -cpp # C pre-processor
eval addflag FFLAGS -pedantic-errors -std=f2008 # folow f08 strictly
eval addflag FFLAGS -J build # where to put mods
......@@ -246,7 +245,6 @@ case "$1,$2,$3" in
fi
;;
ifort*)
eval addflag FFLAGS -autodouble
eval addflag FFLAGS -module build
eval addflag FFLAGS -fpp
eval addflag FFLAGS -stand f03
......@@ -567,7 +565,7 @@ Name: handyG
Description: Numerical evaluation of GPL
Version: 0.1.0
Libs: -L\${libdir} -lhandyg
Cflags: -I\${includedir} -fdefault-real-8
Cflags: -I\${includedir}
EOF
fi
......
......@@ -67,7 +67,7 @@ CONTAINS
case(' ',',',';')
read(line(last:i-1), *) buf(1)
last=i+1
input(nin) = cmplx(buf(1), buf(2))
input(nin) = cmplx(buf(1), buf(2),kind=prec)
nin = nin + 1
buf = 0.
......@@ -77,7 +77,7 @@ CONTAINS
case('I','i','j')
read(line(last:i-1), *) buf(2)
last=i+1
input(nin) = cmplx(buf(1), buf(2))
input(nin) = cmplx(buf(1), buf(2),kind=prec)
nin = nin + 1
buf = 0.
endselect
......
......@@ -4,16 +4,17 @@ MODULE globals
integer, parameter :: prec = selected_real_kind(15,32)
real, parameter :: zero = 1e-15 ! values smaller than this count as zero
real, parameter :: pi = 3.14159265358979323846
real(kind=prec), parameter :: zero = 1e-15_prec ! values smaller than this count as zero
real(kind=prec), parameter :: pi = 3.14159265358979323846_prec
! The following parameters control the accuracy of the evaluation
real(kind=prec), protected :: MPLdel = 1e-15 ! 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 ! when to apply Hoelder convolution?
real(kind=prec), protected :: MPLdel = 1e-15_prec ! 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 /)
! = (/ mmax, n /). At most n polylogs with weight mmax will be cached
complex(kind=prec), parameter :: i_ = (0.,1._prec)
integer :: verb = 0
CONTAINS
......
......@@ -25,7 +25,7 @@ CONTAINS
if (real(y).gt.0) then
GPL_zero_zi = 1.0_prec/factorial(l) * log(real(y)) ** l
else
GPL_zero_zi = 1.0_prec/factorial(l) * (log(-real(y))+cmplx(0,y%i0*pi)) ** l
GPL_zero_zi = 1.0_prec/factorial(l) * (log(-real(y))+i_*(y%i0*pi)) ** l
endif
else
GPL_zero_zi = 1.0_prec/factorial(l) * log(y%c) ** l
......@@ -436,14 +436,14 @@ CONTAINS
! simpler notation for flat evaluation
real(kind=prec) :: g(:)
complex(kind=prec) :: res
res = G_flat(toinum(cmplx(g(1:size(g)-1))), inum(cmplx(g(size(g))),di0))
res = G_flat(toinum(g(1:size(g)-1)), toinum(g(size(g))))
END FUNCTION G_real
FUNCTION G_int(g) result(res)
! simpler notation for flat evaluation
integer:: g(:)
complex(kind=prec) :: res
res = G_flat(toinum(cmplx(g(1:size(g)-1))), inum(cmplx(g(size(g))),di0))
res = G_flat(toinum(real(g(1:size(g)-1),kind=prec)), toinum(real(g(size(g)),kind=prec)))
END FUNCTION G_int
RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res)
......@@ -502,7 +502,7 @@ CONTAINS
FUNCTION G_FLATr(Z_FLAT,Y)
real(kind=prec), intent(in) :: z_flat(:), y
complex(kind=prec) :: g_flatr
g_flatr = G_flat(toinum(cmplx(z_flat)), inum(cmplx(y),di0))
g_flatr = G_flat(toinum(z_flat), toinum(y))
END FUNCTION
FUNCTION G_FLATc(Z_FLAT,Y)
......
......@@ -8,9 +8,9 @@ MODULE maths_functions
module procedure polylog1, polylog2
end interface polylog
real(kind=prec), parameter :: zeta(2:10) = (/1.6449340668482262, 1.2020569031595942, 1.0823232337111381, &
1.03692775514337, 1.0173430619844488, 1.008349277381923, &
1.0040773561979441, 1.0020083928260821, 1.000994575127818/)
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/)
type el
type(inum) :: c
......@@ -26,10 +26,10 @@ CONTAINS
integer :: m
complex(kind=prec) :: x, res
integer :: i
res=0.
res=0._prec
do i=1,PolylogInfinity
if(i**m.lt.0) return ! roll over
if(abs(x**i).lt.1.e-250) return
if(abs(x**i).lt.1.e-250_prec) return
res = res+x**i/i**m
enddo
END FUNCTION naive_polylog
......@@ -131,7 +131,7 @@ CONTAINS
complex(kind=prec) :: x
if(abs(aimag(x)) < zero ) then
res = Li2(real(x))
res = Li2(real(x,kind=prec))
else
res = naive_polylog(2,x)
endif
......@@ -145,47 +145,45 @@ CONTAINS
real (kind=prec):: X,S,A
real (kind=prec):: CA(0:18),HA,ALFAA,BA0,BA1,BA2, YA
real (kind=prec):: CB(0:18),HB,ALFAB,BB0,BB1,BB2, YB
DATA CA(0) / 0.4617293928601208/
DATA CA(1) / 0.4501739958855029/
DATA CA(2) / -0.010912841952292843/
DATA CA(3) / 0.0005932454712725702/
DATA CA(4) / -0.00004479593219266303/
DATA CA(5) / 4.051545785869334e-6/
DATA CA(6) / -4.1095398602619446e-7/
DATA CA(7) / 4.513178777974119e-8/
DATA CA(8) / -5.254661564861129e-9/
DATA CA(9) / 6.398255691618666e-10/
DATA CA(10) / -8.071938105510391e-11/
DATA CA(11) / 1.0480864927082917e-11/
DATA CA(12) / -1.3936328400075057e-12/
DATA CA(13) / 1.8919788723690422e-13/
DATA CA(14) / -2.6097139622039465e-14/
DATA CA(15) / 3.774985548158685e-15/
DATA CA(16) / -5.671361978114946e-16/
DATA CA(17) / 1.1023848202712794e-16/
DATA CA(18) / -5.0940525990875006e-17/
DATA CB(0) / -0.016016180449195803/
DATA CB(1) / -0.5036424400753012/
DATA CB(2) / -0.016150992430500253/
DATA CB(3) / -0.0012440242104245127/
DATA CB(4) / -0.00013757218124463538/
DATA CB(5) / -0.000018563818526041144/
DATA CB(6) / -2.841735345177361e-6/
DATA CB(7) / -4.7459967908588557e-7/
DATA CB(8) / -8.448038544563037e-8/
DATA CB(9) / -1.5787671270014e-8/
DATA CB(10) / -3.0657620579122164e-9/
DATA CB(11) / -6.140791949281482e-10/
DATA CB(12) / -1.2618831590198e-10/
DATA CB(13) / -2.64931268635803e-11/
DATA CB(14) / -5.664711482422879e-12/
DATA CB(15) / -1.2303909436235178e-12/
DATA CB(16) / -2.7089360852246495e-13/
DATA CB(17) / -6.024075373994343e-14/
DATA CB(18) / -1.2894320641440237e-14/
DATA CA(0) / 0.4617293928601208_prec/
DATA CA(1) / 0.4501739958855029_prec/
DATA CA(2) / -0.010912841952292843_prec/
DATA CA(3) / 0.0005932454712725702_prec/
DATA CA(4) / -0.00004479593219266303_prec/
DATA CA(5) / 4.051545785869334e-6_prec/
DATA CA(6) / -4.1095398602619446e-7_prec/
DATA CA(7) / 4.513178777974119e-8_prec/
DATA CA(8) / -5.254661564861129e-9_prec/
DATA CA(9) / 6.398255691618666e-10_prec/
DATA CA(10) / -8.071938105510391e-11_prec/
DATA CA(11) / 1.0480864927082917e-11_prec/
DATA CA(12) / -1.3936328400075057e-12_prec/
DATA CA(13) / 1.8919788723690422e-13_prec/
DATA CA(14) / -2.6097139622039465e-14_prec/
DATA CA(15) / 3.774985548158685e-15_prec/
DATA CA(16) / -5.671361978114946e-16_prec/
DATA CA(17) / 1.1023848202712794e-16_prec/
DATA CA(18) / -5.0940525990875006e-17_prec/
DATA CB(0) / -0.016016180449195803_prec/
DATA CB(1) / -0.5036424400753012_prec/
DATA CB(2) / -0.016150992430500253_prec/
DATA CB(3) / -0.0012440242104245127_prec/
DATA CB(4) / -0.00013757218124463538_prec/
DATA CB(5) / -0.000018563818526041144_prec/
DATA CB(6) / -2.841735345177361e-6_prec/
DATA CB(7) / -4.7459967908588557e-7_prec/
DATA CB(8) / -8.448038544563037e-8_prec/
DATA CB(9) / -1.5787671270014e-8_prec/
DATA CB(10) / -3.0657620579122164e-9_prec/
DATA CB(11) / -6.140791949281482e-10_prec/
DATA CB(12) / -1.2618831590198e-10_prec/
DATA CB(13) / -2.64931268635803e-11_prec/
DATA CB(14) / -5.664711482422879e-12_prec/
DATA CB(15) / -1.2303909436235178e-12_prec/
DATA CB(16) / -2.7089360852246495e-13_prec/
DATA CB(17) / -6.024075373994343e-14_prec/
DATA CB(18) / -1.2894320641440237e-14_prec/
real (kind=prec):: Li3
real (kind=prec), parameter :: zeta2 = 1.6449340668482264365
real (kind=prec), parameter :: zeta3 = 1.2020569031595942854
integer :: i
......@@ -196,16 +194,16 @@ CONTAINS
endif
IF(X > 0.999999_prec) THEN
LI3=zeta3
LI3=zeta(3)
RETURN
ELSE IF( abs(x+1) < zero) THEN
LI3=-0.75_prec*zeta3
LI3=-0.75_prec*zeta(3)
RETURN
END IF
IF(X .LE. -1._prec) THEN
YA=1._prec/x ; YB=0._prec
S=-1._prec
A=-LOG(-X)*(zeta2+LOG(-x)**2/6._prec)
A=-LOG(-X)*(zeta(2)+LOG(-x)**2/6._prec)
ELSE IF(X .LE. 0._prec) THEN
YA=x ; YB=0._prec
S=-1._prec
......@@ -217,15 +215,15 @@ CONTAINS
ELSE IF(X .LE. 1._prec) THEN
YA=(x-1._prec)/x ; YB=1._prec-x
S=1._prec
A=zeta3 + zeta2*Log(x) - (Log(1._prec - X)*Log(X)**2)/2._prec + Log(X)**3/6._prec
A=zeta(3) + zeta(2)*Log(x) - (Log(1._prec - X)*Log(X)**2)/2._prec + Log(X)**3/6._prec
ELSE IF(X .LE. 2._prec) THEN
YA=1._prec - X ; YB=(X-1._prec)/X
S=1._prec
A=zeta3 + zeta2*Log(x) - (Log(X - 1._prec)*Log(X)**2)/2._prec + Log(X)**3/6._prec
A=zeta(3) + zeta(2)*Log(x) - (Log(X - 1._prec)*Log(X)**2)/2._prec + Log(X)**3/6._prec
ELSE
YA=0._prec ; YB=1._prec/X
S=-1._prec
A=2*zeta2*Log(x)-Log(x)**3/6._prec
A=2*zeta(2)*Log(x)-Log(x)**3/6._prec
END IF
......@@ -246,7 +244,7 @@ CONTAINS
complex(kind=prec) :: res
complex(kind=prec) :: x
if(abs(aimag(x)) < zero ) then
res = Li3(real(x))
res = Li3(real(x,kind=prec))
else
res = naive_polylog(3,x)
endif
......@@ -329,12 +327,12 @@ CONTAINS
#endif
if ((m.le.9).and.(abs(x%c-1.).lt.zero)) then
res = zeta(m)
else if ((m.le.9).and.(abs(x%c+1.).lt.zero)) then
res = -(1. - 2.**(1-m))*zeta(m)
else if ((m.le.9).and.(abs(x%c+1._prec).lt.zero)) then
res = -(1._prec - 2._prec**(1-m))*zeta(m)
else if (abs(x) .gt. 1) then
inv = inum(1./x%c, x%i0)
inv = inum(1._prec/x%c, x%i0)
res = (-1)**(m-1)*polylog(m,inv) &
- cmplx(0,2*pi)**m * bernoulli_polynomial(m, 0.5-cmplx(0.,1.)*conjg(log(-x%c))/2/pi) / factorial(m)
- (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*conjg(log(-x%c))/2._prec/pi) / factorial(m)
else if(m == 2) then
res = dilog(x%c)
else if(m == 3) then
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment