test.f90 30.3 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
ulrich_y's avatar
ulrich_y committed
12
  use ttools
Luca's avatar
Luca committed
13
  implicit none
ulrich_y's avatar
ulrich_y committed
14 15
  integer :: i
  character(len=32) :: arg
16

ulrich_y's avatar
ulrich_y committed
17

18
  i = 1
ulrich_y's avatar
ulrich_y committed
19 20 21 22 23 24 25
  do
    call get_command_argument(i, arg)
    if (len_trim(arg) == 0) exit

    ! parse verbosity
    select case(trim(arg))
      case('-verb')
ulrich_y's avatar
ulrich_y committed
26
#ifdef DEBUG
ulrich_y's avatar
ulrich_y committed
27
        verb = readint(trim(arg),i)
ulrich_y's avatar
ulrich_y committed
28
#else
ulrich_y's avatar
ulrich_y committed
29
        call iprint("Argument -verb is not available, compile with --debug", 2)
ulrich_y's avatar
ulrich_y committed
30
#endif
31 32 33 34 35 36 37 38 39 40
      case('-mpl-test')
        tol = zero * 1.e5_prec
        call do_MPL_tests
      case('-gpl-test')
        tol = zero * 1.e5_prec
        call do_GPL_tests
      case('-chen-test')
        tol = 8.0e-7
        tests_successful = tests_successful .and. do_chen_test()

ulrich_y's avatar
ulrich_y committed
41
#if defined(HAVE_GINAC) && defined(HAVE_MM)
42 43 44
      case('-ginac-tests')
        tol = 8.0e-7
        call do_ginac_tests
ulrich_y's avatar
ulrich_y committed
45
      case('-speed-tests')
ulrich_y's avatar
ulrich_y committed
46
        call do_timing_tests(readint(trim(arg),i))
ulrich_y's avatar
ulrich_y committed
47
      case('-hw-tests')
48 49
        tol = 8.0e-7
        call do_high_weight_tests
ulrich_y's avatar
ulrich_y committed
50 51
#else
      case('-ginac-tests', '-speed-tests', '-hw-tests')
ulrich_y's avatar
ulrich_y committed
52
        call iprint("Argument "//trim(arg)//" is not available, compile with --with-ginac --with-mcc", 2)
ulrich_y's avatar
ulrich_y committed
53 54 55
#endif

#ifdef HAVE_GINAC
ulrich_y's avatar
ulrich_y committed
56
      case('-long-test')
57 58
        tol = 8.0e-7
        call do_long_test
ulrich_y's avatar
ulrich_y committed
59 60
#else
      case('-long-test')
ulrich_y's avatar
ulrich_y committed
61
        call iprint("Argument "//trim(arg)//" is not available, compile with --with-ginac", 2)
62 63
#endif
#ifdef DEBUG
ulrich_y's avatar
ulrich_y committed
64
      case('-report')
65 66 67 68 69 70
        verb = 1000
        tol = zero * 1.e5_prec
        call do_MPL_tests
        call do_GPL_tests
        tol = 8.0e-7
        tests_successful = tests_successful .and. do_chen_test()
ulrich_y's avatar
ulrich_y committed
71
#if defined(HAVE_GINAC) && defined(HAVE_MM)
72
        call do_ginac_tests
ulrich_y's avatar
ulrich_y committed
73
#endif
ulrich_y's avatar
ulrich_y committed
74 75 76
#else
      case('-report')
        call iprint("Argument -report is not available, compile with --debug", 2)
77 78
#endif

ulrich_y's avatar
ulrich_y committed
79 80
      case('--help','--h','-h','-help')
        call printhelp
ulrich_y's avatar
ulrich_y committed
81
        stop 0
ulrich_y's avatar
ulrich_y committed
82
      case default
ulrich_y's avatar
ulrich_y committed
83
        call iprint("Unknown argument "//trim(arg)//". Try -h to get help",2)
ulrich_y's avatar
ulrich_y committed
84 85 86 87

    end select
    i = i+1
  end do
ulrich_y's avatar
ulrich_y committed
88

89
  ! call do_shuffle_tests() ! put this somewhere else
Luca's avatar
tests  
Luca committed
90 91 92


  if(tests_successful) then
ulrich_y's avatar
ulrich_y committed
93
    call iprint('All tests passed. ', 0)
Luca's avatar
tests  
Luca committed
94
  else 
ulrich_y's avatar
ulrich_y committed
95
    call iprint('Some tests failed. ',2)
Luca's avatar
tests  
Luca committed
96
  end if
97

98
CONTAINS
99
   
ulrich_y's avatar
ulrich_y committed
100

ulrich_y's avatar
ulrich_y committed
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
  subroutine printhelp
    character(len=32) :: arg
    call get_command_argument(0,arg)

#ifdef DEBUG
    print*,"Usage: "//trim(arg)//" [-verb <n>] [opts]"
#else
    print*,"Usage: "//trim(arg)//" [opts]"
#endif
    print*,"Runs a set of tests for handyG"
    print*,""
    print*,"Possible tests are:"
    print*,"    -mpl-test        performs tests on the series expansion of"
    print*,"                     convergent MPLs"
    print*,"    -gpl-test        tests GPLs and their reduction. This includes"
    print*,"                     real, ieps, and complex arguments"
    print*,"    -chen-test       compares all GPLs needed in [1811.06461] to"
    print*,"                     reference values"
#ifdef HAVE_GINAC
#ifdef HAVE_MM
    print*,"    -ginac-tests     compare all GPLs needed by [1801.01033],"
    print*,"                     [1709.07435], and [1806.08241] to GiNaC"
    print*,"    -speed-tests <n> compare the evaluation speed of the"
    print*,"                     aforementioned GPLs to GiNaC by averaging over"
    print*,"                     <n> evaluations"
    print*,"    -hw-tests        compares `random' GPLs with high weight to GiNaC"
#endif
    print*,"    -long-test       compares many `random' GPLs with weight up to"
    print*,"                     four to GiNaC"
#endif
#ifdef DEBUG
    print*,"    -verb <n>        sets the verbosity level to <n>"
    print*,"    -report          performs a coverage test"
#endif

  end subroutine

Luca's avatar
Luca committed
138
  subroutine do_MPL_tests()
Luca's avatar
Luca committed
139
    complex(kind=prec) :: ref
Luca Naterop's avatar
Luca Naterop committed
140
    print*, 'doing MPL tests...'
Luca's avatar
Luca committed
141
    
ulrich_y's avatar
ulrich_y committed
142 143 144 145 146
    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
147
    
ulrich_y's avatar
ulrich_y committed
148 149
    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
150
    
151 152
    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
153

154
    ref = (-0.0319989639656484125354488353677231989818697987217232462808_prec,0.)
ulrich_y's avatar
ulrich_y committed
155
    call test_one_MPL((/2, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.5')
Luca's avatar
Luca committed
156
  end subroutine do_MPL_tests
Luca's avatar
Luca committed
157

Luca Naterop's avatar
Luca Naterop committed
158
  subroutine do_GPL_tests()
ulrich_y's avatar
ulrich_y committed
159
    complex(kind=prec) :: ref, res
Luca Naterop's avatar
Luca Naterop committed
160
    real(kind=prec) :: z, xchen
Luca Naterop's avatar
Luca Naterop committed
161
    print*, 'doing GPL tests...'
Luca's avatar
Luca committed
162
    
163 164
    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
165
    
166 167 168 169 170
    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
171

Luca Naterop's avatar
Luca Naterop committed
172
    ! requires making convergent
173 174
    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
175

176 177
    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
178 179

    ! requires hoelder convolution
180 181
    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
182 183

    ! here the tests from the mathematica nb start
184 185
    ! --------------------------------------------------------------------------

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

188 189
    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
190

191 192
    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
193

194 195
    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
196

197 198
    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
199
    
200 201
    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
202
    
203 204
    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
205
    
206 207
    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')
208 209 210 211 212


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

213 214 215
    z = 1/100._prec
    ref = (-1.203972804325935992622746217761838502954_prec,0.0)
    call test_one_flat((/ (0._prec,0.), (0.3_prec,0.) /),ref,'4.1')
216

217 218
    ref = (10.6037962209567960211233327771880353832_prec,0.0)
    call test_one_flat( (/ (0._prec,0.), (0._prec,0.), (0.01_prec,0.) /),ref,'4.2')
219

220 221
    ref = (5.0429900651807654463744449762857296353E-4_prec,0.0)
    call test_one_flat((/(100._prec,0.), (1._prec,0.), (0.3_prec,0.) /),ref,'4.3')
222

223 224
    ref = (0.05630861877185110944118569266547272128_prec,0.0)
    call test_one_flat( (/ (1._prec,0.), (0._prec,0.), (0.01_prec,0.) /) ,ref,'4.4')
225

226 227
    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')
228

229 230
    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')
231

232 233
    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')
234

235 236
    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')
237

238 239
    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')
240

241 242
    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')
243

244 245
    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')
246

ulrich_y's avatar
ulrich_y committed
247

ulrich_y's avatar
ulrich_y committed
248
    ! Increase test coverage
249
    ref = (-1.203972804325935992622746217761838502953611_prec, pi)
ulrich_y's avatar
ulrich_y committed
250
    call test_one_flat((/ (0._prec,0.), (-0.3_prec,0.) /),ref,'E.1')
ulrich_y's avatar
ulrich_y committed
251

252
    ref = (-0.539404830685964989184608085018712655769250_prec, 1.0303768265243124637877433270311515319634364_prec)
ulrich_y's avatar
ulrich_y committed
253
    call test_one_flat([(0._prec,0._prec), (0.3_prec, 0.5_prec)],ref,'E.2')
254 255

    ref = (3.79489584700899518570534417045601144619E-5_prec,0.)
ulrich_y's avatar
ulrich_y committed
256
    call iprint('  testing GPL E.3 ...',-1)
257
    res = G((/100._prec, 100._prec, 1._prec, 0._prec, 1._prec/))
ulrich_y's avatar
ulrich_y committed
258 259
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
260
    call iprint('  testing GPL E.4 ...',-1)
ulrich_y's avatar
ulrich_y committed
261 262 263
    res = G((/100,100,1,0,1/))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
264
    call iprint('  testing GPL E.5 ...',-1)
ulrich_y's avatar
ulrich_y committed
265 266 267
    res = G_superflatn(toinum((/100,100,1,0,1/)), 5)
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
268
    call iprint('  testing GPL E.6 ...',-1)
269
    res = G((/100._prec,100._prec,1._prec,0._prec/), 1._prec)
ulrich_y's avatar
ulrich_y committed
270 271
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
272
    call iprint('  testing GPL E.7 ...',-1)
273
    res = G((/(100._prec,0.),(100._prec,0.),(1._prec,0.),(0._prec,0.)/), (1._prec,0.))
ulrich_y's avatar
ulrich_y committed
274 275
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
276
    call iprint('  testing GPL E.8 ...',-1)
277
    res = G( (/1,1,1,1/) , (/ 100._prec, 100._prec, 1._prec, 0._prec /), 1._prec)
ulrich_y's avatar
ulrich_y committed
278 279 280
    call check(res,ref)

    !ref = cmplx(0.01592795952537145)
281
    ref = (0.0485035531168477720320998551314479928648_prec,-2.82326885118931135859594797186107474877_prec)
ulrich_y's avatar
ulrich_y committed
282
    call iprint('  testing GPL E.9 ...',-1)
283
    res = G( (/ 3, 2 /), toinum((/ 1.3_prec, 1.1_prec /)), toinum(4._prec) )
ulrich_y's avatar
ulrich_y committed
284 285
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
286
    call iprint('  testing GPL E.10 ...',-1)
287
    res = G( (/ 3, 2 /), (/ 1.3_prec, 1.1_prec /), 4._prec )
ulrich_y's avatar
ulrich_y committed
288 289
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
290
    call iprint('  testing GPL E.11 ...',-1)
291 292 293
    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
294
    call iprint('  testing GPL E.12 ...',-1)
295 296
    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
297 298
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
299
    call iprint('  testing GPL E.13 ...',-1)
ulrich_y's avatar
ulrich_y committed
300 301 302 303
    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
304
    call iprint('  testing GPL E.14 ...',-1)
ulrich_y's avatar
ulrich_y committed
305 306 307 308
    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
309
    call iprint('  testing GPL E.15 ...',-1)
ulrich_y's avatar
ulrich_y committed
310 311 312 313
    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
314
    call iprint('  testing GPL E.17 ...',-1)
ulrich_y's avatar
ulrich_y committed
315 316 317 318
    ref = (0.190800137777535619036913153766083992418_prec, 0.)
    res = G((/ 6, 1, 1 /))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
319
    call iprint('  testing GPL E.18 ...',-1)
ulrich_y's avatar
ulrich_y committed
320 321 322 323 324
    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
325
    call iprint('  testing GPL E.19a ...',-1)
ulrich_y's avatar
ulrich_y committed
326 327
    res = G((/ inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
328
    call iprint('  testing GPL E.19b ...',-1)
ulrich_y's avatar
ulrich_y committed
329 330
    res = G((/ inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
331
    call iprint('  testing GPL E.19c ...',-1)
ulrich_y's avatar
ulrich_y committed
332 333
    res = G((/ inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
334
    call iprint('  testing GPL E.19d ...',-1)
ulrich_y's avatar
ulrich_y committed
335 336 337
    res = G((/ inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
338
    call iprint('  testing GPL E.19e ...',-1)
ulrich_y's avatar
ulrich_y committed
339 340
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
341
    call iprint('  testing GPL E.19f ...',-1)
ulrich_y's avatar
ulrich_y committed
342 343
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
344
    call iprint('  testing GPL E.19g ...',-1)
ulrich_y's avatar
ulrich_y committed
345 346
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
347
    call iprint('  testing GPL E.19h ...',-1)
ulrich_y's avatar
ulrich_y committed
348 349
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
Luca Naterop's avatar
Luca Naterop committed
350

Luca Naterop's avatar
Luca Naterop committed
351

ulrich_y's avatar
ulrich_y committed
352
    ref = -(2.374395270272480200677499763071638424_prec, - 1.273806204919600530933131685580471698_prec)
ulrich_y's avatar
ulrich_y committed
353
    call iprint('  testing GPL E.20a ...',-1)
ulrich_y's avatar
ulrich_y committed
354 355
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
356
    call iprint('  testing GPL E.20b ...',-1)
ulrich_y's avatar
ulrich_y committed
357 358
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
359
    call iprint('  testing GPL E.20c ...',-1)
ulrich_y's avatar
ulrich_y committed
360 361
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
362
    call iprint('  testing GPL E.20d ...',-1)
ulrich_y's avatar
ulrich_y committed
363 364 365
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
366
    call iprint('  testing GPL E.20e ...',-1)
ulrich_y's avatar
ulrich_y committed
367 368
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
369
    call iprint('  testing GPL E.20f ...',-1)
ulrich_y's avatar
ulrich_y committed
370 371
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
372
    call iprint('  testing GPL E.20g ...',-1)
ulrich_y's avatar
ulrich_y committed
373 374
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
375
    call iprint('  testing GPL E.20h ...',-1)
ulrich_y's avatar
ulrich_y committed
376 377
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
378 379 380 381


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

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

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

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

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

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

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

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

    ref = (-1.11161035333074623447215317094270858897_prec,-0.875225967273157437459269290416981432294_prec)
ulrich_y's avatar
ulrich_y committed
406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431
    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')
432

ulrich_y's avatar
ulrich_y committed
433 434
    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')
435

ulrich_y's avatar
ulrich_y committed
436 437 438 439 440 441

    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)
ulrich_y's avatar
ulrich_y committed
442
    call iprint('  testing GPL F.15 ...',-1)
ulrich_y's avatar
ulrich_y committed
443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459
    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
460 461
  end subroutine do_GPL_tests
  
ulrich_y's avatar
ulrich_y committed
462 463 464



ulrich_y's avatar
ulrich_y committed
465
#ifdef HAVE_GINAC
ulrich_y's avatar
ulrich_y committed
466 467

  function evalt(arr, what)
ulrich_y's avatar
ulrich_y committed
468 469 470
#if KINDREAL==16
    use ieps, only: inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
471
    implicit none
ulrich_y's avatar
ulrich_y committed
472 473
    complex(kind=prec) :: arr(:), evalt
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
474 475 476 477 478
    integer what, i, l

    evalt =0.

    do i=1,size(arr)
479
      if (abs(arr(i)) .gt. 1.e10) then ! isnan?
ulrich_y's avatar
ulrich_y committed
480 481 482 483 484 485 486 487 488 489 490 491
        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
492 493 494
#if KINDREAL==16
      evalt = cmplx(geval(inum2inum(toinum(arr(1:l))),l),kind=prec)
#else
ulrich_y's avatar
ulrich_y committed
495
      evalt = geval(toinum(arr(1:l)),l)
ulrich_y's avatar
ulrich_y committed
496
#endif
ulrich_y's avatar
ulrich_y committed
497
    endif
ulrich_y's avatar
ulrich_y committed
498 499
  end function

ulrich_y's avatar
ulrich_y committed
500 501 502

#ifdef HAVE_MM

ulrich_y's avatar
ulrich_y committed
503
  subroutine perform_ginacv(n, args)
ulrich_y's avatar
ulrich_y committed
504
    use maths_functions, only:clearcache
ulrich_y's avatar
ulrich_y committed
505 506
    complex(kind=prec) :: args(:,:)
    integer i,n
ulrich_y's avatar
ulrich_y committed
507
    call clearcache
ulrich_y's avatar
ulrich_y committed
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524

    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
525 526 527
    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
528

ulrich_y's avatar
ulrich_y committed
529
  end subroutine
ulrich_y's avatar
ulrich_y committed
530

ulrich_y's avatar
ulrich_y committed
531 532 533 534
  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
535 536
    integer(kind=8) cstart, cend, count_rate
    real(kind=prec) :: time(2), ttime(2)
ulrich_y's avatar
ulrich_y committed
537 538 539 540 541
    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
542
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
543 544 545
      do i=1,size(args,3)
        res=evalt(args(j,:,i),0)
      enddo
ulrich_y's avatar
ulrich_y committed
546
      call system_clock(cend, count_rate=count_rate)
547
      time(1) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
548
      if (time(1).lt.zero) print*,j, cend-cstart,count_rate
ulrich_y's avatar
ulrich_y committed
549 550
      ttime(1) = ttime(1) + time(1)

ulrich_y's avatar
ulrich_y committed
551
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
552 553 554
      do i=1,size(args,3)
        res=evalt(args(j,:,i),1)
      enddo
ulrich_y's avatar
ulrich_y committed
555
      call system_clock(cend, count_rate=count_rate)
556
      time(2) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
557 558 559 560 561 562 563 564
      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*,
565
    print*,ttime
ulrich_y's avatar
ulrich_y committed
566 567
    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
568 569 570 571 572 573
    ! 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)
574
    ttime(1) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
575 576 577 578 579 580

    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)
581
    ttime(2) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
582 583 584

    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
585
900 FORMAT(a , 'Function ',i4,'/',i4,' for ',a)
586
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
587 588 589 590
  end subroutine

  subroutine do_timing_tests(n)
    use gtestchen   , only: inichen   =>args
ulrich_y's avatar
ulrich_y committed
591
    use gtestchenff , only: inichenff =>args
ulrich_y's avatar
ulrich_y committed
592 593 594 595 596 597
    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
598
    complex(kind=prec) :: fargs(  540,5,n)
ulrich_y's avatar
ulrich_y committed
599 600 601 602 603 604 605 606 607
    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
608 609
      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
610 611 612

      w = ran2(ranseed) ! 0<w<1
      z = ran2(ranseed) * (sqrt(1-w+w**2)-sqrt(w)) + sqrt(w)
613
      nargs(:,:,i) = inimuonenp(cmplx(w,kind=prec), cmplx(z,kind=prec))
ulrich_y's avatar
ulrich_y committed
614 615 616

      x = ran2(ranseed)
      y = ran2(ranseed)
617
      pargs(:,:,i) = inimuone(cmplx(x,kind=prec), cmplx(y,kind=prec))
ulrich_y's avatar
ulrich_y committed
618 619
    enddo

ulrich_y's avatar
ulrich_y committed
620 621
    cargs(1181,:,:)=1.e15
    fargs(367,:,:)=1.e15
ulrich_y's avatar
ulrich_y committed
622 623 624


    open(unit=9, file="stats.txt")
ulrich_y's avatar
ulrich_y committed
625 626
    write(9,*) "Chen form factor"
    call do_one_speed_test(fargs,9,"Chen FF")
ulrich_y's avatar
ulrich_y committed
627 628 629 630 631 632 633 634 635
    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
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 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728

  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
729
#endif
ulrich_y's avatar
ulrich_y committed
730

ulrich_y's avatar
ulrich_y committed
731
  SUBROUTINE DO_LONG_TEST
ulrich_y's avatar
ulrich_y committed
732 733 734
#if KINDREAL==16
    use ieps,only:inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
735 736 737 738 739
    implicit none
    integer,parameter :: nzero = 10
    integer,parameter :: nieps = 30
    integer,parameter :: ncmpl = 30

740
    integer,parameter :: perweight(4) = (/ 100, 500, 30000, 50000 /)
ulrich_y's avatar
ulrich_y committed
741 742 743 744
    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
745
    integer i, j, w, seed, oldseed
ulrich_y's avatar
ulrich_y committed
746
    integer(kind=1) i0
ulrich_y's avatar
ulrich_y committed
747 748 749
    real(kind=prec) :: v, maxd
    complex(kind=prec) :: ans(2)
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
750
    character,parameter :: cr = achar(13)
ulrich_y's avatar
ulrich_y committed
751
    maxd=0._prec
ulrich_y's avatar
ulrich_y committed
752 753

    seed = 112312
ulrich_y's avatar
ulrich_y committed
754
    tol = 0.01
ulrich_y's avatar
ulrich_y committed
755

ulrich_y's avatar
ulrich_y committed
756
    open(unit=9, action='write', form='unformatted', file="long-test.txt")
ulrich_y's avatar
ulrich_y committed
757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774
    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
775 776
        oldseed = seed
        write( * , 900, advance='no' ) cr, i, perweight(w), w
ulrich_y's avatar
ulrich_y committed
777
        args(1:w) = basis((/ (1+int(size(basis)*ran2(seed)),j=1,w) /))
ulrich_y's avatar
ulrich_y committed
778 779 780 781 782
#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
783
        ans(2) = G(args(1:w),ione)
784 785 786 787 788 789
        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
790
        write(9) w,i,oldseed, abs(ans(1)-ans(2))
791
        flush(9)
ulrich_y's avatar
ulrich_y committed
792
        if(abs(ans(1)-ans(2)) > maxd) maxd = abs(ans(1)-ans(2))
ulrich_y's avatar
ulrich_y committed
793 794
        if(abs(ans(1)-ans(2)) > tol) goto 123
      enddo
ulrich_y's avatar
ulrich_y committed
795 796
      print*,' done. largest=',maxd
      maxd=0.
ulrich_y's avatar
ulrich_y committed
797
    enddo
ulrich_y's avatar
ulrich_y committed
798
    print*,maxd
ulrich_y's avatar
ulrich_y committed
799

ulrich_y's avatar
ulrich_y committed
800 801
    close(9)

ulrich_y's avatar
ulrich_y committed
802 803 804 805 806
    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
807 808 809
    close(9)

900 FORMAT(a , 'Testing ',i5,'/',i5,' GPLs with w=',i1)
ulrich_y's avatar
ulrich_y committed
810 811 812


  END SUBROUTINE
ulrich_y's avatar
ulrich_y committed
813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830
      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
831

832
#endif
ulrich_y's avatar
ulrich_y committed
833

834 835 836
  ! subroutine do_shuffle_tests() 
  !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
  !   complex(kind=prec) :: w(2) = cmplx((/3,4/))
837

838 839
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
840

Luca's avatar
Luca committed
841
END PROGRAM TEST
Luca's avatar
Luca committed
842 843

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