c calculates six dimensional sums for three loop master diagrams. c copyright Richard Easther, Gerry Guralnik and Stephen Hahn (1999) c c You are free to use, copy and modify this code. c However, please ask us before distributing a modified version. c c for more information see our paper, hep-ph/9912255 c c the main block below calls the function "evaluate(gnum)", c which evaluates the F_n c c this code is a little rough, and is not intended to be a finished c product. Hopefully the comments provide some guidance! c c h is the parameter that fixes the accuracy of the Sinc function c representation. c the truncation of the individual sums is determined by the values c of tolrel and tolabs. Making these values smaller includes more terms c in the sums, and minimizes truncation error. For smaller values of h, c these numbers need to be reduced. c real*8 answer,evaluate,hvals,tolvals integer gnum include 'master.h' tolrel = 1d-4 open(20,file='results.dat') c c these loops evaluate the 9 non-trivial F_n for a series of values of c h and tolabs. c c do hvals=1d0,.69999d0,-.05d0 do tolvals=5,5 h=hvals tolabs=10d0**(-tolvals) do gnum = 2,10 answer = evaluate(gnum) write(20,*) gnum,h,TOLABS,answer write(6,*) gnum,h,TOLABS,answer enddo enddo enddo end c --------------------------------------------------------------- c c evaluate calculates F_n (defined in paper); parameter values are c passed via common blocks defined in master.h c real*8 function evaluate(graph) include 'master.h' integer finished,clicker,graph real*8 result,factor,thrown integer mxlocate,kmin,kmax,midpt logical OuterCompQ real*8 diag,aitken,aitken2 c k: indices on sums (ie {k1,k2,k3,...,kn} c direct: direction each sum is master.h (+/- 1) c midpoint: location of maximum term at each level c startvals: value of kN to begin with c inner: boolean for loop c F_1 is trivial; so we check for this in case it gets evaluated by c accident. if(graph.eq.1) then evaluate=0d0 return endif c this is the prefactor in the front of the sum. factor = h**6 c we fill an array with the values of masses for each propagator c in the diagram V_n, and in V_1, before forming F_n = V_1-V_n call massfix(graph) c we precomputed the p(k) and c(k) since their valeus depend only c on k. call initialize(0) c other variables initialized here. finished= 0 clicker =1 result=diag(1) kmin=-10 kmax=10 midpt = startvals(PROPS) c this is the main loop do while(finished.eq.0) c the counter being iterated is the innermost loop, so set depth c to the number of props. Then begin the loop for the innermost sum. c diag( ) computes the general term, as a function of the k(i). c diag(0) assumes that k_1...k_5 have not changed since the last call, c whereas diag(1) "initializes" the function. The call here is designed c to ensure the diag in properly initialized before doing the sum c over k_6. This *is* a little crude! thrown = diag(1) k(PROPS) = startvals(PROPS) c endpoints determines the values of k_6 the sum needs to run between c in order to get the desired level of accuracy. call endpoints(kmin,kmax,midpt,result) c depth fixes the sum we are doing (ie the sum over k(depth) depth = PROPS c direct(i) determines the direction of the i-th sum (since the sums c run in both positive and negative directions) direct(PROPS)=-1 do loop1 = midpt,kmin,-1 k(6)=loop1 partial(loop1,PROPS) = diag(0) running(PROPS) =running(PROPS)+partial(loop1,PROPS) enddo c if possible, the result is "polished" using Aitken extrapolation at c second order. The test here is to make sure we are not trying to c "improve" a result before k_6 has become large enough for the general c term to be dominated by its asymptotic behaviour. if(abs(running(PROPS)/aitken(PROPS)-1d0) & .lt.compabs(PROPS)*100d0)running(PROPS)= aitken2(PROPS) c now we do the sum in the positive k direction. direct(PROPS)=1 do loop1 = midpt+1,kmax k(6)=loop1 partial(loop1,PROPS) = diag(0) running(PROPS) =running(PROPS)+partial(loop1,PROPS) enddo if(abs(running(PROPS)/aitken(PROPS)-1d0) & .lt.compabs(PROPS)*100d0) running(PROPS)=aitken2(PROPS) c end of inner loop. c result keeps a running total of the overall sum. result = result+running(depth) c now we look at the k_1 -- k_5 sums. clicker =1 do while((clicker.eq.1).and.(depth.gt.1)) depth=depth-1 partial(k(depth),depth) = running(depth+1) running(depth) = running(depth) + running(depth+1) running(depth+1) = 0d0 c this produces some output to the screen. It is useful for debugging, c or to keep an eye on progress if you are impatient. if(depth.lt.3) then write(6,*) running(1)*factor write(6,*) k(1),k(2),k(3),k(4),k(5),k(6) write(6,*) running(depth)*factor,aitken2(depth)*factor, & partial(k(depth),depth) write(6,*) endif c if terms no longer contribute to sum, either turn round or finish c this tests to see whether the k(depth) sum can be truncated. c OuterCompQ returns .TRUE. if the summation should continue. if(OuterCompQ(result).or. & (abs(k(depth)-startvals(depth)).lt.4)) then k(depth) = k(depth) + direct(depth) clicker =0 else c check to see whether we have finished the overall sum, or have just c finished the summation for increasing k. if(direct(depth).eq.1) then c we only do the delta^2 improvement if depth > 1, since the c general term for increasing k_1 does not have the right form c to be helped by the extrapolation. if(abs(running(depth)/aitken(depth)-1d0).lt.1d-2)then if(depth.ge.1) running(depth) = aitken2(depth) endif direct(depth)= -1 endpos(depth) = k(depth) k(depth) = startvals(depth)-1 clicker=0 else c this block is used if the current sum (at depth) is completed. if(abs(running(depth)/aitken(depth)-1d0).lt.1d-2) & running(depth) =aitken2(depth) direct(depth)=1 endneg(depth)=k(depth) c this is only true if the whole sum is completed. if(depth.eq.1) then finished=1 else c fixes the starting value of k(depth) be the value that gave c the largest general term on the previous pass. startvals(depth)= mxlocate(startvals(depth)) k(depth)=startvals(depth) endif endif endif enddo enddo c scale final result by h^6 and return! evaluate= running(1) * factor return end c this does a simple aithken extrapolation (see numerical recipes) real*8 function aitken(deep) integer deep include 'master.h' aitken = running(deep)+partial(k(deep),deep)**2 / & ( partial(k(deep)-direct(deep),deep) & - partial(k(deep),deep)) return end c calculates second order aitken extrapolation. real*8 function aitken2(deep) real*8 aitken,num,dem integer deep include 'master.h' num = partial(k(deep),deep) - & & partial(k(deep),deep)**2/(partial(k(deep),deep)- & partial(k(deep)-direct(deep),deep)) + & & partial(k(deep)-direct(deep),deep)**2/ & (partial(k(deep)-direct(deep),deep)- & partial(k(deep)-2*direct(deep),deep)) dem = partial(k(deep),deep) - & partial(k(deep)-direct(deep),deep) - & & partial(k(deep),deep)**2/(partial(k(deep),deep)- & partial(k(deep)-direct(deep),deep)) + & & 2d0* partial(k(deep)-direct(deep),deep)**2/ & (partial(k(deep)-direct(deep),deep)- & partial(k(deep)-2*direct(deep),deep)) - & & partial(k(deep)-2*direct(deep),deep)**2/ & (partial(k(deep)-2*direct(deep),deep)- & partial(k(deep)-3*direct(deep),deep)) if(abs(dem).lt.TINY) then aitken2 = aitken(deep) else aitken2 = running(deep) - & partial(k(deep),deep)**2/(partial(k(deep),deep)- & partial(k(deep)-direct(deep),deep)) - & num*num/dem endif return end c tests to see if an outerloop is completed. c c this is a heuristic (and messy) approach. The tests are the last term c added to the i-th sum is small compared to c c the running total for the i-th sum AND c the running total over the whole sum AND c that the terms are decreasing. c c c moreover, it also checks that the value of k is large enough so the c aitken extrapolation will actually help. c c c OuterCompQ is .TRUE. if the summation is not completed. c c logical function OuterCompQ(result) include 'master.h' real*8 result,hold,aitken hold =aitken(depth) OuterCompQ = & ((partial(k(depth),depth)/result.gt.compabs(depth)).or. & (partial(k(depth),depth).gt. & partial(k(depth)-direct(depth),depth)).or. & (partial(k(depth),depth)/running(depth).gt.comprel(depth))) if((depth.gt.1).OR.(direct(depth).eq.-1)) then if((.not.OuterCompQ).and. & (abs(1d0 -lastaitken(depth)/hold).gt.compaitken(depth))) & OuterCompQ= .TRUE. endif lastaitken(depth)=hold return end c find endpoints and midpoint for innermost summation. subroutine endpoints(kmin,kmax,midpt,result) include 'master.h' integer lowend,midpt,hiend,kmax,kmin,ktest1,ktest2,STEP real*8 result real*8 dkmax,dkmin,dkmid,target real*8 diag parameter(STEP=10) c start by finding maximum value of diag between max and min k(PROPS) = kmax dkmax = diag(1) k(PROPS) = kmin dkmin = diag(1) k(PROPS) = midpt dkmid = diag(1) c make sure we do actually bracket a maximum. do while(dkmid.lt.dkmin) kmin = kmin-STEP k(PROPS) = kmin dkmin = diag(1) enddo do while(dkmid.lt.dkmax) kmax = kmax+STEP k(PROPS) = kmax dkmax = diag(1) enddo c now we know that we have bracketed the maximum, we can find it. hiend = kmax lowend = kmin do while((hiend-lowend).gt.4) k(PROPS) = (midpt+lowend)/2 test = diag(1) if(test.gt.dkmid) then hiend=midpt midpt = (lowend+hiend)/2 else k(PROPS) = (midpt+hiend)/2 test = diag(1) if(test.gt.dkmid) then lowend=midpt midpt= (lowend+hiend)/2 else lowend= (lowend+midpt)/2 hiend= (midpt+hiend)/2 endif endif k(PROPS) = midpt dkmid = diag(1) enddo k(PROPS) = kmax dkmax = diag(1) k(PROPS) = kmin dkmin = diag(1) c The next step is to locate the values of k that correspond to the c minimum terms we want to include. We solve for these by bisection c too. For plumbing reasons it is easier to have the bisection routines c included here, rather than spinning them off as separate subroutines. target= min(result*compabs(PROPS),dkmid*comprel(PROPS)) do while(dkmin.gt.target) kmin = kmin - STEP k(PROPS) = kmin dkmin = diag(1) enddo do while(dkmax.gt.target) kmax = kmax+ STEP k(PROPS) = kmax dkmax = diag(1) enddo ktest2 = midpt ktest1= (kmin+ktest2)/2 do while(abs(ktest1-ktest2).gt.3) k(PROPS) =ktest1 if(diag(1)-target.lt.0d0)then kmin = ktest1 else ktest2=ktest1 endif ktest1= (kmin+ktest2)/2 enddo ktest2 = midpt ktest1= (kmax+ktest2)/2 do while(abs(ktest1-ktest2).gt.3) k(PROPS) =ktest1 if(diag(1)-target.lt.0d0)then kmax = ktest1 else ktest2=ktest1 endif ktest1= (kmax+ktest2)/2 enddo return end c this function finds the value of k_i for which made the largest c contribution to the i-th sum, to use as a starting point for the c next summation. integer function mxlocate(midpt) c c This function locates the maximal value of partial(k_n,depth) c in terms of k. It provides the starting value of k for the c next trip through the loop, thereby optimizing performance c c This is very crude, but it works. c Replace with parabolic version later, for better convergence? c include 'master.h' integer lowend,midpt,hiend lowend=endneg(depth) hiend=endpos(depth) c if endpoints are close, just split the difference and return. if((hiend-lowend).lt.8) then mxlocate= (hiend+lowend)/2 return endif c check to make sure lowend,midpt and hiend really do bracket a maximum c assumes does not make a trip through the while loop if c test is false to begin with. May be compiler dependent. do while((partial(lowend,depth).gt.partial(midpt,depth)).or. & (partial(hiend,depth).gt.partial(midpt,depth))) if(partial(lowend,depth).gt.partial(hiend,depth)) then hiend = midpt midpt = (lowend+hiend)/2 else lowend = midpt midpt = (lowend+hiend)/2 endif if((hiend-lowend).lt.4) then mxlocate=(hiend+lowend)/2 return endif enddo c now we know that partial(midpt,depth) is bigger than c partial(lowend,depth) and partial(hiend,depth) do while((hiend-lowend).gt.8) if(partial((midpt+lowend)/2,depth).gt. & partial(midpt,depth)) then hiend=midpt midpt = (lowend+hiend)/2 else if(partial((midpt+hiend)/2,depth).gt. & partial(midpt,depth)) then lowend=midpt midpt= (lowend+hiend)/2 else lowend= (lowend+midpt)/2 hiend= (midpt+hiend)/2 endif endif enddo mxlocate = midpt return end c this calculates the general term in the sum. real*8 function diag(reset) integer reset real*8 share include 'master.h' c if reset =1, one (or more) of k_1 to k_5 is assumed to have changed. if(reset.eq.1)then pre1 = (p(7,k(1))*p(8,k(1)+k(2))*p(9,k(1)+k(3))* - p(10,k(1)+k(4))*p(11,k(1)+k(5))) pre2 = (p(1,k(1))*p(2,k(1)+k(2))*p(3,k(1)+k(3))* - p(4,k(1)+k(4))*p(5,k(1)+k(5))) term1 = (c(2,k(1)+k(2)) + c(5,k(1)+k(5)))*c(3,k(1)+k(3))* - c(4,k(1)+k(4)) + - c(2,k(1)+k(2))*c(5,k(1)+k(5))*(c(3,k(1)+k(3))+ - c(4,k(1)+k(4))) + - c(1,k(1))*(c(2,k(1)+k(2))+c(3,k(1)+k(3)))* - (c(4,k(1)+k(4))+c(5,k(1)+k(5))) term2 = (c(2,k(1)+k(2))+c(4,k(1)+k(4)))*(c(3,k(1)+k(3))+ - c(5,k(1)+k(5)))+ - c(1,k(1))*(c(2,k(1)+k(2))+c(3,k(1)+k(3))+ - c(4,k(1)+k(4))+c(5,k(1)+k(5))) endif c if reset = 0, the previous terms are not recomputed. share = 1d0/((term1+term2*c(6,k(6)))*(term1+term2*c(6,k(6)))) diag = share*(pre1*p(12,k(6))-pre2*p(6,k(6))) return end c this sets the individual masses. subroutine massfix(graph) include 'master.h' integer graph mass(7)=1d0 mass(8)=0d0 mass(9)=0d0 mass(10)=0d0 mass(11)=0d0 mass(12)=0d0 if(graph.eq.2) then ! V_2A mass(1)=1d0 mass(2)=1d0 mass(3)=0d0 mass(4)=0d0 mass(5)=0d0 mass(6)=0d0 return endif if(graph.eq.3) then ! V_2N mass(1)=1d0 mass(2)=0d0 mass(3)=0d0 mass(4)=0d0 mass(5)=0d0 mass(6)=1d0 return endif if(graph.eq.4) then ! V_3T mass(1)=1d0 mass(2)=1d0 mass(3)=0d0 mass(4)=1d0 mass(5)=0d0 mass(6)=0d0 return endif if(graph.eq.5) then ! V_3S mass(1)=1d0 mass(2)=1d0 mass(3)=1d0 mass(4)=0d0 mass(5)=0d0 mass(6)=0d0 return endif if(graph.eq.6) then ! V_3L mass(1)=1d0 mass(2)=1d0 mass(3)=0d0 mass(4)=0d0 mass(5)=1d0 mass(6)=0d0 return endif if(graph.eq.7) then ! V_4A mass(1)=1d0 mass(2)=0d0 mass(3)=1d0 mass(4)=0d0 mass(5)=1d0 mass(6)=1d0 return endif if(graph.eq.8) then ! V_4N mass(1)=1d0 mass(2)=0d0 mass(3)=1d0 mass(4)=1d0 mass(5)=0d0 mass(6)=1d0 return endif if(graph.eq.9) then ! V_5N mass(1)=1d0 mass(2)=1d0 mass(3)=1d0 mass(4)=1d0 mass(5)=1d0 mass(6)=0d0 return endif if(graph.eq.10) then ! V_6 mass(1)=1d0 mass(2)=1d0 mass(3)=1d0 mass(4)=1d0 mass(5)=1d0 mass(6)=1d0 return endif return end c we initialize the arrays here (particularly the p(k) and c(k) subroutine initialize(dummy) integer dummy,loops1,loops2 include 'master.h' c initialize p's and c's do loops1 = 1,PROPS+6 c c(j) = exp(ind*h) do loops2 = MAXNEG,MAXPOS c(loops1,loops2) = exp(dreal(loops2)*h) if(c(loops1,loops2) .ge.0d0) then p(loops1,loops2) = exp(dreal(loops2)*h - . mass(loops1)*mass(loops1)* exp(dreal(loops2)*h)) else p(loops1,loops2) = 1d0 endif enddo enddo c these numbers fix the truncation of the sums. In general the c innermost sums need to be done more accurately the outer ones, c to prevent a build up of numerical error. do loops1 =1,PROPS compabs(loops1) = tolabs/(40d0**(loops1-1)) compaitken(loops1) = tolrel comprel(loops1) = tolrel enddo comprel(props) = tolrel/50d0 c initialize arrays used by code. do loops1 = 1,MAXPROPS running(loops1)= 0d0 direct(loops1)= 1 startvals(loops1)= -1 k(loops1)= startvals(loops1) enddo end