test.f90 32 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. 
ulrich_y's avatar
ulrich_y committed
16 17
  integer :: i
  character(len=32) :: arg
18

ulrich_y's avatar
ulrich_y committed
19

20
  i = 1
ulrich_y's avatar
ulrich_y committed
21 22 23 24 25 26 27
  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
28
#ifdef DEBUG
ulrich_y's avatar
ulrich_y committed
29
        verb = readint(trim(arg),i)
ulrich_y's avatar
ulrich_y committed
30
#else
ulrich_y's avatar
ulrich_y committed
31
        call errprint("Argument -verb is not available, compile with --debug")
ulrich_y's avatar
ulrich_y committed
32
#endif
33 34 35 36 37 38 39 40 41 42
      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
43
#if defined(HAVE_GINAC) && defined(HAVE_MM)
44 45 46
      case('-ginac-tests')
        tol = 8.0e-7
        call do_ginac_tests
ulrich_y's avatar
ulrich_y committed
47
      case('-speed-tests')
ulrich_y's avatar
ulrich_y committed
48
        call do_timing_tests(readint(trim(arg),i))
ulrich_y's avatar
ulrich_y committed
49
      case('-hw-tests')
50 51
        tol = 8.0e-7
        call do_high_weight_tests
ulrich_y's avatar
ulrich_y committed
52 53 54 55 56 57
#else
      case('-ginac-tests', '-speed-tests', '-hw-tests')
        call errprint("Argument "//trim(arg)//" is not available, compile with --with-ginac --with-mcc")
#endif

#ifdef HAVE_GINAC
ulrich_y's avatar
ulrich_y committed
58
      case('-long-test')
59 60
        tol = 8.0e-7
        call do_long_test
ulrich_y's avatar
ulrich_y committed
61 62 63
#else
      case('-long-test')
        call errprint("Argument "//trim(arg)//" is not available, compile with --with-ginac")
64 65
#endif
#ifdef DEBUG
ulrich_y's avatar
ulrich_y committed
66
      case('-report')
67 68 69 70 71 72
        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
73
#if defined(HAVE_GINAC) && defined(HAVE_MM)
74
        call do_ginac_tests
ulrich_y's avatar
ulrich_y committed
75
#endif
76 77
#endif

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

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

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


  if(tests_successful) then
    print*, 'All tests passed. '
  else 
ulrich_y's avatar
ulrich_y committed
93
    call errprint('Some tests failed. ')
Luca's avatar
tests  
Luca committed
94
  end if
ulrich_y's avatar
ulrich_y committed
95
  stop
96

97
CONTAINS
98
   
ulrich_y's avatar
ulrich_y committed
99 100 101 102 103 104
  subroutine errprint(msg)
    character(len=*) msg
    write(0,*) "Error:", msg
    stop 1
  end subroutine

ulrich_y's avatar
ulrich_y committed
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 138 139 140 141
  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

ulrich_y's avatar
ulrich_y committed
142 143 144 145 146 147 148 149 150 151 152
  function readint(prev,i)
    integer i, readint, st
    character(len=32) :: arg
    character(len=*) :: prev
    i=i+1
    call get_command_argument(i,arg)
    if (len_trim(arg) == 0) call errprint("Argument "//prev//" requires a number")
    read(arg,*,iostat=st) readint
    if (st .ne. 0) call errprint("For argument "//prev//": "//trim(arg)//" is not a number")
  end function

Luca Naterop's avatar
Luca Naterop committed
153 154 155 156
  subroutine check(res, ref)
    complex(kind=prec) :: res, ref
    real(kind=prec) :: delta

Luca Naterop's avatar
Luca Naterop committed
157
    delta = abs(res-ref)
Luca Naterop's avatar
Luca Naterop committed
158
    if(delta < tol) then
Luca Naterop's avatar
Luca Naterop committed
159
      print*, '       ',' passed with delta = ', delta
Luca Naterop's avatar
Luca Naterop committed
160
    else 
161
      print*, '  ',' FAILED with delta = ', delta
162
      tests_successful = .false.
Luca Naterop's avatar
Luca Naterop committed
163 164
    end if
  end subroutine check
165 166 167

  subroutine test_one_MPL(m,x,ref, test_id)
    integer :: m(:)
168
    complex(kind=prec) :: x(:), ref, res
169
    character(len=*) :: test_id
Luca's avatar
Luca committed
170
    
171 172 173 174 175
    print*, '  ', 'testing MPL ', test_id, ' ...'
    res = MPL(m,x)
    call check(res,ref)
  end subroutine test_one_MPL

Luca's avatar
Luca committed
176
  subroutine do_MPL_tests()
Luca's avatar
Luca committed
177
    complex(kind=prec) :: ref
Luca Naterop's avatar
Luca Naterop committed
178
    print*, 'doing MPL tests...'
Luca's avatar
Luca committed
179
    
ulrich_y's avatar
ulrich_y committed
180 181 182 183 184
    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
185
    
ulrich_y's avatar
ulrich_y committed
186 187
    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
188
    
189 190
    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
191

192
    ref = (-0.0319989639656484125354488353677231989818697987217232462808_prec,0.)
ulrich_y's avatar
ulrich_y committed
193
    call test_one_MPL((/2, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.5')
Luca's avatar
Luca committed
194
  end subroutine do_MPL_tests
Luca's avatar
Luca committed
195

Luca's avatar
Luca committed
196
  subroutine test_one_condensed(m,z,y,k,ref,test_id)
197 198 199
    integer :: m(:), k
    complex(kind=prec) :: z(:), y, res, ref
    character(len=*) :: test_id
Luca Naterop's avatar
Luca Naterop committed
200

201
    print*, '  ', 'testing GPL ', test_id, ' ...'
ulrich_y's avatar
ulrich_y committed
202
    res = G_condensed(m,toinum(z),inum(y,di0),k)
Luca Naterop's avatar
Luca Naterop committed
203
    call check(res,ref)
Luca's avatar
Luca committed
204 205
  end subroutine test_one_condensed

206 207
  subroutine test_one_flat(z,ref,test_id)
    complex(kind=prec) :: z(:), res, ref
Luca's avatar
Luca committed
208 209 210
    character(len=*) :: test_id

    print*, '  ', 'testing GPL ', test_id, ' ...'
ulrich_y's avatar
ulrich_y committed
211
    res = G(z)
Luca's avatar
Luca committed
212 213
    call check(res,ref)
  end subroutine test_one_flat
Luca Naterop's avatar
Luca Naterop committed
214 215

  subroutine do_GPL_tests()
ulrich_y's avatar
ulrich_y committed
216
    complex(kind=prec) :: ref, res
Luca Naterop's avatar
Luca Naterop committed
217
    real(kind=prec) :: z, xchen
Luca Naterop's avatar
Luca Naterop committed
218
    print*, 'doing GPL tests...'
Luca's avatar
Luca committed
219
    
220 221
    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
222
    
223 224 225 226 227
    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
228

Luca Naterop's avatar
Luca Naterop committed
229
    ! requires making convergent
230 231
    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
232

233 234
    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
235 236

    ! requires hoelder convolution
237 238
    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
239 240

    ! here the tests from the mathematica nb start
241 242
    ! --------------------------------------------------------------------------

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

245 246
    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
247

248 249
    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
250

251 252
    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
253

254 255
    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
256
    
257 258
    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
259
    
260 261
    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
262
    
263 264
    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')
265 266 267 268 269


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

270 271 272
    z = 1/100._prec
    ref = (-1.203972804325935992622746217761838502954_prec,0.0)
    call test_one_flat((/ (0._prec,0.), (0.3_prec,0.) /),ref,'4.1')
273

274 275
    ref = (10.6037962209567960211233327771880353832_prec,0.0)
    call test_one_flat( (/ (0._prec,0.), (0._prec,0.), (0.01_prec,0.) /),ref,'4.2')
276

277 278
    ref = (5.0429900651807654463744449762857296353E-4_prec,0.0)
    call test_one_flat((/(100._prec,0.), (1._prec,0.), (0.3_prec,0.) /),ref,'4.3')
279

280 281
    ref = (0.05630861877185110944118569266547272128_prec,0.0)
    call test_one_flat( (/ (1._prec,0.), (0._prec,0.), (0.01_prec,0.) /) ,ref,'4.4')
282

283 284
    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')
285

286 287
    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')
288

289 290
    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')
291

292 293
    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')
294

295 296
    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')
297

298 299
    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')
300

301 302
    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')
303

ulrich_y's avatar
ulrich_y committed
304

ulrich_y's avatar
ulrich_y committed
305
    ! Increase test coverage
306
    ref = (-1.203972804325935992622746217761838502953611_prec, pi)
ulrich_y's avatar
ulrich_y committed
307
    call test_one_flat((/ (0._prec,0.), (-0.3_prec,0.) /),ref,'E.1')
ulrich_y's avatar
ulrich_y committed
308

309
    ref = (-0.539404830685964989184608085018712655769250_prec, 1.0303768265243124637877433270311515319634364_prec)
ulrich_y's avatar
ulrich_y committed
310
    call test_one_flat([(0._prec,0._prec), (0.3_prec, 0.5_prec)],ref,'E.2')
311 312

    ref = (3.79489584700899518570534417045601144619E-5_prec,0.)
ulrich_y's avatar
ulrich_y committed
313
    print*, '  ', 'testing GPL ', 'E.3', ' ...'
314
    res = G((/100._prec, 100._prec, 1._prec, 0._prec, 1._prec/))
ulrich_y's avatar
ulrich_y committed
315 316
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
317
    print*, '  ', 'testing GPL ', 'E.4', ' ...'
ulrich_y's avatar
ulrich_y committed
318 319 320
    res = G((/100,100,1,0,1/))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
321
    print*, '  ', 'testing GPL ', 'E.5', ' ...'
ulrich_y's avatar
ulrich_y committed
322 323 324
    res = G_superflatn(toinum((/100,100,1,0,1/)), 5)
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
325
    print*, '  ', 'testing GPL ', 'E.6', ' ...'
326
    res = G((/100._prec,100._prec,1._prec,0._prec/), 1._prec)
ulrich_y's avatar
ulrich_y committed
327 328
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
329
    print*, '  ', 'testing GPL ', 'E.7', ' ...'
330
    res = G((/(100._prec,0.),(100._prec,0.),(1._prec,0.),(0._prec,0.)/), (1._prec,0.))
ulrich_y's avatar
ulrich_y committed
331 332
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
333
    print*, '  ', 'testing GPL ', 'E.8', ' ...'
334
    res = G( (/1,1,1,1/) , (/ 100._prec, 100._prec, 1._prec, 0._prec /), 1._prec)
ulrich_y's avatar
ulrich_y committed
335 336 337
    call check(res,ref)

    !ref = cmplx(0.01592795952537145)
338
    ref = (0.0485035531168477720320998551314479928648_prec,-2.82326885118931135859594797186107474877_prec)
ulrich_y's avatar
ulrich_y committed
339
    print*, '  ', 'testing GPL ', 'E.9', ' ...'
340
    res = G( (/ 3, 2 /), toinum((/ 1.3_prec, 1.1_prec /)), toinum(4._prec) )
ulrich_y's avatar
ulrich_y committed
341 342
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
343
    print*, '  ', 'testing GPL ', 'E.10', ' ...'
344
    res = G( (/ 3, 2 /), (/ 1.3_prec, 1.1_prec /), 4._prec )
ulrich_y's avatar
ulrich_y committed
345 346
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
347
    print*, '  ', 'testing GPL ', 'E.11', ' ...'
348 349 350
    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
351
    print*, '  ', 'testing GPL ', 'E.12', ' ...'
352 353
    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
354 355
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
356
    print*, '  ', 'testing GPL ', 'E.13', ' ...'
ulrich_y's avatar
ulrich_y committed
357 358 359 360
    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
361
    print*, '  ', 'testing GPL ', 'E.14', ' ...'
ulrich_y's avatar
ulrich_y committed
362 363 364 365
    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
366
    print*, '  ', 'testing GPL ', 'E.15', ' ...'
ulrich_y's avatar
ulrich_y committed
367 368 369 370
    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
371
    print*, '  ', 'testing GPL ', 'E.17', ' ...'
ulrich_y's avatar
ulrich_y committed
372 373 374 375
    ref = (0.190800137777535619036913153766083992418_prec, 0.)
    res = G((/ 6, 1, 1 /))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
376
    print*, '  ', 'testing GPL ', 'E.18', ' ...'
ulrich_y's avatar
ulrich_y committed
377 378 379 380 381
    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
382
    print*, '  ', 'testing GPL ', 'E.19a', ' ...'
ulrich_y's avatar
ulrich_y committed
383 384
    res = G((/ inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
385
    print*, '  ', 'testing GPL ', 'E.19b', ' ...'
ulrich_y's avatar
ulrich_y committed
386 387
    res = G((/ inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
388
    print*, '  ', 'testing GPL ', 'E.19c', ' ...'
ulrich_y's avatar
ulrich_y committed
389 390
    res = G((/ inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
391
    print*, '  ', 'testing GPL ', 'E.19d', ' ...'
ulrich_y's avatar
ulrich_y committed
392 393 394
    res = G((/ inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
395
    print*, '  ', 'testing GPL ', 'E.19e', ' ...'
ulrich_y's avatar
ulrich_y committed
396 397
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
398
    print*, '  ', 'testing GPL ', 'E.19f', ' ...'
ulrich_y's avatar
ulrich_y committed
399 400
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
401
    print*, '  ', 'testing GPL ', 'E.19g', ' ...'
ulrich_y's avatar
ulrich_y committed
402 403
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
404
    print*, '  ', 'testing GPL ', 'E.19h', ' ...'
ulrich_y's avatar
ulrich_y committed
405 406
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
Luca Naterop's avatar
Luca Naterop committed
407

Luca Naterop's avatar
Luca Naterop committed
408

ulrich_y's avatar
ulrich_y committed
409
    ref = -(2.374395270272480200677499763071638424_prec, - 1.273806204919600530933131685580471698_prec)
ulrich_y's avatar
ulrich_y committed
410
    print*, '  ', 'testing GPL ', 'E.20a', ' ...'
ulrich_y's avatar
ulrich_y committed
411 412
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
413
    print*, '  ', 'testing GPL ', 'E.20b', ' ...'
ulrich_y's avatar
ulrich_y committed
414 415
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
416
    print*, '  ', 'testing GPL ', 'E.20c', ' ...'
ulrich_y's avatar
ulrich_y committed
417 418
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
419
    print*, '  ', 'testing GPL ', 'E.20d', ' ...'
ulrich_y's avatar
ulrich_y committed
420 421 422
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
423
    print*, '  ', 'testing GPL ', 'E.20e', ' ...'
ulrich_y's avatar
ulrich_y committed
424 425
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
426
    print*, '  ', 'testing GPL ', 'E.20f', ' ...'
ulrich_y's avatar
ulrich_y committed
427 428
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
429
    print*, '  ', 'testing GPL ', 'E.20g', ' ...'
ulrich_y's avatar
ulrich_y committed
430 431
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
432
    print*, '  ', 'testing GPL ', 'E.20h', ' ...'
ulrich_y's avatar
ulrich_y committed
433 434
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
435 436 437 438


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

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

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

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

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

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

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

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

    ref = (-1.11161035333074623447215317094270858897_prec,-0.875225967273157437459269290416981432294_prec)
ulrich_y's avatar
ulrich_y committed
463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488
    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')
489

ulrich_y's avatar
ulrich_y committed
490 491
    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')
492

ulrich_y's avatar
ulrich_y committed
493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516

    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
517 518
  end subroutine do_GPL_tests
  
ulrich_y's avatar
ulrich_y committed
519 520 521



ulrich_y's avatar
ulrich_y committed
522
#ifdef HAVE_GINAC
ulrich_y's avatar
ulrich_y committed
523 524

  function evalt(arr, what)
ulrich_y's avatar
ulrich_y committed
525 526 527
#if KINDREAL==16
    use ieps, only: inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
528
    implicit none
ulrich_y's avatar
ulrich_y committed
529 530
    complex(kind=prec) :: arr(:), evalt
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
531 532 533 534 535
    integer what, i, l

    evalt =0.

    do i=1,size(arr)
536
      if (abs(arr(i)) .gt. 1.e10) then ! isnan?
ulrich_y's avatar
ulrich_y committed
537 538 539 540 541 542 543 544 545 546 547 548
        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
549 550 551
#if KINDREAL==16
      evalt = cmplx(geval(inum2inum(toinum(arr(1:l))),l),kind=prec)
#else
ulrich_y's avatar
ulrich_y committed
552
      evalt = geval(toinum(arr(1:l)),l)
ulrich_y's avatar
ulrich_y committed
553
#endif
ulrich_y's avatar
ulrich_y committed
554
    endif
ulrich_y's avatar
ulrich_y committed
555 556
  end function

ulrich_y's avatar
ulrich_y committed
557 558 559

#ifdef HAVE_MM

ulrich_y's avatar
ulrich_y committed
560
  subroutine perform_ginacv(n, args)
ulrich_y's avatar
ulrich_y committed
561
    use maths_functions, only:clearcache
ulrich_y's avatar
ulrich_y committed
562 563
    complex(kind=prec) :: args(:,:)
    integer i,n
ulrich_y's avatar
ulrich_y committed
564
    call clearcache
ulrich_y's avatar
ulrich_y committed
565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581

    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
582 583 584
    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
585

ulrich_y's avatar
ulrich_y committed
586
  end subroutine
ulrich_y's avatar
ulrich_y committed
587

ulrich_y's avatar
ulrich_y committed
588 589 590 591
  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
592 593
    integer(kind=8) cstart, cend, count_rate
    real(kind=prec) :: time(2), ttime(2)
ulrich_y's avatar
ulrich_y committed
594 595 596 597 598
    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
599
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
600 601 602
      do i=1,size(args,3)
        res=evalt(args(j,:,i),0)
      enddo
ulrich_y's avatar
ulrich_y committed
603
      call system_clock(cend, count_rate=count_rate)
604
      time(1) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
605
      if (time(1).lt.zero) print*,j, cend-cstart,count_rate
ulrich_y's avatar
ulrich_y committed
606 607
      ttime(1) = ttime(1) + time(1)

ulrich_y's avatar
ulrich_y committed
608
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
609 610 611
      do i=1,size(args,3)
        res=evalt(args(j,:,i),1)
      enddo
ulrich_y's avatar
ulrich_y committed
612
      call system_clock(cend, count_rate=count_rate)
613
      time(2) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
614 615 616 617 618 619 620 621
      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*,
622
    print*,ttime
ulrich_y's avatar
ulrich_y committed
623 624
    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
625 626 627 628 629 630
    ! 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)
631
    ttime(1) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
632 633 634 635 636 637

    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)
638
    ttime(2) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
639 640 641

    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
642
900 FORMAT(a , 'Function ',i4,'/',i4,' for ',a)
643
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
644 645 646 647
  end subroutine

  subroutine do_timing_tests(n)
    use gtestchen   , only: inichen   =>args
ulrich_y's avatar
ulrich_y committed
648
    use gtestchenff , only: inichenff =>args
ulrich_y's avatar
ulrich_y committed
649 650 651 652 653 654
    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
655
    complex(kind=prec) :: fargs(  540,5,n)
ulrich_y's avatar
ulrich_y committed
656 657 658 659 660 661 662 663 664
    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
665 666
      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
667 668 669

      w = ran2(ranseed) ! 0<w<1
      z = ran2(ranseed) * (sqrt(1-w+w**2)-sqrt(w)) + sqrt(w)
670
      nargs(:,:,i) = inimuonenp(cmplx(w,kind=prec), cmplx(z,kind=prec))
ulrich_y's avatar
ulrich_y committed
671 672 673

      x = ran2(ranseed)
      y = ran2(ranseed)
674
      pargs(:,:,i) = inimuone(cmplx(x,kind=prec), cmplx(y,kind=prec))
ulrich_y's avatar
ulrich_y committed
675 676
    enddo

ulrich_y's avatar
ulrich_y committed
677 678
    cargs(1181,:,:)=1.e15
    fargs(367,:,:)=1.e15
ulrich_y's avatar
ulrich_y committed
679 680 681


    open(unit=9, file="stats.txt")
ulrich_y's avatar
ulrich_y committed
682 683
    write(9,*) "Chen form factor"
    call do_one_speed_test(fargs,9,"Chen FF")
ulrich_y's avatar
ulrich_y committed
684 685 686 687 688 689 690 691 692
    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
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 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785

  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
786
#endif
ulrich_y's avatar
ulrich_y committed
787

ulrich_y's avatar
ulrich_y committed
788
  SUBROUTINE DO_LONG_TEST
ulrich_y's avatar
ulrich_y committed
789 790 791
#if KINDREAL==16
    use ieps,only:inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
792 793 794 795 796
    implicit none
    integer,parameter :: nzero = 10
    integer,parameter :: nieps = 30
    integer,parameter :: ncmpl = 30

797
    integer,parameter :: perweight(4) = (/ 100, 500, 30000, 50000 /)
ulrich_y's avatar
ulrich_y committed
798 799 800 801
    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
802
    integer i, j, w, seed, oldseed
ulrich_y's avatar
ulrich_y committed
803
    integer(kind=1) i0
ulrich_y's avatar
ulrich_y committed
804 805 806
    real(kind=prec) :: v, maxd
    complex(kind=prec) :: ans(2)
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
807
    character,parameter :: cr = achar(13)
ulrich_y's avatar
ulrich_y committed
808
    maxd=0._prec
ulrich_y's avatar
ulrich_y committed
809 810

    seed = 112312
ulrich_y's avatar
ulrich_y committed
811
    tol = 0.01
ulrich_y's avatar
ulrich_y committed
812

ulrich_y's avatar
ulrich_y committed
813
    open(unit=9, action='write', form='unformatted', file="long-test.txt")
ulrich_y's avatar
ulrich_y committed
814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831
    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
832 833
        oldseed = seed
        write( * , 900, advance='no' ) cr, i, perweight(w), w
ulrich_y's avatar
ulrich_y committed
834
        args(1:w) = basis((/ (1+int(size(basis)*ran2(seed)),j=1,w) /))
ulrich_y's avatar
ulrich_y committed
835 836 837 838 839
#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
840
        ans(2) = G(args(1:w),ione)
841 842 843 844 845 846
        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
847
        write(9) w,i,oldseed, abs(ans(1)-ans(2))
848
        flush(9)
ulrich_y's avatar
ulrich_y committed
849
        if(abs(ans(1)-ans(2)) > maxd) maxd = abs(ans(1)-ans(2))
ulrich_y's avatar
ulrich_y committed
850 851
        if(abs(ans(1)-ans(2)) > tol) goto 123
      enddo
ulrich_y's avatar
ulrich_y committed
852 853
      print*,' done. largest=',maxd
      maxd=0.
ulrich_y's avatar
ulrich_y committed
854
    enddo
ulrich_y's avatar
ulrich_y committed
855
    print*,maxd
ulrich_y's avatar
ulrich_y committed
856

ulrich_y's avatar
ulrich_y committed
857 858
    close(9)

ulrich_y's avatar
ulrich_y committed
859 860 861 862 863
    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
864 865 866
    close(9)

900 FORMAT(a , 'Testing ',i5,'/',i5,' GPLs with w=',i1)
ulrich_y's avatar
ulrich_y committed
867 868 869


  END SUBROUTINE
ulrich_y's avatar
ulrich_y committed
870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887
      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
888

889
#endif
ulrich_y's avatar
ulrich_y committed
890

891 892 893
  ! subroutine do_shuffle_tests() 
  !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
  !   complex(kind=prec) :: w(2) = cmplx((/3,4/))
894

895 896
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
897

Luca's avatar
Luca committed
898
END PROGRAM TEST
Luca's avatar
Luca committed
899 900

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