test.f90 30.4 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
    character(len=40) :: msg
ulrich_y's avatar
ulrich_y committed
508
    call clearcache
ulrich_y's avatar
ulrich_y committed
509 510

    do i=1,size(args,1)
ulrich_y's avatar
ulrich_y committed
511 512
      write(msg,900) n,i
      call iprint(msg, -1)
ulrich_y's avatar
ulrich_y committed
513 514 515 516 517 518 519 520 521 522 523 524 525 526
      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
527 528 529
    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
530

ulrich_y's avatar
ulrich_y committed
531
  end subroutine
ulrich_y's avatar
ulrich_y committed
532

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

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

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

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

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

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

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

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


    open(unit=9, file="stats.txt")
ulrich_y's avatar
ulrich_y committed
627 628
    write(9,*) "Chen form factor"
    call do_one_speed_test(fargs,9,"Chen FF")
ulrich_y's avatar
ulrich_y committed
629 630 631 632 633 634 635 636 637
    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
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 729 730

  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
731
#endif
ulrich_y's avatar
ulrich_y committed
732

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

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

    seed = 112312
ulrich_y's avatar
ulrich_y committed
756
    tol = 0.01
ulrich_y's avatar
ulrich_y committed
757

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

ulrich_y's avatar
ulrich_y committed
802 803
    close(9)

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

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


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

834
#endif
ulrich_y's avatar
ulrich_y committed
835

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

840 841
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
842

Luca's avatar
Luca committed
843
END PROGRAM TEST
Luca's avatar
Luca committed
844 845

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