← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to vendor/CutTools/src/qcdloop/ffcxs4.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/CutTools/src/qcdloop/ffcxs4.f
 
1
*       $Id: ffcxs4.f,v 1.3 1995/10/17 06:55:09 gj Exp $
 
2
*       $Log: ffcxs4.f,v $
 
3
c Revision 1.3  1995/10/17  06:55:09  gj
 
4
c Fixed ieps error in ffdcrr (ffcxs4.f), added real case in ffcrr, debugging
 
5
c info in ffxd0, and warned against remaining errors for del2=0 in ffrot4
 
6
c (ffxd0h.f)
 
7
c
 
8
c Revision 1.2  1995/10/06  09:17:22  gj
 
9
c Found stupid typo in ffxc0p which caused the result to be off by pi^2/3 in
 
10
c some equal-mass cases.  Added checks to ffcxs4.f ffcrr.f.
 
11
c
 
12
*###[ ffcxs4:
 
13
        subroutine ffcxs4(cs3,ipi12,w,y,z,dwy,dwz,dyz,d2yww,d2yzz,
 
14
     +                                  xpi,piDpj,ii,ns,isoort,ier)
 
15
***#[*comment:***********************************************************
 
16
*                                                                       *
 
17
*       Calculate the 8 Spence functions = 4 R's = 2 dR's               *
 
18
*                                                                       *
 
19
*                                                                       *
 
20
*                                                                       *
 
21
***#]*comment:*********************************************************** 
 
22
*  #[ declarations:
 
23
        implicit none
 
24
*
 
25
*       arguments:
 
26
*
 
27
        integer ipi12(4),ii,ns,isoort(4),ier
 
28
        DOUBLE COMPLEX cs3(40)
 
29
        DOUBLE PRECISION w(4),y(4),z(4),dwy(2,2),dwz(2,2),dyz(2,2),
 
30
     +          d2yww,d2yzz,xpi(ns),piDpj(ns,ns),x00(3)
 
31
*
 
32
*       local variables:
 
33
*
 
34
        integer iepz(2),iepw(2)
 
35
        logical ld2yzz,ld2yww
 
36
*       DOUBLE COMPLEX c
 
37
*       DOUBLE PRECISION absc
 
38
*
 
39
*       common blocks
 
40
*
 
41
        include 'ff.h'
 
42
*       absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
43
*  #] declarations: 
 
44
*  #[ groundwork:
 
45
        if ( ltest .and. ns .ne. 6 )
 
46
     +      print *,'ffcxs4: error: only for ns=6, not ',ns
 
47
        if ( isoort(2) .eq. 0 ) then
 
48
            ld2yzz = .FALSE.
 
49
        else
 
50
            ld2yzz = .TRUE.
 
51
        endif
 
52
        if ( isoort(4) .eq. 0 ) then
 
53
            ld2yww = .FALSE.
 
54
        else
 
55
            ld2yww = .TRUE.
 
56
        endif
 
57
        if ( isoort(2) .ne. 0 ) then
 
58
            if ( z(2) .gt. z(1) .eqv. xpi(ii+3) .gt. 0 ) then
 
59
                iepz(1) = +1
 
60
                iepz(2) = -1
 
61
            else
 
62
                iepz(1) = -1
 
63
                iepz(2) = +1
 
64
            endif
 
65
        else
 
66
            print *,'ffcxs4: error: untested algorithm'
 
67
            if ( piDpj(ii,ii+3) .gt. 0 ) then
 
68
                iepz(1) = +1
 
69
            else
 
70
                iepz(1) = -1
 
71
            endif
 
72
        endif
 
73
        if ( isoort(4) .ne. 0 ) then
 
74
            if ( w(2) .gt. w(1) .eqv. xpi(5) .gt. 0 ) then
 
75
                iepw(1) = 1
 
76
                iepw(2) = -1
 
77
            else
 
78
                iepw(1) = -1
 
79
                iepw(2) = 1
 
80
            endif
 
81
        else
 
82
            print *,'ffcxs4: error: untested algorithm'
 
83
            if ( piDpj(2,5) .gt. 0 ) then
 
84
                iepw(1) = +1
 
85
            else
 
86
                iepw(1) = -1
 
87
            endif
 
88
        endif
 
89
*  #] groundwork: 
 
90
*  #[ zm and wp:
 
91
        if ( isoort(4) .eq. 0 ) then
 
92
            if (lwrite) print *,'ffcxs4: to ffcxr(zm)'
 
93
            call ffcxr(cs3(1),ipi12(1),y(2),y(4),z(1),z(3),dyz(2,1),
 
94
     +          ld2yzz,d2yzz,z(2),z(4),.FALSE.,x00,iepz(1),ier)
 
95
        else
 
96
            if (lwrite) print *,'ffcxs4: to ffdcxr(zm,wp)'
 
97
            if ( .not. ( dwz(2,1).eq.0 .and. iepz(1).eq.iepw(2) ) )
 
98
     +          call ffdcxr(cs3( 1),ipi12(1),y(2),y(4),z(1),z(3),
 
99
     +                  z(2),z(4),d2yzz,w(2),w(4),w(1),w(3),d2yww,
 
100
     +                  dyz(2,1),dwy(2,2),dwz(2,1),iepz(1),iepw(2),ier)
 
101
        endif
 
102
*  #] zm and wp: 
 
103
*  #[ zp and wm:
 
104
        if ( isoort(2) .eq. 0 ) then
 
105
            if (lwrite) print *,'ffcxs4: to ffcxr(wm)'
 
106
            call ffcxr(cs3(1),ipi12(1),y(2),y(4),w(1),w(3),-dwy(1,2),
 
107
     +          ld2yww,d2yww,w(2),w(4),.FALSE.,x00,iepw(1),ier)
 
108
        else
 
109
            if (lwrite) print *,'ffcxs4: to ffdcxr(zp,wm)'
 
110
            if ( .not. ( dwz(1,2).eq.0 .and. iepz(2).eq.iepw(1) ) )
 
111
     +          call ffdcxr(cs3(21),ipi12(3),y(2),y(4),z(2),z(4),
 
112
     +                  z(1),z(3),d2yzz,w(1),w(3),w(2),w(4),d2yww,
 
113
     +                  dyz(2,2),dwy(1,2),dwz(1,2),iepz(2),iepw(1),ier)
 
114
        endif
 
115
*  #] zp and wm: 
 
116
*###] ffcxs4: 
 
117
        end
 
118
*###[ ffcs4:
 
119
        subroutine ffcs4(cs3,ipi12,cw,cy,cz,cdwy,cdwz,cdyz,cd2yww,cd2yzz
 
120
     +                  ,cpi,cpiDpj,cp2p,cetami,ii,ns,isoort,ier)
 
121
***#[*comment:***********************************************************
 
122
*                                                                       *
 
123
*       Calculate the 8 Spence functions = 4 R's = 2 dR's               *
 
124
*                                                                       *
 
125
*                                                                       *
 
126
*                                                                       *
 
127
***#]*comment:*********************************************************** 
 
128
*  #[ declarations:
 
129
        implicit none
 
130
*
 
131
*       arguments:
 
132
*
 
133
        integer ipi12(4),ii,ns,isoort(4),ier
 
134
        DOUBLE COMPLEX cs3(40)
 
135
        DOUBLE COMPLEX cw(4),cy(4),cz(4),cdwy(2,2),cdwz(2,2),cdyz(2,2)
 
136
        DOUBLE COMPLEX cd2yww,cd2yzz,cpi(ns),cp2p,cpiDpj(ns,ns),
 
137
     +          cetami(6)
 
138
*
 
139
*       local variables:
 
140
*
 
141
        logical ld2yzz,ld2yww
 
142
        integer i,j,ip,iepz(2),iepw(2),nz(4),nw(4),ntot,i2pi
 
143
        DOUBLE COMPLEX c,cc,clogy,c2y1,cdyw(2,2)
 
144
        DOUBLE COMPLEX zfflo1,zfflog
 
145
        DOUBLE PRECISION absc
 
146
*
 
147
*       common blocks
 
148
*
 
149
        include 'ff.h'
 
150
*
 
151
*       statement function
 
152
*
 
153
        absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
154
*  #] declarations: 
 
155
*  #[ get counters:
 
156
        if ( ltest ) then
 
157
            if ( ns .ne. 6 ) then
 
158
                print *,'ffcs4: error: only for ns=6, not ',ns
 
159
                stop
 
160
            endif
 
161
            do i=1,4
 
162
                if ( ipi12(i).ne.0 ) then
 
163
                    print *,'ffcs4: error: ipi12(',i,') non-zero! ',
 
164
     +                  ipi12(i)
 
165
                endif
 
166
            enddo
 
167
        endif
 
168
        ip = ii+3
 
169
        if ( isoort(2) .eq. 0 ) then
 
170
            ld2yzz = .FALSE.
 
171
        else
 
172
            ld2yzz = .TRUE.
 
173
        endif
 
174
        if ( isoort(4) .eq. 0 ) then
 
175
            ld2yww = .FALSE.
 
176
        else
 
177
            ld2yww = .TRUE.
 
178
        endif
 
179
        call ffieps(iepz,cz,cpi(ip),cpiDpj(ip,ii),isoort)
 
180
        call ffieps(iepw,cw,cp2p,cpiDpj(ip,ii),isoort(3))
 
181
        if ( isoort(4) .eq. 0 ) then
 
182
            print *,'ffcs4: error: case not implemented'
 
183
            ier = ier + 50
 
184
        endif
 
185
*  #] get counters: 
 
186
*  #[ R's:
 
187
        if ( isoort(4) .eq. 0 ) then
 
188
            call ffcrr(cs3(1),ipi12(1),cy(2),cy(4),cz(1),cz(3),cdyz(2,1)
 
189
     +          ,ld2yzz,cd2yzz,cz(2),cz(4),isoort(4),iepz(1),ier)
 
190
        else
 
191
            if (lwrite) print *,'ffcs4: to ffdcrr(zm,wp)'
 
192
            if ( .not. ( cdwz(2,1).eq.0 .and. iepz(1).eq.iepw(2) ) )
 
193
     +      call ffdcrr(cs3( 1),ipi12(1),cy(2),cy(4),cz(1),cz(3),cz(2),
 
194
     +          cz(4),cd2yzz,cw(2),cw(4),cw(1),cw(3),cd2yww,cdyz(2,1),
 
195
     +          cdwy(2,2),cdwz(2,1),isoort(4),iepz(1),iepw(2),ier)
 
196
        endif
 
197
        if ( isoort(2) .eq. 0 ) then
 
198
            call ffcrr(cs3(1),ipi12(1),cy(2),cy(4),cw(1),cw(3),-cdwy(1,2
 
199
     +          ),ld2yww,cd2yww,cw(2),cw(4),isoort(2),iepw(1),ier)
 
200
        else
 
201
            if (lwrite) print *,'ffcs4: to ffdcrr(zp,wm)'
 
202
            if ( .not. ( cdwz(1,2).eq.0 .and. iepz(2).eq.iepw(1) ) )
 
203
     +      call ffdcrr(cs3(21),ipi12(3),cy(2),cy(4),cz(2),cz(4),cz(1),
 
204
     +          cz(3),cd2yzz,cw(1),cw(3),cw(2),cw(4),cd2yww,cdyz(2,2),
 
205
     +          cdwy(1,2),cdwz(1,2),iepz(2),isoort(2),iepw(1),ier)
 
206
        endif
 
207
*  #] R's: 
 
208
*  #[ eta's:
 
209
        if ( DIMAG(cpi(ip)) .eq. 0 ) then
 
210
            call ffgeta(nz,cz,cdyz,cd2yzz,
 
211
     +          cpi(ip),cpiDpj(ii,ip),iepz,isoort,ier)
 
212
            do 120 i=1,2
 
213
                do 110 j=1,2
 
214
                    cdyw(i,j) = cdwy(j,i)
 
215
  110           continue
 
216
  120       continue
 
217
            call ffgeta(nw,cw,cdyw,cd2yww,
 
218
     +          cp2p,cpiDpj(ii,ip),iepw,isoort(3),ier)
 
219
        else
 
220
            print *,'ffcs4: error: not ready for complex D0 yet'
 
221
        endif
 
222
        ntot = nz(1)+nz(2)+nz(3)+nz(4)-nw(1)-nw(2)-nw(3)-nw(4)
 
223
        if ( ntot .ne. 0 ) then
 
224
            i2pi = 0
 
225
            if ( 1/absc(cy(2)) .lt. xloss ) then
 
226
                clogy = zfflo1(1/cy(2),ier)
 
227
            else
 
228
                c = -cy(4)/cy(2)
 
229
                if ( DBLE(c) .gt. -abs(DIMAG(c)) ) then
 
230
                    clogy = zfflog(c,0,c0,ier)
 
231
                else
 
232
*                   take out the factor 2*pi^2
 
233
                    cc = c+1
 
234
                    if ( absc(cc) .lt. xloss ) then
 
235
                        c2y1 = -cd2yzz - cz(1) + cz(4)
 
236
                        if ( absc(c2y1) .lt. xloss*max(absc(cz(1)),
 
237
     +                       absc(cz(4))) ) then
 
238
                            c2y1 = -cd2yzz - cz(2) + cz(3)
 
239
                            if ( lwarn .and. absc(c2y1) .lt. xloss*max(
 
240
     +                          absc(cz(2)),absc(cz(3))) ) then
 
241
                                call ffwarn(134,ier,absc(c2y1),
 
242
     +                                  absc(cy(2)))
 
243
                            endif
 
244
                        endif
 
245
                        if ( lwrite ) then
 
246
                            print *,'1+c         = ',1+c
 
247
                            print *,'-c2y1/cy(2) = ',-c2y1/cy(2)
 
248
                        endif
 
249
                        clogy = zfflo1(-c2y1/cy(2),ier)
 
250
                    else
 
251
                        clogy = zfflog(-c,0,c0,ier)
 
252
                    endif
 
253
                    if ( DIMAG(c) .lt. 0 ) then
 
254
                        i2pi = -1
 
255
                    elseif ( DIMAG(c) .gt. 0 ) then
 
256
                        i2pi = +1
 
257
                    else
 
258
                        call fferr(51,ier)
 
259
                        i2pi = 0
 
260
                    endif
 
261
                    ipi12(2) = ipi12(2) - ntot*24*i2pi
 
262
                endif
 
263
            endif
 
264
            if ( cs3(40) .ne. 0 ) print *,'ffcs4: error: cs3(40) != 0'
 
265
            cs3(40) = ntot*c2ipi*clogy
 
266
        endif
 
267
        if ( lwrite ) then
 
268
            print *,'eta''s:'
 
269
            print *,'nzi  :',nz
 
270
            print *,'nwi  :',nw
 
271
            print *,'total:',ntot*c2ipi*clogy
 
272
            if ( i2pi .ne. 0 ) print *,'     +',-ntot*24*i2pi*pi12
 
273
            print *,'     =',ntot,' *( ',c2ipi*clogy,' + ',24*i2pi*pi12,
 
274
     +          ')'
 
275
        endif
 
276
*  #] eta's: 
 
277
*###] ffcs4:
 
278
        end
 
279
*###[ ffdcxr:
 
280
        subroutine ffdcxr(cs3,ipi12,y,y1,z,z1,zp,zp1,d2yzz,
 
281
     +                  w,w1,wp,wp1,d2yww,dyz,dwy,dwz,iepsz,iepsw,ier)
 
282
***#[*comment:***********************************************************
 
283
*                                                                       *
 
284
*       Calculate                                                       *
 
285
*                                                                       *
 
286
*               R(y,z,iepsz) - R(y,w,iepsw)                             *
 
287
*                                                                       *
 
288
*       Input:                                                          *
 
289
*               a = [yzw]       (real)          see definition          *
 
290
*               a1 = 1 - a      (real)                                  *
 
291
*               dab = a - b     (real)                                  *
 
292
*               ieps[zw]        (integer)       sign of imaginary part  *
 
293
*                                               of argument logarithm   *
 
294
*               cs3(20)         (complex)       assumed zero            *
 
295
*                                                                       *
 
296
*       Output:                                                         *
 
297
*               cs3(20)         (complex)       the results, not added  *
 
298
*               ipi12(2)        (integer)       factors pi^2/12         *
 
299
*                                                                       *
 
300
*       Calls:  ffcxr                                                   *
 
301
*                                                                       *
 
302
***#]*comment:*********************************************************** 
 
303
*  #[ declarations:
 
304
        implicit none
 
305
*
 
306
*       arguments:
 
307
*
 
308
        integer ipi12(2),iepsz,iepsw,ier
 
309
        DOUBLE COMPLEX cs3(20)
 
310
        DOUBLE PRECISION y,z,w,y1,z1,w1,dyz,dwy,dwz,zp,zp1,d2yzz,wp,wp1,
 
311
     +          d2yww
 
312
*
 
313
*       local variables:
 
314
*
 
315
        integer i,ieps,ipi12p(2),ier1,ier2,isign,inorm
 
316
        logical again
 
317
        DOUBLE PRECISION yy,yy1,zz,zz1,dyyzz,xx1,xx1n,term,tot,d2,d3,
 
318
     +          d21,d31,d2n,d3n,d21n1,d31n1,dw,xlogy,x00(3)
 
319
        DOUBLE COMPLEX csum,csum1,csum2,cs3p(20),chulp
 
320
        DOUBLE PRECISION dfflo1
 
321
*
 
322
*       common blocks
 
323
*
 
324
        include 'ff.h'
 
325
*
 
326
*       statement function
 
327
*
 
328
*       absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
329
        inorm = 0
 
330
*  #] declarations: 
 
331
*  #[ groundwork:
 
332
        if ( dwz .eq. 0 .and. iepsz .eq. iepsw ) return
 
333
        if ( dyz .eq. 0 ) then
 
334
            call fferr(75,ier)
 
335
            return
 
336
        endif
 
337
        xx1 = y/dyz
 
338
        dw = dwz/dyz
 
339
        if ( xx1 .le. x05 .or. xx1 .gt. 2 ) then
 
340
            d2 = 1/y
 
341
            dw = dw*y/w
 
342
        else
 
343
            d2 = 1/z1
 
344
        endif
 
345
        again = .FALSE.
 
346
  123   continue
 
347
*  #] groundwork: 
 
348
*  #[ trivial case:
 
349
        if ( dw .eq. 0 ) then
 
350
            if ( lwrite ) print *,'  Trivial case'
 
351
*  #] trivial case: 
 
352
*  #[ normal case:
 
353
        elseif ( abs(dw) .gt. xloss .or. again ) then
 
354
*           nothing's the matter
 
355
            if ( lwrite ) print *,'  Normal case'
 
356
            inorm = 1
 
357
            call ffcxr(cs3( 1),ipi12(1),y,y1,z,z1,dyz,
 
358
     +          .TRUE.,d2yzz,zp,zp1,.FALSE.,x00,iepsz,ier)
 
359
            call ffcxr(cs3(11),ipi12(2),y,y1,w,w1,-dwy,
 
360
     +          .TRUE.,d2yww,wp,wp1,.FALSE.,x00,iepsw,ier)
 
361
            do 10 i=11,20
 
362
   10           cs3(i) = -cs3(i)
 
363
            ipi12(2) = -ipi12(2)
 
364
*  #] normal case: 
 
365
*  #[ only cancellations in w, not in y:
 
366
        elseif ( abs(d2) .gt. xloss ) then
 
367
*           there are no cancellations the other way:
 
368
            if ( lwrite ) print *,'  Cancellations one way, turned Rs'
 
369
            if ( iepsz .ne. iepsw .and. ( y/dyz .gt. 1 .or.-y/dwy .gt.
 
370
     +          1 ) ) then
 
371
                again = .TRUE.
 
372
                if ( lwrite ) then
 
373
                    print *,'ffdcxr: problems with ieps, solvable,'
 
374
                    print *,'        but for the moment just call the ',
 
375
     +                  'normal case'
 
376
                endif
 
377
                again = .TRUE.
 
378
                goto 123
 
379
*               call fferr(21,ier)
 
380
            endif
 
381
            yy = dwy/dwz
 
382
            zz = yy*z/y
 
383
            yy1 = dyz/dwz
 
384
            zz1 = yy1*w/y
 
385
            dyyzz = yy*dyz/y
 
386
            if ( y .lt. 0 ) then
 
387
                ieps = iepsz
 
388
            else
 
389
                ieps = -iepsz
 
390
            endif
 
391
            call ffcxr(cs3( 1),ipi12(1),yy,yy1,zz,zz1,dyyzz,.FALSE.,x0,
 
392
     +                  x0,x0,.FALSE.,x00,2*ieps,ier)
 
393
            zz = yy*z1/y1
 
394
            zz1 = yy1*w1/y1
 
395
            dyyzz = -yy*dyz/y1
 
396
            if ( y1 .gt. 0 ) then
 
397
                ieps = iepsz
 
398
            else
 
399
                ieps = -iepsz
 
400
            endif
 
401
            call ffcxr(cs3(11),ipi12(2),yy,yy1,zz,zz1,dyyzz,.FALSE.,x0,
 
402
     +                  x0,x0,.FALSE.,x00,2*ieps,ier)
 
403
            do 20 i=11,20
 
404
                cs3(i) = -cs3(i)
 
405
   20       continue
 
406
            ipi12(2) = -ipi12(2)
 
407
*  #] only cancellations in w, not in y: 
 
408
*  #[ Hill identity:
 
409
        elseif (  ( 1 .gt. xloss*abs(y) .or. abs(xx1) .gt. xloss )
 
410
     +      .and. ( 1 .gt. xloss*abs(z) .or. abs(z/dyz) .gt. xloss )
 
411
     +      .and. ( 1 .gt. xloss*abs(y) .or. abs(dyz/y) .gt. xloss )
 
412
     +          ) then
 
413
*           do a Hill identity on the y,y-1 direction
 
414
            if ( lwrite ) print *,'  Hill identity to split z,w'
 
415
            yy = -y*w1/dwy
 
416
            yy1 = w*y1/dwy
 
417
            zz = -z*w1/dwz
 
418
            zz1 = w*z1/dwz
 
419
            dyyzz = -w*w1*(dyz/(dwy*dwz))
 
420
            if ( y*dwz .gt. 0 .eqv. (y+dwz) .gt. 0 ) then
 
421
                ieps = 2*iepsw
 
422
            else
 
423
                ieps = -2*iepsw
 
424
            endif
 
425
            call ffcxr(cs3( 1),ipi12(1),yy,yy1,zz,zz1,dyyzz,.FALSE.,x0,
 
426
     +                  x0,x0,.FALSE.,x00,ieps,ier)
 
427
            yy = w1
 
428
            yy1 = w
 
429
            zz = -w1*z/dwz
 
430
            zz1 = w*z1/dwz
 
431
            dyyzz = w*w1/dwz
 
432
            call ffcxr(cs3( 9),ipi12(2),yy,yy1,zz,zz1,dyyzz,.FALSE.,x0,
 
433
     +                  x0,x0,.FALSE.,x00,ieps,ier)
 
434
            do 30 i=9,16
 
435
   30           cs3(i) = -cs3(i)
 
436
            ipi12(2) = -ipi12(2)
 
437
*           the extra logarithms ...
 
438
            if ( 1 .lt. xloss*abs(w) ) then
 
439
                chulp = dfflo1(1/w,ier)
 
440
            elseif ( w1 .lt. 0 .or. w .lt. 0 ) then
 
441
                chulp = log(-w1/w)
 
442
            else
 
443
                chulp = DCMPLX(DBLE(log(w1/w)),DBLE(-iepsw*pi))
 
444
            endif
 
445
            cs3(20) = -DBLE(dfflo1(dwz/dwy,ier))*chulp
 
446
*  #] Hill identity: 
 
447
*  #[ Taylor expansion:
 
448
        elseif ( (w.lt.0..or.w1.lt.0) .and. (z.lt.0..or.z1.lt.0) ) then
 
449
*           do a Taylor expansion
 
450
            if ( abs(xx1) .lt. xloss ) then
 
451
                if ( lwrite ) print *,'ffdcxr: Taylor expansion, normal'
 
452
                d3 = dwz/dwy
 
453
*               isign = 1
 
454
                xx1n = xx1
 
455
                d2n = d2
 
456
                d3n = d3
 
457
                d21 = 1-d2
 
458
                d21n1 = 1
 
459
                d31 = 1-d3
 
460
                d31n1 = 1
 
461
                tot = xx1*d2*d3
 
462
                do 50 i=2,20
 
463
                    xx1n = xx1n*xx1
 
464
                    d21n1 = d21n1*d21
 
465
                    d31n1 = d31n1*d31
 
466
                    d2n = d2n + d2*d21n1
 
467
                    d3n = d3n + d3*d31n1
 
468
                    term = xx1n*d2n*d3n*xn2inv(i)
 
469
                    tot = tot + term
 
470
                    if ( abs(term) .le. precx*abs(tot) ) goto 51
 
471
   50           continue
 
472
                if ( lwarn ) call ffwarn(46,ier,tot,term)
 
473
   51           continue
 
474
*               if ( isign .eq. 1 ) then
 
475
                    cs3(1) = tot
 
476
*               else
 
477
*                   cs3(1) = -tot
 
478
*               endif
 
479
            elseif ( abs(z/dyz) .lt. xloss ) then
 
480
                if ( lwrite ) print *,'  Normal case'
 
481
                inorm = 1
 
482
                call ffcxr(cs3( 1),ipi12(1),y,y1,z,z1,dyz,
 
483
     +               .TRUE.,d2yzz,zp,zp1,.FALSE.,x00,iepsz,ier)
 
484
                call ffcxr(cs3(11),ipi12(2),y,y1,w,w1,-dwy,
 
485
     +               .TRUE.,d2yww,wp,wp1,.FALSE.,x00,iepsw,ier)
 
486
                do 110 i=11,20
 
487
  110               cs3(i) = -cs3(i)
 
488
*               if ( lwrite ) print *,'ffdcxr: Taylor expansion, 1-x'
 
489
*               print *,'NOT YET READY !!!'
 
490
*               ier = ier + 100
 
491
*               yy = y1*dwz/(z1*dwy)
 
492
*               if ( abs(yy) .lt. xloss ) then
 
493
*                   cs3(10) = -dfflo1(1/y,ier)*dfflo1(yy,ier)
 
494
*               else
 
495
*                   yy1 = -w1*dyz/(z1*dwy)
 
496
*                   if ( yy1 .gt. xalogm ) then
 
497
*                       cs3(10) = -dfflo1(1/y,ier)*log(yy1)
 
498
*                   elseif ( yy1 .gt. -xalogm ) then
 
499
*                       if ( lwarn ) call ffwarn(80,ier,yy1,xalogm)
 
500
*                   else
 
501
*                       xlogy = log(-yy1)
 
502
*                       if ( lwarn .and. iepsz.ne.iepsw )
 
503
*     +                         call ffwarn(81,ier,x1,x1)
 
504
*                       if ( (w1+dyz)*dwz*y1*iepsz .lt. 0 ) then
 
505
*                           cs3(10) = -dfflo1(1/y,ier)*DCMPLX(DBLE(xlogy),DBLE(pi))
 
506
*                       else
 
507
*                           cs3(10) = -dfflo1(1/y,ier)*DCMPLX(DBLE(xlogy),DBLE(-pi))
 
508
*                       endif
 
509
*                   endif
 
510
*               endif
 
511
*               cs3(11) = -dfflo1(1/z,ier)*dfflo1(dwz/dwy,ier)
 
512
*               yy = dwz/(w*z1)
 
513
*               if ( abs(yy) .lt. xloss ) then
 
514
*                   cs3(12) = -dfflo1(w/dwy,ier)*dfflo1(yy,ier)
 
515
*               else
 
516
*                   yy1 = z*w1/(w*z1)
 
517
*                   if ( yy1 .gt. xalogm ) then
 
518
*                       cs3(12) = -dfflo1(w/dwy,ier)*log(yy1)
 
519
*                   elseif (  yy .gt. -xalogm ) then
 
520
*                       if ( lwarn ) call ffwarn(80,ier,yy,xalogm)
 
521
*                   else
 
522
*                       xlogy = log(-yy1)
 
523
*                       if ( lwarn .and. iepsz.ne.iepsw )
 
524
*     +                         call ffwarn(81,ier,x1,x1)
 
525
*                       if ( dwz*(dwz+1)*ieps .gt. 0 ) then
 
526
*                           cs3(12) = -dfflo1(w/dwy,ier)*DCMPLX(DBLE(xlogy),DBLE(pi))
 
527
*                       else
 
528
*                           cs3(12) =-dfflo1(w/dwy,ier)*DCMPLX(DBLE(xlogy),DBLE(-pi))
 
529
*                       endif
 
530
*                   endif
 
531
*               endif
 
532
*               isign = -1
 
533
*               xx1 = -z/dyz
 
534
*               d2 = 1/z
 
535
*               d3 = dwz/dwy
 
536
            else
 
537
                if ( lwrite ) print *,'ffdcxr: Taylor expansion, 1/x'
 
538
                call fferr(22,ier)
 
539
                return
 
540
            endif
 
541
        else
 
542
            if ( lwrite ) print *,'Not clear, take normal route'
 
543
            inorm = 1
 
544
            call ffcxr(cs3( 1),ipi12(1),y,y1,z,z1,dyz,.FALSE.,x0,x0,x0,
 
545
     +                          .FALSE.,x00,iepsz,ier)
 
546
            call ffcxr(cs3(11),ipi12(2),y,y1,w,w1,-dwy,.FALSE.,x0,x0,x0,
 
547
     +                          .FALSE.,x00,iepsw,ier)
 
548
            do 40 i=11,20
 
549
   40           cs3(i) = -cs3(i)
 
550
            ipi12(2) = -ipi12(2)
 
551
        endif
 
552
*  #] Taylor expansion: 
 
553
*  #[ debug output:
 
554
        if ( lwrite ) then
 
555
            csum = 0
 
556
            do 900 i=1,20
 
557
                csum = csum + cs3(i)
 
558
                print '(i2,2g16.8)',i,cs3(i)
 
559
  900       continue
 
560
            print '(a)','---------------------------------'
 
561
            print '(2x,2g16.8,2i3)',csum,ipi12
 
562
            print '(a,i3)','ier = ',ier
 
563
            if ( inorm .eq. 0 ) then
 
564
            lwrite = .FALSE.
 
565
            ier1 = 0
 
566
            ier2 = 0
 
567
            do 905 i=1,20
 
568
  905           cs3p(i) = 0
 
569
            ipi12p(1) = 0
 
570
            ipi12p(2) = 0
 
571
            call ffcxr(cs3p( 1),ipi12p(1),y,y1,z,z1,dyz,.FALSE.,x0,x0,
 
572
     +                          x0,.FALSE.,x00,iepsz,ier1)
 
573
            call ffcxr(cs3p(11),ipi12p(2),y,y1,w,w1,-dwy,.FALSE.,x0,x0,
 
574
     +                          x0,.FALSE.,x00,iepsw,ier2)
 
575
            csum1 = 0
 
576
            do 910 i=1,10
 
577
  910           csum1 = csum1 + cs3p(i)
 
578
            csum2 = 0
 
579
            do 920 i=11,20
 
580
  920           csum2 = csum2 - cs3p(i)
 
581
            csum = csum1 + csum2 + (ipi12p(1)-ipi12p(2))*DBLE(pi12)
 
582
            print *,'cmp with:'
 
583
            print '(i2,2g16.8,i3)',1,csum1,ier1
 
584
            print '(i2,2g16.8,i3)',2,csum2,ier2
 
585
            print *,'------------------+'
 
586
            print '(2x,2g16.8,3i3)',csum1+csum2,ipi12p,max(ier1,ier2)
 
587
            print '(2x,2g16.8,3i3)',csum
 
588
            lwrite = .TRUE.
 
589
            endif
 
590
        endif
 
591
*  #] debug output: 
 
592
*###] ffdcxr: 
 
593
        end
 
594
*###[ ffdcrr:
 
595
        subroutine ffdcrr(cs3,ipi12,cy,cy1,cz,cz1,czp,czp1,cd2yzz,cw,cw1
 
596
     +          ,cwp,cwp1,cd2yww,cdyz,cdwy,cdwz,isoort,iepsz,iepsw,ier)
 
597
***#[*comment:***********************************************************
 
598
*                                                                       *
 
599
*       Calculate                                                       *
 
600
*                                                                       *
 
601
*               R(cy,cz,iepsz) - R(cy,cw,iepsw)                         *
 
602
*                                                                       *
 
603
*       Input:                                                          *
 
604
*               a = [yzw]       (real)          see definition          *
 
605
*               a1 = 1 - a      (real)                                  *
 
606
*               dab = a - b     (real)                                  *
 
607
*               ieps[zw]        (integer)       sign of imaginary part  *
 
608
*                                               of argument logarithm   *
 
609
*               cs3(20)         (complex)       assumed zero            *
 
610
*                                                                       *
 
611
*       Output:                                                         *
 
612
*               cs3(20)         (complex)       the results, not added  *
 
613
*               ipi12(2)        (integer)       factors pi^2/12         *
 
614
*                                                                       *
 
615
*       Calls:  ffcrr                                                   *
 
616
*                                                                       *
 
617
***#]*comment:*********************************************************** 
 
618
*  #[ declarations:
 
619
        implicit none
 
620
*
 
621
*       arguments:
 
622
*
 
623
        integer ipi12(2),isoort,iepsz,iepsw,ier
 
624
        DOUBLE COMPLEX cs3(20)
 
625
        DOUBLE COMPLEX cy,cz,czp,cw,cwp,cy1,cz1,czp1,cw1,cwp1,
 
626
     +          cdyz,cdwy,cdwz,cd2yzz,cd2yww
 
627
*
 
628
*       local variables:
 
629
*
 
630
        integer i,ieps,ieps1,ieps2,ipi12p(2),ier1,ier2,isign,inorm,i2pi,
 
631
     +          nffeta,nffet1,n1,n2,n3,n4,n5,n6
 
632
        logical ld2yyz
 
633
        DOUBLE COMPLEX cyy,cyy1,czz,czz1,cdyyzz,chulp,zfflo1,zfflog,
 
634
     +          cc1,cdw,cc1n,cterm,ctot,cd2,cd3,
 
635
     +          cd21,cd31,cd2n,cd3n,cd21n1,cd31n1,
 
636
     +          cc2,cfactz,cfactw,czzp,czzp1,cd2yyz
 
637
        DOUBLE COMPLEX csum,csum1,csum2,cs3p(20),c,check
 
638
        DOUBLE PRECISION absc,xlosn
 
639
*
 
640
*       common blocks
 
641
*
 
642
        include 'ff.h'
 
643
*
 
644
*       statement function
 
645
*
 
646
        absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
647
        inorm = 0
 
648
*  #] declarations: 
 
649
*  #[ check input:
 
650
        if ( ltest ) then
 
651
            xlosn = xloss*DBLE(10)**(-1-mod(ier,50))
 
652
            check = cd2yzz - 2*cy + cz + czp
 
653
            if ( xlosn*absc(check) .gt. precc*max(2*absc(cy),absc(cz),
 
654
     +                  absc(czp)) ) then
 
655
                print *,'ffdcrr: error:  cd2yzz != 2*cy - cz - czp:',
 
656
     +                  cd2yzz,cy,cz,czp,check
 
657
            endif
 
658
            check = cd2yww - 2*cy + cw + cwp
 
659
            if ( xlosn*absc(check) .gt. precc*max(2*absc(cy),absc(cw),
 
660
     +                  absc(cwp)) ) then
 
661
                print *,'ffdcrr: error:  cd2yww != 2*cy - cw - cwp:',
 
662
     +                  cd2yww,cy,cw,cwp,check
 
663
            endif
 
664
        endif
 
665
*  #] check input:
 
666
*  #[ groundwork:
 
667
        if ( cdwz .eq. 0 ) then
 
668
            if ( abs(DIMAG(cz)) .gt. precc*abs(DBLE(cz)) .or.
 
669
     +          iepsz .eq. iepsw ) return
 
670
            if ( DBLE(cz) .ge. 0 .and. DBLE(cz1) .ge. 0 ) return
 
671
            call fferr(76,ier)
 
672
            return
 
673
        endif
 
674
        if ( cdyz .eq. 0 ) then
 
675
            call fferr(77,ier)
 
676
            return
 
677
        endif
 
678
        cc1 = cy/cdyz
 
679
        cdw = cdwz/cdyz
 
680
        if ( DBLE(cc1) .le. x05 .or. abs(cc1-1) .gt. 1 ) then
 
681
            cd2 = 1/cy
 
682
            cdw = cdw*cy/cw
 
683
        else
 
684
            cd2 = 1/cz1
 
685
        endif
 
686
*  #] groundwork:
 
687
*  #[ trivial case:
 
688
        if ( absc(cdw) .eq. 0 ) then
 
689
            if ( lwrite ) print *,'  Trivial case'
 
690
*  #] trivial case: 
 
691
*  #[ normal case:
 
692
*
 
693
*       if no cancellations are expected OR the imaginary signs differ
 
694
*       and are significant
 
695
*
 
696
        elseif ( absc(cdw) .gt. xloss .or. (iepsz.ne.iepsw .and.
 
697
     +          (DBLE(cy/cdyz).gt.1 .or. DBLE(-cy1/cdyz).gt.1) ) ) then
 
698
*           nothing's the matter
 
699
            if ( lwrite ) print *,'ffdcrr: Normal case'
 
700
            inorm = 1
 
701
*           special case to avoid bug found 15-oct=1995
 
702
            if ( iepsz.eq.iepsw ) then
 
703
                if ( DIMAG(cz).eq.0 .and. DIMAG(cz1).eq.0 ) then
 
704
                    print *,'ffdcrr: flipping sign iepsz'
 
705
                    iepsz = -iepsz
 
706
                elseif ( DIMAG(cw).eq.0 .and. DIMAG(cw1).eq.0 ) then
 
707
                    print *,'ffdcrr: flipping sign iepsw'
 
708
                    iepsw = -iepsw
 
709
                else
 
710
                    print *,'ffdcrr: error: missing eta terms!'
 
711
                    ier = ier + 100
 
712
                endif
 
713
            endif
 
714
            call ffcrr(cs3(1),ipi12(1),cy,cy1,cz,cz1,cdyz,.TRUE.,
 
715
     +          cd2yzz,czp,czp1,isoort,iepsz,ier)
 
716
            call ffcrr(cs3(8),ipi12(2),cy,cy1,cw,cw1,-cdwy,.TRUE.,
 
717
     +          cd2yww,cwp,cwp1,isoort,iepsw,ier)
 
718
            do 10 i=8,14
 
719
                cs3(i) = -cs3(i)
 
720
   10       continue
 
721
            ipi12(2) = -ipi12(2)
 
722
*  #] normal case:
 
723
*  #[ only cancellations in cw, not in cy:
 
724
        elseif ( absc(cd2) .gt. xloss ) then
 
725
*           there are no cancellations the other way:
 
726
            if ( lwrite ) print *,'ffdcrr: Cancellations one way, ',
 
727
     +          'turned Rs'
 
728
            cyy = cdwy/cdwz
 
729
            czz = cz*cyy/cy
 
730
            cyy1 = cdyz/cdwz
 
731
            czz1 = cyy1*cw/cy
 
732
            cdyyzz = cdyz*cyy/cy
 
733
            if ( DBLE(cy) .gt. 0 ) then
 
734
                ieps1 = -3*iepsz
 
735
            else
 
736
                ieps1 = +3*iepsz
 
737
            endif
 
738
*           Often 2y-z-z is relevant, but 2*yy-zz-zz is not, solve by
 
739
*           introducing zzp.
 
740
            czzp = czp*cyy/cy
 
741
            cd2yyz = cd2yzz*cyy/cy
 
742
            czzp1 = 1 - czzp
 
743
            if ( absc(czzp1) .lt. xloss ) then
 
744
*               later try more possibilities
 
745
                ld2yyz = .FALSE.
 
746
            else
 
747
                ld2yyz = .TRUE.
 
748
            endif
 
749
            call ffcrr(cs3(1),ipi12(1),cyy,cyy1,czz,czz1,cdyyzz,
 
750
     +          ld2yyz,cd2yyz,czzp,czzp1,isoort,ieps1,ier)
 
751
            czz = cyy*cz1/cy1
 
752
            czz1 = cyy1*cw1/cy1
 
753
            if ( DBLE(-cy1) .gt. 0 ) then
 
754
                ieps2 = -3*iepsz
 
755
            else
 
756
                ieps2 = +3*iepsz
 
757
            endif
 
758
            cdyyzz = -cyy*cdyz/cy1
 
759
            czzp = czp1*cyy/cy1
 
760
            cd2yyz = -cd2yzz*cyy/cy1
 
761
            czzp1 = 1 - czzp
 
762
            if ( absc(czzp1) .lt. xloss ) then
 
763
*               later try more possibilities
 
764
                ld2yyz = .FALSE.
 
765
            else
 
766
                ld2yyz = .TRUE.
 
767
            endif
 
768
            call ffcrr(cs3(8),ipi12(2),cyy,cyy1,czz,czz1,cdyyzz,
 
769
     +          .TRUE.,cd2yyz,czzp,czzp1,isoort,ieps2,ier)
 
770
            do 20 i=8,14
 
771
                cs3(i) = -cs3(i)
 
772
   20       continue
 
773
            ipi12(2) = -ipi12(2)
 
774
*           eta terms (are not calculated in ffcrr as ieps = 3)
 
775
            cfactz = 1/cdyz
 
776
            if ( DIMAG(cz) .eq. 0 ) then
 
777
                if ( DIMAG(cy) .eq. 0 ) then
 
778
                    n1 = 0
 
779
                    n2 = 0
 
780
                else
 
781
                    n1 = nffet1(DCMPLX(DBLE(0),DBLE(iepsz)),cfactz,
 
782
     +                  -cz*cfactz,ier)
 
783
                    n2 = nffet1(DCMPLX(DBLE(0),DBLE(iepsz)),cfactz,
 
784
     +                  cz1*cfactz,ier)
 
785
                endif
 
786
            else
 
787
                n1 = nffeta(-cz,cfactz,ier)
 
788
                n2 = nffeta(cz1,cfactz,ier)
 
789
            endif
 
790
            cfactw = -1/cdwy
 
791
            if ( DIMAG(cw) .eq. 0 ) then
 
792
                if ( DIMAG(cy) .eq. 0 ) then
 
793
                    n4 = 0
 
794
                    n5 = 0
 
795
                else
 
796
                    n4 = nffet1(DCMPLX(DBLE(0),DBLE(iepsw)),cfactw,
 
797
     +                  -cw*cfactw,ier)
 
798
                    n5 = nffet1(DCMPLX(DBLE(0),DBLE(iepsw)),cfactw,
 
799
     +                  cw1*cfactw,ier)
 
800
                endif
 
801
            else
 
802
                n4 = nffeta(-cw,cfactw,ier)
 
803
                n5 = nffeta(cw1,cfactw,ier)
 
804
            endif
 
805
*
 
806
*           we assume that cs3(15-17) are not used, this is always true
 
807
*
 
808
            n3 = 0
 
809
            n6 = 0
 
810
            if ( n1.eq.n4 ) then
 
811
                if ( n1.eq.0 ) then
 
812
*                   nothing to do
 
813
                else
 
814
                    cc1 = cdwz/cdyz
 
815
                    if ( absc(cc1) .lt. xloss ) then
 
816
                        cs3(15) = n1*c2ipi*zfflo1(cc1,ier)
 
817
                    else
 
818
                        cc1 = -cdwy/cdyz
 
819
                        cs3(15) = n1*c2ipi*zfflog(cc1,0,c0,ier)
 
820
                    endif
 
821
                    cc1 = cy*cfactz
 
822
                    cc2 = cy*cfactw
 
823
                    if ( DIMAG(cc1).eq.0 .or. DIMAG(cc2).eq.0 ) then
 
824
                        n3 = 0
 
825
                    else
 
826
                        n3 = nffeta(cc1,1/cc2,ier)
 
827
                    endif
 
828
                    if ( n3.ne.0 ) then
 
829
                        print *,'ffdcrr: error: untested algorithm'
 
830
                        ier = ier + 50
 
831
                        ipi12(1) = ipi12(1) + 4*12*n1*n3
 
832
                    endif
 
833
                endif
 
834
            else
 
835
                cc1 = cy*cfactz
 
836
                cc2 = cy*cfactw
 
837
                cs3(15) = (n1*zfflog(cc1,ieps1,c0,ier) + 
 
838
     +                     n4*zfflog(cc2,ieps1,c0,ier))*c2ipi
 
839
            endif
 
840
            if ( n2.eq.n5 ) then
 
841
                if ( n2.eq.0 ) then
 
842
*                   nothing to do
 
843
                else
 
844
                    cc1 = cdwz/cdyz
 
845
                    if ( absc(cc1) .lt. xloss ) then
 
846
                        cs3(16) = n2*c2ipi*zfflo1(cc1,ier)
 
847
                    else
 
848
                        cc1 = -cdwy/cdyz
 
849
                        cs3(16) = n2*c2ipi*zfflog(cc1,0,c0,ier)
 
850
                    endif
 
851
                    cc1 = -cy1*cfactz
 
852
                    cc2 = -cy1*cfactw
 
853
                    if ( DIMAG(cc1).eq.0 .or. DIMAG(cc2).eq.0 ) then
 
854
                        n6 = 0
 
855
                    else
 
856
                        n6 = nffeta(cc1,1/cc2,ier)
 
857
                    endif
 
858
                    if ( n6.ne.0 ) then
 
859
                        print *,'ffdcrr: error: untested algorithm'
 
860
                        ier = ier + 50
 
861
                        ipi12(2) = ipi12(2) + 4*12*n2*n6
 
862
                    endif
 
863
                endif
 
864
            else
 
865
                cc1 = -cy1*cfactz
 
866
                cc2 = -cy1*cfactw
 
867
                cs3(15) = (n2*zfflog(cc1,ieps2,c0,ier) + 
 
868
     +                     n5*zfflog(cc2,ieps2,c0,ier))*c2ipi
 
869
            endif
 
870
            if ( lwrite ) then
 
871
                print *,'  eta''s z are :',n1,n2,n3
 
872
                print *,'  eta''s w are :',n4,n5,n6
 
873
            endif
 
874
*  #] only cancellations in cw, not in cy: 
 
875
*  #[ Hill identity:
 
876
        elseif (  ( 1.gt.xloss*absc(cy) .or. absc(cc1).gt.xloss )
 
877
     +      .and. ( 1.gt.xloss*absc(cz) .or. absc(cz/cdyz).gt.xloss )
 
878
     +      .and. ( 1.gt.xloss*absc(cy) .or. absc(cdyz/cy).gt.xloss )
 
879
     +          ) then
 
880
*           do a Hill identity on the cy,cy-1 direction
 
881
            if ( lwrite ) print *,'ffdcrr: Hill identity to split cz,cw'
 
882
            cyy = -cy*cw1/cdwy
 
883
            cyy1 = cw*cy1/cdwy
 
884
            czz = -cz*cw1/cdwz
 
885
            czz1 = cw*cz1/cdwz
 
886
            cdyyzz = -cw*cw1*(cdyz/(cdwy*cdwz))
 
887
            ieps = -2*iepsz
 
888
            call ffcrr(cs3(1),ipi12(1),cyy,cyy1,czz,czz1,cdyyzz,
 
889
     +          .FALSE.,c0,c0,c0,isoort,ieps,ier)
 
890
            cyy = cw1
 
891
            cyy1 = cw
 
892
            czz = -cw1*cz/cdwz
 
893
            czz1 = cw*cz1/cdwz
 
894
            cdyyzz = cw*cw1/cdwz
 
895
            call ffcrr(cs3(8),ipi12(2),cyy,cyy1,czz,czz1,cdyyzz,
 
896
     +          .FALSE.,c0,c0,c0,isoort,0,ier)
 
897
            do 30 i=8,14
 
898
   30           cs3(i) = -cs3(i)
 
899
            ipi12(2) = -ipi12(2)
 
900
*           the extra logarithms ...
 
901
            if ( 1 .lt. xloss*absc(cw) ) then
 
902
                chulp = zfflo1(1/cw,ier)
 
903
            else
 
904
                chulp = zfflog(-cw1/cw,0,c0,ier)
 
905
            endif
 
906
            cs3(15) = -zfflo1(cdwz/cdwy,ier)*chulp
 
907
*  #] Hill identity: 
 
908
*  #[ Taylor expansion:
 
909
        else
 
910
*           Do a Taylor expansion
 
911
            if ( absc(cc1) .lt. xloss ) then
 
912
                if ( lwrite ) print *,'ffdcrr: Taylor expansion, normal'
 
913
                cd3 = cdwz/cdwy
 
914
*               isign = 1
 
915
                cc1n = cc1
 
916
                cd2n = cd2
 
917
                cd3n = cd3
 
918
                cd21 = 1-cd2
 
919
                cd21n1 = 1
 
920
                cd31 = 1-cd3
 
921
                cd31n1 = 1
 
922
                ctot = cc1*cd2*cd3
 
923
                do 50 i=2,20
 
924
                    cc1n = cc1n*cc1
 
925
                    cd21n1 = cd21n1*cd21
 
926
                    cd31n1 = cd31n1*cd31
 
927
                    cd2n = cd2n + cd2*cd21n1
 
928
                    cd3n = cd3n + cd3*cd31n1
 
929
                    cterm = cc1n*cd2n*cd3n*DBLE(xn2inv(i))
 
930
                    ctot = ctot + cterm
 
931
                    if ( absc(cterm) .lt. precc*absc(ctot) ) goto 51
 
932
   50           continue
 
933
                if ( lwarn ) call ffwarn(45,ier,absc(ctot),absc(cterm))
 
934
   51           continue
 
935
*               if ( isign .eq. 1 ) then
 
936
                    cs3(1) = ctot
 
937
*               else
 
938
*                   cs3(1) = -ctot
 
939
*               endif
 
940
            elseif ( absc(cz/cdyz) .lt. xloss ) then
 
941
                if ( lwrite ) print *,'ffdcrr: Normal case'
 
942
                inorm = 1
 
943
                call ffcrr(cs3(1),ipi12(1),cy,cy1,cz,cz1,cdyz,.TRUE.,
 
944
     +                  cd2yzz,czp,czp1,isoort,iepsz,ier)
 
945
                call ffcrr(cs3(8),ipi12(2),cy,cy1,cw,cw1,-cdwy,.TRUE.,
 
946
     +                  cd2yww,cwp,cwp1,isoort,iepsw,ier)
 
947
                do 110 i=8,14
 
948
  110                   cs3(i) = -cs3(i)
 
949
                ipi12(2) = -ipi12(2)
 
950
*               if ( lwrite ) print *,'ffdcrr: Taylor expansion, 1-x'
 
951
*               print *,'NOT YET READY !!'
 
952
*               ier = ier + 100
 
953
*               cyy = cy1*cdwz/(cz1*cdwy)
 
954
*               if ( absc(cyy) .lt. xloss ) then
 
955
*                   cs3(10) = -zfflo1(1/cy,ier)*zfflo1(cyy,ier)
 
956
*               else
 
957
*                   cyy1 = -cw1*cdyz/(cz1*cdwy)
 
958
*                   cs3(10) = -zfflo1(1/cy,ier)*zfflog(cyy1,0,cy,ier)
 
959
*               endif
 
960
*               cs3(11) = -zfflo1(1/cz,ier)*zfflo1(cdwz/cdwy,ier)
 
961
*               cyy = cdwz/(cw*cz1)
 
962
*               if ( absc(cyy) .lt. xloss ) then
 
963
*                   cs3(12) = -zfflo1(cw/cdwy,ier)*zfflo1(cyy,ier)
 
964
*               else
 
965
*                   cyy1 = cz*cw1/(cw*cz1)
 
966
*                   cs3(12) = -zfflo1(cw/cdwy,ier)*zfflog(cyy1,0,c0,ier)
 
967
*               endif
 
968
*               isign = -1
 
969
*               cc1 = -cz/cdyz
 
970
*               cd2 = 1/cz
 
971
*               cd3 = cdwz/cdwy
 
972
            else
 
973
                if ( lwrite ) print *,'ffdcrr: Taylor expansion, 1/x'
 
974
                call fferr(20,ier)
 
975
                return
 
976
            endif
 
977
        endif
 
978
*  #] Taylor expansion: 
 
979
*  #[ debug output:
 
980
        if ( lwrite ) then
 
981
            csum = 0
 
982
            do 900 i=1,20
 
983
                csum = csum + cs3(i)
 
984
                print '(i2,2g16.8)',i,cs3(i)
 
985
  900       continue
 
986
            print '(a)','---------------------------------'
 
987
            print '(2x,2g16.8,2i3)',csum,ipi12
 
988
            print '(a,2g16.8)','= ',csum+(ipi12(1)+ipi12(2))*DBLE(pi12)
 
989
            print '(a,i3)','ier = ',ier
 
990
            if ( inorm .eq. 0 ) then
 
991
                lwrite = .FALSE.
 
992
                ier1 = 0
 
993
                ier2 = 0
 
994
                do 905 i=1,14
 
995
  905                   cs3p(i) = 0
 
996
                ipi12p(1) = 0
 
997
                ipi12p(2) = 0
 
998
                call ffcrr(cs3p(1),ipi12p(1),cy,cy1,cz,cz1,cdyz,
 
999
     +                  .TRUE.,cd2yzz,czp,czp1,isoort,iepsz,ier1)
 
1000
                call ffcrr(cs3p(8),ipi12p(2),cy,cy1,cw,cw1,-cdwy,
 
1001
     +                  .TRUE.,cd2yww,cwp,cwp1,isoort,iepsw,ier2)
 
1002
                csum1 = 0
 
1003
                do 910 i=1,7
 
1004
  910                   csum1 = csum1 + cs3p(i)
 
1005
                csum2 = 0
 
1006
                do 920 i=8,14
 
1007
  920                   csum2 = csum2 - cs3p(i)
 
1008
                print *,'cmp with:'
 
1009
                print '(i2,2g16.8,i2)',1,csum1,ier1
 
1010
                print '(i2,2g16.8,i2)',2,csum2,ier2
 
1011
                print *,'------------------+'
 
1012
                print '(2x,2g16.8,3i3)',csum1+csum2,ipi12p,
 
1013
     +                  max(ier1,ier2)
 
1014
                print '(a,2g16.8,3i3)','= ',csum1+csum2+
 
1015
     +                  (ipi12p(1)-ipi12p(2))*DBLE(pi12)
 
1016
                lwrite = .TRUE.
 
1017
            endif
 
1018
        endif
 
1019
*  #] debug output: 
 
1020
*###] ffdcrr:
 
1021
        end