test.f90 28.7 KB
Newer Older
Luca's avatar
Luca committed
1

Luca's avatar
Luca committed
2
! These tests assume that GPLInfinity = 30
Luca's avatar
Luca committed
3

Luca's avatar
Luca committed
4
PROGRAM TEST
5
  use globals
6
  use utils
7
  use shuffle
8
  use maths_functions
Luca's avatar
Luca committed
9 10
  use mpl_module
  use gpl_module
ulrich_y's avatar
ulrich_y committed
11
  use chenreftest, only: do_chen_test
Luca's avatar
Luca committed
12
  implicit none
Luca's avatar
Luca committed
13
  
14
  real(kind=prec) :: tol = 8.0e-7
Luca's avatar
Luca committed
15
  logical :: tests_successful = .true. 
16

ulrich_y's avatar
ulrich_y committed
17
#ifdef DEBUG
Luca's avatar
tests  
Luca committed
18
  call parse_cmd_args()
ulrich_y's avatar
ulrich_y committed
19
#endif
ulrich_y's avatar
ulrich_y committed
20

ulrich_y's avatar
ulrich_y committed
21
  tol = zero * 1.e5_prec
Luca's avatar
tests  
Luca committed
22 23
  call do_MPL_tests() 
  call do_GPL_tests()
24
  ! call do_shuffle_tests() ! put this somewhere else
ulrich_y's avatar
ulrich_y committed
25
  tol = 8.0e-7
26
  tests_successful = tests_successful .and. do_chen_test()
Luca's avatar
tests  
Luca committed
27

ulrich_y's avatar
ulrich_y committed
28
#ifdef HAVE_GINAC
ulrich_y's avatar
ulrich_y committed
29 30
#ifdef HAVE_MM
  call do_ginac_tests
ulrich_y's avatar
ulrich_y committed
31
#ifndef NOSPEED
ulrich_y's avatar
ulrich_y committed
32
  call do_timing_tests(5)
ulrich_y's avatar
ulrich_y committed
33
  call do_high_weight_tests
ulrich_y's avatar
ulrich_y committed
34
#endif
ulrich_y's avatar
ulrich_y committed
35 36 37
#ifdef LONG_TEST
  call do_long_test
#endif
ulrich_y's avatar
ulrich_y committed
38
#endif
ulrich_y's avatar
ulrich_y committed
39
#endif
Luca's avatar
tests  
Luca committed
40 41 42 43 44 45 46

  if(tests_successful) then
    print*, 'All tests passed. '
  else 
    print*, 'Some tests failed. '
    stop 1
  end if
ulrich_y's avatar
ulrich_y committed
47
  stop
48

49
CONTAINS
50
   
Luca Naterop's avatar
Luca Naterop committed
51 52 53 54
  subroutine check(res, ref)
    complex(kind=prec) :: res, ref
    real(kind=prec) :: delta

Luca Naterop's avatar
Luca Naterop committed
55
    delta = abs(res-ref)
Luca Naterop's avatar
Luca Naterop committed
56
    if(delta < tol) then
Luca Naterop's avatar
Luca Naterop committed
57
      print*, '       ',' passed with delta = ', delta
Luca Naterop's avatar
Luca Naterop committed
58
    else 
59
      print*, '  ',' FAILED with delta = ', delta
60
      tests_successful = .false.
Luca Naterop's avatar
Luca Naterop committed
61 62
    end if
  end subroutine check
63 64 65

  subroutine test_one_MPL(m,x,ref, test_id)
    integer :: m(:)
66
    complex(kind=prec) :: x(:), ref, res
67
    character(len=*) :: test_id
Luca's avatar
Luca committed
68
    
69 70 71 72 73
    print*, '  ', 'testing MPL ', test_id, ' ...'
    res = MPL(m,x)
    call check(res,ref)
  end subroutine test_one_MPL

Luca's avatar
Luca committed
74
  subroutine do_MPL_tests()
Luca's avatar
Luca committed
75
    complex(kind=prec) :: ref
Luca Naterop's avatar
Luca Naterop committed
76
    print*, 'doing MPL tests...'
Luca's avatar
Luca committed
77
    
ulrich_y's avatar
ulrich_y committed
78 79 80 81 82
    ref = (0.0226974036143118403367048874190211851348_prec,0.)
    call test_one_MPL((/ 1,1 /),cmplx([ 119._prec/377._prec,  35._prec / 102._prec ],kind=prec),ref, '1.1')

    ref = (2.31346156303083313941538625891818239608127E-4_prec,0.)
    call test_one_MPL((/ 1,1 /),cmplx([0.03_prec, 100._prec-30._prec*sqrt(11._prec) ], kind=prec),ref, '1.2')
Luca's avatar
Luca committed
83
    
ulrich_y's avatar
ulrich_y committed
84 85
    ref = (2.34617486908178608355021488772249510275409004880846912378356e-5_prec,0.)
    call test_one_MPL((/ 2,1,2 /),cmplx([0.03_prec, 100._prec-30._prec*sqrt(11._prec), 250._prec/9._prec + sqrt(6875._prec)/3._prec ], kind=prec),ref, '1.3')
Luca's avatar
Luca committed
86
    
87 88
    ref = (-0.06565799418838388900030511939327151809905663800733251335678_prec,0.)
    call test_one_MPL((/1, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.4')
ulrich_y's avatar
ulrich_y committed
89

90
    ref = (-0.0319989639656484125354488353677231989818697987217232462808_prec,0.)
ulrich_y's avatar
ulrich_y committed
91
    call test_one_MPL((/2, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.5')
Luca's avatar
Luca committed
92
  end subroutine do_MPL_tests
Luca's avatar
Luca committed
93

Luca's avatar
Luca committed
94
  subroutine test_one_condensed(m,z,y,k,ref,test_id)
95 96 97
    integer :: m(:), k
    complex(kind=prec) :: z(:), y, res, ref
    character(len=*) :: test_id
Luca Naterop's avatar
Luca Naterop committed
98

99
    print*, '  ', 'testing GPL ', test_id, ' ...'
ulrich_y's avatar
ulrich_y committed
100
    res = G_condensed(m,toinum(z),inum(y,di0),k)
Luca Naterop's avatar
Luca Naterop committed
101
    call check(res,ref)
Luca's avatar
Luca committed
102 103
  end subroutine test_one_condensed

104 105
  subroutine test_one_flat(z,ref,test_id)
    complex(kind=prec) :: z(:), res, ref
Luca's avatar
Luca committed
106 107 108
    character(len=*) :: test_id

    print*, '  ', 'testing GPL ', test_id, ' ...'
ulrich_y's avatar
ulrich_y committed
109
    res = G(z)
Luca's avatar
Luca committed
110 111
    call check(res,ref)
  end subroutine test_one_flat
Luca Naterop's avatar
Luca Naterop committed
112 113

  subroutine do_GPL_tests()
ulrich_y's avatar
ulrich_y committed
114
    complex(kind=prec) :: ref, res
Luca Naterop's avatar
Luca Naterop committed
115
    real(kind=prec) :: z, xchen
Luca Naterop's avatar
Luca Naterop committed
116
    print*, 'doing GPL tests...'
Luca's avatar
Luca committed
117
    
118 119
    ref = (0.0819393734128676670196050250748145117414_prec,0.)
    call test_one_condensed((/ 1,1 /),(/ (1.3_prec,0.), (1.1_prec,0.) /),(0.4_prec,0.),2,ref,'2.1')
Luca's avatar
Luca committed
120
    
121 122 123 124 125
    ref = (0.0159279595253714815409530431586662402173_prec,0.)
    call test_one_condensed((/ 3,2 /),(/ (1.3_prec,0.), (1.1_prec,0.) /),(0.4_prec,0.),2,ref,'2.2')

    ref = (0.0020332632172573960663410568998748341439_prec,0.)
    call test_one_condensed((/ 4 /),(/ (0._prec,0.) /),(1.6_prec,0.),1,ref,'2.3')
Luca's avatar
Luca committed
126

Luca Naterop's avatar
Luca Naterop committed
127
    ! requires making convergent
128 129
    ref = (0.0959304167763934268889451723647597269243_prec,-0.882935179519785044294098901083860819793_prec)
    call test_one_flat(cmplx([0,1,3,2],kind=prec),ref,'2.4')
Luca Naterop's avatar
Luca Naterop committed
130

131 132
    ref = (0.00994794778992862285848386211876302459484_prec,0.0)
    call test_one_flat([ (0._prec,0.), (0._prec,0.), cmplx(-10._prec/3._prec,kind=prec), cmplx(-10._prec/3._prec,kind=prec), (1._prec,0.) ],ref,'2.5')
Luca Naterop's avatar
Luca Naterop committed
133 134

    ! requires hoelder convolution
135 136
    ref = (-0.01270994282825128082040127705282437854604_prec,0.0)
    call test_one_flat([ (0._prec,0.), cmplx(10._prec/3._prec,kind=prec), (1._prec,0.), cmplx(10._prec/3._prec,kind=prec), (1._prec,0.) ],ref,'2.6')
Luca Naterop's avatar
Luca Naterop committed
137 138

    ! here the tests from the mathematica nb start
139 140
    ! --------------------------------------------------------------------------

141
    z = 1._prec/200; xchen = 0.3_prec;
Luca Naterop's avatar
Luca Naterop committed
142

143 144
    ref = (-0.00501254182354428204309373895836778138658_prec,0.0)
    call test_one_flat(cmplx([1./z,1.0_prec],kind=prec),ref,'3.1')
Luca Naterop's avatar
Luca Naterop committed
145

146 147
    ref = (-0.001501126126267145650881556104825717433567_prec,0.0)
    call test_one_flat(cmplx([1./(xchen*z),1.0_prec],kind=prec),ref,'3.2')
Luca Naterop's avatar
Luca Naterop committed
148

149 150
    ref = (-7.50286081781062398521044689345496255877E-4_prec,0.0)
    call test_one_flat(cmplx([(1+sqrt(1-z**2))/(xchen*z),1.0_prec],kind=prec),ref,'3.3')
Luca Naterop's avatar
Luca Naterop committed
151

152 153
    ref = (0.0074335969894769006767497160582533005912_prec,0.0)
    call test_one_flat(cmplx([-1._prec/xchen,-1._prec/xchen,1._prec,1._prec,1._prec],kind=prec),ref,'3.4')
Luca Naterop's avatar
Luca Naterop committed
154
    
155 156
    ref = (-8.40378597494480117619893114520705717404E-6_prec,0.0)
    call test_one_flat(cmplx([-1._prec/xchen,0._prec,-1._prec/xchen,1._prec/(xchen*z),1.0_prec],kind=prec),ref,'3.5')
Luca Naterop's avatar
Luca Naterop committed
157
    
158 159
    ref = (0.492575584745015188956099120118429131062_prec,2.63892140547433046711653044177665570864_prec)
    call test_one_flat(cmplx([-1._prec,-1._prec,z,z,1._prec],kind=prec),ref,'3.6')
Luca Naterop's avatar
Luca Naterop committed
160
    
161 162
    ref = (-0.001531771317885848316868878120481266022195_prec,-4.50039113668843745227527662921010599E-4_prec)
    call test_one_flat(cmplx([0._prec,0._prec,(1._prec-sqrt(1-z**2))/(xchen*z), 1._prec/(xchen*z),1._prec],kind=prec),ref,'3.7')
163 164 165 166 167


    ! here the chen integral tests start
    ! ----------------------------------------------------------------------

168 169 170
    z = 1/100._prec
    ref = (-1.203972804325935992622746217761838502954_prec,0.0)
    call test_one_flat((/ (0._prec,0.), (0.3_prec,0.) /),ref,'4.1')
171

172 173
    ref = (10.6037962209567960211233327771880353832_prec,0.0)
    call test_one_flat( (/ (0._prec,0.), (0._prec,0.), (0.01_prec,0.) /),ref,'4.2')
174

175 176
    ref = (5.0429900651807654463744449762857296353E-4_prec,0.0)
    call test_one_flat((/(100._prec,0.), (1._prec,0.), (0.3_prec,0.) /),ref,'4.3')
177

178 179
    ref = (0.05630861877185110944118569266547272128_prec,0.0)
    call test_one_flat( (/ (1._prec,0.), (0._prec,0.), (0.01_prec,0.) /) ,ref,'4.4')
180

181 182
    ref = (0.0403289515087272334386799760345612761274_prec,0.29227096475688417627233177065883250904_prec)
    call test_one_flat([cmplx(z,sqrt(1._prec-z**2),kind=prec), (0.3_prec,0._prec)],ref,'4.5')
183

184 185
    ref = (2.52106450983403521007132254962389653479E-5_prec,0.)
    call test_one_flat(cmplx([1._prec/z,  (1+sqrt(1-z**2))/z, 1._prec], kind=prec),ref,'4.6')
186

187 188
    ref = (0.0556470547634805907802456701296765199105_prec,0.)
    call test_one_flat((/ (-1._prec,0.), (0.01_prec,0.), (0._prec,0.), (0.01_prec,0.) /),ref,'4.7')
189

190 191
    ref = (3.79489584700899518570534417045601144619E-5_prec,0.)
    call test_one_flat((/ (100._prec,0.), (100._prec,0.), (1._prec,0.), (0._prec,0.), (1._prec,0.) /),ref,'4.8')
192

193 194
    ref = (7.57428443370789293292663421167820669347E-5_prec,0.)
    call test_one_flat((/ (100._prec,0.), (1._prec,0.), (100._prec,0.), (0._prec,0.), (1._prec,0.) /),ref,'4.9')
195

196 197
    ref = (1.88019115860102242751826218248903485831_prec,2.50943413859806280406885100241705929527_prec)
    call test_one_flat((/ (0.01_prec,0.), (-1.0_prec,0.), (0.01_prec,0.), (1._prec,0.) /),ref,'4.10')
198

199 200
    ref = (-0.0125391083150538283259376060156453327489_prec,-0.0154142501684376766660782277248386124052_prec)
    call test_one_flat(cmplx([z, (1+sqrt(1-z**2))/z, z, 1._prec],kind=prec),ref,'4.11')
201

ulrich_y's avatar
ulrich_y committed
202

ulrich_y's avatar
ulrich_y committed
203
    ! Increase test coverage
204
    ref = (-1.203972804325935992622746217761838502953611_prec, pi)
ulrich_y's avatar
ulrich_y committed
205
    call test_one_flat((/ (0._prec,0.), (-0.3_prec,0.) /),ref,'E.1')
ulrich_y's avatar
ulrich_y committed
206

207
    ref = (-0.539404830685964989184608085018712655769250_prec, 1.0303768265243124637877433270311515319634364_prec)
ulrich_y's avatar
ulrich_y committed
208
    call test_one_flat([(0._prec,0._prec), (0.3_prec, 0.5_prec)],ref,'E.2')
209 210

    ref = (3.79489584700899518570534417045601144619E-5_prec,0.)
ulrich_y's avatar
ulrich_y committed
211
    print*, '  ', 'testing GPL ', 'E.3', ' ...'
212
    res = G((/100._prec, 100._prec, 1._prec, 0._prec, 1._prec/))
ulrich_y's avatar
ulrich_y committed
213 214
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
215
    print*, '  ', 'testing GPL ', 'E.4', ' ...'
ulrich_y's avatar
ulrich_y committed
216 217 218
    res = G((/100,100,1,0,1/))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
219
    print*, '  ', 'testing GPL ', 'E.5', ' ...'
ulrich_y's avatar
ulrich_y committed
220 221 222
    res = G_superflatn(toinum((/100,100,1,0,1/)), 5)
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
223
    print*, '  ', 'testing GPL ', 'E.6', ' ...'
224
    res = G((/100._prec,100._prec,1._prec,0._prec/), 1._prec)
ulrich_y's avatar
ulrich_y committed
225 226
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
227
    print*, '  ', 'testing GPL ', 'E.7', ' ...'
228
    res = G((/(100._prec,0.),(100._prec,0.),(1._prec,0.),(0._prec,0.)/), (1._prec,0.))
ulrich_y's avatar
ulrich_y committed
229 230
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
231
    print*, '  ', 'testing GPL ', 'E.8', ' ...'
232
    res = G( (/1,1,1,1/) , (/ 100._prec, 100._prec, 1._prec, 0._prec /), 1._prec)
ulrich_y's avatar
ulrich_y committed
233 234 235
    call check(res,ref)

    !ref = cmplx(0.01592795952537145)
236
    ref = (0.0485035531168477720320998551314479928648_prec,-2.82326885118931135859594797186107474877_prec)
ulrich_y's avatar
ulrich_y committed
237
    print*, '  ', 'testing GPL ', 'E.9', ' ...'
238
    res = G( (/ 3, 2 /), toinum((/ 1.3_prec, 1.1_prec /)), toinum(4._prec) )
ulrich_y's avatar
ulrich_y committed
239 240
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
241
    print*, '  ', 'testing GPL ', 'E.10', ' ...'
242
    res = G( (/ 3, 2 /), (/ 1.3_prec, 1.1_prec /), 4._prec )
ulrich_y's avatar
ulrich_y committed
243 244
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
245
    print*, '  ', 'testing GPL ', 'E.11', ' ...'
246 247 248
    res = G( (/ 3, 2 /), (/ (1.3_prec,0.), (1.1_prec,0.) /), (4._prec,0.) )
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
249
    print*, '  ', 'testing GPL ', 'E.12', ' ...'
250 251
    ref = (0.107961231635965271810236308594279564021882373573100109007684973130605145011141_prec, 0.)
    res = G((/ 0._prec, 1.3_prec, 0._prec, 0._prec, 4._prec, 1.1_prec /))
ulrich_y's avatar
ulrich_y committed
252 253
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
254
    print*, '  ', 'testing GPL ', 'E.13', ' ...'
ulrich_y's avatar
ulrich_y committed
255 256 257 258
    ref = (-1.229846637025798984989235266565630713100_prec,-1.05705870670612913651142424798975000765_prec)
    res = G([inum(2._prec, +1), inum(7._prec, +1)], inum(5._prec, +1))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
259
    print*, '  ', 'testing GPL ', 'E.14', ' ...'
ulrich_y's avatar
ulrich_y committed
260 261 262 263
    ref = (-1.229846637025798984989235266565630713100_prec,+1.05705870670612913651142424798975000765_prec)
    res = G([inum(2._prec, -1), inum(7._prec, +1)], inum(5._prec, +1))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
264
    print*, '  ', 'testing GPL ', 'E.15', ' ...'
ulrich_y's avatar
ulrich_y committed
265 266 267 268
    ref = (0.2982254208688088675254638762780704094718_prec,0.)
    res = G([inum(2._prec, -1), inum(7._prec, +1)], inum(-5._prec, -1))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
269
    print*, '  ', 'testing GPL ', 'E.17', ' ...'
ulrich_y's avatar
ulrich_y committed
270 271 272 273
    ref = (0.190800137777535619036913153766083992418_prec, 0.)
    res = G((/ 6, 1, 1 /))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
274
    print*, '  ', 'testing GPL ', 'E.18', ' ...'
ulrich_y's avatar
ulrich_y committed
275 276 277 278 279
    ref = (0.058192342415778512650117048874978455691_prec, 0.)
    res = G((/ 6, 1, -1 /))
    call check(res,ref)

    ref = cmplx(log(0.5_prec), pi,kind=prec)
ulrich_y's avatar
ulrich_y committed
280
    print*, '  ', 'testing GPL ', 'E.19a', ' ...'
ulrich_y's avatar
ulrich_y committed
281 282
    res = G((/ inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
283
    print*, '  ', 'testing GPL ', 'E.19b', ' ...'
ulrich_y's avatar
ulrich_y committed
284 285
    res = G((/ inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
286
    print*, '  ', 'testing GPL ', 'E.19c', ' ...'
ulrich_y's avatar
ulrich_y committed
287 288
    res = G((/ inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
289
    print*, '  ', 'testing GPL ', 'E.19d', ' ...'
ulrich_y's avatar
ulrich_y committed
290 291 292
    res = G((/ inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
293
    print*, '  ', 'testing GPL ', 'E.19e', ' ...'
ulrich_y's avatar
ulrich_y committed
294 295
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
296
    print*, '  ', 'testing GPL ', 'E.19f', ' ...'
ulrich_y's avatar
ulrich_y committed
297 298
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
299
    print*, '  ', 'testing GPL ', 'E.19g', ' ...'
ulrich_y's avatar
ulrich_y committed
300 301
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
302
    print*, '  ', 'testing GPL ', 'E.19h', ' ...'
ulrich_y's avatar
ulrich_y committed
303 304
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
Luca Naterop's avatar
Luca Naterop committed
305

Luca Naterop's avatar
Luca Naterop committed
306

ulrich_y's avatar
ulrich_y committed
307
    ref = -(2.374395270272480200677499763071638424_prec, - 1.273806204919600530933131685580471698_prec)
ulrich_y's avatar
ulrich_y committed
308
    print*, '  ', 'testing GPL ', 'E.20a', ' ...'
ulrich_y's avatar
ulrich_y committed
309 310
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
311
    print*, '  ', 'testing GPL ', 'E.20b', ' ...'
ulrich_y's avatar
ulrich_y committed
312 313
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
314
    print*, '  ', 'testing GPL ', 'E.20c', ' ...'
ulrich_y's avatar
ulrich_y committed
315 316
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
317
    print*, '  ', 'testing GPL ', 'E.20d', ' ...'
ulrich_y's avatar
ulrich_y committed
318 319 320
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
321
    print*, '  ', 'testing GPL ', 'E.20e', ' ...'
ulrich_y's avatar
ulrich_y committed
322 323
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
324
    print*, '  ', 'testing GPL ', 'E.20f', ' ...'
ulrich_y's avatar
ulrich_y committed
325 326
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
327
    print*, '  ', 'testing GPL ', 'E.20g', ' ...'
ulrich_y's avatar
ulrich_y committed
328 329
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
330
    print*, '  ', 'testing GPL ', 'E.20h', ' ...'
ulrich_y's avatar
ulrich_y committed
331 332
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
333 334 335 336


    ! Thanks to Roman Zwicky and Ben Pullin for these tests
    ref = (0.392707112217551702879328061598355445_prec, - 1.274969448494380061814943180491890080_prec)
ulrich_y's avatar
ulrich_y committed
337
    call test_one_flat((/ (0._prec,0.), (1._prec,0.), (0._prec, 1.5_prec) /),ref,'F.1')
338 339

    ref = (2.09624324167194065961839363660174566_prec, 0.62644605179087454819421905065313983_prec)
ulrich_y's avatar
ulrich_y committed
340
    call test_one_flat( (/ (1._prec,0.), (-1._prec,0.), (2.5_prec, -2.4_prec) /) , ref, 'F.2')
341 342

    ref = conjg(ref)
ulrich_y's avatar
ulrich_y committed
343
    call test_one_flat( (/ (1._prec,0.), (-1._prec,0.), (2.5_prec, +2.4_prec) /) , ref, 'F.3')
344 345

    ref = (0.6874619495289224183167486785286777066_prec, -1.8261934916106546308783576818077830287_prec)
ulrich_y's avatar
ulrich_y committed
346
    call test_one_flat( (/ (-1._prec,0.), (1._prec,0.), (2.5_prec, +2.4_prec) /) , ref, 'F.4')
347 348

    ref = (0.320016770069038023050391296154549_prec , 0.064263450879286017225353602844105_prec)
ulrich_y's avatar
ulrich_y committed
349
    call test_one_flat( (/ (0._prec,0.), (1._prec,0.), (-1._prec,0.), (0.3_prec, -1.2_prec) /) , ref, 'F.5')
350 351

    ref = conjg(ref)
ulrich_y's avatar
ulrich_y committed
352
    call test_one_flat( (/ (0._prec,0.), (1._prec,0.), (-1._prec,0.), (0.3_prec, +1.2_prec) /) , ref, 'F.6')
353 354

    ref = (0.8382358254435272068734922352_prec, - 0.3702062785327198487992149341_prec)
ulrich_y's avatar
ulrich_y committed
355
    call test_one_flat( (/ (0._prec,0.), (1._prec,0.), (-1._prec,0.), (1.1_prec,2._prec) /) , ref, 'F.7')
356 357

    ref = (0.185156872427485220072774923047908301422_prec,-0.249989197161429773744237773045427322847_prec)
ulrich_y's avatar
ulrich_y committed
358
    call test_one_flat( (/ (1._prec,0.), (-1._prec,0.), (2.3_prec,-1._prec), (1.1_prec,-0.1_prec) /) , ref, 'F.8')
359 360

    ref = (-1.11161035333074623447215317094270858897_prec,-0.875225967273157437459269290416981432294_prec)
ulrich_y's avatar
ulrich_y committed
361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386
    call test_one_flat( (/ (1._prec,0.), (1.1_prec,-0.1_prec), (1._prec,0.), (2.3_prec,-1.2_prec), (3.4_prec,9.7_prec) /) , ref, 'F.9')

    ref = (-0.1109203962012759970680116512959180346924_prec,0.114779444698808046532782358409267479097_prec)
    call test_one_flat( (/ (0._prec, 0.) , (-0.1_prec, - 0.3_prec), (1._prec,0.), (0.5_prec, 0._prec) /) , ref, 'F.10')

    ref = (-0.108930339927322068475622176344923224046_prec,0.1122893584726586134006078500824852031925_prec)
    call test_one_flat( (/ (-0.005_prec, 0.) , (-0.1_prec, - 0.3_prec), (1._prec,0.), (0.5_prec, 0._prec)/) , ref, 'F.11')

    ref = (-1.00642162671475190543135817048062920163_prec,0.75972801214517468120067717498723544972_prec)
    call test_one_flat( (/ (-0.005_prec, 0.) , (-0.1_prec, - 0.3_prec), (1._prec,0.), (1.5_prec, 0._prec)/) , ref, 'F.12')

    ref = (-0.764543307991553066235740993950606803628_prec, 0.54861507469010653128533842086292369104_prec)
    call test_one_flat( (/ (-0.1_prec, - 0.3_prec), (1._prec,0.), (-0.005_prec, 0.) , (0.5_prec, 0._prec)/) , ref, 'F.13')

    ref = (4.52004218277686921073832347986485146573_prec,-1.812384329980915889835338149845575221705_prec)
    call test_one_flat( (/ (1.005_prec,0.), (0._prec,0.), (1.1_prec,0.3_prec), (1._prec,0.) /) , ref, 'F.14')

    ref = (1.6706939608077788063504210821631008732_prec,-0.631820518538505729899318406678510198266_prec)
    call test_one_flat( (/ (1.1_prec,0.), (0._prec,0.), (1.1_prec,0.3_prec), (1._prec,0.) /) , ref, 'F.10')

    ref = (0.198693810038777868937967610639933071867_prec,+0.72807783717224553483870820314753060868_prec)
    call test_one_flat( (/ (0._prec,0.), (1.0_prec,0._prec), (1.1_prec,0.), (1.1_prec,+0.3_prec) /) , ref, 'F.11')


    ref = (1.6706939608077788063504210821631008732_prec,-0.631820518538505729899318406678510198266_prec)
    call test_one_flat( (/ (1.1_prec,0.), (0._prec,0._prec), (1.1_prec,+0.3_prec), (1._prec,0.) /) , ref, 'F.12')
387

ulrich_y's avatar
ulrich_y committed
388 389
    ref = (0.1810339553848393655117844582129810543006_prec,+0.82543851493400141056822586864697825094_prec)
    call test_one_flat( (/ (0._prec,0.), (1.0_prec,0._prec), (1._prec,0.), (1.1_prec,+0.3_prec) /) , ref, 'F.13')
390

ulrich_y's avatar
ulrich_y committed
391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414

    ref = (-3.028608056170828558746818459366566689807_prec,-0.52999686156911157999452896083600882192_prec)
    call test_one_flat( (/ (1._prec,0._prec), (0.0_prec,0._prec), (0._prec,1._prec), (1.1_prec,0.0_prec) /) , ref, 'F.14')

    ! Here the branch cut matters, in Mathematica this is entered as G(1_-, 0, -I, -1.1)
    ref = (0.0538179677556747874824671812943783830465_prec,0.340519960719077282936661573786600733497_prec)
    print*, '  ', 'testing GPL ', 'F.15', ' ...'
    res = G((/ inum(1._prec,-1), izero, inum((0._prec,-1._prec), -1) /), inum(-1.1_prec, +1))
    call check(res,ref)

    ref = (0.0075181720252360930529649927100295277234_prec,0.009279206321129807482628652573319115651_prec)
    call test_one_flat( (/ (0._prec,0.), (0._prec,0.), (1.0_prec,0._prec), (0.1_prec,-0.1_prec), (0.11_prec,0._prec) /) , ref, 'F.16')

    ref = (0.0135493735561310633082925053080300174678_prec,+0.01851696361979639851211857931403696163_prec)
    call test_one_flat( (/ (0._prec,0.), (0._prec,0.), (1.0_prec,0._prec), (0.1_prec,-0.1_prec), (0.15_prec,0._prec) /) , ref, 'F.17')

    ref = (-2.41176903997917647290181947143140323970805e-4_prec,0.00129690184788114442363777169655122)
    call test_one_flat( (/ (0._prec,0.), (1._prec,0.), (0.5_prec,-0.1_prec), (-1.3_prec,10.2_prec), (0.4_prec,0._prec) /) , ref, 'F.18')

    ref = (-0.003113755848404291707649322614093090379815_prec,-4.6914973893242872720303973859400498154E-4_prec)
    call test_one_flat( (/ (0._prec,0.), (1._prec,0.), (0.1_prec,-0.1_prec), (-1.3_prec,10.2_prec), (0.4_prec,0._prec) /) , ref, 'F.19')



Luca Naterop's avatar
Luca Naterop committed
415 416
  end subroutine do_GPL_tests
  
ulrich_y's avatar
ulrich_y committed
417 418 419



ulrich_y's avatar
ulrich_y committed
420
#ifdef HAVE_GINAC
421
#ifdef HAVE_MM
ulrich_y's avatar
ulrich_y committed
422 423

  function evalt(arr, what)
ulrich_y's avatar
ulrich_y committed
424 425 426
#if KINDREAL==16
    use ieps, only: inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
427
    implicit none
ulrich_y's avatar
ulrich_y committed
428 429
    complex(kind=prec) :: arr(:), evalt
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
430 431 432 433 434
    integer what, i, l

    evalt =0.

    do i=1,size(arr)
435
      if (abs(arr(i)) .gt. 1.e10) then ! isnan?
ulrich_y's avatar
ulrich_y committed
436 437 438 439 440 441 442 443 444 445 446 447
        l = i-1
        goto 123
      endif
    enddo
    l = size(arr)

123 continue
    if (l==0) return

    if (what .eq. 0) then
      evalt = G(arr(1:l))
    elseif (what.eq.1) then
ulrich_y's avatar
ulrich_y committed
448 449 450
#if KINDREAL==16
      evalt = cmplx(geval(inum2inum(toinum(arr(1:l))),l),kind=prec)
#else
ulrich_y's avatar
ulrich_y committed
451
      evalt = geval(toinum(arr(1:l)),l)
ulrich_y's avatar
ulrich_y committed
452
#endif
ulrich_y's avatar
ulrich_y committed
453
    endif
ulrich_y's avatar
ulrich_y committed
454 455 456
  end function

  subroutine perform_ginacv(n, args)
ulrich_y's avatar
ulrich_y committed
457
    use maths_functions, only:clearcache
ulrich_y's avatar
ulrich_y committed
458 459
    complex(kind=prec) :: args(:,:)
    integer i,n
ulrich_y's avatar
ulrich_y committed
460
    call clearcache
ulrich_y's avatar
ulrich_y committed
461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477

    do i=1,size(args,1)
      write(*,900) n,i
      call check(evalt(args(i,:),0),evalt(args(i,:),1))
    enddo

900 format("   testing GPL ",I1,".",I4.4," ...")

  end subroutine

  subroutine do_ginac_tests
    use maths_functions, only:clearcache
    use gtestchen   , only: inichen   =>args
    use gtestmuone  , only: inimuone  =>args
    use gtestmuonenp, only: inimuonenp=>args
    implicit none
    tol = 6.0e-6
478 479 480
    call perform_ginacv( 6, inichen   ((0.3_prec,0.), (0.1_prec,0.)) )
    call perform_ginacv( 7, inimuone  ((0.5_prec,0.), (0.6_prec,0.)) )
    call perform_ginacv( 8, inimuonenp((0.3_prec,0.), (0.6_prec,0.)) )
ulrich_y's avatar
ulrich_y committed
481

ulrich_y's avatar
ulrich_y committed
482
  end subroutine
ulrich_y's avatar
ulrich_y committed
483

ulrich_y's avatar
ulrich_y committed
484
#ifndef NOSPEED
ulrich_y's avatar
ulrich_y committed
485 486 487 488
  subroutine do_one_speed_test(args, u, msg)
    use maths_functions, only:clearcache
    implicit none
    complex(kind=prec) :: args(:,:,:),res
ulrich_y's avatar
ulrich_y committed
489 490
    integer(kind=8) cstart, cend, count_rate
    real(kind=prec) :: time(2), ttime(2)
ulrich_y's avatar
ulrich_y committed
491 492 493 494 495
    integer i,j, u
    character,parameter :: cr = achar(13)
    character(len=*) msg
    do j=1,size(args,1)
      ! try function a bunch of times
ulrich_y's avatar
ulrich_y committed
496
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
497 498 499
      do i=1,size(args,3)
        res=evalt(args(j,:,i),0)
      enddo
ulrich_y's avatar
ulrich_y committed
500
      call system_clock(cend, count_rate=count_rate)
501
      time(1) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
502
      if (time(1).lt.zero) print*,j, cend-cstart,count_rate
ulrich_y's avatar
ulrich_y committed
503 504
      ttime(1) = ttime(1) + time(1)

ulrich_y's avatar
ulrich_y committed
505
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
506 507 508
      do i=1,size(args,3)
        res=evalt(args(j,:,i),1)
      enddo
ulrich_y's avatar
ulrich_y committed
509
      call system_clock(cend, count_rate=count_rate)
510
      time(2) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
511 512 513 514 515 516 517 518
      if (time(2).lt.zero) print*,j
      ttime(2) = ttime(2) + time(2)

      write(u,*) time

      write( * , 900, advance='no' ) cr, j, size(args,1), msg
    enddo
    print*,
519
    print*,ttime
ulrich_y's avatar
ulrich_y committed
520 521
    write(*,901) msg, size(args,1)/ttime(2)/1000., size(args,1)/ttime(1)/1000., int(ttime(2)/ttime(1))

ulrich_y's avatar
ulrich_y committed
522 523 524 525 526 527
    ! Lets to another, fair comparison
    call system_clock(cstart, count_rate=count_rate)
    do j=1,size(args,1)
      res=evalt(args(j,:,1),0)
    enddo
    call system_clock(cend, count_rate=count_rate)
528
    ttime(1) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
529 530 531 532 533 534

    call system_clock(cstart, count_rate=count_rate)
    do j=1,size(args,1)
      res=evalt(args(j,:,1),1)
    enddo
    call system_clock(cend, count_rate=count_rate)
535
    ttime(2) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
536 537 538

    write(*,901) msg, size(args,1)/ttime(2)/1000., size(args,1)/ttime(1)/1000., int(ttime(2)/ttime(1))

ulrich_y's avatar
ulrich_y committed
539
900 FORMAT(a , 'Function ',i4,'/',i4,' for ',a)
540
901 format('Evaluating ',A,' using GiNaC at ',F9.2,'kG/s and handyG at ',F9.2,'kG/s (',I3,'x)')
ulrich_y's avatar
ulrich_y committed
541 542 543 544
  end subroutine

  subroutine do_timing_tests(n)
    use gtestchen   , only: inichen   =>args
ulrich_y's avatar
ulrich_y committed
545
    use gtestchenff , only: inichenff =>args
ulrich_y's avatar
ulrich_y committed
546 547 548 549 550 551
    use gtestmuone  , only: inimuone  =>args
    use gtestmuonenp, only: inimuonenp=>args
    implicit none
    integer, intent(in) :: n
    integer i
    complex(kind=prec) :: cargs( 1399,5,n)
ulrich_y's avatar
ulrich_y committed
552
    complex(kind=prec) :: fargs(  540,5,n)
ulrich_y's avatar
ulrich_y committed
553 554 555 556 557 558 559 560 561
    complex(kind=prec) :: pargs(  198,5,n)
    complex(kind=prec) :: nargs( 1733,5,n)
    real(kind=prec) :: z, x, y, w
    integer ranseed
    ranseed = 233123

    do i=1,n
      z = ran2(ranseed) / 2.
      x = ran2(ranseed)*(1-z) + z
562 563
      cargs(:,:,i) = inichen  (cmplx(x,kind=prec), cmplx(z,kind=prec))
      fargs(:,:,i) = inichenff(cmplx(x,kind=prec), cmplx(z,kind=prec))
ulrich_y's avatar
ulrich_y committed
564 565 566

      w = ran2(ranseed) ! 0<w<1
      z = ran2(ranseed) * (sqrt(1-w+w**2)-sqrt(w)) + sqrt(w)
567
      nargs(:,:,i) = inimuonenp(cmplx(w,kind=prec), cmplx(z,kind=prec))
ulrich_y's avatar
ulrich_y committed
568 569 570

      x = ran2(ranseed)
      y = ran2(ranseed)
571
      pargs(:,:,i) = inimuone(cmplx(x,kind=prec), cmplx(y,kind=prec))
ulrich_y's avatar
ulrich_y committed
572 573
    enddo

ulrich_y's avatar
ulrich_y committed
574 575
    cargs(1181,:,:)=1.e15
    fargs(367,:,:)=1.e15
ulrich_y's avatar
ulrich_y committed
576 577 578


    open(unit=9, file="stats.txt")
ulrich_y's avatar
ulrich_y committed
579 580
    write(9,*) "Chen form factor"
    call do_one_speed_test(fargs,9,"Chen FF")
ulrich_y's avatar
ulrich_y committed
581 582 583 584 585 586 587 588 589
    write(9,*) "Chen"
    call do_one_speed_test(cargs,9,"Chen")
    write(9,*) "MUonE-planar"
    call do_one_speed_test(pargs,9,"Muone planar")
    write(9,*) "MUonE-non-planar"
    call do_one_speed_test(nargs,9,"Muone non planar")
    close(unit=9)

  end subroutine
ulrich_y's avatar
ulrich_y committed
590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682

  subroutine fyshuffle(list, seed)
    real(kind=prec) :: list(:), tmp
    integer j,i, seed

    do i=size(list),2,-1
      j = int(ran2(seed) * i) + 1
      tmp = list(j)
      list(j) = list(i)
      list(i) = tmp
    enddo
  end subroutine

  function test_one_high_weight(length, fzero, funiq, alphabet, seed, what) result(res)
    integer length,seed, i, what
    integer nzero, nuniq, nsmall, c
    real(kind=prec) :: fzero, funiq, fsmall
    real(kind=prec) :: arguments(length), alphabet(10)
    complex(kind=prec) :: Geval, res
    nzero  = nint(fzero*length)
    nuniq  = nint(       funiq*length)
    nsmall = nint(fsmall*funiq*length)

    if (nzero+nuniq > length) nuniq = length - nzero

    !do i=1,nsmall
    !  alphabet(i) = 0.8_prec * ran2(seed)
    !enddo
    !do i=nsmall+1,nuniq
    !  alphabet(i) = 1.5_prec/(0.15_prec + ran2(seed))
    !enddo

    arguments(1:nzero) = 0._prec
    do i=1,nuniq
      arguments(nzero+i) = alphabet(i)
    enddo
    do i=nzero+nuniq,length
      c = int(ran2(seed)*nuniq)+1
      arguments(i) = alphabet(c)
    enddo

    call fyshuffle(arguments, seed)
    do while(arguments(length)<zero)
      call fyshuffle(arguments, seed)
    enddo
    write(9,*)"args",arguments
    if (what.eq.1) then
      res = G(arguments, 1._prec)
    elseif (what.eq.2) then
      res = Geval(cmplx([arguments, 1._prec],kind=prec),length+1)
    endif
  end function


  subroutine do_high_weight_tests
    integer, parameter :: ntests=10
    real(kind=prec), parameter :: fzero = 0.28
    real(kind=prec), parameter :: funiq = 0.63
    real(kind=prec), parameter :: fsmal = 0.42
    integer(kind=8) cstart, cend, count_rate
    complex(kind=prec) :: res(ntests,2)
    real(kind=prec) :: alphabet(10)

    integer i,j
    integer seed, seedold
    seed = 123123

    open(unit=9, file="stats-hi.txt")
    alphabet = (/ 0.3, 0.6, 1.8, 5.3, 3.6, 8.3, 0.1, 0.4, 0.2, 10.2 /)
    do i=2,7
      seedold=seed
      do j=1,ntests
        call system_clock(cstart, count_rate=count_rate)
        res(j,1) = test_one_high_weight(i, fzero, funiq, alphabet, seed,1)
        call system_clock(cend, count_rate=count_rate)
        write(9,*) "time",i,j,real(cend-cstart)/real(count_rate,kind=prec)!/ntests
      enddo

      seed=seedold
      do j=1,ntests
        call system_clock(cstart, count_rate=count_rate)
        res(j,2) = test_one_high_weight(i, fzero, funiq, alphabet, seed,2)
        call system_clock(cend, count_rate=count_rate)
        write(9,*)"ginac",i,j,real(cend-cstart)/real(count_rate,kind=prec)!/ntests
      enddo

      write(9,*)"del",abs(res(:,1)-res(:,2)) / tol
      write(9,*)
      write(9,*)
    enddo
    close(unit=9)
  end subroutine

ulrich_y's avatar
ulrich_y committed
683
#endif
ulrich_y's avatar
ulrich_y committed
684
#ifdef LONG_TEST
ulrich_y's avatar
ulrich_y committed
685

ulrich_y's avatar
ulrich_y committed
686
  SUBROUTINE DO_LONG_TEST
ulrich_y's avatar
ulrich_y committed
687 688 689
#if KINDREAL==16
    use ieps,only:inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
690 691 692 693 694
    implicit none
    integer,parameter :: nzero = 10
    integer,parameter :: nieps = 30
    integer,parameter :: ncmpl = 30

695
    integer,parameter :: perweight(4) = (/ 100, 500, 30000, 50000 /)
ulrich_y's avatar
ulrich_y committed
696 697 698 699
    real(kind=prec), parameter :: rrange = 1.5
    type(inum), dimension(nzero+nieps+ncmpl) :: basis
    type(inum), dimension(size(perweight)) :: args
    type(inum), parameter :: ione = inum((1._prec,0._prec), di0)
ulrich_y's avatar
ulrich_y committed
700
    integer i, j, w, seed, oldseed
ulrich_y's avatar
ulrich_y committed
701
    integer(kind=1) i0
ulrich_y's avatar
ulrich_y committed
702 703 704
    real(kind=prec) :: v, maxd
    complex(kind=prec) :: ans(2)
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
705
    character,parameter :: cr = achar(13)
ulrich_y's avatar
ulrich_y committed
706
    maxd=0._prec
ulrich_y's avatar
ulrich_y committed
707 708

    seed = 112312
ulrich_y's avatar
ulrich_y committed
709
    tol = 0.01
ulrich_y's avatar
ulrich_y committed
710

ulrich_y's avatar
ulrich_y committed
711
    open(unit=9, action='write', form='unformatted', file="long-test.txt")
ulrich_y's avatar
ulrich_y committed
712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729
    basis(1:nzero) = izero

    do i=nzero+1,nzero+nieps
      if (ran2(seed).gt.0.5) then
        i0 = +1_1
      else
        i0 = -1_1
      endif
      v = 2*rrange*(ran2(seed) - 0.5)
      basis(i) = inum(cmplx(v,kind=prec), i0)
    enddo
    do i=nzero+nieps+1,nzero+nieps+ncmpl
      v = 2*rrange*(ran2(seed) - 0.5)
      basis(i) = toinum(v*exp(i_*2*pi*ran2(seed)))
    enddo

    do w=1,size(perweight)
      do i=1,perweight(w)
ulrich_y's avatar
ulrich_y committed
730 731
        oldseed = seed
        write( * , 900, advance='no' ) cr, i, perweight(w), w
ulrich_y's avatar
ulrich_y committed
732
        args(1:w) = basis((/ (1+int(size(basis)*ran2(seed)),j=1,w) /))
ulrich_y's avatar
ulrich_y committed
733 734 735 736 737
#if KINDREAL==16
        ans(1) = cmplx(geval(inum2inum([args(1:w),ione]), w+1),kind=prec)
#else
        ans(1) = cmplx(geval([args(1:w),ione], w+1),kind=prec)
#endif
ulrich_y's avatar
ulrich_y committed
738
        ans(2) = G(args(1:w),ione)
739 740 741 742 743 744
        if ((abs(ans(1)) .gt. 1.e10).or.(abs(ans(2)) .gt. 1.e10)) then
          print*," can't deal with",args," (seed was",oldseed,')'
          print*,"GiNaC : ",ans(1)
          print*,"handyG: ",ans(2)
          cycle
        endif
ulrich_y's avatar
ulrich_y committed
745
        write(9) w,i,oldseed, abs(ans(1)-ans(2))
746
        flush(9)
ulrich_y's avatar
ulrich_y committed
747
        if(abs(ans(1)-ans(2)) > maxd) maxd = abs(ans(1)-ans(2))
ulrich_y's avatar
ulrich_y committed
748 749
        if(abs(ans(1)-ans(2)) > tol) goto 123
      enddo
ulrich_y's avatar
ulrich_y committed
750 751
      print*,' done. largest=',maxd
      maxd=0.
ulrich_y's avatar
ulrich_y committed
752
    enddo
ulrich_y's avatar
ulrich_y committed
753
    print*,maxd
ulrich_y's avatar
ulrich_y committed
754

ulrich_y's avatar
ulrich_y committed
755 756
    close(9)

ulrich_y's avatar
ulrich_y committed
757 758 759 760 761
    return
123 continue
    print*,"Failed with delta",abs(ans(1)-ans(2))
    print*,"Offending G was",args(1:w)
    print*,ans
ulrich_y's avatar
ulrich_y committed
762 763 764
    close(9)

900 FORMAT(a , 'Testing ',i5,'/',i5,' GPLs with w=',i1)
ulrich_y's avatar
ulrich_y committed
765 766 767


  END SUBROUTINE
ulrich_y's avatar
ulrich_y committed
768
#endif
ulrich_y's avatar
ulrich_y committed
769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786
      FUNCTION RAN2(randy)

            ! This is the usual "random"

      implicit none
      real(kind=prec) :: MINV,RAN2
      integer m,a,Qran,r,hi,lo,randy
      PARAMETER(M=2147483647,A=16807,Qran=127773,R=2836)
      PARAMETER(MINV=0.46566128752458e-09)
      HI = RANDY/Qran
      LO = MOD(RANDY,Qran)
      RANDY = A*LO - R*HI
      IF(RANDY.LE.0) RANDY = RANDY + M
      RAN2 = RANDY*MINV
      END FUNCTION RAN2



ulrich_y's avatar
ulrich_y committed
787

788
#endif
ulrich_y's avatar
ulrich_y committed
789
#endif
790 791 792
  ! subroutine do_shuffle_tests() 
  !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
  !   complex(kind=prec) :: w(2) = cmplx((/3,4/))
793

794 795
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
796

Luca's avatar
Luca committed
797
END PROGRAM TEST
Luca's avatar
Luca committed
798 799

! In terminal kann man den exit code bekommen via echo $?