test.f90 33 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 iprint("Argument -verb is not available, compile with --debug", 2)
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
#else
      case('-ginac-tests', '-speed-tests', '-hw-tests')
ulrich_y's avatar
ulrich_y committed
54
        call iprint("Argument "//trim(arg)//" is not available, compile with --with-ginac --with-mcc", 2)
ulrich_y's avatar
ulrich_y committed
55 56 57
#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
#else
      case('-long-test')
ulrich_y's avatar
ulrich_y committed
63
        call iprint("Argument "//trim(arg)//" is not available, compile with --with-ginac", 2)
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
ulrich_y's avatar
ulrich_y committed
76 77 78
#else
      case('-report')
        call iprint("Argument -report is not available, compile with --debug", 2)
79 80
#endif

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

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

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


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

100
CONTAINS
101
   
ulrich_y's avatar
ulrich_y committed
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
  subroutine iprint(imsg, typ)
    character(len=*) imsg
    character(len=100) msg
    integer :: typ
    character(len=5),parameter :: red   = char(27)//'[31m'
    character(len=5),parameter :: green = char(27)//'[32m'
    character(len=5),parameter :: orange= char(27)//'[33m'
    character(len=4),parameter :: norm  = char(27)//'[0m'
    character       ,parameter :: cr    = achar(13)
    integer, save :: prevlen, prevtype
    character(len=100), save :: prevmsg

    if (prevtype == -1) then
      msg = prevmsg(1:prevlen) // imsg
    else
      msg = imsg
    endif

    select case(typ)
      case(0)
        print*,green //'[PASS]'//norm//' '//trim(msg)
      case(1)
        print*, red  //'[FAIL]'//norm//' '//trim(msg)
      case(2)
        print*, red  //'[FATL]'//norm//' '//trim(msg)
        stop 1
      case(3)
        print*,orange//'[WARN]'//norm//' '//trim(msg)
      case(-1)
        write(*,'(a)',advance='no')' [    ]'//' '//trim(msg)//cr
    end select
    prevtype = typ
    prevlen  = len_trim(msg)
    prevmsg  = msg
ulrich_y's avatar
ulrich_y committed
136 137
  end subroutine

ulrich_y's avatar
ulrich_y committed
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
  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
175 176 177 178 179 180
  function readint(prev,i)
    integer i, readint, st
    character(len=32) :: arg
    character(len=*) :: prev
    i=i+1
    call get_command_argument(i,arg)
ulrich_y's avatar
ulrich_y committed
181
    if (len_trim(arg) == 0) call iprint("Argument "//prev//" requires a number",2)
ulrich_y's avatar
ulrich_y committed
182
    read(arg,*,iostat=st) readint
ulrich_y's avatar
ulrich_y committed
183
    if (st .ne. 0) call iprint("For argument "//prev//": "//trim(arg)//" is not a number",2)
ulrich_y's avatar
ulrich_y committed
184 185
  end function

Luca Naterop's avatar
Luca Naterop committed
186 187 188
  subroutine check(res, ref)
    complex(kind=prec) :: res, ref
    real(kind=prec) :: delta
ulrich_y's avatar
ulrich_y committed
189
    character(len=40) :: msg
Luca Naterop's avatar
Luca Naterop committed
190

Luca Naterop's avatar
Luca Naterop committed
191
    delta = abs(res-ref)
Luca Naterop's avatar
Luca Naterop committed
192
    if(delta < tol) then
ulrich_y's avatar
ulrich_y committed
193 194
      write(msg, 900) delta
      call iprint(trim(msg), 0)
Luca Naterop's avatar
Luca Naterop committed
195
    else 
ulrich_y's avatar
ulrich_y committed
196 197
      write(msg, 900) delta
      call iprint(trim(msg), 1)
198
      tests_successful = .false.
Luca Naterop's avatar
Luca Naterop committed
199
    end if
ulrich_y's avatar
ulrich_y committed
200 201

900 format(" with delta = ",ES10.3)
Luca Naterop's avatar
Luca Naterop committed
202
  end subroutine check
203 204 205

  subroutine test_one_MPL(m,x,ref, test_id)
    integer :: m(:)
206
    complex(kind=prec) :: x(:), ref, res
207
    character(len=*) :: test_id
Luca's avatar
Luca committed
208
    
ulrich_y's avatar
ulrich_y committed
209
    call iprint('  testing MPL '//test_id//' ...',-1)
210 211 212 213
    res = MPL(m,x)
    call check(res,ref)
  end subroutine test_one_MPL

Luca's avatar
Luca committed
214
  subroutine do_MPL_tests()
Luca's avatar
Luca committed
215
    complex(kind=prec) :: ref
Luca Naterop's avatar
Luca Naterop committed
216
    print*, 'doing MPL tests...'
Luca's avatar
Luca committed
217
    
ulrich_y's avatar
ulrich_y committed
218 219 220 221 222
    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
223
    
ulrich_y's avatar
ulrich_y committed
224 225
    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
226
    
227 228
    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
229

230
    ref = (-0.0319989639656484125354488353677231989818697987217232462808_prec,0.)
ulrich_y's avatar
ulrich_y committed
231
    call test_one_MPL((/2, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.5')
Luca's avatar
Luca committed
232
  end subroutine do_MPL_tests
Luca's avatar
Luca committed
233

Luca's avatar
Luca committed
234
  subroutine test_one_condensed(m,z,y,k,ref,test_id)
235 236 237
    integer :: m(:), k
    complex(kind=prec) :: z(:), y, res, ref
    character(len=*) :: test_id
Luca Naterop's avatar
Luca Naterop committed
238

ulrich_y's avatar
ulrich_y committed
239
    call iprint('  testing GPL '//test_id//' ...',-1)
ulrich_y's avatar
ulrich_y committed
240
    res = G_condensed(m,toinum(z),inum(y,di0),k)
Luca Naterop's avatar
Luca Naterop committed
241
    call check(res,ref)
Luca's avatar
Luca committed
242 243
  end subroutine test_one_condensed

244 245
  subroutine test_one_flat(z,ref,test_id)
    complex(kind=prec) :: z(:), res, ref
Luca's avatar
Luca committed
246 247
    character(len=*) :: test_id

ulrich_y's avatar
ulrich_y committed
248
    call iprint('  testing GPL '//test_id//' ...',-1)
ulrich_y's avatar
ulrich_y committed
249
    res = G(z)
Luca's avatar
Luca committed
250 251
    call check(res,ref)
  end subroutine test_one_flat
Luca Naterop's avatar
Luca Naterop committed
252 253

  subroutine do_GPL_tests()
ulrich_y's avatar
ulrich_y committed
254
    complex(kind=prec) :: ref, res
Luca Naterop's avatar
Luca Naterop committed
255
    real(kind=prec) :: z, xchen
Luca Naterop's avatar
Luca Naterop committed
256
    print*, 'doing GPL tests...'
Luca's avatar
Luca committed
257
    
258 259
    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
260
    
261 262 263 264 265
    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
266

Luca Naterop's avatar
Luca Naterop committed
267
    ! requires making convergent
268 269
    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
270

271 272
    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
273 274

    ! requires hoelder convolution
275 276
    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
277 278

    ! here the tests from the mathematica nb start
279 280
    ! --------------------------------------------------------------------------

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

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

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

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

292 293
    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
294
    
295 296
    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
297
    
298 299
    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
300
    
301 302
    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')
303 304 305 306 307


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

308 309 310
    z = 1/100._prec
    ref = (-1.203972804325935992622746217761838502954_prec,0.0)
    call test_one_flat((/ (0._prec,0.), (0.3_prec,0.) /),ref,'4.1')
311

312 313
    ref = (10.6037962209567960211233327771880353832_prec,0.0)
    call test_one_flat( (/ (0._prec,0.), (0._prec,0.), (0.01_prec,0.) /),ref,'4.2')
314

315 316
    ref = (5.0429900651807654463744449762857296353E-4_prec,0.0)
    call test_one_flat((/(100._prec,0.), (1._prec,0.), (0.3_prec,0.) /),ref,'4.3')
317

318 319
    ref = (0.05630861877185110944118569266547272128_prec,0.0)
    call test_one_flat( (/ (1._prec,0.), (0._prec,0.), (0.01_prec,0.) /) ,ref,'4.4')
320

321 322
    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')
323

324 325
    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')
326

327 328
    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')
329

330 331
    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')
332

333 334
    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')
335

336 337
    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')
338

339 340
    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')
341

ulrich_y's avatar
ulrich_y committed
342

ulrich_y's avatar
ulrich_y committed
343
    ! Increase test coverage
344
    ref = (-1.203972804325935992622746217761838502953611_prec, pi)
ulrich_y's avatar
ulrich_y committed
345
    call test_one_flat((/ (0._prec,0.), (-0.3_prec,0.) /),ref,'E.1')
ulrich_y's avatar
ulrich_y committed
346

347
    ref = (-0.539404830685964989184608085018712655769250_prec, 1.0303768265243124637877433270311515319634364_prec)
ulrich_y's avatar
ulrich_y committed
348
    call test_one_flat([(0._prec,0._prec), (0.3_prec, 0.5_prec)],ref,'E.2')
349 350

    ref = (3.79489584700899518570534417045601144619E-5_prec,0.)
ulrich_y's avatar
ulrich_y committed
351
    call iprint('  testing GPL E.3 ...',-1)
352
    res = G((/100._prec, 100._prec, 1._prec, 0._prec, 1._prec/))
ulrich_y's avatar
ulrich_y committed
353 354
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
355
    call iprint('  testing GPL E.4 ...',-1)
ulrich_y's avatar
ulrich_y committed
356 357 358
    res = G((/100,100,1,0,1/))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
359
    call iprint('  testing GPL E.5 ...',-1)
ulrich_y's avatar
ulrich_y committed
360 361 362
    res = G_superflatn(toinum((/100,100,1,0,1/)), 5)
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
363
    call iprint('  testing GPL E.6 ...',-1)
364
    res = G((/100._prec,100._prec,1._prec,0._prec/), 1._prec)
ulrich_y's avatar
ulrich_y committed
365 366
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
367
    call iprint('  testing GPL E.7 ...',-1)
368
    res = G((/(100._prec,0.),(100._prec,0.),(1._prec,0.),(0._prec,0.)/), (1._prec,0.))
ulrich_y's avatar
ulrich_y committed
369 370
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
371
    call iprint('  testing GPL E.8 ...',-1)
372
    res = G( (/1,1,1,1/) , (/ 100._prec, 100._prec, 1._prec, 0._prec /), 1._prec)
ulrich_y's avatar
ulrich_y committed
373 374 375
    call check(res,ref)

    !ref = cmplx(0.01592795952537145)
376
    ref = (0.0485035531168477720320998551314479928648_prec,-2.82326885118931135859594797186107474877_prec)
ulrich_y's avatar
ulrich_y committed
377
    call iprint('  testing GPL E.9 ...',-1)
378
    res = G( (/ 3, 2 /), toinum((/ 1.3_prec, 1.1_prec /)), toinum(4._prec) )
ulrich_y's avatar
ulrich_y committed
379 380
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
381
    call iprint('  testing GPL E.10 ...',-1)
382
    res = G( (/ 3, 2 /), (/ 1.3_prec, 1.1_prec /), 4._prec )
ulrich_y's avatar
ulrich_y committed
383 384
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
385
    call iprint('  testing GPL E.11 ...',-1)
386 387 388
    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
389
    call iprint('  testing GPL E.12 ...',-1)
390 391
    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
392 393
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
394
    call iprint('  testing GPL E.13 ...',-1)
ulrich_y's avatar
ulrich_y committed
395 396 397 398
    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
399
    call iprint('  testing GPL E.14 ...',-1)
ulrich_y's avatar
ulrich_y committed
400 401 402 403
    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
404
    call iprint('  testing GPL E.15 ...',-1)
ulrich_y's avatar
ulrich_y committed
405 406 407 408
    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
409
    call iprint('  testing GPL E.17 ...',-1)
ulrich_y's avatar
ulrich_y committed
410 411 412 413
    ref = (0.190800137777535619036913153766083992418_prec, 0.)
    res = G((/ 6, 1, 1 /))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
414
    call iprint('  testing GPL E.18 ...',-1)
ulrich_y's avatar
ulrich_y committed
415 416 417 418 419
    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
420
    call iprint('  testing GPL E.19a ...',-1)
ulrich_y's avatar
ulrich_y committed
421 422
    res = G((/ inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
423
    call iprint('  testing GPL E.19b ...',-1)
ulrich_y's avatar
ulrich_y committed
424 425
    res = G((/ inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
426
    call iprint('  testing GPL E.19c ...',-1)
ulrich_y's avatar
ulrich_y committed
427 428
    res = G((/ inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
429
    call iprint('  testing GPL E.19d ...',-1)
ulrich_y's avatar
ulrich_y committed
430 431 432
    res = G((/ inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
433
    call iprint('  testing GPL E.19e ...',-1)
ulrich_y's avatar
ulrich_y committed
434 435
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
436
    call iprint('  testing GPL E.19f ...',-1)
ulrich_y's avatar
ulrich_y committed
437 438
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
439
    call iprint('  testing GPL E.19g ...',-1)
ulrich_y's avatar
ulrich_y committed
440 441
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
442
    call iprint('  testing GPL E.19h ...',-1)
ulrich_y's avatar
ulrich_y committed
443 444
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
Luca Naterop's avatar
Luca Naterop committed
445

Luca Naterop's avatar
Luca Naterop committed
446

ulrich_y's avatar
ulrich_y committed
447
    ref = -(2.374395270272480200677499763071638424_prec, - 1.273806204919600530933131685580471698_prec)
ulrich_y's avatar
ulrich_y committed
448
    call iprint('  testing GPL E.20a ...',-1)
ulrich_y's avatar
ulrich_y committed
449 450
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
451
    call iprint('  testing GPL E.20b ...',-1)
ulrich_y's avatar
ulrich_y committed
452 453
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
454
    call iprint('  testing GPL E.20c ...',-1)
ulrich_y's avatar
ulrich_y committed
455 456
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
457
    call iprint('  testing GPL E.20d ...',-1)
ulrich_y's avatar
ulrich_y committed
458 459 460
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
461
    call iprint('  testing GPL E.20e ...',-1)
ulrich_y's avatar
ulrich_y committed
462 463
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
464
    call iprint('  testing GPL E.20f ...',-1)
ulrich_y's avatar
ulrich_y committed
465 466
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
467
    call iprint('  testing GPL E.20g ...',-1)
ulrich_y's avatar
ulrich_y committed
468 469
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
470
    call iprint('  testing GPL E.20h ...',-1)
ulrich_y's avatar
ulrich_y committed
471 472
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
473 474 475 476


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

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

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

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

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

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

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

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

    ref = (-1.11161035333074623447215317094270858897_prec,-0.875225967273157437459269290416981432294_prec)
ulrich_y's avatar
ulrich_y committed
501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526
    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')
527

ulrich_y's avatar
ulrich_y committed
528 529
    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')
530

ulrich_y's avatar
ulrich_y committed
531 532 533 534 535 536

    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
537
    call iprint('  testing GPL F.15 ...',-1)
ulrich_y's avatar
ulrich_y committed
538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554
    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
555 556
  end subroutine do_GPL_tests
  
ulrich_y's avatar
ulrich_y committed
557 558 559



ulrich_y's avatar
ulrich_y committed
560
#ifdef HAVE_GINAC
ulrich_y's avatar
ulrich_y committed
561 562

  function evalt(arr, what)
ulrich_y's avatar
ulrich_y committed
563 564 565
#if KINDREAL==16
    use ieps, only: inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
566
    implicit none
ulrich_y's avatar
ulrich_y committed
567 568
    complex(kind=prec) :: arr(:), evalt
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
569 570 571 572 573
    integer what, i, l

    evalt =0.

    do i=1,size(arr)
574
      if (abs(arr(i)) .gt. 1.e10) then ! isnan?
ulrich_y's avatar
ulrich_y committed
575 576 577 578 579 580 581 582 583 584 585 586
        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
587 588 589
#if KINDREAL==16
      evalt = cmplx(geval(inum2inum(toinum(arr(1:l))),l),kind=prec)
#else
ulrich_y's avatar
ulrich_y committed
590
      evalt = geval(toinum(arr(1:l)),l)
ulrich_y's avatar
ulrich_y committed
591
#endif
ulrich_y's avatar
ulrich_y committed
592
    endif
ulrich_y's avatar
ulrich_y committed
593 594
  end function

ulrich_y's avatar
ulrich_y committed
595 596 597

#ifdef HAVE_MM

ulrich_y's avatar
ulrich_y committed
598
  subroutine perform_ginacv(n, args)
ulrich_y's avatar
ulrich_y committed
599
    use maths_functions, only:clearcache
ulrich_y's avatar
ulrich_y committed
600 601
    complex(kind=prec) :: args(:,:)
    integer i,n
ulrich_y's avatar
ulrich_y committed
602
    call clearcache
ulrich_y's avatar
ulrich_y committed
603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619

    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
620 621 622
    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
623

ulrich_y's avatar
ulrich_y committed
624
  end subroutine
ulrich_y's avatar
ulrich_y committed
625

ulrich_y's avatar
ulrich_y committed
626 627 628 629
  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
630 631
    integer(kind=8) cstart, cend, count_rate
    real(kind=prec) :: time(2), ttime(2)
ulrich_y's avatar
ulrich_y committed
632 633 634 635 636
    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
637
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
638 639 640
      do i=1,size(args,3)
        res=evalt(args(j,:,i),0)
      enddo
ulrich_y's avatar
ulrich_y committed
641
      call system_clock(cend, count_rate=count_rate)
642
      time(1) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
643
      if (time(1).lt.zero) print*,j, cend-cstart,count_rate
ulrich_y's avatar
ulrich_y committed
644 645
      ttime(1) = ttime(1) + time(1)

ulrich_y's avatar
ulrich_y committed
646
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
647 648 649
      do i=1,size(args,3)
        res=evalt(args(j,:,i),1)
      enddo
ulrich_y's avatar
ulrich_y committed
650
      call system_clock(cend, count_rate=count_rate)
651
      time(2) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
652 653 654 655 656 657 658 659
      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*,
660
    print*,ttime
ulrich_y's avatar
ulrich_y committed
661 662
    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
663 664 665 666 667 668
    ! 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)
669
    ttime(1) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
670 671 672 673 674 675

    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)
676
    ttime(2) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
677 678 679

    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
680
900 FORMAT(a , 'Function ',i4,'/',i4,' for ',a)
681
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
682 683 684 685
  end subroutine

  subroutine do_timing_tests(n)
    use gtestchen   , only: inichen   =>args
ulrich_y's avatar
ulrich_y committed
686
    use gtestchenff , only: inichenff =>args
ulrich_y's avatar
ulrich_y committed
687 688 689 690 691 692
    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
693
    complex(kind=prec) :: fargs(  540,5,n)
ulrich_y's avatar
ulrich_y committed
694 695 696 697 698 699 700 701 702
    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
703 704
      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
705 706 707

      w = ran2(ranseed) ! 0<w<1
      z = ran2(ranseed) * (sqrt(1-w+w**2)-sqrt(w)) + sqrt(w)
708
      nargs(:,:,i) = inimuonenp(cmplx(w,kind=prec), cmplx(z,kind=prec))
ulrich_y's avatar
ulrich_y committed
709 710 711

      x = ran2(ranseed)
      y = ran2(ranseed)
712
      pargs(:,:,i) = inimuone(cmplx(x,kind=prec), cmplx(y,kind=prec))
ulrich_y's avatar
ulrich_y committed
713 714
    enddo

ulrich_y's avatar
ulrich_y committed
715 716
    cargs(1181,:,:)=1.e15
    fargs(367,:,:)=1.e15
ulrich_y's avatar
ulrich_y committed
717 718 719


    open(unit=9, file="stats.txt")
ulrich_y's avatar
ulrich_y committed
720 721
    write(9,*) "Chen form factor"
    call do_one_speed_test(fargs,9,"Chen FF")
ulrich_y's avatar
ulrich_y committed
722 723 724 725 726 727 728 729 730
    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
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 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823

  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
824
#endif
ulrich_y's avatar
ulrich_y committed
825

ulrich_y's avatar
ulrich_y committed
826
  SUBROUTINE DO_LONG_TEST
ulrich_y's avatar
ulrich_y committed
827 828 829
#if KINDREAL==16
    use ieps,only:inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
830 831 832 833 834
    implicit none
    integer,parameter :: nzero = 10
    integer,parameter :: nieps = 30
    integer,parameter :: ncmpl = 30

835
    integer,parameter :: perweight(4) = (/ 100, 500, 30000, 50000 /)
ulrich_y's avatar
ulrich_y committed
836 837 838 839
    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
840
    integer i, j, w, seed, oldseed
ulrich_y's avatar
ulrich_y committed
841
    integer(kind=1) i0
ulrich_y's avatar
ulrich_y committed
842 843 844
    real(kind=prec) :: v, maxd
    complex(kind=prec) :: ans(2)
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
845
    character,parameter :: cr = achar(13)
ulrich_y's avatar
ulrich_y committed
846
    maxd=0._prec
ulrich_y's avatar
ulrich_y committed
847 848

    seed = 112312
ulrich_y's avatar
ulrich_y committed
849
    tol = 0.01
ulrich_y's avatar
ulrich_y committed
850

ulrich_y's avatar
ulrich_y committed
851
    open(unit=9, action='write', form='unformatted', file="long-test.txt")
ulrich_y's avatar
ulrich_y committed
852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869
    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
870 871
        oldseed = seed
        write( * , 900, advance='no' ) cr, i, perweight(w), w
ulrich_y's avatar
ulrich_y committed
872
        args(1:w) = basis((/ (1+int(size(basis)*ran2(seed)),j=1,w) /))
ulrich_y's avatar
ulrich_y committed
873 874 875 876 877
#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
878
        ans(2) = G(args(1:w),ione)
879 880 881 882 883 884
        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
885
        write(9) w,i,oldseed, abs(ans(1)-ans(2))
886
        flush(9)
ulrich_y's avatar
ulrich_y committed
887
        if(abs(ans(1)-ans(2)) > maxd) maxd = abs(ans(1)-ans(2))
ulrich_y's avatar
ulrich_y committed
888 889
        if(abs(ans(1)-ans(2)) > tol) goto 123
      enddo
ulrich_y's avatar
ulrich_y committed
890 891
      print*,' done. largest=',maxd
      maxd=0.
ulrich_y's avatar
ulrich_y committed
892
    enddo
ulrich_y's avatar
ulrich_y committed
893
    print*,maxd
ulrich_y's avatar
ulrich_y committed
894

ulrich_y's avatar
ulrich_y committed
895 896
    close(9)

ulrich_y's avatar
ulrich_y committed
897 898 899 900 901
    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
902 903 904
    close(9)

900 FORMAT(a , 'Testing ',i5,'/',i5,' GPLs with w=',i1)
ulrich_y's avatar
ulrich_y committed
905 906 907


  END SUBROUTINE
ulrich_y's avatar
ulrich_y committed
908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925
      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
926

927
#endif
ulrich_y's avatar
ulrich_y committed
928

929 930 931
  ! subroutine do_shuffle_tests() 
  !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
  !   complex(kind=prec) :: w(2) = cmplx((/3,4/))
932

933 934
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
935

Luca's avatar
Luca committed
936
END PROGRAM TEST
Luca's avatar
Luca committed
937 938

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