...
 
Commits (43)
......@@ -14,6 +14,11 @@ cache:
paths:
- makefile
- build/
- checks/test-chen.f90
- checks/test-chenff.f90
- checks/test-muone.f90
- checks/test-muoneNP.f90
key: ${CI_COMMIT_REF_SLUG}
stages:
- configure
......
......@@ -78,9 +78,17 @@ The compilation process creates the following results
difference between two successive terms at which the series
expansion (6) is truncated.
* `integer LiInf = 1000`
* `real(kind=prec) :: Lidel = 1e-15`
difference $`|\delta_i|`$ between two successive terms at which the
series expansion for
```math
{\rm Li}_n(x) = \sum_{i=1}^\infty \frac{x^i}{i^n}
= \sum_{i=1}^\infty \delta_i
```
is truncated.
number of terms in the expansion of classical polylogarithms.
* `real(kind=prec) :: hCircle = 1.1`
......
......@@ -2,24 +2,13 @@ module chenreftest
contains
function test(z,ref,test_id)
use globals, only: prec
use gpl_module
use ttools
implicit none
complex(kind=prec) :: z(:), res, ref
complex(kind=prec) :: z(:), ref
character(len=*) :: test_id
real(kind=prec) :: delta
logical test
call test_one_flat(z,ref,test_id, test)
print*, ' ', 'testing GPL ', test_id, ' ...'
res = G(z)
delta = abs(res-ref)
if(delta < 8e-7) then
print*, ' ',' passed with delta = ', delta
test = .true.
else
print*, ' ',' FAILED with delta = ', delta
test = .false.
endif
end function
function do_chen_test() result(success)
use globals, only: prec
......
This diff is collapsed.
MODULE TTOOLS
use globals
use mpl_module
use gpl_module
real(kind=prec) :: tol = 8.0e-7
logical :: tests_successful = .true.
contains
subroutine iprint(imsg, typ)
character(len=*) imsg
character(len=200) msg
integer :: typ
character(len=5),parameter :: red = char(27)//'[31m'
character(len=5),parameter :: green = char(27)//'[32m'
character(len=5),parameter :: orange= char(27)//'[33m'
character(len=4),parameter :: norm = char(27)//'[0m'
character ,parameter :: cr = achar(13)
integer, save :: prevlen, prevtype
character(len=200), save :: prevmsg
if ( (prevtype == -1).and.(typ .ne. -1) ) then
msg = prevmsg(1:prevlen) // trim(imsg)
else
msg = trim(imsg)
endif
select case(typ)
case(0)
print*,green //'[PASS]'//norm//' '//trim(msg)
case(1)
print*, red //'[FAIL]'//norm//' '//trim(msg)
case(2)
print*, red //'[FATL]'//norm//' '//trim(msg)
stop 1
case(3)
print*,orange//'[WARN]'//norm//' '//trim(msg)
case(4)
print*,orange//'[INFO]'//norm//' '//trim(msg)
case(-1)
write(*,'(a)',advance='no')' [ ]'//' '//trim(msg)//cr
end select
prevtype = typ
prevlen = len_trim(msg)
prevmsg = msg
end subroutine
function readint(prev,i)
integer i, readint, st
character(len=32) :: arg
character(len=*) :: prev
i=i+1
call get_command_argument(i,arg)
if (len_trim(arg) == 0) call iprint("Argument "//prev//" requires a number",2)
read(arg,*,iostat=st) readint
if (st .ne. 0) call iprint("For argument "//prev//": "//trim(arg)//" is not a number",2)
end function
subroutine check(res, ref, ans, ttol)
complex(kind=prec) :: res, ref
real(kind=prec) :: delta, mytol
real(kind=prec), optional :: ttol
character(len=40) :: msg
logical, optional :: ans
mytol = tol
if (present(ttol)) mytol = ttol
delta = abs(res-ref)
if(delta < mytol) then
write(msg, 900) delta
call iprint(trim(msg), 0)
if (present(ans)) ans = .true.
else
write(msg, 900) delta
call iprint(trim(msg), 1)
tests_successful = .false.
if (present(ans)) ans = .false.
end if
900 format(" with delta = ",ES10.3)
end subroutine check
subroutine test_one_MPL(m,x,ref, test_id)
integer :: m(:)
complex(kind=prec) :: x(:), ref, res
character(len=*) :: test_id
call iprint(' testing MPL '//test_id//' ...',-1)
res = MPL(m,x)
call check(res,ref)
end subroutine test_one_MPL
subroutine test_one_condensed(m,z,y,k,ref,test_id)
integer :: m(:), k
complex(kind=prec) :: z(:), y, res, ref
character(len=*) :: test_id
call iprint(' testing GPL '//test_id//' ...',-1)
res = G_condensed(m,toinum(z),inum(y,di0),k)
call check(res,ref)
end subroutine test_one_condensed
subroutine test_one_flat(z,ref,test_id, ans)
complex(kind=prec) :: z(:), res, ref
character(len=*) :: test_id
logical, optional :: ans
call iprint(' testing GPL '//test_id//' ...',-1)
res = G(z)
if (present(ans)) then
call check(res,ref,ans)
else
call check(res,ref)
endif
end subroutine test_one_flat
END MODULE TTOOLS
This diff is collapsed.
......@@ -5,21 +5,33 @@ using namespace GiNaC;
#include <iostream>
typedef struct {double r,i;} complex_t;
typedef struct {complex_t c; signed char i0;} inum_t;
extern "C"{
complex_t geval_(complex_t * z, int* n);
complex_t geval_(inum_t * z, int* n);
};
complex_t geval_(complex_t * z, int* n) {
complex_t geval_(inum_t * z, int* n) {
cln::cl_inhibit_floating_point_underflow = true;
lst w;
lst w,s;
for(long i=0;i<(*n)-1;i++)
{
w.append((z->r)+(z->i)*I);
ex zz;
if (abs(z->c.i) < 1e-15)
w.append((z->c.r));
else
w.append((z->c.r)+(z->c.i)*I);
s.append(z->i0);
z++;
}
ex ans = G(w,z->r).evalf();
return {
.r = ex_to<numeric>(evalf(real_part(ans))).to_double(),
.i = ex_to<numeric>(evalf(imag_part(ans))).to_double()
};
try {
ex ans = G(w,s,z->c.r).evalf();
return {
.r = ex_to<numeric>(evalf(real_part(ans))).to_double(),
.i = ex_to<numeric>(evalf(imag_part(ans))).to_double()
};
} catch (...) {
std::cout << "Caught!!\n";
return {.r = 1.e11, .i = 1.e11};
}
}
......@@ -15,8 +15,8 @@ MODULE globals
real(kind=prec), parameter :: pi = 3.1415926535897932384626433832795028841971693993751_prec
! The following parameters control the accuracy of the evaluation
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 :: MPLdelta = zero ! if the MPL sum changes less then del it is truncated.
real(kind=prec), protected :: Lidelta = zero ! like MPLdelta but 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
......@@ -46,11 +46,10 @@ CONTAINS
END SUBROUTINE parse_cmd_args
#endif
SUBROUTINE SET_OPTIONS(mpldel, liinf, hcircle)
real(kind=prec), optional :: hcircle, mpldel
integer, optional :: liinf
if (present(mpldel)) MPLdel = mpldel
if (present(liinf)) PolyLogInfinity = liinf
SUBROUTINE SET_OPTIONS(mpldel, lidel, hcircle)
real(kind=prec), optional :: hcircle, mpldel, lidel
if (present(MPLdel)) MPLdelta = mpldel
if (present( Lidel)) LiDelta = lidel
if (present(hcircle)) HoelderCircle = hcircle
END SUBROUTINE
......
......@@ -52,9 +52,10 @@ CONTAINS
if(.not. present(y)) print*, 'G(', abs(z_flat), ')'
END SUBROUTINE print_G
RECURSIVE FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p) result(res)
RECURSIVE FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p, srs) result(res)
! here what is passed is not the full a vector, only a1, ..., ak without the trailing zeroes
integer :: m, i, j, n
integer(1) :: srs
type(inum) :: a(:), y2, s(m), p(:)
complex(kind=prec) :: res
type(inum) :: alpha(product((/(i,i=1,size(a)+size(s))/))/ &
......@@ -70,7 +71,7 @@ CONTAINS
print*, 'PI with p=',real(p),'i=',m,'g =',real([zeroes(m-1),y2])
end if
#endif
res = G_flat(a,y2)*pending_integral(p,m,[zeroes(m-1),y2])
res = G_flat(a,y2)*pending_integral(p,m,[zeroes(m-1),y2], srs)
#ifdef DEBUG
if(verb >= 50) print*, 'also mapping to'
#endif
......@@ -81,16 +82,17 @@ CONTAINS
if(verb >= 50) print*, 'PI with p=',real(p),'i=',n,'g =',&
real([alpha(j,1:n-1),alpha(j,n+1:size(alpha,2)),y2])
#endif
res = res - pending_integral(p, n, [alpha(j,1:n-1),alpha(j,n+1:size(alpha,2)),y2])
res = res - pending_integral(p, n, [alpha(j,1:n-1),alpha(j,n+1:size(alpha,2)),y2], srs)
end do
END FUNCTION remove_sr_from_last_place_in_PI
RECURSIVE FUNCTION pending_integral(p,i,g) result(res)
RECURSIVE FUNCTION pending_integral(p,i,g,srs) result(res)
! evaluates a pending integral by reducing it to simpler ones and g functions
complex(kind=prec) :: res
type(inum) :: p(:), g(:)
type(inum) :: y1, y2, b(size(p)-1), a(size(g)-1)
integer :: i, m
integer(1):: srs
res = 0
#ifdef DEBUG
......@@ -126,11 +128,8 @@ CONTAINS
#endif
!res = pending_integral(p,2,[sub_ieps(g(1))]) - pending_integral(p,2,[cmplx(0.0)]) &
! + G_flat(p(2:size(p)), p(1)) * log(-sub_ieps(g(1)))
if(abs(aimag(p(1))).gt.zero) then
p(1)%i0 = int(sign(1._prec, aimag(p(1))),1)
endif
res = pending_integral(p,2,[inum( g(1)%c,-g(1)%i0 ) ]) - pending_integral(p,2,[izero]) &
+ G_flat(p(2:size(p)), p(1)) * (log(g(1)%c) + p(1)%i0 * pi * i_)
res = pending_integral(p,2,[inum( g(1)%c,-srs ) ], srs) - pending_integral(p,2,[izero], srs) &
+ G_flat(p(2:size(p)), p(1)) * (log(g(1)%c) + srs * pi * i_)
return
end if
......@@ -151,9 +150,9 @@ CONTAINS
print*, 'PI with p=',tocmplx([p,izero]),'i=',m-1,'g =',tocmplx(zeroes(0))
end if
#endif
res = -zeta(m)*pending_integral(p,0,zeroes(0)) &
+ pending_integral([y2,izero],m-1,[zeroes(m-2),y2])*pending_integral(p,0,zeroes(0)) &
- pending_integral([p,izero] ,m-1,[zeroes(m-2),y2])
res = -zeta(m)*pending_integral(p,0,zeroes(0),srs) &
+ pending_integral([y2,izero],m-1,[zeroes(m-2),y2],srs)*pending_integral(p,0,zeroes(0),srs) &
- pending_integral([p,izero] ,m-1,[zeroes(m-2),y2],srs)
return
end if
......@@ -174,10 +173,10 @@ CONTAINS
end if
#endif
res = pending_integral(p,0,zeroes(0)) * G_flat([izero,a],y2) &
+ pending_integral([p,y2],0,zeroes(0)) * G_flat(a,y2) &
+ pending_integral([p,a(1)],1,[a(2:size(a)),y2]) &
- pending_integral([p,a(1)],0,zeroes(0)) * G_flat(a,y2)
res = pending_integral(p,0,zeroes(0),srs) * G_flat([izero,a],y2) &
+ pending_integral([p,y2],0,zeroes(0),srs) * G_flat(a,y2) &
+ pending_integral([p,a(1)],1,[a(2:size(a)),y2],srs) &
- pending_integral([p,a(1)],0,zeroes(0),srs) * G_flat(a,y2)
return
end if
......@@ -187,7 +186,7 @@ CONTAINS
if(verb >= 30) print*, 's_r at the end under PI, need to shuffle'
#endif
m = find_amount_trailing_zeros(a) + 1
res = remove_sr_from_last_place_in_PI(a(1:size(a)-(m-1)), y2, m, p)
res = remove_sr_from_last_place_in_PI(a(1:size(a)-(m-1)), y2, m, p, srs)
return
end if
......@@ -196,11 +195,11 @@ CONTAINS
if(verb >= 30) print*, 's_r in the middle under PI'
#endif
res = +pending_integral(p,1,zeroes(0)) * G_flat([a(1:i-1),izero,a(i:size(a))],y2) &
- pending_integral([p,a(i-1)],i-1,[a(1:i-2),a(i:size(a)),y2]) &
+ pending_integral([p,a(i-1)],1,zeroes(0)) * G_flat(a,y2) &
+ pending_integral([p,a(i)], i, [a(1:i-1), a(i+1:size(a)),y2]) &
- pending_integral([p,a(i)],1,zeroes(0)) * G_flat(a,y2)
res = +pending_integral(p,1,zeroes(0),srs) * G_flat([a(1:i-1),izero,a(i:size(a))],y2) &
- pending_integral([p,a(i-1)],i-1,[a(1:i-2),a(i:size(a)),y2],srs) &
+ pending_integral([p,a(i-1)],1,zeroes(0),srs) * G_flat(a,y2) &
+ pending_integral([p,a(i)], i, [a(1:i-1), a(i+1:size(a)),y2],srs) &
- pending_integral([p,a(i)],1,zeroes(0),srs) * G_flat(a,y2)
END FUNCTION pending_integral
RECURSIVE FUNCTION remove_sr_from_last_place_in_G(a,y2,m,sr) result(res)
......@@ -254,10 +253,9 @@ CONTAINS
print*, '--------------------------------------------------'
end if
#endif
res = G_flat([izero, a(i+1:size(a))], y2) &
+ G_flat([y2], sr) * G_flat(a(i+1:size(a)), y2) &
+ pending_integral([sr, a(i+1)], i, [a(i+2:size(a)), y2]) &
+ pending_integral([sr, a(i+1)], i, [a(i+2:size(a)), y2],sr%i0) &
- G_flat([a(i+1)],sr) * G_flat(a(i+1:size(a)), y2)
return
end if
......@@ -279,9 +277,9 @@ CONTAINS
#endif
res = G_flat([a(1:i-1),izero,a(i+1:size(a))],y2) &
- pending_integral([sr,a(i-1)], i-1, [a(1:i-2),a(i+1:size(a)),y2]) &
- pending_integral([sr,a(i-1)], i-1, [a(1:i-2),a(i+1:size(a)),y2],sr%i0) &
+ G_flat([a(i-1)],sr) * G_flat([a(1:i-1),a(i+1:size(a))],y2) &
+ pending_integral([sr,a(i+1)], i, [a(1:i-1),a(i+2:size(a)),y2]) &
+ pending_integral([sr,a(i+1)], i, [a(1:i-1),a(i+2:size(a)),y2],sr%i0) &
- G_flat([a(i+1)],sr) * G_flat([a(1:i-1),a(i+1:size(a))],y2)
END FUNCTION make_convergent
......@@ -398,6 +396,9 @@ CONTAINS
! ieps, which is what GiNaC does (l. 1013)
do j=1,size(z_flat)
znorm(j) = inum(z_flat(j)%c/y%c, z_flat(j)%i0)
if (abs(aimag(znorm(j)))>zero) then
znorm(j)%i0 = int(sign(1._prec, aimag(znorm(j))),1)
endif
enddo
res = G_flat(znorm,inum((1.,0.), y%i0))
return
......@@ -418,6 +419,9 @@ CONTAINS
! ieps, which is what GiNaC does (l. 1013)
do j=1,size(z_flat)
znorm(j) = inum(z_flat(j)%c/y%c, z_flat(j)%i0)
if (abs(aimag(znorm(j)))>zero) then
znorm(j)%i0 = int(sign(1._prec, aimag(znorm(j))),1)
endif
enddo
res = improve_convergence(znorm)
return
......@@ -443,7 +447,7 @@ CONTAINS
FUNCTION G_superflat(g) result(res)
! simpler notation for flat evaluation
complex(kind=prec) :: g(:), res
res = G_flat(toinum(g(1:size(g)-1)), inum(g(size(g)),di0))
res = G_flat(toinum(g(1:size(g)-1)), toinum(g(size(g))))
END FUNCTION G_superflat
FUNCTION G_real(g) result(res)
......@@ -522,7 +526,7 @@ CONTAINS
FUNCTION G_FLATc(Z_FLAT,Y)
complex(kind=prec), intent(in) :: z_flat(:), y
complex(kind=prec) :: g_flatc
g_flatc = G_flat(toinum(z_flat), inum(y,di0))
g_flatc = G_flat(toinum(z_flat), toinum(y))
END FUNCTION
......
......@@ -18,7 +18,7 @@ MODULE ieps
end interface abs
interface toinum
module procedure toinum_cmplx, toinum_real, toinum_reals, toinum_int
module procedure toinum_cmplxs, toinum_cmplx, toinum_real, toinum_reals, toinum_int
end interface toinum
interface tocmplx
module procedure tocmplxv, tocmplxs
......@@ -29,6 +29,20 @@ MODULE ieps
interface aimag
module procedure imags, imagv
end interface aimag
#if KINDREAL==16
#ifdef HAVE_GINAC
type inumD
complex(kind=8) :: c
integer(1) :: i0
end type inumD
interface inum2inum
module procedure inum2inumS, inum2inumV
end interface inum2inum
#endif
#endif
CONTAINS
FUNCTION ABSINUM(n1)
implicit none
......@@ -44,6 +58,21 @@ CONTAINS
absinumv = abs(n1%c)
END FUNCTION ABSINUMV
FUNCTION TOINUM_cmplxs(z, s)
complex(kind=prec) :: z
type(inum) :: toinum_cmplxs
integer(1),optional :: s
integer(1) :: ss
if (present(s)) then
ss = s
else
ss = di0
endif
toinum_cmplxs = inum(z, ss)
if (abs(aimag(z))>zero) then
toinum_cmplxs%i0 = int(sign(1._prec, aimag(z)),1)
endif
END FUNCTION TOINUM_cmplxs
FUNCTION TOINUM_cmplx(z, s)
complex(kind=prec) :: z(:)
type(inum) :: toinum_cmplx(size(z))
......@@ -57,6 +86,9 @@ CONTAINS
endif
do i=1,size(z)
toinum_cmplx(i) = inum(z(i), ss)
if (abs(aimag(z(i)))>zero) then
toinum_cmplx(i)%i0 = int(sign(1._prec, aimag(z(i))),1)
endif
enddo
END FUNCTION TOINUM_cmplx
......@@ -141,4 +173,23 @@ CONTAINS
imags = aimag(z%c)
END FUNCTION
#if KINDREAL==16
#ifdef HAVE_GINAC
FUNCTION INUM2INUMS(i)
type(inum ) :: i
type(inumD) :: inum2inums
inum2inums = inumD( cmplx(i%c, kind=8), i%i0)
END FUNCTION INUM2INUMS
FUNCTION INUM2INUMV(i)
type(inum ) :: i(:)
type(inumD) :: inum2inumv(size(i))
integer j
do j=1,size(i)
inum2inumv(j) = inumD( cmplx(i(j)%c, kind=8), i(j)%i0)
enddo
END FUNCTION INUM2INUMV
#endif
#endif
END MODULE IEPS
......@@ -11,7 +11,11 @@
:Evaluate:
args2r[a_]:=Re[N[a/.SubPlus|SubMinus->Identity]];
args2i[a_]:=Im[N[a/.SubPlus|SubMinus->Identity]];
args2e[a_]:=Switch[Head[#], SubPlus, 1, SubMinus, -1, _, 1]& /@ a;
args2e[a_]:=Switch[Head[#],
SubPlus, 1,
SubMinus, -1,
Complex, Sign[Im[#]],
_, 1]& /@ a;
:Begin:
:Function: gpl
:Pattern: G[a__]/;And @@ (NumberQ /@ ({a} /. SubPlus | SubMinus -> Identity))
......
......@@ -13,6 +13,12 @@ MODULE maths_functions
1.0369277551433699263313654864570341680570809195019_prec, 1.0173430619844491397145179297909205279018174900329_prec, &
1.0083492773819228268397975498497967595998635605652_prec, 1.0040773561979443393786852385086524652589607906499_prec, &
1.0020083928260822144178527692324120604856058513949_prec, 1.0009945751278180853371459589003190170060195315645_prec /)
real(kind=prec), parameter :: DirichletBeta(2:10) = (/ 0.91596559417721901505460351493238411077414937428167_prec, &
0.96894614625936938048363484584691860006954026768391_prec, 0.98894455174110533610842263322837782131586088706273_prec, &
0.99615782807708806400631936863097528151139552938826_prec, 0.99868522221843813544160078786020654967836454612651_prec, &
0.99955450789053990949634654989905898300218848194998_prec, 0.99984999024682965633806705924046378147600743300743_prec, &
0.99994968418722008982135887329384752737274799691796_prec, 0.99998316402619687740554072995833414145685781649717_prec /)
type el
complex(kind=prec) :: c
complex(kind=prec) ans
......@@ -25,14 +31,19 @@ CONTAINS
FUNCTION naive_polylog(m,x) result(res)
! Computes the classical polylogarithm Li_m(x) using series representation up to order n
integer :: m
complex(kind=prec) :: x, res
complex(kind=prec) :: x, res, del
integer(kind=ikin) :: i
res=0._prec
do i=1,PolylogInfinity
i = 1
del = 1._prec
do while (abs(del) > zero)
if(i**m.lt.0) return ! roll over
if(abs(x**i).lt.1.e-250_prec) return
res = res+x**i/i**m
enddo
del = x**i/i**m
res = res+del
i = i+1
end do
END FUNCTION naive_polylog
FUNCTION Li2(x)
......@@ -394,6 +405,19 @@ CONTAINS
END FUNCTION
FUNCTION MYLOG(x)
complex(kind=prec) :: x, mylog
if (abs(aimag(x)) < zero) then
if (real(x) > 0) then
mylog = cmplx(log(real(+x,kind=prec)), kind=prec)
else
mylog = cmplx(log(real(-x,kind=prec)), pi, kind=prec)
endif
else
mylog = log(x)
endif
END FUNCTION
RECURSIVE FUNCTION polylog1(m,x) result(res)
! computes the polylog
......@@ -420,10 +444,14 @@ CONTAINS
res = zeta(m)
else if ((m.le.9).and.(abs(x+1._prec).lt.zero)) then
res = -(1._prec - 2._prec**(1-m))*zeta(m)
else if ((m.le.9).and.(abs(x-i_).lt.zero)) then
res = -(0.5_prec**m - 0.5_prec**(2*m-1)) * zeta(m) + i_*dirichletbeta(m)
else if ((m.le.9).and.(abs(x+i_).lt.zero)) then
res = -(0.5_prec**m - 0.5_prec**(2*m-1)) * zeta(m) - i_*dirichletbeta(m)
else if (abs(x) .gt. 1) then
inv = 1._prec/x
res = (-1)**(m-1)*polylog(m,inv) &
- (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*log(-x)/2._prec/pi) / factorial(m)
- (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*mylog(-x)/2._prec/pi) / factorial(m)
else if(m == 2) then
res = dilog(x)
else if(m == 3) then
......@@ -479,7 +507,7 @@ CONTAINS
plog1 = plog1 - b%i0*i_*pi
endif
else
plog1 = log(1.-a%c/b%c)
plog1 = mylog(1.-a%c/b%c)
endif
END FUNCTION
......
......@@ -105,7 +105,7 @@ CONTAINS
enddo
if (mod(q,2_ikin) .eq. 1) then
if (abs(t(1)-res).lt.MPLdel/100.) exit
if (abs(t(1)-res).lt.MPLdelta/100.) exit
endif
enddo
res = t(1)
......