← Back to branch summary

~chirality-flow/chiralityflow/ChiralityFlowMG

« back to all changes in this revision

Viewing changes to vendor/IREGI/src/funlib.f90

  • Committer: andrew.lifson at lu
  • Date: 2021-09-02 13:57:34 UTC
  • Revision ID: andrew.lifson@thep.lu.se-20210902135734-4eybgli0iljkax9b
added fresh copy of MG5_aMC_v3.2.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/IREGI/src/funlib.f90
 
1
MODULE funlib
 
2
  USE global
 
3
  IMPLICIT NONE
 
4
CONTAINS
 
5
  RECURSIVE FUNCTION factorial(i) RESULT(fac)
 
6
    IMPLICIT NONE
 
7
    INTEGER::fac
 
8
    INTEGER,INTENT(IN)::i
 
9
    INTEGER::init=0,j
 
10
    LOGICAL,DIMENSION(12)::lfacsave
 
11
    INTEGER,DIMENSION(12)::facsave
 
12
    SAVE facsave,init,lfacsave
 
13
    IF(init.EQ.0)THEN
 
14
       lfacsave(1:12)=.FALSE.
 
15
       init=1
 
16
    ENDIF
 
17
    IF(i.LT.0)THEN
 
18
       WRITE(*,*)"ERROR: i<0 in factorial with i=",i
 
19
       STOP
 
20
    ENDIF
 
21
    IF(i.GT.12)THEN
 
22
       WRITE(*,*)"ERROR: i > 12, please take long type (KIND=LINT) integer"
 
23
       STOP
 
24
    ENDIF
 
25
    IF(i.EQ.0)THEN
 
26
       fac=1
 
27
    ELSE
 
28
       IF(lfacsave(i))THEN
 
29
          fac=facsave(i)
 
30
       ELSE
 
31
          fac=i*factorial(i-1)
 
32
          facsave(i)=fac
 
33
          lfacsave(i)=.TRUE.
 
34
       ENDIF
 
35
    ENDIF
 
36
  END FUNCTION factorial
 
37
 
 
38
  SUBROUTINE PCL2PE(NLOOPLINE,PCL,PE)
 
39
    ! transfer the loop momenta q,q+p1,q+p1+p2,q+p1+p2+p3,... 
 
40
    ! to p1,p2,p3,p4,... 
 
41
    IMPLICIT NONE
 
42
    INTEGER,INTENT(IN)::NLOOPLINE
 
43
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(IN)::PCL
 
44
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(OUT)::PE
 
45
    INTEGER::i
 
46
    DO i=1,NLOOPLINE-1
 
47
       PE(i,0:3)=PCL(i+1,0:3)-PCL(i,0:3)
 
48
    ENDDO
 
49
    PE(NLOOPLINE,0:3)=PCL(1,0:3)-PCL(NLOOPLINE,0:3)
 
50
    RETURN
 
51
  END SUBROUTINE PCL2PE
 
52
 
 
53
  RECURSIVE SUBROUTINE calc_all_integers(n,ntot,i,sol)
 
54
    ! finds all the non-negative solutions to x1,...,xn
 
55
    ! that x1+x2+...+xn==i
 
56
    ! the number of solutions should be ntot=C(i+n-1)^(n-1)=(i+n-1)!/(n-1)!/i!
 
57
    ! it can be recycled for all phase space point
 
58
    IMPLICIT NONE
 
59
    INTEGER,INTENT(IN)::n,ntot,i
 
60
    INTEGER,DIMENSION(ntot,n),INTENT(OUT)::sol
 
61
    INTEGER::j,k,k0
 
62
    IF(i.EQ.0)THEN
 
63
       sol(1,1:n)=0
 
64
       RETURN
 
65
    ENDIF
 
66
    IF(n.EQ.1)THEN
 
67
       sol(1,1)=i
 
68
       RETURN
 
69
    ENDIF
 
70
    k0=0
 
71
    DO j=0,i
 
72
       k=xiarray_arg1(i-j,n-1)
 
73
       sol(k0+1:k0+k,1)=j
 
74
       CALL calc_all_integers(n-1,k,i-j,sol(k0+1:k0+k,2:n))
 
75
       k0=k0+k
 
76
    ENDDO
 
77
    RETURN
 
78
  END SUBROUTINE calc_all_integers
 
79
 
 
80
  SUBROUTINE all_Integers(n,ntot,i,sol,factor)
 
81
    ! finds all the non-negative solutions to x1,...,xn 
 
82
    ! that x1+x2+...+xn==i 
 
83
    ! the number of solutions should be ntot=C(i+n-1)^(n-1)=(i+n-1)!/(n-1)!/i!
 
84
    ! it can be recycled for all phase space point 
 
85
    IMPLICIT NONE
 
86
    INTEGER,INTENT(IN)::n,ntot,i
 
87
    INTEGER,DIMENSION(ntot,n),INTENT(OUT)::sol
 
88
    REAL(KIND(1d0)),DIMENSION(ntot),INTENT(OUT)::factor
 
89
    INTEGER::ifirst=0,j,jk,k,ntemptot
 
90
    INTEGER::maxxiarray,inum,nnum,maxxiarray1,maxxiarray2
 
91
    INTEGER::maxfactor_xiarray
 
92
    SAVE ifirst
 
93
    ! calculate xiarray_i_n first
 
94
    IF(ifirst.EQ.0)THEN
 
95
       maxxiarray=(MAXRANK_IREGI+1)*(MAXNLOOP_IREGI-1)
 
96
       IF(.NOT.ALLOCATED(xiarray))THEN
 
97
          ALLOCATE(xiarray(maxxiarray))
 
98
       ENDIF
 
99
       maxfactor_xiarray=(MAXRANK_IREGI+1)*(MAXNLOOP_IREGI-2)
 
100
       IF(.NOT.ALLOCATED(factor_xiarray))THEN
 
101
          ALLOCATE(factor_xiarray(MAXRANK_IREGI+2:maxfactor_xiarray+MAXRANK_IREGI+1))
 
102
       ENDIF
 
103
       ! x1+x2+...+xn==i
 
104
       ! C(i+n-1)^(n-1)=(i+n-1)!/(n-1)!/i!
 
105
       DO j=1,maxxiarray
 
106
          inum=MOD(j-1,MAXRANK_IREGI+1) ! i
 
107
          nnum=(j-1)/(MAXRANK_IREGI+1)+1 ! n
 
108
          maxxiarray1=xiarray_arg1(inum,nnum)
 
109
          maxxiarray2=nnum
 
110
          IF(.NOT.ALLOCATED(xiarray(j)%xiarray_i_n))THEN
 
111
             ALLOCATE(xiarray(j)%xiarray_i_n(maxxiarray1,maxxiarray2))
 
112
          ENDIF
 
113
          IF(nnum.EQ.1)THEN
 
114
             xiarray(j)%xiarray_i_n(1,1)=j-1
 
115
             CYCLE
 
116
          ENDIF
 
117
          CALL calc_all_integers(nnum,maxxiarray1,inum,&
 
118
               xiarray(j)%xiarray_i_n(1:maxxiarray1,1:maxxiarray2))
 
119
          IF(.NOT.ALLOCATED(factor_xiarray(j)%factor_xiarray_i_n))THEN
 
120
             ALLOCATE(factor_xiarray(j)%factor_xiarray_i_n(maxxiarray1))
 
121
          ENDIF
 
122
          DO jk=1,maxxiarray1
 
123
             factor_xiarray(j)%factor_xiarray_i_n(jk)=DBLE(factorial(inum))
 
124
             DO k=1,nnum
 
125
                factor_xiarray(j)%factor_xiarray_i_n(jk)=factor_xiarray(j)%factor_xiarray_i_n(jk)/&
 
126
                     DBLE(factorial(xiarray(j)%xiarray_i_n(jk,k)))
 
127
             ENDDO
 
128
          ENDDO
 
129
          ntot_xiarray(inum,nnum)=maxxiarray1
 
130
       ENDDO
 
131
       ifirst=1
 
132
    ENDIF
 
133
    IF(n.GT.MAXNLOOP_IREGI-1.OR.n.LT.1)THEN
 
134
       WRITE(*,100)"ERROR: n is out of range 1<=n<=",MAXNLOOP_IREGI-1," in all_integers"
 
135
       STOP
 
136
    ENDIF
 
137
    IF(i.GT.MAXRANK_IREGI.OR.i.LT.0)THEN
 
138
       WRITE(*,100)"ERROR: r is out of range 0<=r<=",MAXRANK_IREGI," in all_integers"
 
139
       STOP
 
140
    ENDIF
 
141
    IF(n.EQ.1.AND.ntot.NE.1)THEN
 
142
       WRITE(*,*)"ERROR:ntot should be 1 when n=1 in all_integers"
 
143
       STOP
 
144
    ENDIF
 
145
    IF(n.GE.2.AND.n.LE.MAXNLOOP_IREGI-1)THEN
 
146
       ! Make it work in MadLoop, otherwise it is wrong
 
147
       IF(ntot.NE.ntot_xiarray(i,n))THEN
 
148
          WRITE(*,*)"ERROR: ntot is not correct in all_integers"
 
149
          STOP
 
150
       ENDIF
 
151
    ENDIF
 
152
    j=(n-1)*(MAXRANK_IREGI+1)+i+1
 
153
    SELECT CASE(n)
 
154
       CASE(1)
 
155
          sol(1:ntot,1:n)=xiarray(j)%xiarray_i_n(1:1,1:n)
 
156
          factor(1)=1d0
 
157
       CASE DEFAULT
 
158
          sol(1:ntot,1:n)=xiarray(j)%xiarray_i_n(1:ntot,1:n)
 
159
          factor(1:ntot)=factor_xiarray(j)%factor_xiarray_i_n(1:ntot)
 
160
    END SELECT
 
161
    RETURN
 
162
100 FORMAT(2X,A31,I2,A16)
 
163
  END SUBROUTINE all_Integers
 
164
 
 
165
  SUBROUTINE calc_factorial_pair
 
166
    IMPLICIT NONE
 
167
    INTEGER::i,j
 
168
    factorial_pair(1:MAXINDICES_IREGI,0)=1d0
 
169
    DO i=1,MAXINDICES_IREGI
 
170
       DO j=1,MAXRANK_IREGI
 
171
          factorial_pair(i,j)=DBLE(factorial(i+j-1))&
 
172
               /DBLE(factorial(i-1))/DBLE(factorial(j))
 
173
       ENDDO
 
174
    ENDDO
 
175
    RETURN
 
176
  END SUBROUTINE calc_factorial_pair
 
177
 
 
178
  FUNCTION number_coefs_for_rank(rank) RESULT(res)
 
179
    IMPLICIT NONE
 
180
    INTEGER,INTENT(IN)::rank
 
181
    INTEGER::res,i
 
182
    IF(rank.LT.0)res=0
 
183
    res=0
 
184
    DO i=0,rank
 
185
       res=res+((3+i)*(2+i)*(1+i))/6
 
186
    ENDDO
 
187
    RETURN
 
188
  END FUNCTION number_coefs_for_rank
 
189
 
 
190
  SUBROUTINE timestamp
 
191
    IMPLICIT NONE    
 
192
    CHARACTER(len=8)::ampm
 
193
    INTEGER::d
 
194
    INTEGER::h
 
195
    INTEGER::m,n
 
196
    INTEGER::mm
 
197
    CHARACTER( len = 9 ),PARAMETER,DIMENSION(12) :: month = (/ &
 
198
         'January  ', 'February ', 'March    ', 'April    ', &
 
199
         'May      ', 'June     ', 'July     ', 'August   ', &
 
200
         'September', 'October  ', 'November ', 'December ' /)
 
201
    INTEGER::dn
 
202
    INTEGER::s
 
203
    INTEGER,DIMENSION(8)::values
 
204
    INTEGER::y
 
205
    CALL date_and_time(values=values)
 
206
    y = values(1)
 
207
    m = values(2)
 
208
    d = values(3)
 
209
    h = values(5)
 
210
    n = values(6)
 
211
    s = values(7)
 
212
    mm = values(8)
 
213
    IF( h.LT.12 )THEN
 
214
       ampm = 'AM'
 
215
    ELSEIF( h.EQ.12 )THEN
 
216
       IF( n.EQ.0 .AND. s.EQ.0 )THEN
 
217
          ampm = 'Noon'
 
218
       ELSE
 
219
          ampm = 'PM'
 
220
       ENDIF
 
221
    ELSE
 
222
       h = h - 12
 
223
       IF( h.LT.12 )THEN
 
224
          ampm = 'PM'
 
225
       ELSEIF( h.EQ.12 )THEN
 
226
          IF(n.EQ.0.AND.s.EQ.0)THEN
 
227
             ampm = 'Midnight'
 
228
          ELSE
 
229
             ampm = 'AM'
 
230
          ENDIF
 
231
       ENDIF
 
232
    ENDIF
 
233
    WRITE( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
 
234
         d, TRIM( month(m) ), y, h, ':', n, ':', s, '.', mm, TRIM( ampm )  
 
235
    RETURN
 
236
  END SUBROUTINE timestamp
 
237
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
238
!                                                                               !
 
239
! FOLLOWING SUBROUTINES/FUNCTIONS ARE ONLY FOR THE REORDERING THE TIR COEFS     !
 
240
!                                                                               !
 
241
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
242
  SUBROUTINE SORT_IREGICOEFS(RANK,NLOOPCOEFS,OLDCOEFS,NEWCOEFS)
 
243
    !
 
244
    ! CONVERT THE OUTPUT OF IREGI TO THAT OF (NEW) MADLOOP 
 
245
    !
 
246
    ! THE NEW OUTPUT OF COEFS FROM MADLOOP IS
 
247
    ! RANK=0: (,)
 
248
    ! RANK=1: (0,),(1,),(2,),(3,)
 
249
    ! RANK=2: (0,0),(0,1),(1,1),(0,2),(1,2),(2,2),(0,3),(1,3),(2,3),(3,3)
 
250
    ! ...
 
251
    ! THE OLD OUTPUT OF COEFS FROM MADLOOP IS
 
252
    ! RANK=0: (,)
 
253
    ! RANK=1: (0,),(1,),(2,),(3,)
 
254
    ! RANK=2: (0,0),(0,1),(0,2),(0,3),(1,1),(2,1),(3,1),(2,2),(2,3),(3,3)
 
255
    ! ...
 
256
    !
 
257
    ! ARGUMENTS
 
258
    IMPLICIT NONE
 
259
    INTEGER,INTENT(IN)::RANK,NLOOPCOEFS
 
260
    COMPLEX(KIND(1d0)),DIMENSION(0:NLOOPCOEFS-1,3),INTENT(IN)::OLDCOEFS
 
261
    COMPLEX(KIND(1d0)),DIMENSION(0:NLOOPCOEFS-1,3),INTENT(OUT)::NEWCOEFS
 
262
    !
 
263
    ! LOCAL VARIABLES
 
264
    !
 
265
    INTEGER::I
 
266
    ! MAX RANK SET AS 10
 
267
    ! Sum((3+r)*(2+r)*(1+r)/6,{r,0,10})=1001
 
268
    INTEGER,PARAMETER::LOOPMAXCOEFS_IREGI=1001
 
269
    INTEGER,DIMENSION(0:LOOPMAXCOEFS_IREGI-1)::POS
 
270
    SAVE POS
 
271
    LOGICAL::INIT=.TRUE.
 
272
    SAVE INIT
 
273
    ! ----------
 
274
    ! BEGIN CODE
 
275
    ! ----------
 
276
            
 
277
    IF(INIT)THEN
 
278
       IF(NLOOPCOEFS.GT.LOOPMAXCOEFS_IREGI)&
 
279
            STOP "ERROR:LOOPMAXCOEFS_IREGI IS TOO SMALL!!!"
 
280
       INIT=.FALSE.
 
281
       ! ASSIGN THE POSITION OF POS FOR SWAP
 
282
       CALL ASSIGN_PJPOS(POS)
 
283
    ENDIF
 
284
    
 
285
    DO I=0,NLOOPCOEFS-1
 
286
    !   NEWCOEFS(I,1:3)=OLDCOEFS(POS(I),1:3)
 
287
       NEWCOEFS(POS(I),1:3)=OLDCOEFS(I,1:3)
 
288
    ENDDO
 
289
    
 
290
  END SUBROUTINE SORT_IREGICOEFS
 
291
 
 
292
  SUBROUTINE ASSIGN_PJPOS(POS)
 
293
    !
 
294
    ! ASSIGN THE POSITION OF POS FOR SWAP
 
295
    !
 
296
    !
 
297
    IMPLICIT NONE
 
298
    !
 
299
    ! CONSTANS
 
300
    !
 
301
    ! MAX RANK SET AS 10
 
302
    ! Sum((3+r)*(2+r)*(1+r)/6,{r,0,10})=1001 
 
303
    INTEGER,PARAMETER::LOOPMAXCOEFS_IREGI=1001,MAXRANK=10
 
304
    ! 
 
305
    ! ARGUMENTS
 
306
    ! 
 
307
    INTEGER,DIMENSION(0:LOOPMAXCOEFS_IREGI-1),INTENT(OUT)::POS
 
308
    !
 
309
    ! LOCAL VARIABLES
 
310
    !
 
311
    INTEGER::I,J,K,SHIFT,DN
 
312
    INTEGER,DIMENSION(MAXRANK)::POSINDEX,PJPOSINDEX
 
313
    ! ----------
 
314
    ! BEGIN CODE
 
315
    ! ----------
 
316
    POS(0)=0
 
317
    DO I=1,4
 
318
       POS(I)=I
 
319
    ENDDO
 
320
    SHIFT=4
 
321
    DO J=2,MAXRANK
 
322
       DN=(J+3)*(J+2)*(J+1)/6
 
323
       POSINDEX(1:MAXRANK)=0
 
324
       PJPOSINDEX(1:MAXRANK)=0
 
325
       DO I=1,DN
 
326
          IF(I.GT.1)CALL NEXTINDEX(J,POSINDEX)
 
327
          CALL CONVERT_PJPOSINDEX(J,POSINDEX,PJPOSINDEX)
 
328
          K=DN-QPOLYPOS(J,PJPOSINDEX)+1+SHIFT
 
329
          !POS(K)=I+SHIFT
 
330
          POS(I+SHIFT)=K
 
331
       ENDDO
 
332
       SHIFT=SHIFT+DN
 
333
    ENDDO
 
334
    
 
335
  END SUBROUTINE ASSIGN_PJPOS
 
336
 
 
337
  SUBROUTINE NEXTINDEX(RANK,POSINDEX)
 
338
    !
 
339
    ! CALL FOR THE NEXT INDEX
 
340
    !
 
341
    IMPLICIT NONE
 
342
    !
 
343
    ! CONSTANTS
 
344
    !
 
345
    INTEGER,PARAMETER::MAXRANK=10
 
346
    !
 
347
    ! ARGUMENTS
 
348
    !
 
349
    INTEGER,INTENT(IN)::RANK
 
350
    INTEGER,DIMENSION(MAXRANK),INTENT(INOUT)::POSINDEX
 
351
    !
 
352
    ! LOCAL VARIABLES
 
353
    !
 
354
    INTEGER::I
 
355
    ! ----------
 
356
    ! BEGIN CODE
 
357
    ! ----------
 
358
    DO I=1,RANK
 
359
       POSINDEX(I)=POSINDEX(I)+1
 
360
       IF(POSINDEX(I).GT.3)THEN
 
361
          POSINDEX(I)=0
 
362
          IF(I.EQ.RANK)THEN
 
363
             RETURN
 
364
          ENDIF
 
365
       ELSE
 
366
          IF(I.GT.1)THEN
 
367
             POSINDEX(1:I-1)=POSINDEX(I)
 
368
          ENDIF
 
369
          RETURN
 
370
       ENDIF
 
371
    ENDDO
 
372
    
 
373
  END SUBROUTINE NEXTINDEX
 
374
 
 
375
  SUBROUTINE CONVERT_PJPOSINDEX(RANK,POSINDEX,PJPOSINDEX)
 
376
    !
 
377
    ! CONVERT POSINDEX TO PJPOSINDEX
 
378
    !
 
379
    IMPLICIT NONE
 
380
    !
 
381
    ! CONSTANTS
 
382
    !
 
383
    INTEGER,PARAMETER::MAXRANK=10
 
384
    !
 
385
    ! ARGUMENTS
 
386
    !
 
387
    INTEGER,INTENT(IN)::RANK
 
388
    INTEGER,DIMENSION(MAXRANK),INTENT(IN)::POSINDEX
 
389
    INTEGER,DIMENSION(MAXRANK),INTENT(OUT)::PJPOSINDEX
 
390
    !
 
391
    ! LOCAL VARIABLES
 
392
    !
 
393
    INTEGER::I
 
394
    ! ----------
 
395
    ! BEGIN CODE
 
396
    ! ----------
 
397
    DO I=1,RANK
 
398
       PJPOSINDEX(RANK+1-I)=3-POSINDEX(I)
 
399
    ENDDO
 
400
    RETURN
 
401
  END SUBROUTINE CONVERT_PJPOSINDEX
 
402
 
 
403
  FUNCTION QPOLYPOS(RANK,POSINDEX)
 
404
    !
 
405
    ! COMPUTATION THE RELATIVE POSITION OF INDEX WITH RANK
 
406
    ! IN THE *OLD* MADLOOP CONVENTION
 
407
    !
 
408
    IMPLICIT NONE
 
409
    !
 
410
    ! CONSTANTS
 
411
    !
 
412
    INTEGER,PARAMETER::MAXRANK=10
 
413
    !
 
414
    ! ARGUMENTS
 
415
    !
 
416
    INTEGER,INTENT(IN)::RANK
 
417
    INTEGER,DIMENSION(MAXRANK),INTENT(IN)::POSINDEX
 
418
    INTEGER::QPOLYPOS
 
419
    !
 
420
    ! LOCAL VARIABLES
 
421
    !
 
422
    INTEGER::I,J,IMIN
 
423
    ! ----------
 
424
    ! BEGIN CODE
 
425
    ! ----------
 
426
 
 
427
    IF(RANK.EQ.0)THEN
 
428
       QPOLYPOS=1
 
429
       RETURN
 
430
    ENDIF
 
431
    
 
432
    IF(RANK.EQ.1)THEN
 
433
       QPOLYPOS=POSINDEX(1)+1
 
434
       RETURN
 
435
    ENDIF
 
436
    
 
437
    QPOLYPOS=POSINDEX(1)-POSINDEX(2)+1
 
438
    DO I=2,RANK
 
439
       IF(I.EQ.RANK)THEN
 
440
          IMIN=0
 
441
       ELSE
 
442
          IMIN=POSINDEX(I+1)
 
443
       ENDIF
 
444
       DO J=IMIN,POSINDEX(I)-1
 
445
          QPOLYPOS=QPOLYPOS+QPOLYNUMBER(J,I-1)
 
446
       ENDDO
 
447
    ENDDO
 
448
    RETURN
 
449
  END FUNCTION QPOLYPOS
 
450
 
 
451
  FUNCTION NEWQPOLYPOS(RANK,POSINDEX)
 
452
    !
 
453
    ! COMPUTATION THE RELATIVE POSITION OF INDEX WITH RANK
 
454
    ! IN THE *NEW* MADLOOP CONVENTION
 
455
    !
 
456
    IMPLICIT NONE
 
457
    !
 
458
    ! CONSTANTS
 
459
    !
 
460
    INTEGER,PARAMETER::MAXRANK=10
 
461
    !
 
462
    ! ARGUMENTS
 
463
    !
 
464
    INTEGER,INTENT(IN)::RANK
 
465
    INTEGER,DIMENSION(MAXRANK),INTENT(IN)::POSINDEX
 
466
    INTEGER::NEWQPOLYPOS
 
467
    !
 
468
    ! LOCAL VARIABLES
 
469
    !
 
470
    INTEGER::I,J,IMIN
 
471
    ! ----------
 
472
    ! BEGIN CODE
 
473
    ! ----------                                                                                                                                                                                            
 
474
 
 
475
    IF(RANK.EQ.0)THEN
 
476
       NEWQPOLYPOS=1
 
477
       RETURN
 
478
    ENDIF
 
479
 
 
480
    IF(RANK.EQ.1)THEN
 
481
       NEWQPOLYPOS=POSINDEX(1)+1
 
482
       RETURN
 
483
    ENDIF
 
484
 
 
485
    NEWQPOLYPOS=1
 
486
    DO I=1,RANK
 
487
       IF(POSINDEX(RANK-I+1).EQ.0)THEN
 
488
          IMIN=0
 
489
       ELSE
 
490
          ! Eq.(3.4.6) in Valentin's PhD Thesis
 
491
          IMIN=factorial(POSINDEX(RANK-I+1)+I-1)/factorial(I)&
 
492
               /factorial(POSINDEX(RANK-I+1)-1)
 
493
       ENDIF
 
494
       NEWQPOLYPOS=NEWQPOLYPOS+IMIN
 
495
    ENDDO
 
496
    RETURN
 
497
  END FUNCTION NEWQPOLYPOS
 
498
 
 
499
  FUNCTION QPOLYNUMBER(I,RANK)
 
500
    !
 
501
    ! THE INDEPENDENT NUMBER OF Q POLY WITH \MU=I,...,3 AND RANK
 
502
    !
 
503
    IMPLICIT NONE
 
504
    !
 
505
    ! CONSTANTS
 
506
    !
 
507
    !
 
508
    ! ARGUMENTS
 
509
    !
 
510
    INTEGER,INTENT(IN)::I,RANK
 
511
    INTEGER::QPOLYNUMBER
 
512
    !
 
513
    ! LOCAL VARIABLES
 
514
    !
 
515
    ! ----------
 
516
    ! BEGIN CODE
 
517
    ! ----------
 
518
    SELECT CASE(I)
 
519
    CASE(0)
 
520
       QPOLYNUMBER=(3+RANK)*(2+RANK)*(1+RANK)/6
 
521
    CASE(1)
 
522
       QPOLYNUMBER=(2+RANK)*(1+RANK)/2
 
523
    CASE(2)
 
524
       QPOLYNUMBER=(1+RANK)
 
525
    CASE(3)
 
526
       QPOLYNUMBER=1
 
527
    CASE DEFAULT
 
528
       STOP 'I must be >= 0 and <=3 in QPOLYNUMBER.'
 
529
    END SELECT
 
530
    RETURN
 
531
  END FUNCTION QPOLYNUMBER
 
532
 
 
533
  FUNCTION POS2RANK(ipos)
 
534
    IMPLICIT NONE
 
535
    INTEGER,INTENT(IN)::ipos
 
536
    INTEGER::POS2RANK
 
537
    LOGICAL::INIT=.TRUE.
 
538
    ! NOW ITS MAX RANK IS 10
 
539
    INTEGER,PARAMETER::MAXRANK=10
 
540
    INTEGER,DIMENSION(0:MAXRANK)::IRANGE
 
541
    INTEGER::I
 
542
    SAVE IRANGE,INIT
 
543
    IF(INIT)THEN
 
544
       IRANGE(0)=0
 
545
       DO I=1,MAXRANK
 
546
          IRANGE(I)=IRANGE(I-1)+(I+3)*(I+2)*(I+1)/6
 
547
       ENDDO
 
548
       INIT=.FALSE.
 
549
    ENDIF
 
550
    IF(ipos.EQ.0)THEN
 
551
       POS2RANK=0
 
552
       RETURN
 
553
    ENDIF
 
554
    DO I=1,MAXRANK
 
555
       IF(ipos.LE.IRANGE(I))THEN
 
556
          POS2RANK=I
 
557
          RETURN
 
558
       ENDIF
 
559
    ENDDO
 
560
    WRITE(*,*)"ERROR in POS2RANK,ipos=",ipos
 
561
    STOP
 
562
  END FUNCTION POS2RANK
 
563
 
 
564
  FUNCTION RANK_NUM(MAXRANK)
 
565
    IMPLICIT NONE
 
566
    INTEGER::RANK_NUM
 
567
    INTEGER,INTENT(IN)::MAXRANK
 
568
    INTEGER::I
 
569
    RANK_NUM=0
 
570
    IF(MAXRANK.LT.0)RETURN
 
571
    DO I=0,MAXRANK
 
572
       RANK_NUM=RANK_NUM+(I+3)*(I+2)*(I+1)/6
 
573
    ENDDO
 
574
    RETURN
 
575
  END FUNCTION RANK_NUM
 
576
 
 
577
  SUBROUTINE SHIFT_MOM(MOM,NCOEFS,TICOEFS_SAVE,TICOEFS)
 
578
    ! SHIFT MOMENTUM MOM
 
579
    IMPLICIT NONE
 
580
    !
 
581
    ! CONSTANS
 
582
    !
 
583
    ! MAX RANK SET AS 10
 
584
    ! Sum((3+r)*(2+r)*(1+r)/6,{r,0,10})=1001 --> LOOPMAXCOEFS_IREGI
 
585
    ! MAX OF (i1+1)*(i2+1)*(i3+1)*(i4+1) with i1+i2+i3+i4<=MAXRANK = 144 --> MAXSPLIT
 
586
    INTEGER,PARAMETER::LOOPMAXCOEFS_IREGI=1001,MAXRANK=10,MAXSPLIT=144
 
587
    REAL(KIND(1d0)),DIMENSION(0:3),INTENT(IN)::MOM
 
588
    INTEGER,INTENT(IN)::NCOEFS
 
589
    COMPLEX(KIND(1d0)),DIMENSION(0:NCOEFS-1,1:4),INTENT(IN)::TICOEFS_SAVE
 
590
    COMPLEX(KIND(1d0)),DIMENSION(0:NCOEFS-1,1:4),INTENT(OUT)::TICOEFS
 
591
    LOGICAL::INIT=.TRUE.
 
592
    INTEGER,DIMENSION(0:LOOPMAXCOEFS_IREGI-1)::NTERM
 
593
    REAL(KIND(1d0)),DIMENSION(0:LOOPMAXCOEFS_IREGI-1,MAXSPLIT)::SPLIT_FACTOR
 
594
    ! i=0-3 -> number of lorentz indices =i
 
595
    ! i=4 -> the position of TICOEFS that should be multiplied
 
596
    ! e.g. (q+MOM)^mu1(q+MOM)^mu2...(q+MOM)^muk
 
597
    ! it can be splitted into q^mu1...q^muk+q^mu1...q^mui Mom^mui+1 ...q^muk ...
 
598
    ! for some specific term configurateion mu1=0,mu2=0,...muk=3
 
599
    ! NTERM stores the total number of terms
 
600
    ! SPLIT_FACTOR is the factor for each term
 
601
    ! SPLIT_INFO(J,K,0-3) is the number of 0,1,2,3 with MOM^mu
 
602
    ! SPLIT_INFO(J,K,4) is the corresponding q^mu (remaining lorentz indices)
 
603
    INTEGER,DIMENSION(0:LOOPMAXCOEFS_IREGI-1,MAXSPLIT,0:4)::SPLIT_INFO  
 
604
    SAVE INIT,NTERM,SPLIT_FACTOR,SPLIT_INFO
 
605
    INTEGER::I,J,K,L,NINDEX,r,incr
 
606
    INTEGER::i1,i2,i3,i4,nsplit,itit,rte,ND
 
607
    REAL(KIND(1d0))::temp
 
608
    INTEGER,DIMENSION(MAXRANK)::POSINDEX
 
609
    IF(INIT)THEN
 
610
       incr=0
 
611
       ND=0
 
612
       DO r=0,MAXRANK
 
613
          DO I=r,0,-1
 
614
             DO J=r-I,0,-1
 
615
                DO K=r-I-J,0,-1
 
616
                   L=r-I-J-K
 
617
                   ! exhaust all of the subsets
 
618
                   nsplit=0
 
619
                   NTERM(incr)=(I+1)*(J+1)*(K+1)*(L+1)
 
620
                   DO i1=0,I
 
621
                      !DO itit=1,i1
 
622
                      !   POSINDEX(itit)=0
 
623
                      !ENDDO
 
624
                      DO i2=0,J
 
625
                         !DO itit=i1+1,i1+i2
 
626
                         !   POSINDEX(itit)=1
 
627
                         !ENDDO
 
628
                         DO i3=0,K
 
629
                            !DO itit=i1+i2+1,i1+i2+i3
 
630
                            !   POSINDEX(itit)=2
 
631
                            !ENDDO
 
632
                            DO i4=0,L
 
633
                               !DO itit=i1+i2+i3+1,i1+i2+i3+i4
 
634
                               !   POSINDEX(itit)=3
 
635
                               !ENDDO
 
636
                               nsplit=nsplit+1
 
637
                               rte=i1+i2+i3+i4
 
638
                               DO itit=1,i1
 
639
                                  POSINDEX(rte-itit+1)=0
 
640
                               ENDDO
 
641
                               DO itit=i1+1,i1+i2
 
642
                                  POSINDEX(rte-itit+1)=1
 
643
                               ENDDO
 
644
                               DO itit=i1+i2+1,i1+i2+i3
 
645
                                  POSINDEX(rte-itit+1)=2
 
646
                               ENDDO
 
647
                               DO itit=i1+i2+i3+1,i1+i2+i3+i4
 
648
                                  POSINDEX(rte-itit+1)=3
 
649
                               ENDDO
 
650
                               SPLIT_INFO(incr,nsplit,4)=QPOLYPOS(rte,POSINDEX)+RANK_NUM(rte-1)-1
 
651
                               SPLIT_INFO(incr,nsplit,0)=I-i1
 
652
                               SPLIT_INFO(incr,nsplit,1)=J-i2
 
653
                               SPLIT_INFO(incr,nsplit,2)=K-i3
 
654
                               SPLIT_INFO(incr,nsplit,3)=L-i4
 
655
                               SPLIT_FACTOR(incr,nsplit)=DBLE(factorial(I))&
 
656
                                    /DBLE(factorial(i1)*factorial(I-i1))
 
657
                               SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)*&
 
658
                                    DBLE(factorial(J))/DBLE(factorial(i2)*factorial(J-i2))
 
659
                               SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)*&
 
660
                                    DBLE(factorial(K))/DBLE(factorial(i3)*factorial(K-i3))
 
661
                               SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)*&
 
662
                                    DBLE(factorial(L))/DBLE(factorial(i4)*factorial(L-i4))
 
663
                               IF(.NOT.ML5_CONVENTION)THEN
 
664
                                  SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)*&
 
665
                                       DBLE(factorial(r))/DBLE(factorial(I)*factorial(J))&
 
666
                                       /DBLE(factorial(K)*factorial(L))
 
667
                                  SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)/&
 
668
                                       DBLE(factorial(rte))*DBLE(factorial(i1))*DBLE(factorial(i2))*&
 
669
                                       DBLE(factorial(i3))*DBLE(factorial(i4))
 
670
                               ENDIF
 
671
                            ENDDO
 
672
                         ENDDO
 
673
                      ENDDO
 
674
                   ENDDO
 
675
                   incr=incr+1
 
676
                ENDDO
 
677
             ENDDO
 
678
          ENDDO
 
679
          ND=ND+(3+r)*(2+r)*(1+r)/6
 
680
       ENDDO
 
681
       INIT=.FALSE.
 
682
    ENDIF
 
683
    DO I=0,NCOEFS-1
 
684
       TICOEFS(I,1:4)=DCMPLX(0d0)
 
685
       DO J=1,NTERM(I)
 
686
          temp=SPLIT_FACTOR(I,J)
 
687
          DO K=0,3
 
688
             NINDEX=SPLIT_INFO(I,J,K)
 
689
             IF(NINDEX.GT.0)THEN
 
690
                temp=temp*MOM(K)**NINDEX
 
691
             ENDIF
 
692
          ENDDO
 
693
          TICOEFS(I,1:4)=TICOEFS(I,1:4)&
 
694
               +TICOEFS_SAVE(SPLIT_INFO(I,J,4),1:4)*temp
 
695
       ENDDO
 
696
    ENDDO
 
697
    RETURN
 
698
  END SUBROUTINE SHIFT_MOM
 
699
 
 
700
  FUNCTION xiarray_arg1(i,n)
 
701
    IMPLICIT NONE
 
702
    INTEGER,INTENT(IN)::i,n
 
703
    INTEGER::xiarray_arg1
 
704
    INTEGER::imax,imin,j
 
705
    !INTEGER(KIND=LINT)::itemp
 
706
    IF(i+n-1.GT.12)THEN
 
707
       imax=MAX(n-1,i)
 
708
       imin=MIN(n-1,i)
 
709
       xiarray_arg1=1
 
710
       DO j=imax+1,i+n-1
 
711
          xiarray_arg1=xiarray_arg1*j
 
712
       ENDDO
 
713
       IF(imin.LE.12)THEN
 
714
          xiarray_arg1=xiarray_arg1/factorial(imin)
 
715
       ELSE
 
716
          DO j=1,imin
 
717
             xiarray_arg1=xiarray_arg1/j
 
718
          ENDDO
 
719
       ENDIF
 
720
    ELSE
 
721
       xiarray_arg1=factorial(i+n-1)/factorial(n-1)/factorial(i)
 
722
    ENDIF
 
723
    RETURN
 
724
  END FUNCTION xiarray_arg1
 
725
END MODULE FUNLIB