JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
reductionTN.F90
Go to the documentation of this file.
1 !!
2 !! File reductionTN.F90 is part of COLLIER
3 !! - A Complex One-Loop Library In Extended Regularizations
4 !!
5 !! Copyright (C) 2015, 2016 Ansgar Denner, Stefan Dittmaier, Lars Hofer
6 !!
7 !! COLLIER is licenced under the GNU GPL version 3, see COPYING for details.
8 !!
9 
10 !#define TNtest
11 
12 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13 !
14 ! ************************
15 ! * module reductionTN *
16 ! * by Lars Hofer *
17 ! ************************
18 !
19 ! functions and subroutines:
20 ! CalcTN
21 !
22 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23 
24 
25 
26 
27 
28 module reductiontn
29 
30  use reductionefg
31 
32  implicit none
33 
34 ! integer, allocatable :: AddToCInd(:,:,:,:), DropCInd(:,:,:,:)
35 
36 contains
37 
38 
39  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
40  ! subroutine CalcTN(TN,TNuv,MomInv,masses2,rmax,id,TNerr,TNerr2)
41  !
42  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
43 
44  subroutine calctn(TN,TNuv,MomInv,masses2,N,rmax,id,TNerr,TNerr2)
45 
46  integer, intent(in) :: N,rmax,id
47  double complex, intent(in) :: MomInv(BinomTable(2,N)), masses2(0:N-1)
48  double complex, intent(inout) :: TN(NCoefs(rmax,N))
49  double complex, intent(inout) :: TNuv(NCoefs(rmax,N))
50  double precision, intent(out) :: TNerr(0:rmax),TNerr2(0:rmax)
51  double complex, allocatable :: TNaux(:,:,:),TNuvaux(:,:,:)
52  double complex :: B(0:rmax,0:rmax), Buv(0:rmax,0:rmax)
53  double complex :: C(0:rmax,0:rmax,0:rmax), Cuv(0:rmax,0:rmax,0:rmax)
54  double complex :: D(0:rmax,0:rmax,0:rmax,0:rmax), Duv(0:rmax,0:rmax,0:rmax,0:rmax)
55  double complex :: E(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
56  double complex :: Euv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
57  double complex :: x(BinomTable(2,N)+N)
58  double complex, allocatable :: fct(:)
59  integer :: r,n0,n1,n2,n3,n4,i,rank,bino,bino_m2,cnt
60  logical :: nocalc,wrica
61 
62 
63  select case (n)
64 
65  case(1)
66  call calca(tn,tnuv,masses2(0),rmax,tnerr)
67  tnerr2=tnerr
68  return
69 
70  case(2)
71  call calcb(b,buv,mominv(1),masses2(0),masses2(1),rmax,id,tnerr)
72 
73  cnt = 1
74  do r=0,rmax
75  do n0=r/2,0,-1
76  n1 = r-2*n0
77  tn(cnt) = b(n0,n1)
78  cnt = cnt+1
79  end do
80  end do
81  cnt = 1
82  do r=0,rmax
83  do n0=r/2,0,-1
84  n1 = r-2*n0
85  tnuv(cnt) = buv(n0,n1)
86  cnt = cnt+1
87  end do
88  end do
89  tnerr2=tnerr
90  return
91 
92  case(3)
93  call calcc(c,cuv,mominv(1),mominv(2),mominv(3), &
94  masses2(0),masses2(1),masses2(2),rmax,id,tnerr,tnerr2)
95 
96  cnt = 1
97  do r=0,rmax
98  do n0=r/2,0,-1
99  do n1=r-2*n0,0,-1
100  n2 = r-2*n0-n1
101  tn(cnt) = c(n0,n1,n2)
102  cnt = cnt+1
103  end do
104  end do
105  end do
106  tnuv = 0d0
107  do r=2,rmax
108  cnt = ncoefs(r-1,3)+1
109  do n0=r/2,1,-1
110  do n1=r-2*n0,0,-1
111  n2 = r-2*n0-n1
112  tnuv(cnt) = cuv(n0,n1,n2)
113  cnt = cnt+1
114  end do
115  end do
116  end do
117 
118 ! write(*,*) 'CalcTN 3 TNerr ',TNerr
119 ! write(*,*) 'CalcTN 3 TNerr2',TNerr2
120 
121  return
122 
123 
124 
125  case(4)
126  call calcd(d,duv,mominv(1),mominv(2),mominv(3),mominv(4),mominv(5),mominv(6), &
127  masses2(0),masses2(1),masses2(2),masses2(3),rmax,id,tnerr,tnerr2)
128 
129  cnt = 1
130  do r=0,rmax
131  do n0=r/2,0,-1
132  do n1=r-2*n0,0,-1
133  do n2=r-2*n0-n1,0,-1
134  n3 = r-2*n0-n1-n2
135  tn(cnt) = d(n0,n1,n2,n3)
136  cnt = cnt+1
137  end do
138  end do
139  end do
140  end do
141  tnuv = 0d0
142  do r=4,rmax
143  cnt = ncoefs(r-1,4)+1
144  do n0=r/2,2,-1
145  do n1=r-2*n0,0,-1
146  do n2=r-2*n0-n1,0,-1
147  n3 = r-2*n0-n1-n2
148  tnuv(cnt) = duv(n0,n1,n2,n3)
149  cnt = cnt+1
150  end do
151  end do
152  end do
153  end do
154 
155 ! write(*,*) 'CalcTN 4 TNerr ',TNerr
156 ! write(*,*) 'CalcTN 4 TNerr2',TNerr2
157 
158  return
159 
160  case(5)
161  call calce(e,euv,mominv(1),mominv(2),mominv(3),mominv(4),mominv(5), &
162  mominv(6),mominv(7),mominv(8),mominv(9),mominv(10), &
163  masses2(0),masses2(1),masses2(2),masses2(3),masses2(4) &
164  ,rmax,id,tnerr,tnerr2)
165 
166  cnt = 1
167  do r=0,rmax
168  do n0=r/2,0,-1
169  do n1=r-2*n0,0,-1
170  do n2=r-2*n0-n1,0,-1
171  do n3=r-2*n0-n1-n2,0,-1
172  n4 = r-2*n0-n1-n2-n3
173  tn(cnt) = e(n0,n1,n2,n3,n4)
174  cnt = cnt+1
175  end do
176  end do
177  end do
178  end do
179  end do
180  tnuv = 0d0
181  do r=6,rmax
182  cnt = ncoefs(r-1,5)+1
183  do n0=r/2,3,-1
184  do n1=r-2*n0,0,-1
185  do n2=r-2*n0-n1,0,-1
186  do n3=r-2*n0-n1-n2,0,-1
187  n4 = r-2*n0-n1-n2-n3
188  tnuv(cnt) = euv(n0,n1,n2,n3,n4)
189  cnt = cnt+1
190  end do
191  end do
192  end do
193  end do
194  end do
195  return
196 
197  case default
198  allocate(tnaux(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax))
199  allocate(tnuvaux(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax))
200  call calctnint(tnaux,tnuvaux,mominv,masses2,n,rmax,id,tnerr,tnerr2)
201 
202  cnt = 1
203  do r=0,rmax
204  do n0=r/2,0,-1
205  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
206  tn(cnt) = tnaux(i,n0,r)
207  cnt = cnt+1
208  end do
209  end do
210  end do
211  tnuv = 0d0
212  do r=2*n-4,rmax
213  cnt = ncoefs(r-1,n)+1
214  do n0=r/2,n-2,-1
215  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
216  tnuv(cnt) = tnuvaux(i,n0,r)
217  cnt = cnt+1
218  end do
219  end do
220  end do
221  return
222 
223  end select
224 
225 
226  end subroutine calctn
227 
228 
229 
230 
231 
232  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
233  ! subroutine CalcTNint(TN,TNuv,MomInv,masses2,N,rmax,id,TNerr,TNerr2)
234  !
235  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
236 
237  recursive subroutine calctnint(TN,TNuv,MomInv,masses2,N,rmax,id,TNerr,TNerr2)
238 
239  integer, intent(in) :: n,rmax,id
240  double complex, intent(in) :: mominv(binomtable(2,n)), masses2(0:n-1)
241  double complex, intent(inout) :: tn(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax)
242  double complex, intent(inout) :: tnuv(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax)
243  double precision, intent(out) :: tnerr(0:rmax),tnerr2(0:rmax)
244  double complex, allocatable :: tnaux(:,:,:), tnuvaux(:,:,:)
245  double complex :: e(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
246  double complex :: euv(0:rmax/2,0:rmax,0:rmax,0:rmax,0:rmax)
247  double precision, allocatable :: tnerraux(:),tnerr2aux(:)
248  double complex :: x(binomtable(2,n)+n)
249  double complex, allocatable :: fct(:)
250  integer :: r,n0,n1,n2,n3,n4,i,rank,bino,bino_m2,cnt
251  logical :: nocalc,wrica,noten
252 
253 
254 
255  if (n.eq.5) then
256  call calce(e,euv,mominv(1),mominv(2),mominv(3),mominv(4),mominv(5), &
257  mominv(6),mominv(7),mominv(8),mominv(9),mominv(10), &
258  masses2(0),masses2(1),masses2(2),masses2(3),masses2(4),rmax,id,tnerr,tnerr2)
259 
260  do r=0,rmax
261  do n0=r/2,0,-1
262  cnt = 1
263  do n1=r-2*n0,0,-1
264  do n2=r-2*n0-n1,0,-1
265  do n3=r-2*n0-n1-n2,0,-1
266  n4 = r-2*n0-n1-n2-n3
267  tn(cnt,n0,r) = e(n0,n1,n2,n3,n4)
268  cnt = cnt+1
269  end do
270  end do
271  end do
272  end do
273  end do
274  do r=6,rmax
275  do n0=r/2,3,-1
276  cnt = 1
277  do n1=r-2*n0,0,-1
278  do n2=r-2*n0-n1,0,-1
279  do n3=r-2*n0-n1-n2,0,-1
280  n4 = r-2*n0-n1-n2-n3
281  tnuv(cnt,n0,r) = euv(n0,n1,n2,n3,n4)
282  cnt = cnt+1
283  end do
284  end do
285  end do
286  end do
287  end do
288  return
289  end if
290 
291  if ((use_cache_system).and.(tencache.gt.n)) then
292  if ((ncache.gt.0).and.(ncache.le.ncache_max)) then
293  bino = binomtable(2,n)
294  x(1:bino) = mominv
295  x(bino+1:bino+n) = masses2
296  rank = rmax
297 
298  if(rmax.ge.2*n-4) then
299  allocate(fct(ncoefs(rmax,n)+ncoefs(rmax-2*n+4,n)+2*(rmax+1)))
300  call readcache(fct,ncoefs(rmax,n)+ncoefs(rmax-2*n+4,n)+2*(rmax+1),x,bino+n,1,id,n,rank,nocalc,wrica)
301  else
302  allocate(fct(ncoefs(rmax,n)+2*(rmax+1)))
303  call readcache(fct,ncoefs(rmax,n)+2*(rmax+1),x,bino+n,1,id,n,rank,nocalc,wrica)
304  end if
305 
306  if(nocalc)then
307  cnt = 0
308  do r=0,rmax
309  do n0=r/2,0,-1
310  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
311  cnt = cnt+1
312  tn(i,n0,r) = fct(cnt)
313  end do
314  end do
315  do n0=r/2,n-2,-1
316  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
317  cnt = cnt+1
318  tnuv(i,n0,r) = fct(cnt)
319  end do
320  end do
321  cnt = cnt+1
322  tnerr(r) = real(fct(cnt))
323  cnt = cnt+1
324  tnerr2(r) = real(fct(cnt))
325  end do
326  return
327  endif
328 
329 
330  if(rank.eq.rmax) then
331 
332  call calctnred(tn,tnuv,mominv,masses2,n,rank,id,tnerr,tnerr2)
333 
334  if (wrica) then
335  cnt = 0
336  do r=0,rank
337  do n0=r/2,0,-1
338  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
339 
340  cnt = cnt+1
341  fct(cnt) = tn(i,n0,r)
342 
343  end do
344  end do
345  do n0=r/2,n-2,-1
346  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
347 
348  cnt = cnt+1
349  fct(cnt) = tnuv(i,n0,r)
350 
351  end do
352  end do
353  cnt = cnt+1
354  fct(cnt) = tnerr(r)
355  cnt = cnt+1
356  fct(cnt) = tnerr2(r)
357  end do
358 
359  if(rank.ge.2*n-4) then
360  call writecache(fct,ncoefs(rank,n)+ncoefs(rank-2*n+4,n)+2*(rank+1),id,n,rank)
361  else
362  call writecache(fct,ncoefs(rank,n)+2*(rank+1),id,n,rank)
363  end if
364 
365  end if
366 
367  return
368 
369 
370  else
371  allocate(tnaux(binomtable(rank,n+rank-2),0:rank/2,0:rank))
372  allocate(tnuvaux(binomtable(rank,n+rank-2),0:rank/2,0:rank))
373  allocate(tnerraux(0:rank))
374  allocate(tnerr2aux(0:rank))
375 
376  call calctnred(tnaux,tnuvaux,mominv,masses2,n,rank,id,tnerraux,tnerr2aux)
377 
378  if (wrica) then
379  cnt = 0
380  deallocate(fct)
381  if(rank.ge.2*n-4) then
382  allocate(fct(ncoefs(rank,n)+ncoefs(rank-2*n+4,n)+2*(rank+1)))
383  else
384  allocate(fct(ncoefs(rank,n)+2*(rank+1)))
385  end if
386  do r=0,rank
387  do n0=r/2,0,-1
388  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
389 
390  cnt = cnt+1
391  fct(cnt) = tnaux(i,n0,r)
392 
393  end do
394  end do
395  do n0=r/2,n-2,-1
396  do i=1,binomtable(r-2*n0,n+r-2*n0-2)
397 
398  cnt = cnt+1
399  fct(cnt) = tnuvaux(i,n0,r)
400 
401  end do
402  end do
403  cnt = cnt+1
404  fct(cnt) = tnerraux(r)
405  cnt = cnt+1
406  fct(cnt) = tnerr2aux(r)
407  end do
408 
409  if(rank.ge.2*n-4) then
410  call writecache(fct,ncoefs(rank,n)+ncoefs(rank-2*n+4,n)+2*(rank+1),id,n,rank)
411  else
412  call writecache(fct,ncoefs(rank,n)+2*(rank+1),id,n,rank)
413  end if
414 
415  end if
416 
417  tn = tnaux(1:binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax)
418  tnuv = tnuvaux(1:binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax)
419  tnerr = tnerraux(0:rmax)
420  tnerr2 = tnerr2aux(0:rmax)
421  return
422 
423  end if
424  end if
425  end if
426 
427  call calctnred(tn,tnuv,mominv,masses2,n,rmax,id,tnerr,tnerr2)
428 
429 
430  end subroutine calctnint
431 
432 
433 
434 
435 
436  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
437  ! subroutine CalcTNred(TN,TNuv,MomInv,masses2,N,rmax,id)
438  !
439  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
440 
441  recursive subroutine calctnred(TN,TNuv,MomInv,masses2,N,rmax,id,TNerr,TNerr2)
442 
443  integer, intent(in) :: n,rmax,id
444  double complex, intent(in) :: mominv(binomtable(2,n)), masses2(0:n-1)
445  double precision, intent(out) :: tnerr(0:rmax),tnerr2(0:rmax)
446  double complex :: q10,q21,q32,q43,q54,q20,q31,q42,q53,q50,q30,q41,q52,q40,q51
447  double complex :: mm02,mm12,mm22,mm32,mm42,mm52
448  double complex :: mx(0:5,0:5), mxinv(0:5,0:5),mx0k(5,5),mx0kinv(5,5),ff(n-1)
449  double complex :: det,newdet,tnaux,mxinvs,mx0kinvs(5)
450  double complex :: zmx0kinv(5,5),z(5,5),zmx0kinvs(5)
451  double precision :: maxz,maxzmx0kinv(5),maxzmx0kinvs
452  double complex, intent(inout) :: tn(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax)
453  double complex, intent(inout) :: tnuv(binomtable(rmax,n+rmax-2),0:rmax/2,0:rmax)
454  double complex, allocatable :: tnm1_0(:,:,:), tnm1uv_0(:,:,:), tnm1_i(:,:,:,:)
455  double complex, allocatable :: tnm1_0aux(:,:,:), tnm1uv_0aux(:,:,:), tnm1uv_i(:,:,:,:)
456  double complex :: s(5), elimminf2_coli,chdet,gram(1:5,1:5),gramdet
457  double precision :: tnm1err(0:5,0:rmax),tnm1err2(0:5,0:rmax)
458  integer :: r,n0,n1,n2,n3,n4,n5,n6,k,i,j,m,nid(0:5),r0,bin,kbest,nm1,inds(n-1)
459  integer :: bino,bino_0,bino_i,bino_m1,bino_m2,cnt,iaux,rmax_m1,shift,cind,rbcd
460  integer :: mia1,mia2,mia3
461  logical :: errorwriteflag
462  character(len=*),parameter :: fmt10 = "(A18,'(',d25.18,' , ',d25.18,' )')"
463 
464 #ifdef TNtest
465  character(len=*),parameter :: fmt999 = "(' gm(',i1,',',i1,') = ',d23.16,' , ',d23.16)"
466 #endif
467 
468 ! write(*,*) 'CalcTNred in N',N,rmax
469 ! write(*,*) 'CalcTNred in s',MomInv
470 ! write(*,*) 'CalcTNred in m',masses2
471 
472  ! allocation of TN functions
473  nm1 = n-1
474  rmax_m1 = max(rmax-1,0)
475  bino_0 = binomtable(rmax_m1,n+rmax_m1-2)
476  bino_i = binomtable(rmax_m1,nm1+rmax_m1-2)
477  allocate(tnm1_0(bino_0,0:rmax_m1/2,0:rmax_m1))
478  allocate(tnm1_0aux(bino_i,0:rmax_m1/2,0:rmax_m1))
479  allocate(tnm1_i(bino_i,0:rmax_m1/2,0:rmax_m1,5))
480  allocate(tnm1uv_0(bino_0,0:rmax_m1/2,0:rmax_m1))
481  allocate(tnm1uv_0aux(bino_i,0:rmax_m1/2,0:rmax_m1))
482  allocate(tnm1uv_i(bino_i,0:rmax_m1/2,0:rmax_m1,5))
483 
484  ! determine binaries for TN-coefficients
485  k=0
486  bin = 1
487  do while (k.le.5)
488  if (mod(id/bin,2).eq.0) then
489  nid(k) = id+bin
490  k = k+1
491  end if
492  bin = 2*bin
493  end do
494 
495  call calctnint(tnm1_0aux,tnm1uv_0aux,submominv(n,0,mominv),submasses(n,0,masses2), &
496  nm1,rmax_m1,nid(0),tnerr=tnm1err(0,0:rmax_m1), &
497  tnerr2=tnm1err2(0,0:rmax_m1))
498  do k=1,5
499  call calctnint(tnm1_i(:,:,:,k),tnm1uv_i(:,:,:,k),submominv(n,k,mominv), &
500  submasses(n,k,masses2),nm1,rmax_m1,nid(k), &
501  tnerr=tnm1err(k,0:rmax_m1),tnerr2=tnm1err2(k,0:rmax_m1))
502  end do
503 
504  ! shift of integration momentum in TN\{0}
505  do n0=0,rmax_m1/2
506  tnm1_0(1,n0,2*n0) = tnm1_0aux(1,n0,2*n0)
507  do r=1,rmax_m1-2*n0
508  mia1 = binomtable(r-1,n+r-3)+1
509  mia2 = binomtable(r,n+r-3)
510  mia3 = binomtable(r,n+r-2)
511  tnm1_0(mia1:mia3,n0,r+2*n0) = tnm1_0aux(1:mia2,n0,r+2*n0)
512  end do
513  end do
514  do r=1,rmax_m1
515  do n0=0,(r-1)/2
516  do i=binomtable(r-2*n0-1,n+r-2*n0-3),1,-1
517  tnm1_0(i,n0,r) = -tnm1_0(i,n0,r-1)
518  do j=2,n-1
519  tnm1_0(i,n0,r) = tnm1_0(i,n0,r) &
520  - tnm1_0(addtocind(j,i,r-2*n0-1,n-1),n0,r)
521  end do
522  end do
523  end do
524  end do
525  do n0=n-3,rmax_m1/2
526  tnm1uv_0(1,n0,2*n0) = tnm1uv_0aux(1,n0,2*n0)
527  do r=1,rmax_m1-2*n0
528  mia1 = binomtable(r-1,n+r-3)+1
529  mia2 = binomtable(r,n+r-3)
530  mia3 = binomtable(r,n+r-2)
531  tnm1uv_0(mia1:mia3,n0,r+2*n0) = tnm1uv_0aux(1:mia2,n0,r+2*n0)
532  end do
533  end do
534  do r=2*n-6,rmax_m1
535  do n0=n-3,(r-1)/2
536  do i=binomtable(r-2*n0-1,n+r-2*n0-3),1,-1
537  tnm1uv_0(i,n0,r) = -tnm1uv_0(i,n0,r-1)
538  do j=2,n-1
539  tnm1uv_0(i,n0,r) = tnm1uv_0(i,n0,r) &
540  - tnm1uv_0(addtocind(j,i,r-2*n0-1,n-1),n0,r)
541  end do
542  end do
543  end do
544  end do
545 
546 
547  ! determine inverse modified Caley matrix
548  mm02 = elimminf2_coli(masses2(0))
549  mm12 = elimminf2_coli(masses2(1))
550  mm22 = elimminf2_coli(masses2(2))
551  mm32 = elimminf2_coli(masses2(3))
552  mm42 = elimminf2_coli(masses2(4))
553  mm52 = elimminf2_coli(masses2(5))
554  q10 = elimminf2_coli(mominv(1))
555  q21 = elimminf2_coli(mominv(2))
556  q32 = elimminf2_coli(mominv(3))
557  q43 = elimminf2_coli(mominv(4))
558  q54 = elimminf2_coli(mominv(5))
559  if (n.le.9) then
560  q50 = elimminf2_coli(mominv((n-6)*n+6))
561  else
562  q50 = elimminf2_coli(mominv(4*n+1))
563  end if
564  q20 = elimminf2_coli(mominv(n+1))
565  q31 = elimminf2_coli(mominv(n+2))
566  q42 = elimminf2_coli(mominv(n+3))
567  q53 = elimminf2_coli(mominv(n+4))
568  if (n.le.7) then
569  q40 = elimminf2_coli(mominv((n-5)*n+5))
570  q51 = elimminf2_coli(mominv((n-5)*n+6))
571  else
572  q40 = elimminf2_coli(mominv(3*n+1))
573  q51 = elimminf2_coli(mominv(3*n+2))
574  end if
575  q30 = elimminf2_coli(mominv(2*n+1))
576  q41 = elimminf2_coli(mominv(2*n+2))
577  q52 = elimminf2_coli(mominv(2*n+3))
578  ! write(*,*) 'q50', q50
579  ! write(*,*) 'q20', q20
580  ! write(*,*) 'q31', q31
581  ! write(*,*) 'q42', q42
582  ! write(*,*) 'q53', q53
583  ! write(*,*) 'q40', q40
584  ! write(*,*) 'q51', q51
585  ! write(*,*) 'q30', q30
586  ! write(*,*) 'q41', q41
587  ! write(*,*) 'q52', q52
588 
589  cnt = 1
590  do i=1,n-1
591  ff(i) = elimminf2_coli(mominv(cnt)) + mm02 &
592  - elimminf2_coli(masses2(i))
593  cnt = cnt+n-i
594  end do
595 
596  mx(0,0) = 2d0*mm02
597  mx(1,0) = q10 - mm12 + mm02
598  mx(2,0) = q20 - mm22 + mm02
599  mx(3,0) = q30 - mm32 + mm02
600  mx(4,0) = q40 - mm42 + mm02
601  mx(5,0) = q50 - mm52 + mm02
602  mx(0,1) = mx(1,0)
603  mx(1,1) = 2d0*q10
604  mx(2,1) = q10+q20-q21
605  mx(3,1) = q10+q30-q31
606  mx(4,1) = q10+q40-q41
607  mx(5,1) = q10+q50-q51
608  mx(0,2) = mx(2,0)
609  mx(1,2) = mx(2,1)
610  mx(2,2) = 2d0*q20
611  mx(3,2) = q20+q30-q32
612  mx(4,2) = q20+q40-q42
613  mx(5,2) = q20+q50-q52
614  mx(0,3) = mx(3,0)
615  mx(1,3) = mx(3,1)
616  mx(2,3) = mx(3,2)
617  mx(3,3) = 2d0*q30
618  mx(4,3) = q30+q40-q43
619  mx(5,3) = q30+q50-q53
620  mx(0,4) = mx(4,0)
621  mx(1,4) = mx(4,1)
622  mx(2,4) = mx(4,2)
623  mx(3,4) = mx(4,3)
624  mx(4,4) = 2d0*q40
625  mx(5,4) = q40+q50-q54
626  mx(0,5) = mx(5,0)
627  mx(1,5) = mx(5,1)
628  mx(2,5) = mx(5,2)
629  mx(3,5) = mx(5,3)
630  mx(4,5) = mx(5,4)
631  mx(5,5) = 2d0*q50
632 
633  call chinv(6,mx,mxinv)
634 
635  ! determine X_(0,5)
636  do j=1,5
637  do i=1,5
638  mx0k(i,j) = mx(i,j-1)
639  end do
640  end do
641 
642  det = chdet(5,mx0k)
643  kbest = 5
644 
645 #ifdef TNtest
646  write(*,*) 'det',5,det
647 #endif
648 
649  do j=5,2,-1
650  do i=1,5
651  mx0k(i,j) = mx(i,j)
652  end do
653 
654  newdet = chdet(5,mx0k)
655  if (abs(newdet).gt.abs(det)) then
656  kbest = j-1
657 #ifdef TNtest
658  write(*,*) 'det',j-1,newdet
659 #endif
660  det = newdet
661  end if
662 
663  end do
664 
665 #ifdef TNtest
666 ! kbest=5
667  write(*,*) 'kbest',kbest
668 #endif
669 
670 
671 
672 
673  do i=1,5
674  mx0k(i,1) = mx(i,1)
675  mx0k(i,kbest) = mx(i,0)
676  end do
677 
678  call chinv(5,mx0k,mx0kinv)
679  do i=1,5
680  mx0kinv(kbest,i) = 0d0
681  end do
682 
683  mxinvs = sum(mxinv(0:5,0))
684  do i=1,5
685  mx0kinvs(i) = sum(mx0kinv(i,1:5))
686  end do
687 
688  ! for alternative error estimate
689  z(1:5,1:5) = mx(1:5,1:5)
690 
691  zmx0kinv = matmul(z,mx0kinv)
692 
693  do i=1,5
694  maxzmx0kinv(i) = maxval(abs(zmx0kinv(1:5,i)))
695  zmx0kinvs(i) = sum(zmx0kinv(i,1:5))
696  end do
697 
698  maxzmx0kinvs = maxval(abs(zmx0kinvs(1:5)))
699 
700  maxz = maxval(abs(z))
701 
702  ! calculation of UV-divergent parts
703  ! TN UV-finite for n0<N-2
704  ! TNuv(:,0:min(rmax/2,N-3),:) = 0d0
705 
706  ! PV reduction (5.10)
707  do r=2*n-4,rmax
708  do n0=n-2,r/2
709  do cind=1,binomtable(r-2*n0,n+r-2*n0-2)
710 
711  tnuv(cind,n0,r) = tnm1uv_0(cind,n0-1,r-2) + 2*mm02*tnuv(cind,n0-1,r-2)
712  do i=1,nm1
713  tnuv(cind,n0,r) = tnuv(cind,n0,r) &
714  + ff(i)*tnuv(addtocind(i,cind,r-2*n0,n-1),n0-1,r-1)
715  end do
716  tnuv(cind,n0,r) = tnuv(cind,n0,r)/(2*(r+3-n))
717 
718  end do
719  end do
720  end do
721 
722  ! scalar coefficient
723  tn = 0d0
724  tn(1,0,0) = -mxinv(0,0)*tnm1_0(1,0,0)
725  do k=1,5
726  tn(1,0,0) = tn(1,0,0) &
727  + mxinv(k,0)*(tnm1_i(1,0,0,k)-tnm1_0(1,0,0))
728  end do
729 
730 ! TNerr(0) = max( maxval(abs(mxinv(0:4,0)))*TNm1err(0,0), &
731  tnerr(0) = max(abs(mxinvs)*tnm1err(0,0), &
732  abs(mxinv(1,0))*tnm1err(1,0) , &
733  abs(mxinv(2,0))*tnm1err(2,0) , &
734  abs(mxinv(3,0))*tnm1err(3,0) , &
735  abs(mxinv(4,0))*tnm1err(4,0) , &
736  abs(mxinv(5,0))*tnm1err(5,0) )
737 
738  tnerr2(0) = max(abs(mxinvs)*tnm1err2(0,0), &
739  abs(mxinv(1,0))*tnm1err2(1,0) , &
740  abs(mxinv(2,0))*tnm1err2(2,0) , &
741  abs(mxinv(3,0))*tnm1err2(3,0) , &
742  abs(mxinv(4,0))*tnm1err2(4,0) , &
743  abs(mxinv(5,0))*tnm1err2(5,0) )
744 
745  if (rmax.eq.0) return
746 
747 
748  ! formula (7.13) extended to arbitrary N
749  do r=0,rmax_m1
750  do n0=0,max((r-n+6)/2,0)
751  bino_m1 = binomtable(r-2*n0,n+r-2*n0-2)
752  do i=1,bino_m1
753  inds = calccindarr(nm1,r-2*n0,i)
754 
755  do m=1,5
756  s(m) = -tnm1_0(i,n0,r)
757  end do
758 
759 #ifdef TNtest
760  write(*,*) 'Tnred s(m)',s
761 #endif
762 
763  if (inds(1).eq.0) then
764  s(1) = s(1) + tnm1_i(dropcind(1,i,r-2*n0,nm1),n0,r,1)
765  end if
766  if (inds(2).eq.0) then
767  s(2) = s(2) + tnm1_i(dropcind(2,i,r-2*n0,nm1),n0,r,2)
768  end if
769  if (inds(3).eq.0) then
770  s(3) = s(3) + tnm1_i(dropcind(3,i,r-2*n0,nm1),n0,r,3)
771  end if
772  if (inds(4).eq.0) then
773  s(4) = s(4) + tnm1_i(dropcind(4,i,r-2*n0,nm1),n0,r,4)
774  end if
775  if (inds(5).eq.0) then
776  s(5) = s(5) + tnm1_i(dropcind(5,i,r-2*n0,nm1),n0,r,5)
777  end if
778 
779 #ifdef TNtest
780  write(*,*) 'Tnred s(m)',s
781 #endif
782 
783  do k=1,5
784  tnaux = mx0kinv(k,1)*s(1)+mx0kinv(k,2)*s(2) &
785  + mx0kinv(k,3)*s(3)+mx0kinv(k,4)*s(4)+mx0kinv(k,5)*s(5)
786  iaux = addtocind(k,i,r,n-1)
787  tn(iaux,n0,r+1) = tn(iaux,n0,r+1) + (inds(k)+1)*tnaux/(r+1)
788 #ifdef TNtest
789  write(*,*) 'Tnred comb',k,mx0kinv(k,1)*s(1),mx0kinv(k,2)*s(2) &
790  , mx0kinv(k,3)*s(3),mx0kinv(k,4)*s(4),mx0kinv(k,5)*s(5), tnaux
791 
792  if(abs(tnaux).ne.0d0) then
793  write(*,*) 'Tnred cancel', k, max(abs(mx0kinv(k,1)*s(1)),abs(mx0kinv(k,2)*s(2)) &
794  , abs(mx0kinv(k,3)*s(3)),abs(mx0kinv(k,4)*s(4)),abs(mx0kinv(k,5)*s(5))) &
795  /abs(tnaux)
796  endif
797 #endif
798 
799  end do
800 
801 
802 
803  end do
804  end do
805 
806  if (r.le.rmax-1) then
807 
808  tnerr(r+1) = max( maxval(abs(mx0kinvs(1:5)))*tnm1err(0,r), &
809  maxval(abs(mx0kinv(1:5,1)))*tnm1err(1,r) , &
810  maxval(abs(mx0kinv(1:5,2)))*tnm1err(2,r) , &
811  maxval(abs(mx0kinv(1:5,3)))*tnm1err(3,r) , &
812  maxval(abs(mx0kinv(1:5,4)))*tnm1err(4,r) , &
813  maxval(abs(mx0kinv(1:5,5)))*tnm1err(5,r) )
814 
815  tnerr2(r+1) = max( abs(maxzmx0kinvs)*tnm1err2(0,r), &
816  abs(maxzmx0kinv(1))*tnm1err2(1,r) , &
817  abs(maxzmx0kinv(2))*tnm1err2(2,r) , &
818  abs(maxzmx0kinv(3))*tnm1err2(3,r) , &
819  abs(maxzmx0kinv(4))*tnm1err2(4,r) , &
820  abs(maxzmx0kinv(5))*tnm1err2(5,r) )/maxz
821 
822 
823  if (mode_coli.lt.1) then
824  gram= mx(1:5,1:5)
825 
826  gramdet = chdet(5,gram)
827 
828 #ifdef TNtest
829  write(*,*)
830  do i = 1,5
831  do j = 1,5
832  write(*,fmt999) i,j,gram(i,j)
833  enddo
834  enddo
835  write(*,*) 'TNred relGramdet=',gramdet/det
836 #endif
837 
838  if (max(abs(tn(1,0,0)),maxval(abs(tn(1:6,0,1))))*abs(gramdet/det).gt. &
839  tnerr(r+1)) then
840 
841  tnerr(r+1)= max(tnerr(r+1), &
842  max(abs(tn(1,0,0)),maxval(abs(tn(1:6,0,1))))*abs(gramdet/det) )
843  tnerr2(r+1)= max(tnerr2(r+1), &
844  max(abs(tn(1,0,0)),maxval(abs(tn(1:6,0,1))))*abs(gramdet/det) )
845 
846  if (abs(gramdet/det).gt.reqacc_coli) then
847  call seterrflag_coli(-6)
848  call errout_coli('CalcTNred', &
849  'input momenta inconsistent! (not 4-dimensional)', &
850  errorwriteflag)
851  if (errorwriteflag) then
852  write(nerrout_coli,fmt10) ' CalcTNred: q10 = ',q10
853  write(nerrout_coli,fmt10) ' CalcTNred: q21 = ',q21
854  write(nerrout_coli,fmt10) ' CalcTNred: q32 = ',q32
855  write(nerrout_coli,fmt10) ' CalcTNred: q43 = ',q43
856  write(nerrout_coli,fmt10) ' CalcTNred: q54 = ',q54
857  write(nerrout_coli,fmt10) ' CalcTNred: q50 = ',q50
858  write(nerrout_coli,fmt10) ' CalcTNred: q20 = ',q10
859  write(nerrout_coli,fmt10) ' CalcTNred: q31 = ',q31
860  write(nerrout_coli,fmt10) ' CalcTNred: q42 = ',q42
861  write(nerrout_coli,fmt10) ' CalcTNred: q53 = ',q53
862  write(nerrout_coli,fmt10) ' CalcTNred: q40 = ',q40
863  write(nerrout_coli,fmt10) ' CalcTNred: q51 = ',q51
864  write(nerrout_coli,fmt10) ' CalcTNred: q30 = ',q30
865  write(nerrout_coli,fmt10) ' CalcTNred: q41 = ',q41
866  write(nerrout_coli,fmt10) ' CalcTNred: q52 = ',q52
867  write(nerrout_coli,fmt10) ' CalcTNred: mm02 = ',mm02
868  write(nerrout_coli,fmt10) ' CalcTNred: mm12 = ',mm12
869  write(nerrout_coli,fmt10) ' CalcTNred: mm22 = ',mm22
870  write(nerrout_coli,fmt10) ' CalcTNred: mm32 = ',mm32
871  write(nerrout_coli,fmt10) ' CalcTNred: mm42 = ',mm42
872  write(nerrout_coli,fmt10) ' CalcTNred: mm52 = ',mm52
873  write(nerrout_coli,fmt10) ' CalcTNred: gram = ',gramdet/det
874  end if
875  end if
876 
877  end if
878 
879  end if
880  end if
881 
882 
883 #ifdef TNtest
884  write(*,*) 'Tnred mn err',tnm1err(1:5,r)
885  write(*,*) 'Tnred coef err', (maxval(abs(mx0kinv(1:5,k))),k=1,5)
886  write(*,*) 'Tnred err',tnerr(r+1)
887 #endif
888 
889  end do
890 
891 
892  end subroutine calctnred
893 
894 
895 
896 
897 
898  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
899  ! function SubMasses
900  !
901  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
902 
903  function submasses(N,k,masses)
904 
905  integer :: n,k,i
906  double complex :: masses(0:n-1), submasses(0:n-2)
907 
908  if ((k.gt.n-1).or.(k.lt.0)) then
909  write(nerrout_coli,*) 'SubMasses:'
910  write(nerrout_coli,*) 'inkonsistent argument k', k
911  write(nerrout_coli,*) 0, '<= k <=', n-1, 'required!'
912  stop
913  end if
914 
915  do i=0,n-2
916  if (i.lt.k) then
917  submasses(i) = masses(i)
918  else
919  submasses(i) = masses(i+1)
920  end if
921  end do
922 
923  end function
924 
925 
926 
927 
928 
929  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
930  ! function SubMomVec
931  !
932  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
933 
934  function submomvec(N,k,MomVec)
935 
936  integer :: n,k,i
937  double complex :: momvec(0:3,n-1), submomvec(0:3,n-2)
938 
939  if ((k.gt.n-1).or.(k.lt.0)) then
940  write(nerrout_coli,*) 'SubMomVec:'
941  write(nerrout_coli,*) 'inkonsistent argument k', k
942  write(nerrout_coli,*) 0, '<= k <=', n-1, 'required!'
943  stop
944  end if
945 
946  if (k.eq.0) then
947  do i=1,n-2
948  submomvec(0:3,i) = momvec(0:3,i+1)-momvec(0:3,1)
949  end do
950 
951  else
952  do i=1,n-2
953  if (i.lt.k) then
954  submomvec(0:3,i) = momvec(0:3,i)
955  else
956  submomvec(0:3,i) = momvec(0:3,i+1)
957  end if
958  end do
959 
960  end if
961 
962  end function
963 
964 
965 
966 
967 
968  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
969  ! function SubMomInv
970  ! ******************************************
971  ! * rewritten by Ansgar Denner 25.03.15 *
972  ! ******************************************
973  !
974  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
975 
976  function submominv(N,k,MomInv)
977 
978  integer, intent(in) :: n,k
979  integer :: i,cnt,moms,limit
980  double complex :: mominv(binomtable(2,n)), submominv(binomtable(2,n-1))
981 
982  ! write(*,*) 'SubMomIv', N,k
983  if ((k.gt.n-1).or.(k.lt.0)) then
984  write(nerrout_coli,*) 'SubMomInv:'
985  write(nerrout_coli,*) 'inkonsistent argument k', k
986  write(nerrout_coli,*) 0, '<= k <=', n-1, 'required!'
987  stop
988  end if
989 
990  if (n.eq.2) return
991 
992  cnt = 1
993  do moms=1,(n-3)/2
994  do i=1,n
995  if (i.ne.k+1) then
996  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
997  submominv(cnt) = mominv(n*moms+i)
998  else
999  submominv(cnt) = mominv(n*(moms-1)+i)
1000  end if
1001 
1002  cnt = cnt+1
1003  end if
1004  end do
1005  end do
1006 
1007  if (mod(n,2).eq.1) then
1008  moms = (n-1)/2
1009  if(k.le.(n-1)/2) then
1010  limit = (n+1)/2
1011  else
1012  limit = (n-1)/2
1013  end if
1014  do i=1,limit
1015  if (i.ne.k+1) then
1016  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
1017  submominv(cnt) = mominv(n*moms+i-(n-1)/2)
1018  else
1019  submominv(cnt) = mominv(n*(moms-1)+i)
1020  end if
1021 
1022  cnt = cnt+1
1023  end if
1024  end do
1025 
1026  else
1027  do i=1,n
1028  if (i.ne.k+1) then
1029  moms = (n/2-1)
1030  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
1031  if (i.gt.n/2) then
1032  submominv(cnt) = mominv(n*moms-n/2+i)
1033  else
1034  submominv(cnt) = mominv(n*moms+i)
1035  endif
1036  else
1037  submominv(cnt) = mominv(n*(moms-1)+i)
1038  end if
1039 
1040  cnt = cnt+1
1041  end if
1042  end do
1043 
1044  end if
1045 
1046 
1047  end function
1048 
1049 
1050 
1051 
1052 
1053 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1054 ! ! function CalcCIndArr
1055 ! !
1056 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1057 !
1058 ! function CalcCIndArr(Nm1,r,cind) result(arr)
1059 !
1060 ! integer :: Nm1, r, cind, arr(Nm1), i
1061 ! integer :: combis(r)
1062 !
1063 ! if (r.ge.1) then
1064 ! combis = IndCombisEq(1:r,cind,r,Nm1)
1065 ! end if
1066 ! arr = 0
1067 ! do i=1,r
1068 ! arr(combis(i)) = arr(combis(i))+1
1069 ! end do
1070 !
1071 ! end function CalcCIndArr
1072 
1073 
1074 
1075 
1076 
1077 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1078 ! ! subroutine SetDropCInd(Nm1,r)
1079 ! !
1080 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1081 !
1082 ! subroutine SetDropCInd(Nm1,r)
1083 !
1084 ! integer :: Nm1,r,i,j,k
1085 !
1086 ! if (allocated(DropCInd)) then
1087 ! deallocate(DropCInd)
1088 ! end if
1089 ! allocate(DropCInd(Nm1,BinomTable(r,r+Nm1-1),0:r,Nm1))
1090 !
1091 ! do i=1,Nm1
1092 ! DropCInd(1:i,1,0,i) = 1
1093 ! do j=1,r
1094 ! do k=1,i
1095 ! DropCind(k,1:BinomTable(j,j+i-1),j,i) = CalcDropCInd(i,j,k)
1096 ! end do
1097 ! end do
1098 ! end do
1099 !
1100 ! end subroutine SetDropCInd
1101 
1102 
1103 
1104 
1105 
1106 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1107 ! ! function CalcDropCInd(Nm1,r,k)
1108 ! !
1109 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1110 !
1111 ! function CalcDropCInd(Nm1,r,k) result(dcinds)
1112 !
1113 ! integer :: Nm1,r,k,dcinds(BinomTable(r,r+Nm1-1))
1114 ! integer :: cnt,arr(Nm1),i
1115 !
1116 ! cnt = 1
1117 ! do i=1,BinomTable(r,r+Nm1-1)
1118 ! arr = CalcCIndArr(Nm1,r,i)
1119 ! if (arr(k).eq.0) then
1120 ! dcinds(i) = cnt
1121 ! cnt = cnt+1
1122 ! else
1123 ! dcinds(i) = 0
1124 ! end if
1125 ! end do
1126 !
1127 ! end function CalcDropCInd
1128 
1129 
1130 
1131 
1132 
1133 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1134 ! ! subroutine SetAddToCInd(Nm1,r)
1135 ! !
1136 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1137 !
1138 ! subroutine SetAddToCInd(Nm1,r)
1139 !
1140 ! integer, intent(in) :: Nm1,r
1141 ! integer :: i,j,k,m
1142 !
1143 ! if (allocated(AddToCInd)) then
1144 ! deallocate(AddToCInd)
1145 ! end if
1146 ! allocate(AddToCInd(Nm1,BinomTable(r,Nm1+r-1),0:r,Nm1))
1147 ! AddToCInd = 0
1148 !
1149 ! do i=1,Nm1
1150 ! do j=0,r
1151 ! do k=1,BinomTable(j,i+j-1)
1152 ! do m=1,i
1153 ! AddToCInd(m,k,j,i) = CalcAddToCInd(i,j,k,m)
1154 ! end do
1155 ! end do
1156 ! end do
1157 ! end do
1158 !
1159 ! end subroutine SetAddToCInd
1160 
1161 
1162 
1163 
1164 
1165 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1166 ! ! function CalcAddToCInd
1167 ! !
1168 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1169 !
1170 ! function CalcAddToCInd(Nm1,r,cind,k) result(newcind)
1171 !
1172 ! integer :: Nm1,r,k,cind,newcind,i
1173 ! integer :: indArr(r),newindArr(r+1)
1174 ! logical :: kflag
1175 !
1176 ! if ((r.eq.0).and.(cind.eq.1)) then
1177 ! newcind = k
1178 ! return
1179 ! end if
1180 !
1181 ! indArr = IndCombisEq(1:r,cind,r,Nm1)
1182 !
1183 ! kflag = .true.
1184 ! do i=1,r
1185 ! if (indArr(i).ge.k) then
1186 ! newindArr(i+1) = indArr(i)
1187 ! if (kflag) then
1188 ! newindArr(i) = k
1189 ! kflag = .false.
1190 ! end if
1191 ! else
1192 ! newindArr(i) = indArr(i)
1193 ! end if
1194 ! end do
1195 ! if (kflag) then
1196 ! newindArr(r+1) = k
1197 ! end if
1198 !
1199 ! newcind = CalcPosIndCombisEq(Nm1,r+1,newindArr)
1200 !
1201 ! end function CalcAddToCInd
1202 
1203 
1204 end module reductiontn
submasses
double complex function, dimension(0:n-2) submasses(N, k, masses)
Definition: democache.f90:468
reductiontn::calctn
subroutine calctn(TN, TNuv, MomInv, masses2, N, rmax, id, TNerr, TNerr2)
Definition: reductionTN.F90:45
reductiontn::calctnred
recursive subroutine calctnred(TN, TNuv, MomInv, masses2, N, rmax, id, TNerr, TNerr2)
Definition: reductionTN.F90:442
reductiontn
Definition: reductionTN.F90:28
submominv
double complex function, dimension((n-1) *(n-2)/2) submominv(N, k, MomInv)
Definition: democache.f90:495
reductionefg::calce
subroutine calce(E, Euv, p10, p21, p32, p43, p40, p20, p31, p42, p30, p41, m02, m12, m22, m32, m42, rmax, id, Eerr, Eerr2)
Definition: reductionEFG.F90:43
reductiontn::calctnint
recursive subroutine calctnint(TN, TNuv, MomInv, masses2, N, rmax, id, TNerr, TNerr2)
Definition: reductionTN.F90:238
reductionefg
Definition: reductionEFG.F90:25
reductiontn::submomvec
double complex function, dimension(0:3, n-2) submomvec(N, k, MomVec)
Definition: reductionTN.F90:935