Commit 646afc83 authored by Yannick Ulrich's avatar Yannick Ulrich

Merge branch 'bugfix-polylog-conjg'

parents 87817de8 09592838
This diff is collapsed.
......@@ -534,10 +534,6 @@ CONF_LD=${LD:-$CONF_FC}
if $CONF_QUAD ; then
if $HAVE_GINAC ; then
echo "GiNaC testing is not supported for quad-precision!" 1>&3
exit 1
fi
echo -n "does $CONF_FC support quad-precision... " 1>&3
rm -fr $test*
tee $test.f90 << _EOF_ 1>&2
......@@ -568,7 +564,7 @@ if $HAVE_GINAC ; then
eval addflag CXXFLAGS "-std=c++11"
if [[ ! -z "$CONF_PKGCONFIG" ]]; then
echo -n "Does pkg-config know about GiNaC... " 1>&3
if $CONF_PKGCONFIG --exists ginac ; then
if $CONF_PKGCONFIG --atleast-version=1.7.4 ginac ; then
echo "yes" 1>&3
eval addflag CXXFLAGS `$CONF_PKGCONFIG --cflags ginac`
eval addflag LFLAGS `$CONF_PKGCONFIG --libs ginac`
......@@ -579,24 +575,63 @@ if $HAVE_GINAC ; then
FOUND_GINAC=false
fi
fi
if ! $FOUNDGINAC ; then
FOUNDGINAC=`findlib ginac GINAC && findlib cln CLN`
if ! $FOUND_GINAC ; then
findlib ginac GINAC && findlib cln CLN && FOUND_GINAC=true
eval addflag LFLAGS "-L`dirname $CONF_GINAC`"
eval addflag LFLAGS "-L`dirname $CONF_CLN`"
eval addflag LFLAGS "-lginac -lcln"
fi
if $FOUNDGINAC; then
if $FOUND_GINAC; then
CONF_LD=${LD:-$CONF_CXX}
echo -n "Checking if GiNaC works... " 1>&3
tee $test.ginac.cpp << _EOF_ 1>&2
#include <ginac/ginac.h>
#include <iostream>
int main() {
#include <cln/cln.h>
typedef struct {double r,i;} complex_t;
typedef struct {complex_t c; signed char i0;} inum_t;
complex_t geval_(inum_t * z, int* n) {
cln::cl_inhibit_floating_point_underflow = true;
GiNaC::lst w,s;
for(long i=0;i<(*n)-1;i++)
{
GiNaC::ex zz;
if (abs(z->c.i) < 1e-15)
w.append((z->c.r));
else
w.append((z->c.r)+(z->c.i)*GiNaC::I);
s.append(z->i0);
z++;
}
GiNaC::ex ans = GiNaC::G(w,s,z->c.r).evalf();
return {
.r = GiNaC::ex_to<GiNaC::numeric>(GiNaC::evalf(GiNaC::real_part(ans))).to_double(),
.i = GiNaC::ex_to<GiNaC::numeric>(GiNaC::evalf(GiNaC::imag_part(ans))).to_double()
};
}
int main() {
// Test 1
GiNaC::ex ans = GiNaC::G(0.3,0.5,0.7);
ans -= 0.2876820724517808812+3.1415926535897932385*GiNaC::I;
ans = GiNaC::abs(ans);
if (ans<1e-15)
// Test 2
inum_t x[] = {
{ .c = {.r = 0.73, .i=0}, .i0 = -1 },
{ .c = {.r = 0.00, .i=0}, .i0 = +1 },
{ .c = {.r = 1.00, .i=0}, .i0 = +1 }
};
int n=3;
complex_t cans = geval_(&x[0],&n);
float del = (2.2982889094420660-cans.r)*(2.2982889094420660-cans.r) + (0.98869296399417417-cans.i)*(0.98869296399417417-cans.i);
if ( (ans<1e-15) && (del<1e-15) )
return 0;
else
return 1;
......
......@@ -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,7 +15,7 @@ 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.
real(kind=prec), protected :: MPLdelta = 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 /)
......@@ -49,7 +49,7 @@ CONTAINS
SUBROUTINE SET_OPTIONS(mpldel, liinf, hcircle)
real(kind=prec), optional :: hcircle, mpldel
integer, optional :: liinf
if (present(mpldel)) MPLdel = mpldel
if (present(MPLdel)) MPLdelta = mpldel
if (present(liinf)) PolyLogInfinity = liinf
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
......@@ -394,6 +400,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 +439,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 +502,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)
......
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