Commit 14e345d6 authored by ulrich_y's avatar ulrich_y

Keep track of complex sign of sr

parent b5751dad
......@@ -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
......
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