Commit d921d35f authored by Luca's avatar Luca
Browse files

naive polylog, medium naive dilog, evaluate depth 1 G-functions

parent cbffd701
......@@ -3,8 +3,9 @@ MODULE globals
implicit none
integer, parameter :: prec = selected_real_kind(15,32)
real, parameter :: zero = 1e-15 ! values smaller than this count as zero
integer, parameter :: GPLInfinity = 30 ! the default outermost expansion order for MPLs
real, parameter :: epsilon = 1e-15 ! used for the small imaginary part
real, parameter :: zero = 1e-15 ! values smaller than this count as zero
real, parameter :: pi = 3.14159265358979323846
END MODULE globals
......@@ -2,6 +2,7 @@
MODULE gpl_module
use globals
use utils
use maths_functions
use mpl_module
implicit none
......@@ -63,22 +64,22 @@ CONTAINS
RECURSIVE FUNCTION pending_integral(p,i,g) result(res)
! evaluates a pending integral by reducing it to simpler ones and g functions
complex(kind=prec) :: p(:), g(:), res
integer :: i
res = 0
if(i == size(g)+1) then
res = G_flat([p(2:size(p)),g], p(1))
return
end if
! if depth one and m = 1
if(size(g) == 1) then
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)))
return
end if
! evaluates a pending integral by reducing it to simpler ones and g functions
complex(kind=prec) :: p(:), g(:), res
integer :: i
res = 0
if(i == size(g)+1) then
res = G_flat([p(2:size(p)),g], p(1))
return
end if
! if depth one and m = 1
if(size(g) == 1) then
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)))
return
end if
END FUNCTION pending_integral
......@@ -97,8 +98,6 @@ CONTAINS
- G_flat([a(i+1)], sr) * G_flat(a(i+1:size(a)), y2)
return
end if
END FUNCTION reduce_to_convergent
RECURSIVE FUNCTION G_flat(z_flat,y) result(res)
......@@ -111,31 +110,31 @@ CONTAINS
call print_G(z_flat,y)
! is just a logarithm?
if(size(z_flat) == 1) then
print*, 'is just a logarithm'
if(abs(z_flat(1)) <= zero) then
res = log(y)
return
end if
res = log((z_flat(1) - y)/z_flat(1))
return
end if
! is it a polylog?
m_prime = get_condensed_m(z_flat)
m_1 = m_prime(1)
is_depth_one = (count((m_prime>0)) == 1)
if(is_depth_one) then
! case m >= 2, other already handled above
print*, 'is just a polylog'
res = -polylog(m_1,y/z_flat(m_1))
return
end if
! need make convergent?
if(.not. is_convergent(z_flat,y)) then
print*, 'need to make convergent'
m_prime = get_condensed_m(z_flat)
m_1 = m_prime(1)
is_depth_one = (count((m_prime>0)) == 1)
if(is_depth_one) then
! case m = 1
if(size(z_flat) == 1) then
print*, 'case m = 1'
res = log(y-z_flat(1)) - log(-z_flat(1))
return
end if
! case m >= 2
print*, 'case m > 1'
print*, 'm = ', m_1
res = -zeta(m_1) &
+ pending_integral([y, cmplx(0.0)], m_1-1, [zero_array(m_1-2), y]) &
- pending_integral([z_flat(m_1), cmplx(0.0)], m_1-1, [zero_array(m_1-2), y])
return
end if
res = reduce_to_convergent(z_flat, y)
return
end if
......
......@@ -7,8 +7,7 @@ FFLAGS=-fdefault-real-8 -cpp
LD=gfortran
OBJ= globals.o utils.o shuffle.o mpl_module.o gpl_module.o
OBJ= globals.o utils.o shuffle.o maths_functions.o mpl_module.o gpl_module.o
# Rules for main fortran files
%.o: %.f90
......
MODULE maths_functions
use globals
use utils
implicit none
! integer, parameter :: prec = selected_real_kind(15,32)
! integer, parameter :: GPLInfinity = 30 ! the default outermost expansion order for MPLs
! real(kind=prec), parameter :: epsilon = 1e-15 ! used for the small imaginary part
! real(kind=prec), parameter :: zero = 1e-15 ! values smaller than this count as zero
! real(kind=prec), parameter :: pi = 3.14159265358979323846_prec
CONTAINS
FUNCTION naive_polylog(m,x,n_passed) result(res)
! Computes the classical polylogarithm Li_m(x) using series representation up to order n
integer :: m
integer, optional :: n_passed
complex(kind=prec) :: x, res
integer :: i,n
integer, allocatable :: j(:)
n = merge(n_passed,GPLInfinity,present(n_passed))
j = (/(i, i=1,n,1)/)
res = sum(x**j / j**m)
END FUNCTION naive_polylog
FUNCTION Li2(x)
!! 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):: C(0:18),H,ALFA,B0,B1,B2,LI2_OLD
real (kind=prec):: Li2
integer :: i
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/
DATA C( 2) /-0.0185884366501460_prec/
DATA C( 3) / 0.0014575108406227_prec/
DATA C( 4) /-0.0001430418444234_prec/
DATA C( 5) / 0.0000158841554188_prec/
DATA C( 6) /-0.0000019078495939_prec/
DATA C( 7) / 0.0000002419518085_prec/
DATA C( 8) /-0.0000000319334127_prec/
DATA C( 9) / 0.0000000043454506_prec/
DATA C(10) /-0.0000000006057848_prec/
DATA C(11) / 0.0000000000861210_prec/
DATA C(12) /-0.0000000000124433_prec/
DATA C(13) / 0.0000000000018226_prec/
DATA C(14) /-0.0000000000002701_prec/
DATA C(15) / 0.0000000000000404_prec/
DATA C(16) /-0.0000000000000061_prec/
DATA C(17) / 0.0000000000000009_prec/
DATA C(18) /-0.0000000000000001_prec/
if(X > 1.00000000001_prec) then
print*, 'crashes because Li called with bad arguments'
elseif(X > 1.0_prec) then
X = 1._prec
endif
IF(X > 0.999999_prec) THEN
LI2_OLD=PI6
Li2 = Real(LI2_OLD,prec)
RETURN
ELSE IF(X .EQ. MONE) THEN
LI2_OLD=MALF*PI6
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)
ELSE IF(T .LT. MONE) THEN
Y=MONE-T
S=MONE
A=LOG(-T)
A=-PI6+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))
ELSE IF(T .LT. ZERO) THEN
Y=-T/(ONE+T)
S=MONE
A=HALF*LOG(ONE+T)**2
ELSE IF(T .LE. ONE) THEN
Y=T
S=ONE
A=ZERO
ELSE
Y=ONE/T
S=MONE
A=PI6+HALF*LOG(T)**2
END IF
H=Y+Y-ONE
ALFA=H+H
B1=ZERO
B2=ZERO
DO I = 18,0,-1
B0=C(I)+ALFA*B1-B2
B2=B1
B1=B0
ENDDO
LI2_OLD=-(S*(B0-H*B2)+A)
! Artificial conversion
Li2 = Real(LI2_OLD,prec)
END FUNCTION Li2
FUNCTION dilog_in_unit_circle(x) result(res)
! evaluates for any argument x in unit circle
complex(kind=prec) :: x, res
res = naive_polylog(2,x)
END FUNCTION dilog_in_unit_circle
FUNCTION dilog(x) result(res)
! evaluates dilog for any argument
complex(kind=prec) :: res
complex(kind=prec) :: x
if(abs(x) <= 1.0) then
res = naive_polylog(2,x)
else
res = -dilog_in_unit_circle(1/x) - pi**2/6 - log(dcmplx( add_ieps(-x) ))**2 / 2
end if
END FUNCTION dilog
FUNCTION polylog(m,x) result(res)
! computes the polylog for now naively (except for dilog half-naively)
integer :: m
complex(kind=prec) :: x,res
print*, 'called polylog(',m,',',x,')'
if(m == 2) then
res = dilog(x)
else
res = naive_polylog(m,x)
end if
END FUNCTION polylog
END MODULE maths_functions
! PROGRAM test
! use maths_functions
! implicit none
! complex(kind=prec) :: res
! res = dilog(dcmplx(0.4d0))
! print*, res
! END PROGRAM test
......@@ -6,27 +6,17 @@ MODULE mpl_module
CONTAINS
FUNCTION dilog(x,n)
! Computes the dilog Li_2(x) using the series representation up to order n
integer :: n
complex(kind=prec) :: x, dilog
integer :: i
integer :: j(n)
j = (/(i, i=1,n,1)/)
dilog = sum(x**j / j**2)
END FUNCTION dilog
FUNCTION polylog(m,x,n)
FUNCTION polylog_series(m,x,n_passed) result(res)
! Computes the classical polylogarithm Li_m(x) using series representation up to order n
integer :: m
complex(kind=prec) :: x, polylog
integer, optional :: n_passed
complex(kind=prec) :: x, res
integer :: i,n
integer, allocatable :: j(:)
n = merge(n_passed,GPLInfinity,present(n_passed))
j = (/(i, i=1,n,1)/)
polylog = sum(x**j / j**m)
END FUNCTION polylog
res = sum(x**j / j**m)
END FUNCTION polylog_series
FUNCTION MPL_converges(m,x)
! checks if the MPL converges
......@@ -57,7 +47,7 @@ CONTAINS
if(size(m) == 1) then
! base case
res = polylog(m(1),x(1), n)
res = polylog_series(m(1),x(1), n)
else
! recursion step
res = 0
......
This diff is collapsed.
MODULE GLOBAL_DEF
implicit none
!! ----------
!! parameters
!! ----------
integer, parameter :: prec = selected_real_kind(15,32)
real (kind=prec), parameter :: cw = 0.876613
real (kind=prec), parameter :: sw = 0.481196
real (kind=prec), parameter :: pi = 3.14159265358979323846_prec
real (kind=prec), parameter :: z3 = 1.20205690315959428540_prec
real (kind=prec), parameter :: log2 = 0.693147180559945309417_prec
real (kind=prec), parameter :: conv = 3.893850E+8 ! convert GeV to pb
real (kind=prec), parameter :: xsc = 0._prec ! FDH=>1 vs HV=>0
real (kind=prec), parameter :: Nc = 3._prec
real (kind=prec), parameter :: Tf = 0.5_prec
real (kind=prec), parameter :: Cf = (Nc**2-1)/(2*Nc)
real (kind=prec), parameter :: Nh = 1._prec
real (kind=prec), parameter :: Nf = 5._prec
complex (kind=prec), parameter :: imag = (0.0_prec,1.0_prec)
real (kind=prec), parameter :: zero = 1.0E-50_prec
real (kind=prec), parameter :: alpha_ew = 0.03394_prec
! real (kind=prec), parameter :: alpha = 1/127.9_prec
! real (kind=prec), parameter :: alpha = 1./137.0359997_prec
! real (kind=prec), parameter :: GF = 1.16637E-11_prec ! MeV^-2
real (kind=prec), parameter :: GF = 1._prec
real (kind=prec), parameter :: alpha = 1._prec
real (kind=prec), parameter :: Mmu = 105.658372_prec ! MeV
real (kind=prec), parameter :: Mel = 0.51099893_prec ! MeV
! real (kind=prec), parameter :: Mel = 10._prec ! MeV
real (kind=prec), parameter :: Mtau = 1776.82_prec ! MeV
real (kind=prec), parameter :: xi_sep = 1.0E-10_prec
real (kind=prec), parameter :: del_sep = 1.0E-10_prec
! character (len=3), parameter :: cgamma = "exp"
character (len=3), parameter :: cgamma = "gam"
integer print_ok, throw_away
!! ---------
!! variables
!! ---------
integer :: ran_seed = 1
real (kind=prec) :: p1(4),p2(4),p3(4),p4(4),p5(4),p6(4),p7(4), &
p8(4),p9(4), pol1(4)
real (kind=prec) :: mu, musq, delcut, xinormcut
real (kind=prec) :: xinormcut1, xinormcut2
character (len=8) :: flavour
real (kind=prec) :: Mm ! MeV
real (kind=prec) :: Me ! MeV
contains
SUBROUTINE CRASH(function_name)
character(len=*) :: function_name
write(6,*) "Program crashes because of a call to the function ", &
function_name
stop
END SUBROUTINE CRASH
END MODULE GLOBAL_DEF
......@@ -6,6 +6,7 @@ PROGRAM TEST
use globals
use utils
use shuffle
use maths_functions
use mpl_module
use gpl_module
implicit none
......@@ -19,13 +20,12 @@ PROGRAM TEST
! call do_GPL_tests()
! call do_shuffle_tests() ! put this somewhere else
! res = G_flat(cmplx((/0.0,10.0/)),cmplx(20.0))
! print*, res
res = G_flat(cmplx((/0.0,1.0/)),cmplx(2.0))
print*, res
do i = 1,8
res = zeta(i)
print*, 'zeta(',i,') =', res
end do
! res = polylog(2,cmplx(2.0))
! print*, log((-2.0,-0.0000000000000001))
! print*, log(add_ieps((-2.0,-0.0000000000000001)))
! if(tests_successful) then
! print*, 'All tests passed. '
......
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