Commit bdc545ff by ulrich_y

### Undid some stuff related to ieps

parent 522fa0b7
 ... ... @@ -18,7 +18,7 @@ CONTAINS integer :: l type(inum) :: y complex(kind=prec) :: GPL_zero_zi GPL_zero_zi = 1.0_prec/factorial(l) * log(y) ** l GPL_zero_zi = 1.0_prec/factorial(l) * log(y%c) ** l END FUNCTION GPL_zero_zi FUNCTION is_convergent(z,y) ... ... @@ -103,7 +103,7 @@ CONTAINS !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))) res = pending_integral(p,2,[g(1)]) - pending_integral(p,2,[izero]) & + G_flat(p(2:size(p)), p(1)) * log(neg(g(1))) + G_flat(p(2:size(p)), p(1)) * (log(-g(1)%c)+cmplx(0,2*pi)) return end if ... ... @@ -248,16 +248,19 @@ CONTAINS ! improves the convergence by applying the Hoelder convolution to G(z1,...zk,1) type(inum) :: z(:),oneminusz(size(z)) complex(kind=prec) :: res type(inum), parameter :: p = inum(2.0,+1) complex(kind=prec), parameter :: p = 2.0 integer :: k, j if(verb >= 30) print*, 'requires Hoelder convolution' oneminusz = ione-z !TODO ieps?!?? do j=1,size(z) oneminusz(j) = inum(1.-z(j)%c,-z(j)%i0) enddo k = size(z) res = G_flat(z,ione/p) ! first term of the sum res = res + (-1)**k * G_flat(oneminusz(k:1:-1), ione-ione/p) res = G_flat(z,inum(1./p,di0)) ! first term of the sum res = res + (-1)**k * G_flat(oneminusz(k:1:-1), inum(1.-1/p,di0)) do j = 1,k-1 res = res + (-1)**j * G_flat(oneminusz(j:1:-1),ione-ione/p) * G_flat(z(j+1:k),ione/p) res = res + (-1)**j * G_flat(oneminusz(j:1:-1),inum(1.-1/p,di0)) * G_flat(z(j+1:k),inum(1./p,di0)) end do END FUNCTION improve_convergence ... ... @@ -274,7 +277,7 @@ CONTAINS if(size(z_flat) == 1) then if( abs(z_flat(1) - y) <= zero ) then if( abs(z_flat(1)%c - y%c) <= zero ) then res = 0 return end if ... ... @@ -289,16 +292,16 @@ CONTAINS ! is just a logarithm? if(all(abs(z_flat) < zero)) then if(verb >= 70) print*, 'all z are zero' res = log(y)**size(z_flat) / factorial(size(z_flat)) res = log(y%c)**size(z_flat) / factorial(size(z_flat)) return end if if(size(z_flat) == 1) then if(verb >= 70) print*, 'is just a logarithm' if(abs(z_flat(1)) <= zero) then res = log(y) res = log(y%c) return end if res = log((z_flat(1) - y)/z_flat(1)) res = plog1(y,z_flat(1)) ! log((z_flat(1) - y)/z_flat(1)) return end if ... ... @@ -309,7 +312,7 @@ CONTAINS if(is_depth_one) then ! case m >= 2, other already handled above if(verb >= 70) print*, 'is just a polylog' res = -polylog(m_1,y/z_flat(m_1)) res = -polylog(m_1, y, z_flat(m_1))!-polylog(m_1,y/z_flat(m_1)) return end if ... ... @@ -324,7 +327,7 @@ CONTAINS if(verb >= 30) print*, 'need to remove trailing zeroes' allocate(s(j,j)) s = shuffle_with_zero(z_flat(1:j-1)) res = log(y)*G_flat(z_flat(1:size(z_flat)-1),y) res = log(y%c)*G_flat(z_flat(1:size(z_flat)-1),y) do i = 1,size(s,1) res = res - G_flat([s(i,:),z_flat(j),zeroes(kminusj-1)], y) end do ... ... @@ -342,7 +345,8 @@ CONTAINS ! requires Hoelder convolution? if( any(1.0 <= abs(z_flat%c/y%c) .and. abs(z_flat%c/y%c) <= 1.1) ) then res = improve_convergence(z_flat/y) !TODO res = improve_convergence(toinum(z_flat%c/y%c)) return end if ... ... @@ -390,8 +394,8 @@ CONTAINS ! assumes zero arguments expressed through the m's integer :: m(:), k, i type(inum) :: z(:), y, x(k), z_flat(sum(m)) complex(kind=prec) :: res type(inum) :: z(:), y, z_flat(sum(m)) complex(kind=prec) :: res, x(k) ! print*, 'called G_condensed with args' ! print*, 'm = ', m ... ... @@ -418,15 +422,15 @@ CONTAINS res = G_flat(z_flat,y) return end if x(1) = y/z(1) !TODO is that okay? x(1) = y%c/z(1)%c do i = 2,k x(i) = z(i-1)/z(i) x(i) = z(i-1)%c/z(i)%c end do ! print*, 'computed using Li with ' ! print*, 'm = ', m ! print*, 'x = ', x res = (-1)**k * MPL(m,x%c) res = (-1)**k * MPL(m,x) END FUNCTION G_condensed ... ...
 ... ... @@ -15,27 +15,9 @@ MODULE ieps type(inum), parameter :: marker=inum(0.,5) interface operator (*) module procedure multinumss, multinumvs end interface operator (*) interface operator (+) module procedure addinumss, addinumvs end interface operator (+) interface operator (-) module procedure subinumss,subinumvs,subinumsv end interface operator (-) interface operator (**) module procedure powinum end interface operator (**) interface operator (/) module procedure divint, divinumss, divinumvs end interface operator (/) interface abs module procedure absinum, absinumv end interface abs interface log module procedure loginum end interface log interface toinum module procedure toinum_cmplx, toinum_real, toinum_int ... ... @@ -47,79 +29,6 @@ MODULE ieps module procedure realis, realiv end interface real CONTAINS FUNCTION NEG(n1) implicit none type(inum), intent(in) :: n1 type(inum) :: neg neg = inum(-n1%c,-n1%i0) END FUNCTION FUNCTION MULTINUMSS(n1, n2) implicit none type(inum), intent(in) :: n1, n2 type(inum) :: multinumss multinumss = inum( n1%c*n2%c, int(sign(1._prec,real(n1%c)*n2%i0 + real(n2%c)*n1%i0)) ) END FUNCTION MULTINUMSS FUNCTION MULTINUMVS(n1, n2) implicit none type(inum), intent(in) :: n1(:), n2 type(inum) :: multinumvs(size(n1)) integer i do i = 1,size(n1) multinumvs(i) = inum( n1(i)%c*n2%c, int(sign(1._prec,real(n1(i)%c)*n2%i0 + real(n2%c)*n1(i)%i0)) ) enddo END FUNCTION MULTINUMVS FUNCTION ADDINUMSS(n1, n2) implicit none type(inum), intent(in) :: n1, n2 type(inum) :: addinumss !TODO: what *is* the sum? addinumss = inum(n1%c + n2%c, n1%i0 ) END FUNCTION ADDINUMSS FUNCTION ADDINUMVS(n1, n2) implicit none type(inum), intent(in) :: n1(:), n2 type(inum) :: addinumvs(size(n1)) !TODO: what *is* the sum? integer i do i = 1,size(n1) addinumvs(i) = inum(n1(i)%c + n2%c, n1(i)%i0 ) enddo END FUNCTION ADDINUMVS FUNCTION SUBINUMSS(n1, n2) implicit none type(inum), intent(in) :: n1, n2 type(inum) :: subinumss !TODO: what *is* the sum? subinumss = inum(n1%c - n2%c, n1%i0 ) END FUNCTION SUBINUMSS FUNCTION SUBINUMVS(n1, n2) implicit none type(inum), intent(in) :: n1(:), n2 type(inum) :: subinumvs(size(n1)) !TODO: what *is* the sum? integer i do i = 1,size(n1) subinumvs(i) = inum(n1(i)%c - n2%c, n1(i)%i0 ) enddo END FUNCTION SUBINUMvs FUNCTION SUBINUMSV(n2, n1) implicit none type(inum), intent(in) :: n1(:), n2 type(inum) :: subinumsv(size(n1)) !TODO: what *is* the sum? integer i do i = 1,size(n1) subinumsv(i) = inum(n2%c - n1(i)%c, n1(i)%i0 ) enddo END FUNCTION SUBINUMSV FUNCTION ABSINUM(n1) implicit none type(inum), intent(in) :: n1 ... ... @@ -134,59 +43,6 @@ CONTAINS absinumv = abs(n1%c) END FUNCTION ABSINUMV FUNCTION POWINUM(n1, m) implicit none type(inum), intent(in) :: n1 integer, intent(in) :: m type(inum) :: powinum if (aimag(n1%c)
 ... ... @@ -3,6 +3,9 @@ MODULE maths_functions use globals use utils implicit none interface polylog module procedure polylog1, polylog2 end interface polylog CONTAINS FUNCTION zeta(n) ... ... @@ -297,11 +300,11 @@ CONTAINS END FUNCTION RECURSIVE FUNCTION polylog(m,x) result(res) RECURSIVE FUNCTION polylog1(m,x) result(res) ! computes the polylog integer :: m type(inum) :: x type(inum) :: x, inv complex(kind=prec) :: res if(verb >= 70) print*, 'called polylog(',m,',',x%c,x%i0,')' ... ... @@ -312,8 +315,9 @@ CONTAINS res = -(1. - 2.**(1-m))*zeta(m) return else if (abs(x) .gt. 1) then res = (-1)**(m-1)*polylog(m,ione/x) & - cmplx(0,2*pi)**m * bernoulli_polynomial(m, 0.5-cmplx(0.,1.)*log(neg(x))/2/pi) / factorial(m) inv = inum(1./x%c, x%i0) res = (-1)**(m-1)*polylog(m,inv) & - cmplx(0,2*pi)**m * bernoulli_polynomial(m, 0.5-cmplx(0.,1.)*conjg(log(-x%c))/2/pi) / factorial(m) return endif ... ... @@ -324,7 +328,26 @@ CONTAINS else res = naive_polylog(m,x%c) end if END FUNCTION polylog END FUNCTION polylog1 RECURSIVE FUNCTION polylog2(m,x,y) result(res) type(inum) :: x, y integer m complex(kind=prec) :: res res=polylog1(m,inum(x%c/y%c,di0)) END FUNCTION POLYLOG2 FUNCTION PLOG1(a,b) ! calculates log(1-a/b) implicit none type(inum) :: a,b complex(kind=prec) plog1 plog1 = log(1.-a%c/b%c) END FUNCTION END MODULE maths_functions ... ...
 !!!!!!!!!!!!!!!!!!!!!! MODULE FUNCTIONS !!!!!!!!!!!!!!!!!!!!!! USE GLOBAL_DEF IMPLICIT NONE TYPE MLM ! MassLess Momentum real(kind=prec), dimension(4) :: momentum complex(kind=prec), dimension(2) :: sp ! <+| complex(kind=prec), dimension(2) :: sm ! <-| complex(kind=prec), dimension(2) :: ep ! |+> complex(kind=prec), dimension(2) :: em ! |-> END TYPE MLM TYPE PARTON character(len=1) :: flavour ! choices for flavour: ! "j": light quarks or gluons ! "b": bottom ! "B": antibottom ! "F": fat, bottom + antibottom real (kind=prec), dimension(4) :: momentum END TYPE PARTON INTERFACE OPERATOR(+) module procedure add_parton END INTERFACE INTERFACE OPERATOR(-) module procedure minus_mlm END INTERFACE INTERFACE A module procedure Aversion2, Aversion3, Aversion4, Aversion5 END INTERFACE INTERFACE B module procedure Bversion2, Bversion3, Bversion4, Bversion5 END INTERFACE INTERFACE S module procedure Sversion1, Sversion2, Sversion3, Sversion4 END INTERFACE INTERFACE SQ module procedure SQversion1 END INTERFACE INTERFACE EIK module procedure eik00, eik0M, eikM0, eikMM END INTERFACE INTERFACE IREG module procedure ireg_00, ireg_0M, ireg_M0, ireg_MM, ireg_SELF END INTERFACE INTERFACE ISIN module procedure isin_00, isin_0m, isin_m0, isin_MM, isin_SELF END INTERFACE INTERFACE ILIN module procedure ilin_MM, ilin_SELF END INTERFACE INTERFACE DISCB module procedure discbSMN, discbSMM END INTERFACE INTERFACE SCALARC0IR6 module procedure scalarc0ir6SMM, scalarc0ir6SMN END INTERFACE CONTAINS FUNCTION MAKE_MLM(qq) real(kind=prec), intent(in) :: qq(4) type(mlm) :: make_mlm complex(kind=prec) :: sq if(qq(4) > 0._prec) then sq = sqrt(qq(4) + qq(1) + Tiny(qq(4))) else sq = imag*sqrt(-qq(4) - qq(1) + Tiny(qq(4))) endif make_mlm%momentum = qq make_mlm%sp = (/ (-imag*qq(2)-qq(3))/sq, sq /) make_mlm%sm = (/ (-imag*qq(2)+qq(3))/sq, sq /) make_mlm%ep = (/ sq, (imag*qq(2)-qq(3))/sq /) make_mlm%em = (/ sq, (imag*qq(2)+qq(3))/sq /) END FUNCTION MAKE_MLM FUNCTION MINUS_MLM(ki) type(mlm) :: minus_mlm type(mlm), intent(in) :: ki complex(kind=prec) :: factor minus_mlm%momentum = -ki%momentum if(ki%momentum(4) > 0._prec) then factor = imag else factor = -imag endif minus_mlm%sp = factor*ki%sp minus_mlm%sm = factor*ki%sm minus_mlm%ep = factor*ki%ep minus_mlm%em = factor*ki%em END FUNCTION MINUS_MLM FUNCTION ADD_PARTON(par1,par2) type(parton) :: add_parton type(parton), intent(in) :: par1, par2 character(len=1) :: fl1,fl2 fl1 = par1%flavour fl2 = par2%flavour select case(fl1) case("j") add_parton%flavour = fl2 case("b") select case(fl2) case("j") add_parton%flavour = fl1 case("B") add_parton%flavour = "F" case default call crash("ADD_PARTON_1") end select case("B") select case(fl2) case("j") add_parton%flavour = fl1 case("b") add_parton%flavour = "F" case default call crash("ADD_PARTON_2") end select case("F") select case(fl2) case("j") add_parton%flavour = fl1 case default call crash("ADD_PARTON_3") end select case default call crash("ADD_PARTON_0") end select add_parton%momentum = par1%momentum + par2%momentum END FUNCTION ADD_PARTON FUNCTION Sversion1(k1,k2) ! S(k1,k2) = 2 k1.k2 type(mlm), intent(in) :: k1,k2 real (kind=prec) :: Sversion1,dot_dot,q1(4),q2(4) q1 = k1%momentum q2 = k2%momentum dot_dot = q1(4)*q2(4) - q1(1)*q2(1) - q1(2)*q2(2) - q1(3)*q2(3) Sversion1 = 2*dot_dot END FUNCTION Sversion1 FUNCTION Sversion2(q1,q2) ! S(q1,q2) = 2 q1.q2 real (kind=prec), intent(in) :: q1(4),q2(4) real (kind=prec) :: Sversion2,dot_dot dot_dot = q1(4)*q2(4) - q1(1)*q2(1) - q1(2)*q2(2) - q1(3)*q2(3) Sversion2 = 2*dot_dot END FUNCTION Sversion2 FUNCTION Sversion3(k1,q2) ! S(q1,q2) = 2 q1.q2 type(mlm), intent(in) :: k1 real (kind=prec), intent(in) :: q2(4) real (kind=prec) :: q1(4) real (kind=prec) :: Sversion3,dot_dot q1 = k1%momentum dot_dot = q1(4)*q2(4) - q1(1)*q2(1) - q1(2)*q2(2) - q1(3)*q2(3) Sversion3 = 2*dot_dot END FUNCTION Sversion3 FUNCTION Sversion4(q1,k2) ! S(q1,q2) = 2 q1.q2 type(mlm), intent(in) :: k2 real (kind=prec), intent(in) :: q1(4) real (kind=prec) :: q2(4) real (kind=prec) :: Sversion4,dot_dot q2 = k2%momentum dot_dot = q1(4)*q2(4) - q1(1)*q2(1) - q1(2)*q2(2) - q1(3)*q2(3) Sversion4 = 2*dot_dot END FUNCTION Sversion4 FUNCTION SQversion1(q1) ! SQ(q1) = (q1)^2 real (kind=prec), intent(in) :: q1(4) real (kind=prec) :: SQversion1 SQversion1 = q1(4)**2 - q1(1)**2 - q1(2)**2 - q1(3)**2 END FUNCTION