1
* $Id: ffinit.f,v 1.9 1996/04/26 10:39:03 gj Exp $
4
***#[*comment:***********************************************************
5
* calculate a lot of commonly-used constants in the common block *
6
* /ffcnst/. also set the precision, maximum loss of digits and *
7
* the minimum value the logarithm accepts in /prec/. *
8
***#]*comment:***********************************************************
11
integer i,j,init,ioldp(13,12),isgrop(10,12),ji
13
DOUBLE PRECISION s,sold
16
DOUBLE PRECISION delta
19
data ioldp/1,2,3,4, 5,6,7,8,9,10, 11,12,13,
20
+ 4,1,2,3, 8,5,6,7,10,9, 11,13,12,
21
+ 3,4,1,2, 7,8,5,6,9,10, 11,12,13,
22
+ 2,3,4,1, 6,7,8,5,10,9, 11,13,12,
23
+ 4,2,3,1, 10,6,9,8,7,5, 12,11,13,
24
+ 1,3,2,4, 9,6,10,8,5,7, 12,11,13,
25
+ 1,2,4,3, 5,10,7,9,8,6, 13,12,11,
26
+ 1,4,3,2, 8,7,6,5,9,10, 11,13,12,
27
+ 3,4,2,1, 7,10,5,9,6,8, 13,12,11,
28
+ 2,3,1,4, 6,9,8,10,5,7, 12,13,11,
29
+ 4,2,1,3, 10,5,9,7,8,6, 13,11,12,
30
+ 1,3,4,2, 9,7,10,5,8,6, 13,11,12/
32
+ +1,+1,+1,+1, +1,+1,+1,+1, +1,+1,
33
+ +1,+1,+1,+1, +1,+1,+1,+1, -1,+1,
34
+ +1,+1,+1,+1, +1,+1,+1,+1, -1,-1,
35
+ +1,+1,+1,+1, +1,+1,+1,+1, +1,-1,
36
+ +1,+1,+1,+1, -1,+1,+1,-1, +1,-1,
37
+ +1,+1,+1,+1, -1,-1,+1,+1, -1,+1,
38
+ +1,+1,+1,+1, +1,+1,-1,+1, +1,+1,
39
+ +1,+1,+1,+1, -1,-1,-1,-1, +1,-1,
40
+ +1,+1,+1,+1, -1,+1,+1,+1, -1,-1,
41
+ +1,+1,+1,+1, +1,+1,+1,-1, +1,-1,
42
+ +1,+1,+1,+1, -1,+1,+1,-1, -1,-1,
43
+ +1,+1,+1,+1, -1,-1,+1,+1, -1,-1/
46
* check whether tehre is anything to do
47
if ( init .ne. 0 ) return
49
print *,'===================================================='
50
print *,' FF 2.0, a package to evaluate one-loop integrals'
51
print *,'written by G. J. van Oldenborgh, NIKHEF-H, Amsterdam'
52
print *,'===================================================='
53
print *,'for the algorithms used see preprint NIKHEF-H 89/17,'
54
print *,'''New Algorithms for One-loop Integrals'', by G.J. van'
55
print *,'Oldenborgh and J.A.M. Vermaseren, published in '
56
print *,'Zeitschrift fuer Physik C46(1990)425.'
57
print *,'===================================================='
63
* the loss of accuracy in any single subtraction at which
64
* (timeconsuming) corrective action is to be taken is
68
* the precision to which real calculations are done is
75
if ( s .eq. sold ) goto 2
80
* (take three bits for safety)
81
if ( lwrite ) print *,'ffini: precx = ',precx
83
* the precision to which complex calculations are done is
89
cs = exp(log(DCMPLX(1+precc,x0)))
90
if ( DBLE(cs) .eq. sold ) goto 4
95
* (take three bits for safety)
96
if ( lwrite ) print *,'ffini: precc = ',precc
98
* for efficiency tke them equal if they are not too different
100
if ( precx/precc .lt. 4 .and. precx/precc .gt. .25 ) then
101
precx = max(precc,precx)
102
precc = max(precc,precx)
105
* and the minimum value the logarithm accepts without complaining
106
* about arguments zero is (DOUBLE PRECISION cq DOUBLE COMPLEX)
112
if ( 2*abs(s) .ne. xalogm ) goto 6
116
if ( xalogm.eq.0 ) xalogm = 1d-308
117
if ( lwrite ) print *,'ffini: xalogm = ',xalogm
119
xclogm = abs(DCMPLX(s))
122
if ( 2*abs(DCMPLX(s)) .ne. xclogm ) goto 8
123
xclogm = abs(DCMPLX(s))
126
if ( xclogm.eq.0 ) xclogm = 1d-308
127
if ( lwrite ) print *,'ffini: xclogm = ',xclogm
129
* These values are for Absoft, Apollo fortran (68000):
132
* These values are for VAX g_float
135
* These values are for Gould fort (because of div_zz)
138
xalog2 = sqrt(xalogm)
139
xclog2 = sqrt(xclogm)
143
* calculate the coefficients of the series expansion
144
* li2(x) = sum bn*z^n/(n+1)!, z = -log(1-x), bn are the
145
* bernouilli numbers (zero for odd n>1).
147
bf(1) = - 1.D+0/4.D+0
148
bf(2) = + 1.D+0/36.D+0
149
bf(3) = - 1.D+0/36.D+2
150
bf(4) = + 1.D+0/21168.D+1
151
bf(5) = - 1.D+0/108864.D+2
152
bf(6) = + 1.D+0/52690176.D+1
153
bf(7) = - 691.D+0/16999766784.D+3
154
bf(8) = + 1.D+0/1120863744.D+3
155
bf(9) = - 3617.D+0/18140058832896.D+4
156
bf(10) = + 43867.D+0/97072790126247936.D+3
157
bf(11) = - 174611.D+0/168600109166641152.D+5
158
bf(12) = + 77683.D+0/32432530090601152512.D+4
159
bf(13) = - 236364091.D+0/4234560341829359173632.D+7
160
bf(14) = + 657931.D+0/5025632054039239458816.D+6
161
bf(15) = - 3392780147.D+0/109890470493622010006470656.D+7
162
bf(16)=+172.3168255201D+0/2355349904102724211909.3102313472D+6
163
bf(17)=-770.9321041217D+0/4428491985594062112714.2791446528D+8
164
bf(18)=( 0.4157635644614046176D-28)
165
bf(19)=(-0.9962148488284986022D-30)
166
bf(20)=( 0.2394034424896265390D-31)
168
* inverses of integers:
175
* inverses of faculties of integers:
179
xinfac(i) = xinfac(i-1)/i
182
* inx: p(inx(i,j)) = isgn(i,j)*(s(i)-s(j))
218
iold(j,i) = ioldp(j,i)
221
isgrot(j,i) = isgrop(j,i)
249
* isgn5 is not yet used.
296
if ( ji.gt.+3 ) ji = ji - 6
297
if ( ji.lt.-3 ) ji = ji + 6
300
elseif ( abs(ji).eq.3 ) then
306
elseif ( ji.gt.0 ) then
308
elseif ( ji.lt.0 ) then
311
print *,'ffini: internal error in isgn6'
318
* #[ defaults for flags:
321
* the debugging flags.
332
* Specify which root to take in cases were two are possible
333
* it may be advantageous to change this to -1 (debugging hook)
338
* the cutoff has to be initialized because of the memory mechansim
342
* the scheme used for the complex scalar functions:
344
* nschem = 1: do not use the complex mass at all
345
* 2: only use the complex mass in linearly divergent terms
346
* 3: also use the complex mass in divergent logs UNDEFINED
347
* 4: use the complex mass in the C0 if there are
349
* 5: include the almost-divergent threshold terms from
351
* 6: include the (s-m^2)*log(s-m^2) threshold terms from
352
* (m1+m2),m1,m2) vertices
353
* 7: full complex computation
354
* (only in the ffz... functions):
355
* onshel = .FALSE.: use the offshell p^2 everywhere
356
* .TRUE.: use the onshell p^2 except in complex parts
361
* the precision wanted in the complex D0 (and hence E0) when
362
* nschem=7, these are calculated via Taylor exoansion in the real
363
* one and hence expensive.
367
* in some schemes, for onshel=.FALSE.,
368
* when |p^2-Re(m^2)| < nwidth*|Im(m^2)| special action is taken
372
* a flag to indicate the validity of differences smuggled to the
373
* IR routines in the C0 (ff internal only)
377
* #] defaults for flags:
382
***#[*comment:***********************************************************
383
* check a lot of commonly-used constants in the common block *
385
***#]*comment:***********************************************************
393
* calculate the coefficients of the series expansion
394
* li2(x) = sum bn*z^n/(n+1)!, z = -log(1-x), bn are the
395
* bernouilli numbers (zero for odd n>1).
397
if ( bf(1) .ne. - 1.D+0/4.D+0 )
398
+ print *,'ffexi: error: bf(1) is corrupted'
399
if ( bf(2) .ne. + 1.D+0/36.D+0 )
400
+ print *,'ffexi: error: bf(2) is corrupted'
401
if ( bf(3) .ne. - 1.D+0/36.D+2 )
402
+ print *,'ffexi: error: bf(3) is corrupted'
403
if ( bf(4) .ne. + 1.D+0/21168.D+1 )
404
+ print *,'ffexi: error: bf(4) is corrupted'
405
if ( bf(5) .ne. - 1.D+0/108864.D+2 )
406
+ print *,'ffexi: error: bf(5) is corrupted'
407
if ( bf(6) .ne. + 1.D+0/52690176.D+1 )
408
+ print *,'ffexi: error: bf(6) is corrupted'
409
if ( bf(7) .ne. - 691.D+0/16999766784.D+3 )
410
+ print *,'ffexi: error: bf(7) is corrupted'
411
if ( bf(8) .ne. + 1.D+0/1120863744.D+3 )
412
+ print *,'ffexi: error: bf(8) is corrupted'
413
if ( bf(9) .ne. - 3617.D+0/18140058832896.D+4 )
414
+ print *,'ffexi: error: bf(9) is corrupted'
415
if ( bf(10) .ne. + 43867.D+0/97072790126247936.D+3 )
416
+ print *,'ffexi: error: bf(10) is corrupted'
417
if ( bf(11) .ne. - 174611.D+0/168600109166641152.D+5 )
418
+ print *,'ffexi: error: bf(11) is corrupted'
419
if ( bf(12) .ne. + 77683.D+0/32432530090601152512.D+4 )
420
+ print *,'ffexi: error: bf(12) is corrupted'
421
if ( bf(13) .ne. - 236364091.D+0/4234560341829359173632.D+7 )
422
+ print *,'ffexi: error: bf(13) is corrupted'
423
if ( bf(14) .ne. + 657931.D+0/5025632054039239458816.D+6 )
424
+ print *,'ffexi: error: bf(14) is corrupted'
425
if ( bf(15) .ne. -3392780147.D+0/109890470493622010006470656.D+7
426
+ ) print *,'ffexi: error: bf(15) is corrupted'
427
if ( bf(16).ne.+172.3168255201D+0/2355349904102724211909.3102313
429
+ print *,'ffexi: error: bf(16) is corrupted'
430
if ( bf(17).ne.-770.9321041217D+0/4428491985594062112714.2791446
432
+ print *,'ffexi: error: bf(17) is corrupted'
433
if ( bf(18).ne.( 0.4157635644614046176D-28) )
434
+ print *,'ffexi: error: bf(18) is corrupted'
435
if ( bf(19).ne.(-0.9962148488284986022D-30) )
436
+ print *,'ffexi: error: bf(19) is corrupted'
437
if ( bf(20).ne.( 0.2394034424896265390D-31) )
438
+ print *,'ffexi: error: bf(20) is corrupted'
440
* inverses of integers:
443
if ( abs(xninv(i)-x1/i) .gt. precx*xninv(i) ) print *,
444
+ 'ffexi: error: xninv(',i,') is not 1/',i,': ',
445
+ xninv(i),xninv(i)-x1/i
449
* #[ print summary of errors and warning:
452
* #] print summary of errors and warning:
456
subroutine fferr(nerr,ierr)
457
***#[*comment:***********************************************************
459
* generates an errormessage #nerr with severity 2 *
460
* nerr=999 gives a frequency listing of all errors *
462
***#]*comment:***********************************************************
467
integer nerr,ierr,ifile
468
character*80 error(nmax),error1
470
integer noccur(nmax),init,i,ier,inone,nnerr,nomore
471
save error,noccur,init,locwrt,nomore
479
if ( init.eq.0 ) then
483
+ 'fferr: error: illegal value for ierr'
485
call ffopen(ifile,'fferr.dat',ier)
486
if ( ier .ne. 0 ) goto 100
488
read(ifile,'(a)')error1
489
read(ifile,'(a)')error1
491
read(ifile,'(i4,a80)',end=110,err=110)ier,error1
492
if ( ier .lt. 1 .or. ier .gt. nmax ) then
493
print '(a,i3)','fferr: error: wild error number ',
495
print '(a,a)','>>> ',error1
503
+ 'fferr: warning cannot open fferr.dat with error texts'
509
if ( nerr .eq. 999 ) then
510
* print out total numbers...
512
print '(a)','total number of errors and warnings'
513
print '(a)','==================================='
516
if ( noccur(i) .gt. 0 ) then
517
print '(a,i8,a,i3,a,a)','fferr: ',noccur(i),
518
+ ' times ',i,': ',error(i)
523
if ( inone.eq.1 ) print '(a)','fferr: no errors'
525
call ffwarn(999,ierr,x1,x1)
527
print '(a)','the warning system has been disabled'
534
if ( nerr .lt. 1 .or. nerr .gt. nmax ) then
539
noccur(nnerr) = noccur(nnerr) + 1
542
if ( nevent .eq. nomore ) return
545
print '(a,i6,a,i6,a,i8)','fferr: id nr ',id,'/',idsub,
546
+ ', event nr ',nevent
547
print '(a,i6,a,a)','error nr',nerr,': ',error(nnerr)
550
if ( nerr .eq. 100 ) then
551
* we found a matherror - discard all errors from now till next
560
subroutine ffwarn(nerr,ierr,som,xmax)
561
***#[*comment:***********************************************************
563
* The warning routine. A warning is aloss of precision greater *
564
* than xloss (which is default set in ffini), whenever in a *
565
* subtraction the result is smaller than xloss*max(operands) this *
566
* routine is called. Now the strategy is to remember these *
567
* warnings until a 998 message is obtained; then all warnings of *
568
* the previous event are printed. The rationale is that one *
569
* makes this call if too much preciasion is lost only. *
570
* nerr=999 gives a frequency listing of all warnings *
572
* Input: nerr integer the id of the warning message, see the *
573
* file ffwarn.dat or 998 or 999 *
574
* ierr integer the usual error flag: number of digits *
576
* som real the result of the addition *
577
* xmax real the largest operand *
579
* Output: ierr integer is raised by the number of digits lost *
580
* the tolerated loss of xloss *
582
* NOTE: This routine needs a file ffwarn.dat with the warning *
583
* texts, it is very system dependent where to pick it up *
584
* set the PATH variable to your own taste. *
586
***#]*comment:***********************************************************
595
DOUBLE PRECISION som,xmax
600
parameter (memmax = 1000)
601
character*80 warn(nmax),warn1
602
integer noccur(nmax),init,i,ier,inone,nnerr,ilost,
603
+ nermem(memmax),losmem(memmax),idmem(memmax),
604
+ idsmem(memmax),laseve,imem,ifile
605
DOUBLE PRECISION xlosti(nmax),xlost
606
save warn,noccur,init,xlosti,nermem,losmem,idmem,idsmem,
617
if ( init.eq.0 .and. nerr.ne.999 ) then
621
+ 'ffwarn: warning: illegal value for ierr'
624
call ffopen(ifile,'ffwarn.dat',ier)
625
if ( ier.ne.0 ) goto 100
627
read(ifile,'(a)')warn1
628
read(ifile,'(a)')warn1
630
read(ifile,'(i4,a80)',end=110,err=110)ier,warn1
631
if ( warn1.eq.' ' ) goto 90
632
if ( ier.lt.1 .or. ier.gt.nmax ) then
633
print '(a,i3)','ffwarn: error: wild warning number '
635
print '(a,a)','>>> ',warn1
643
+ 'ffwarn: warning cannot open ffwarn.dat with warning texts'
651
if ( nerr.eq.999 ) then
652
* print out total numbers...
655
if ( noccur(i) .gt. 0 ) then
656
print '(a,i8,a,i3,a,a)','ffwarn: ',noccur(i),
657
+ ' times ',i,': ',warn(i)
659
+ ' (lost at most a factor ',xlosti(i),')'
665
if ( inone.eq.1 ) print '(a)','ffwarn: no warnings'
670
if ( nerr .eq. 998 ) then
671
if ( nevent .ne. laseve ) return
673
if ( nermem(i).ne.0 ) then
674
print '(a,i6,a,i6,a,i8)','ffwarn: id nr ',idmem(i),
675
+ '/',idsmem(i),', event nr ',nevent
676
print '(a,i6,a,a)','warning nr ',nermem(i),': ',
678
print '(a,i3,a)',' (lost ',losmem(i),' digits)'
685
* #[ collect warnings:
689
if ( nerr .lt. 1 .or. nerr .gt. nmax ) then
697
noccur(nnerr) = noccur(nnerr) + 1
698
if ( som .ne. 0 ) then
699
xlost = abs(xmax/som)
700
elseif ( xmax .ne. 0 ) then
705
xlosti(nnerr) = max(xlosti(nnerr),xlost)
706
if ( xlost*xloss .gt. xalogm ) then
707
ilost = 1 + int(abs(log10(xlost*xloss)))
713
* nice place to stop when debugging
715
if ( ilost.ge.10 ) then
716
ilost = ilost + 1 - init
721
if ( laseve .ne. nevent ) then
725
if ( imem .le. memmax ) then
733
* print directly if lwrite TRUE
735
if ( awrite .or. lwrite ) then
737
print '(a,i6,a,i6,a,i8)','ffwarn: id nr ',idmem(imem),'/',
738
+ idsmem(imem),', event nr ',nevent
739
print '(a,i6,a,a)','warning nr ',nermem(imem),': ',
741
print '(a,i3,a)',' (lost ',losmem(imem),' digits)'
743
* #] collect warnings:
747
subroutine ffopen(ifile,name,ier)
749
* opens a data file and returns the unit number.
759
character*128 path,fullname
765
inquire(ifile,opened=lopen)
766
if ( .not.lopen ) goto 20
770
* Adjust PATH to suit your own directory structure
771
* I could use a getenv() here, but that may not work
773
* VMS users: use something like the following lines instead
774
* fullname = 'USR$LOCAL[GEERT]'//name
775
* open(ifile,file=fullname,status='OLD',READONLY,err=100)
777
* first try - my home directory
778
path = '/user/gj/lib/'
779
fullname = path(1:index(path,' ')-1)//name
780
open(ifile,file=fullname,status='OLD',err=30)
783
* second try - the system directory
784
path = '/usr/local/ff/'
785
fullname = path(1:index(path,' ')-1)//name
786
open(ifile,file=fullname,status='OLD',err=40)
788
* file could not be found
790
print *,'ffopen: error: could not open ',fullname
791
print *,' adjust path in ffopen (ffinit.f)'
796
DOUBLE PRECISION function ffbnd(n1,n2,array)
797
*************************************************************************
799
* calculate bound = (precx*|a(n1)/a(n1+n2)|^(1/n2) which is the *
800
* maximum value of x in a series expansion sum_(i=n1)^(n1+n2) *
801
* a(i)*x(i) to give a result of accuracy precx (actually of |next *
804
*************************************************************************
807
DOUBLE PRECISION array(n1+n2)
809
if ( array(n1+n2) .eq. 0 ) then
810
print *,'ffbnd: fatal: array not intialized; did you call ',
814
ffbnd = (precx*abs(array(n1)/array(n1+n2)))**(1/DBLE(n2))
818
DOUBLE PRECISION function ffbndc(n1,n2,carray)
819
*************************************************************************
821
* calculate bound = (precc*|a(n1)/a(n1+n2)|^(1/n2) which is the *
822
* maximum value of x in a series expansion sum_(i=n1)^(n1+n2) *
823
* a(i)*x(i) to give a result of accuracy precc (actually of |next *
826
*************************************************************************
829
DOUBLE COMPLEX carray(n1+n2)
831
if ( carray(n1+n2) .eq. 0 ) then
832
print *,'ffbnd: fatal: array not intialized; did you call ',
836
ffbndc = (precc*abs(carray(n1)/carray(n1+n2)))**(1/DBLE(n2))
840
subroutine ffroot(xm,xp,a,b,c,d,ier)
841
***#[*comment:***********************************************************
843
* Calculate the roots of the equation *
844
* a*x^2 - 2*b*x + c = 0 *
846
* x = (b +/- d )/a xp*xm = c/a *
848
***#]*comment:***********************************************************
855
DOUBLE PRECISION xm,xp,a,b,c,d
859
DOUBLE PRECISION s1,s2,s3,rloss
868
if ( b.gt.0 .eqv. d.gt.0 ) then
877
* if ( lwrite ) print *,'ffroot: a,b,c,d = ',a,b,c,d
883
elseif ( b .gt. 0 .eqv. d .gt. 0 ) then
893
rloss = xloss*DBLE(10)**(-2-mod(ier,50))
894
if ( xm .ne. 0 ) then
898
if ( rloss*abs(s1-s2+s3) .gt. precx*max(abs(s1),abs(s2),
900
print *,'ffroot: error: xm not root! ',s1,s2,s3,
904
if ( xp .ne. 0 ) then
908
if ( rloss*abs(s1-s2+s3) .gt. precx*max(abs(s1),abs(s2),
910
print *,'ffroot: error: xp not root! ',s1,s2,s3,
919
subroutine ffcoot(xm,xp,a,b,c,d,ier)
920
***#[*comment:***********************************************************
922
* Calculate the roots of the equation *
923
* a*x^2 - 2*b*x + c = 0 *
925
* x = (b +/- d )/a xp*xm = c/a *
927
***#]*comment:***********************************************************
934
DOUBLE COMPLEX xm,xp,a,b,c,d
938
DOUBLE COMPLEX s1,s2,s3,cc
939
DOUBLE PRECISION absc,rloss
947
absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
952
if ( DBLE(b).gt.0 .eqv. DBLE(d).gt.0 ) then
961
* if ( lwrite ) print *,'ffroot: a,b,c,d = ',a,b,c,d
968
elseif ( absc(cc) .gt. xloss*absc(d) ) then
978
rloss = xloss**2*DBLE(10)**(-mod(ier,50))
979
if ( absc(xm) .gt. xclogm ) then
984
if ( rloss*absc(cc).gt.precc*max(absc(s1),absc(
985
+ s2),absc(s3)) ) print *,
986
+ 'ffcoot: error: xm not root! ',s1,s2,s3,s1-s2+s3
988
if ( absc(xp) .gt. xclogm ) then
993
if ( rloss*absc(cc).gt.precc*max(absc(s1),absc(
994
+ s2),absc(s3)) ) print *,
995
+ 'ffcoot: error: xp not root! ',s1,s2,s3,s1-s2+s3
1002
subroutine ffxhck(xpi,dpipj,ns,ier)
1003
***#[*comment:***********************************************************
1005
* check whether the differences dpipj are compatible with xpi *
1007
***#]*comment:***********************************************************
1011
DOUBLE PRECISION xpi(ns),dpipj(ns,ns)
1013
DOUBLE PRECISION xheck,rloss
1017
if ( ier.lt.0 ) then
1018
print *,'ffxhck: error: ier < 0 ',ier
1021
rloss = xloss**2*DBLE(10)**(-mod(ier,50))
1024
xheck = dpipj(j,i) - xpi(j) + xpi(i)
1025
if ( rloss*abs(xheck) .gt. precx*max(abs(dpipj(j,i)),
1026
+ abs(xpi(j)),abs(xpi(i))) ) then
1027
print *,'ffxhck: error: dpipj(',j,i,') <> xpi(',j,
1028
+ ') - xpi(',i,'):',dpipj(j,i),xpi(j),xpi(i),
1030
if ( lwrite ) ier = ier + 100
1038
subroutine ffchck(cpi,cdpipj,ns,ier)
1039
***#[*comment:***********************************************************
1041
* check whether the differences cdpipj are compatible with cpi *
1043
***#]*comment:***********************************************************
1047
DOUBLE COMPLEX cpi(ns),cdpipj(ns,ns),c
1049
DOUBLE COMPLEX check
1050
DOUBLE PRECISION absc,rloss
1052
absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
1055
if ( ier.lt.0 ) then
1056
print *,'ffchck: error: ier < 0 ',ier
1059
rloss = xloss**2*DBLE(10)**(-mod(ier,50))
1062
check = cdpipj(j,i) - cpi(j) + cpi(i)
1063
if ( rloss*absc(check) .gt. precc*max(absc(
1064
+ cdpipj(j,i)),absc(cpi(j)),absc(cpi(i))) ) then
1065
print *,'ffchck: error: cdpipj(',j,i,') <> cpi(',j,
1066
+ ') - cpi(',i,'):',cdpipj(j,i),cpi(j),cpi(i),
1068
if ( lwrite ) ier = ier + 100
1076
integer function nffeta(ca,cb,ier)
1077
***#[*comment:***********************************************************
1080
* eta(a,b)/(2*i*pi) = ( thIm(-a)*thIm(-b)*thIm(a*b) *
1081
* - thIm(a)*thIm(b)*thIm(-a*b) ) *
1083
* with thIm(a) = theta(Im(a)) *
1084
***#]*comment:***********************************************************
1088
DOUBLE COMPLEX ca,cb
1089
DOUBLE PRECISION a,b,ab,rab
1095
if ( a*b .lt. 0 ) then
1099
rab = DBLE(ca)*DBLE(cb) - a*b
1100
ab = DBLE(ca)*b + a*DBLE(cb)
1101
if ( abs(ab) .lt. precc*abs(DBLE(ca)*b) ) then
1103
if ( lwrite ) print *,'a,b = ',ca,cb,
1104
+ ' (no precision left in DIMAG(ab)=',ab,')'
1106
if ( a .lt. 0 .and. b .lt. 0 .and. ab .gt. 0 ) then
1108
elseif ( a .gt. 0 .and. b .gt. 0 .and. ab .lt. 0 ) then
1110
elseif ( a .eq. 0 .and. DBLE(ca) .le. 0 .or.
1111
+ b .eq. 0 .and. DBLE(cb) .le. 0 .or.
1112
+ ab .eq. 0 .and. rab .le. 0 ) then
1114
if ( ltest .or. lwrite ) print *,'a,b = ',ca,cb
1123
integer function nffet1(ca,cb,cc,ier)
1124
***#[*comment:***********************************************************
1125
* calculates the same eta with three input variables *
1127
* et1(a,b)/(2*i*pi) = ( thIm(-a)*thIm(-b)*thIm(c) *
1128
* - thIm(a)*thIm(b)*thIm(-c) ) *
1130
* with thIm(a) = theta(Im(a)) *
1131
***#]*comment:***********************************************************
1135
DOUBLE COMPLEX ca,cb,cc,c
1136
DOUBLE PRECISION a,b,ab,abp,absc
1138
absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
1141
if ( ltest .and. DIMAG(ca)*DIMAG(cb) .gt. 0 .and. DBLE(ca)*DBLE(
1145
if ( xloss*abs(abp) .lt. precc*absc(ca)*absc(cb) )
1147
if ( ab .gt. 0 .and. abp .lt. 0 .or. ab .lt. 0 .and. abp
1149
print *,'nffet1: error: sgn im(ca*cb) != sgn im(cc): ',
1157
if ( a .gt. 0 .neqv. b .gt. 0 ) then
1162
if ( a .lt. 0 .and. b .lt. 0 .and. ab .gt. 0 ) then
1164
elseif ( a .gt. 0 .and. b .gt. 0 .and. ab .lt. 0 ) then
1166
elseif ( a .eq. 0 .and. DBLE(ca) .le. 0 .or.
1167
+ b .eq. 0 .and. DBLE(cb) .le. 0 .or.
1168
+ ab .eq. 0 .and. DBLE(cc) .le. 0 ) then
1170
if ( ltest.or.lwrite ) print *,'a,b,ab = ',ca,cb,cc
1179
subroutine ffcayl(cs,z,coeff,n,ier)
1180
***#[*comment:***********************************************************
1182
* Do a Taylor expansion in z with real coefficients coeff(i) *
1184
* Input: z complex *
1188
* Output cs complex \sum_{i=1} z^i coeff(i) *
1190
***#]*comment:***********************************************************
1197
DOUBLE PRECISION coeff(n)
1203
DOUBLE PRECISION absc
1204
DOUBLE COMPLEX c,zi,csi
1210
* statement function
1212
absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
1216
cs = z*DBLE(coeff(1))
1217
if ( absc(z) .lt. precc ) return
1221
csi = zi*DBLE(coeff(i))
1223
if ( absc(csi) .lt. precc*absc(cs) ) goto 20
1225
call ffwarn(9,ier,precc,absc(csi))
1231
subroutine fftayl(s,z,coeff,n,ier)
1232
***#[*comment:***********************************************************
1234
* Do a Taylor expansion in z with real coefficients coeff(i) *
1240
* Output cs real \sum_{i=1} z^i coeff(i) *
1242
***#]*comment:***********************************************************
1249
DOUBLE PRECISION coeff(n),z,s
1254
DOUBLE PRECISION zi,si
1263
if ( abs(z) .lt. precx ) return
1269
if ( abs(si) .lt. precx*abs(s) ) goto 20
1271
call ffwarn(9,ier,precx,si)
added
removed
vendor/IREGI/src/qcdloop/ff/ffinit.f
