Skip to content

GitLab

H
handyG
Commit 6161ad7a authored by Luca Naterop's avatar Luca Naterop
Browse files
Options

do not compare complex numbers, more tests, catch G(1,1)

parent 24510811
Hide whitespace changes
Inline Side-by-side
Showing with 95 additions and 48 deletions
.gitignore
View file @ 6161ad7a
......@@ -13,4 +13,5 @@
*.o
*.mod
*.smod
.DS_Store
\ No newline at end of file
.DS_Store
*.a
\ No newline at end of file
makefile
View file @ 6161ad7a
......@@ -10,7 +10,7 @@ MODE=DEBUG
FC=gfortran
AR= ar rcs
FFLAGS=-fdefault-real-8 -cpp -pedantic-errors -std=f2008
FFLAGS+= -Werror -Wall -Wno-maybe-uninitialized -Wno-uninitialized -Wno-compare-reals
FFLAGS+= -Werror -Wall -Wno-maybe-uninitialized -Wno-uninitialized
ifeq ($(MODE),RELEASE)
FFLAGS += -O3 -funroll-loops -Wtaps
......
src/eval.f90
View file @ 6161ad7a
......@@ -3,6 +3,7 @@ PROGRAM eval
use globals
use gpl_module
use utils
use maths_functions
use shuffle
implicit none
......@@ -11,7 +12,7 @@ PROGRAM eval
call parse_cmd_args()
res = GPL([2,1,3,4])
res = GPL([0.1,-1.,1.])
print*, res
! res = pending_integral(cmplx([1,3]),2,cmplx([2,4]))
......
src/gpl_module.f90
View file @ 6161ad7a
......@@ -32,7 +32,7 @@ CONTAINS
if(all(abs(y) <= abs(z))) then
if(m(1) == 1) then
GPL_has_convergent_series = (y/z(1) /= 1)
GPL_has_convergent_series = .true. !(abs((y/z(1) - 1)) < zero)
else
GPL_has_convergent_series = .true.
end if
......@@ -63,7 +63,7 @@ CONTAINS
SUBROUTINE print_G(z_flat, y)
complex(kind=prec) :: z_flat(:)
complex(kind=prec), optional :: y
if(present(y)) print*, 'G(', abs(z_flat), abs(y), ')'
if(present(y)) print*, 'G(', real(z_flat), real(y), ')'
if(.not. present(y)) print*, 'G(', abs(z_flat), ')'
END SUBROUTINE print_G
......@@ -80,14 +80,14 @@ CONTAINS
if(verb >= 50) then
print*, 'mapping to '
call print_G(a,y2)
print*, 'PI with p=',abs(p),'i=',m,'g =',abs([zeroes(m-1),y2])
print*, 'PI with p=',real(p),'i=',m,'g =',real([zeroes(m-1),y2])
end if
res = GPL(a,y2)*pending_integral(p,m,[zeroes(m-1),y2])
if(verb >= 50) print*, 'also mapping to'
do j = 2,size(alpha, 1)
! find location of s_r
n = find_first_true(alpha(j,:) == 42e50)
if(verb >= 50) print*, 'PI with p=',abs(p),'i=',n,'g =',abs([alpha(j,1:n-1),alpha(j,n+1:size(alpha,2)),y2])
n = find_first_true(abs(alpha(j,:) - 42e50) < zero)
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])
res = res - pending_integral(p, n, [alpha(j,1:n-1),alpha(j,n+1:size(alpha,2)),y2])
end do
END FUNCTION remove_sr_from_last_place_in_PI
......@@ -99,7 +99,7 @@ CONTAINS
integer :: i, m
res = 0
if(verb >= 30) print*, 'evaluating PI with p=',abs(p),'i=',abs(i),'g =',abs(g)
if(verb >= 30) print*, 'evaluating PI with p=',real(p),'i=',real(i),'g =',real(g)
y1 = p(1)
b = p(2:size(p))
......@@ -137,9 +137,9 @@ CONTAINS
if(verb >= 50) then
print*, 'map to'
print*, 'zeta(',m,')'
print*, 'PI with p=',abs(p),'i=',0,'g =',zeroes(0)
print*, 'PI with p=',abs([y2,cmplx(0.0)]),'i=',m-1,'g =',[zeroes(m-2),y2]
print*, 'PI with p=',abs(p),'i=',0,'g =',zeroes(0)
print*, 'PI with p=',real(p),'i=',0,'g =',zeroes(0)
print*, 'PI with p=',real([y2,cmplx(0.0)]),'i=',m-1,'g =',[zeroes(m-2),y2]
print*, 'PI with p=',real(p),'i=',0,'g =',zeroes(0)
print*, 'PI with p=',[p,cmplx(0.0)],'i=',m-1,'g =',zeroes(0)
end if
res = -zeta(m)*pending_integral(p,0,zeroes(0)) &
......@@ -154,12 +154,12 @@ CONTAINS
if(verb >= 50) then
print*, 'map to (using 68)'
print*, 'PI with p=',abs(p),'i=',0,'g =',zeroes(0)
print*, 'PI with p=',real(p),'i=',0,'g =',zeroes(0)
call print_G([cmplx(0.0),a],y2)
print*, 'PI with p=',abs([p,y2]),'i=',0,'g =',zeroes(0)
print*, 'PI with p=',real([p,y2]),'i=',0,'g =',zeroes(0)
call print_G(a,y2)
print*, 'PI with p=',abs([p,a(1)]),'i=',1,'g =',abs([a(2:size(a)),y2])
print*, 'PI with p=',abs([p,a(1)]),'i=',0,'g =',zeroes(0)
print*, 'PI with p=',real([p,a(1)]),'i=',1,'g =',real([a(2:size(a)),y2])
print*, 'PI with p=',real([p,a(1)]),'i=',0,'g =',zeroes(0)
call print_G(a,y2)
end if
......@@ -188,7 +188,7 @@ CONTAINS
- pending_integral([p,a(i)],1,zeroes(0)) * G_flat(a,y2)
END FUNCTION pending_integral
FUNCTION remove_sr_from_last_place_in_G(a,y2,m,sr) result(res)
RECURSIVE FUNCTION remove_sr_from_last_place_in_G(a,y2,m,sr) result(res)
complex(kind=prec) :: a(:), sr, res,y2
integer :: m,i,j
complex(kind=prec) :: alpha(product((/(i,i=1,size(a)+m)/))/ &
......@@ -201,7 +201,7 @@ CONTAINS
end do
END FUNCTION remove_sr_from_last_place_in_G
FUNCTION make_convergent(a,y2) result(res)
RECURSIVE FUNCTION make_convergent(a,y2) result(res)
! goes from G-functions to pending integrals and simpler expressions according to (62),(64),(67) and (68)
complex(kind=prec) :: a(:), y2, res, sr
......@@ -228,7 +228,7 @@ CONTAINS
call print_G([cmplx(0), a(i+1:size(a))], y2)
call print_G([y2], sr)
call print_G(a(i+1:size(a)), y2)
print*, 'PI with p=',abs([sr, a(i+1)]),'i=',i,'g =', abs([a(i+2:size(a)), y2])
print*, 'PI with p=',real([sr, a(i+1)]),'i=',i,'g =', real([a(i+2:size(a)), y2])
call print_G([a(i+1)], sr)
call print_G(a(i+1:size(a)), y2)
print*, '--------------------------------------------------'
......@@ -246,10 +246,10 @@ CONTAINS
print*, '--------------------------------------------------'
print*, 'sr in the middle, map to: '
call print_G([a(1:i-1),cmplx(0),a(i+1:size(a))],y2)
print*, 'PI with p=',abs([sr,a(i-1)]),'i=', i-1,'g =', abs([a(1:i-2),a(i+1:size(a)),y2])
print*, 'PI with p=',real([sr,a(i-1)]),'i=', i-1,'g =', real([a(1:i-2),a(i+1:size(a)),y2])
call print_G([a(i-1)],sr)
call print_G([a(1:i-1),a(i+1:size(a))],y2)
print*, 'and PI with p=',abs([sr,a(i+1)]),'i=',i,'g =', abs([a(1:i-1),a(i+2:size(a)),y2])
print*, 'and PI with p=',real([sr,a(i+1)]),'i=',i,'g =', real([a(1:i-1),a(i+2:size(a)),y2])
call print_G([a(i+1)],sr)
call print_G([a(1:i-1),a(i+1:size(a))],y2)
print*, '--------------------------------------------------'
......@@ -276,7 +276,6 @@ CONTAINS
do j = 1,k-1
res = res + (-1)**j * G_flat(oneminusz(j:1:-1),1.0-1.0/p) * G_flat(z(j+1:k),1.0/p)
end do
END FUNCTION improve_convergence
RECURSIVE FUNCTION G_flat(z_flat,y) result(res)
......@@ -289,6 +288,15 @@ CONTAINS
if(verb >= 50) call print_G(z_flat,y)
! catch G(1,1)
if(size(z_flat) == 1) then
if( abs(z_flat(1) - y) <= zero ) then
print*, 'catch G(1,1)'
res = 0
return
end if
end if
! add small imaginary part if not there
! do i = 1,size(z_flat)
! if(abs(aimag(z_flat(i))) < 1e-25) z_flat(i) = add_ieps(z_flat(i))
......@@ -350,7 +358,7 @@ CONTAINS
end if
! requires Hoelder convolution?
if( any(1.0 <= abs(z_flat/y) .and. abs(z_flat/y) <= 2.0) ) then
if( any(1.0 <= abs(z_flat/y) .and. abs(z_flat/y) <= 1.1) ) then
res = improve_convergence(z_flat/y)
return
end if
......
src/maths_functions.f90
View file @ 6161ad7a
......@@ -61,7 +61,7 @@ CONTAINS
LI2_OLD=PI6
Li2 = Real(LI2_OLD,prec)
RETURN
ELSE IF(X .EQ. MONE) THEN
ELSE IF(abs(x-MONE) < zero) THEN
LI2_OLD=MALF*PI6
RETURN
END IF
......@@ -120,7 +120,7 @@ CONTAINS
res = naive_polylog(2,x)
endif
else
res = -dilog(1/x) - pi**2/6 - log(add_ieps(-x))**2 / 2
res = -dilog(1/x) - (pi**2) /6 - log(add_ieps(-x))**2 / 2
end if
END FUNCTION dilog
......@@ -185,7 +185,7 @@ CONTAINS
IF(X > 0.999999_prec) THEN
LI3=zeta3
RETURN
ELSE IF(X .EQ. -1._prec) THEN
ELSE IF( abs(x+1) < zero) THEN
LI3=-0.75_prec*zeta3
RETURN
END IF
......
src/mpl_module.f90
View file @ 6161ad7a
......@@ -25,7 +25,7 @@ CONTAINS
logical :: MPL_converges
MPL_converges = .false.
if(abs(product(x)) < 1) then
if(m(1) /= 1 .or. x(1) /= 1) then
if(m(1) /= 1 .or. abs(x(1) - 1) < zero) then
MPL_converges = .true.
end if
end if
......
src/test.f90
View file @ 6161ad7a
......@@ -15,10 +15,9 @@ PROGRAM TEST
logical :: tests_successful = .true.
call parse_cmd_args()
call do_MPL_tests()
call do_GPL_tests()
call do_shuffle_tests() ! put this somewhere else
! call do_shuffle_tests() ! put this somewhere else
if(tests_successful) then
......@@ -88,65 +87,103 @@ CONTAINS
subroutine do_GPL_tests()
complex(kind=prec) :: ref
complex(kind=prec), parameter :: epsilon = 1E-14
real(kind=prec) :: z, xchen
print*, 'doing GPL tests...'
ref = cmplx(0.0819393734128676)
call test_one_condensed((/ 1,1 /),cmplx((/ 1.3_prec, 1.1_prec /)),cmplx(0.4),2,ref,'2.1')
call test_one_condensed((/ 1,1 /),cmplx((/ 1.3, 1.1 /)),cmplx(0.4),2,ref,'2.1')
ref = cmplx(0.01592795952537145)
call test_one_condensed((/ 3,2 /),cmplx((/ 1.3_prec, 1.1_prec /)),cmplx(0.4),2,ref,'2.2')
call test_one_condensed((/ 3,2 /),cmplx((/ 1.3, 1.1 /)),cmplx(0.4),2,ref,'2.2')
ref = cmplx(0.0020332632172573974)
call test_one_condensed((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3')
! requires making convergent
ref = cmplx((0.09593041677639341, -0.8829351795197851))
ref = (0.09593041677639341, -0.8829351795197851)
call test_one_flat(cmplx([0,1,3,2]),ref,'2.4')
ref = cmplx((0.009947947789928621,0.0))
ref = (0.009947947789928621,0.0)
call test_one_flat(cmplx([0.0, 0.0, -3.3333333333333335, -3.3333333333333335, 1.0]),ref,'2.5')
! requires hoelder convolution
ref = cmplx((-0.012709942828250949,0.0))
ref = (-0.012709942828250949,0.0)
call test_one_flat(cmplx([0.0, 3.3333333333333335, 1.0, 3.3333333333333335, 1.0]),ref,'2.6')
! here the tests from the mathematica nb start
! --------------------------------------------------------------------------
z = 1./200; xchen = 0.3;
ref = cmplx((-0.0050125418235441935,0.0))
ref = (-0.0050125418235441935,0.0)
call test_one_flat(cmplx([1./z,1.0]),ref,'3.1')
ref = cmplx((-0.0015011261262671913,0.0))
ref = (-0.0015011261262671913,0.0)
call test_one_flat(cmplx([1./(xchen*z),1.0]),ref,'3.2')
ref = cmplx((-0.0007502860817810596,0.0))
ref = (-0.0007502860817810596,0.0)
call test_one_flat(cmplx([(1+sqrt(1-z**2))/(xchen*z),1.0]),ref,'3.3')
ref = cmplx((0.0074335969894765335,0.0))
ref = (0.0074335969894765335,0.0)
call test_one_flat(cmplx([-1./xchen,-1./xchen,1.,1.,1.0]),ref,'3.4')
ref = cmplx((-8.403785974849544e-6,0.0))
ref = (-8.403785974849544e-6,0.0)
call test_one_flat(cmplx([-1./xchen,0.,-1./xchen,1./(xchen*z),1.0]),ref,'3.5')
! ref = cmplx((0.4925755847450199,2.6389214054743295))
! call test_one_flat(cmplx([-1.,-1.,z,z,1.]),ref,'3.6')
ref = cmplx((-0.0015317713178859967,-0.00045003911367000565))
call test_one_flat(cmplx([0.,0.,(1-sqrt(1-z**2))/(xchen*z), 1./(xchen*z),1.]),ref,'3.7')
! ref = cmplx((-0.0015317713178859967,-0.00045003911367000565))
! call test_one_flat(cmplx([0.,0.,(1-sqrt(1-z**2))/(xchen*z), 1./(xchen*z),1.]),ref,'3.7')
! here the chen integral tests start
! ----------------------------------------------------------------------
ref = (-1.2039728043259361,0.0)
call test_one_flat(cmplx([0., 0.3]),ref,'4.1')
ref = (10.603796220956799,0.0)
call test_one_flat(cmplx([0., 0., 0.01]),ref,'4.2')
ref = (0.0005042990065180765,0.0)
call test_one_flat(cmplx([100., 1., 0.3]),ref,'4.3')
ref = (0.05630861877185141,0.0)
call test_one_flat(cmplx([1., 0., 0.01]),ref,'4.4')
ref = (0.04032895150872735,0.2922709647568842)
call test_one_flat([(0.01,0.9999499987499375), (0.3,0)],ref,'4.5')
ref = cmplx(0.000025210645098340354)
call test_one_flat(cmplx([100., 199.99499987499374, 1.]),ref,'4.6')
ref = (0.0556470547632135,0.)
call test_one_flat(cmplx([-1., 0.01, 0., 0.01]),ref,'4.7')
ref = (0.00003794895846844715,0.)
call test_one_flat(cmplx([100., 100., 1., 0., 1.]),ref,'4.8')
ref = (0.00007574284433216497,0.)
call test_one_flat(cmplx([100., 1., 100., 0., 1.]),ref,'4.9')
ref = (1.8801911586012443,2.509434138598062)
call test_one_flat(cmplx([0.01, -1.0, 0.01, 1.]),ref,'4.10')
ref = (-0.012539108315054982, -0.015414250168437678)
call test_one_flat(cmplx([0.01, 199.99499987499374, 0.01, 1.]),ref,'4.10')
! (1+sqrt(1-z**2))/(xchen*z) -1./xchen
end subroutine do_GPL_tests
subroutine do_shuffle_tests()
complex(kind=prec) :: v(2) = cmplx((/1,2/))
complex(kind=prec) :: w(2) = cmplx((/3,4/))
! subroutine do_shuffle_tests()
! complex(kind=prec) :: v(2) = cmplx((/1,2/))
! complex(kind=prec) :: w(2) = cmplx((/3,4/))
call print_matrix(shuffle_product(v,w))
end subroutine do_shuffle_tests
! call print_matrix(shuffle_product(v,w))
! end subroutine do_shuffle_tests
END PROGRAM TEST
......
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