test.f90 30.2 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 30 31 32
        i=i+1
        call get_command_argument(i,arg)
        read(arg,*) verb               ! str to int
#else
ulrich_y's avatar
ulrich_y committed
33
        call errprint("Argument -verb is not available, compile with --debug")
ulrich_y's avatar
ulrich_y committed
34
#endif
35 36 37 38 39 40 41 42 43 44
      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
45
#if defined(HAVE_GINAC) && defined(HAVE_MM)
46 47 48 49 50 51 52 53 54 55 56
      case('-ginac-tests')
        tol = 8.0e-7
        call do_ginac_tests
      case ('-speed-tests')
        i=i+1
        call get_command_argument(i,arg)
        read(arg,*) verb
        call do_timing_tests(i)
      case ('-hw-tests')
        tol = 8.0e-7
        call do_high_weight_tests
ulrich_y's avatar
ulrich_y committed
57 58 59 60 61 62
#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
63 64 65
      case ('-long-test')
        tol = 8.0e-7
        call do_long_test
ulrich_y's avatar
ulrich_y committed
66 67 68
#else
      case('-long-test')
        call errprint("Argument "//trim(arg)//" is not available, compile with --with-ginac")
69 70 71 72 73 74 75 76 77 78
#endif
      case ('-report')
#ifdef DEBUG
        verb = 1000
#endif
        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
79
#if defined(HAVE_GINAC) && defined(HAVE_MM)
80 81 82
        call do_ginac_tests
#endif

ulrich_y's avatar
ulrich_y committed
83 84 85 86 87 88
      case default
        call errprint("Unknown argument "//trim(arg))

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

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


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

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

Luca Naterop's avatar
Luca Naterop committed
108 109 110 111
  subroutine check(res, ref)
    complex(kind=prec) :: res, ref
    real(kind=prec) :: delta

Luca Naterop's avatar
Luca Naterop committed
112
    delta = abs(res-ref)
Luca Naterop's avatar
Luca Naterop committed
113
    if(delta < tol) then
Luca Naterop's avatar
Luca Naterop committed
114
      print*, '       ',' passed with delta = ', delta
Luca Naterop's avatar
Luca Naterop committed
115
    else 
116
      print*, '  ',' FAILED with delta = ', delta
117
      tests_successful = .false.
Luca Naterop's avatar
Luca Naterop committed
118 119
    end if
  end subroutine check
120 121 122

  subroutine test_one_MPL(m,x,ref, test_id)
    integer :: m(:)
123
    complex(kind=prec) :: x(:), ref, res
124
    character(len=*) :: test_id
Luca's avatar
Luca committed
125
    
126 127 128 129 130
    print*, '  ', 'testing MPL ', test_id, ' ...'
    res = MPL(m,x)
    call check(res,ref)
  end subroutine test_one_MPL

Luca's avatar
Luca committed
131
  subroutine do_MPL_tests()
Luca's avatar
Luca committed
132
    complex(kind=prec) :: ref
Luca Naterop's avatar
Luca Naterop committed
133
    print*, 'doing MPL tests...'
Luca's avatar
Luca committed
134
    
ulrich_y's avatar
ulrich_y committed
135 136 137 138 139
    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
140
    
ulrich_y's avatar
ulrich_y committed
141 142
    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
143
    
144 145
    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
146

147
    ref = (-0.0319989639656484125354488353677231989818697987217232462808_prec,0.)
ulrich_y's avatar
ulrich_y committed
148
    call test_one_MPL((/2, 1/), (/(-0.25_prec,0.),(-2._prec,0.) /), ref, '1.5')
Luca's avatar
Luca committed
149
  end subroutine do_MPL_tests
Luca's avatar
Luca committed
150

Luca's avatar
Luca committed
151
  subroutine test_one_condensed(m,z,y,k,ref,test_id)
152 153 154
    integer :: m(:), k
    complex(kind=prec) :: z(:), y, res, ref
    character(len=*) :: test_id
Luca Naterop's avatar
Luca Naterop committed
155

156
    print*, '  ', 'testing GPL ', test_id, ' ...'
ulrich_y's avatar
ulrich_y committed
157
    res = G_condensed(m,toinum(z),inum(y,di0),k)
Luca Naterop's avatar
Luca Naterop committed
158
    call check(res,ref)
Luca's avatar
Luca committed
159 160
  end subroutine test_one_condensed

161 162
  subroutine test_one_flat(z,ref,test_id)
    complex(kind=prec) :: z(:), res, ref
Luca's avatar
Luca committed
163 164 165
    character(len=*) :: test_id

    print*, '  ', 'testing GPL ', test_id, ' ...'
ulrich_y's avatar
ulrich_y committed
166
    res = G(z)
Luca's avatar
Luca committed
167 168
    call check(res,ref)
  end subroutine test_one_flat
Luca Naterop's avatar
Luca Naterop committed
169 170

  subroutine do_GPL_tests()
ulrich_y's avatar
ulrich_y committed
171
    complex(kind=prec) :: ref, res
Luca Naterop's avatar
Luca Naterop committed
172
    real(kind=prec) :: z, xchen
Luca Naterop's avatar
Luca Naterop committed
173
    print*, 'doing GPL tests...'
Luca's avatar
Luca committed
174
    
175 176
    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
177
    
178 179 180 181 182
    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
183

Luca Naterop's avatar
Luca Naterop committed
184
    ! requires making convergent
185 186
    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
187

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

    ! requires hoelder convolution
192 193
    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
194 195

    ! here the tests from the mathematica nb start
196 197
    ! --------------------------------------------------------------------------

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

200 201
    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
202

203 204
    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
205

206 207
    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
208

209 210
    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
211
    
212 213
    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
214
    
215 216
    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
217
    
218 219
    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')
220 221 222 223 224


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

225 226 227
    z = 1/100._prec
    ref = (-1.203972804325935992622746217761838502954_prec,0.0)
    call test_one_flat((/ (0._prec,0.), (0.3_prec,0.) /),ref,'4.1')
228

229 230
    ref = (10.6037962209567960211233327771880353832_prec,0.0)
    call test_one_flat( (/ (0._prec,0.), (0._prec,0.), (0.01_prec,0.) /),ref,'4.2')
231

232 233
    ref = (5.0429900651807654463744449762857296353E-4_prec,0.0)
    call test_one_flat((/(100._prec,0.), (1._prec,0.), (0.3_prec,0.) /),ref,'4.3')
234

235 236
    ref = (0.05630861877185110944118569266547272128_prec,0.0)
    call test_one_flat( (/ (1._prec,0.), (0._prec,0.), (0.01_prec,0.) /) ,ref,'4.4')
237

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

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

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

247 248
    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')
249

250 251
    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')
252

253 254
    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')
255

256 257
    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')
258

ulrich_y's avatar
ulrich_y committed
259

ulrich_y's avatar
ulrich_y committed
260
    ! Increase test coverage
261
    ref = (-1.203972804325935992622746217761838502953611_prec, pi)
ulrich_y's avatar
ulrich_y committed
262
    call test_one_flat((/ (0._prec,0.), (-0.3_prec,0.) /),ref,'E.1')
ulrich_y's avatar
ulrich_y committed
263

264
    ref = (-0.539404830685964989184608085018712655769250_prec, 1.0303768265243124637877433270311515319634364_prec)
ulrich_y's avatar
ulrich_y committed
265
    call test_one_flat([(0._prec,0._prec), (0.3_prec, 0.5_prec)],ref,'E.2')
266 267

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

ulrich_y's avatar
ulrich_y committed
272
    print*, '  ', 'testing GPL ', 'E.4', ' ...'
ulrich_y's avatar
ulrich_y committed
273 274 275
    res = G((/100,100,1,0,1/))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
276
    print*, '  ', 'testing GPL ', 'E.5', ' ...'
ulrich_y's avatar
ulrich_y committed
277 278 279
    res = G_superflatn(toinum((/100,100,1,0,1/)), 5)
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
280
    print*, '  ', 'testing GPL ', 'E.6', ' ...'
281
    res = G((/100._prec,100._prec,1._prec,0._prec/), 1._prec)
ulrich_y's avatar
ulrich_y committed
282 283
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
284
    print*, '  ', 'testing GPL ', 'E.7', ' ...'
285
    res = G((/(100._prec,0.),(100._prec,0.),(1._prec,0.),(0._prec,0.)/), (1._prec,0.))
ulrich_y's avatar
ulrich_y committed
286 287
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
288
    print*, '  ', 'testing GPL ', 'E.8', ' ...'
289
    res = G( (/1,1,1,1/) , (/ 100._prec, 100._prec, 1._prec, 0._prec /), 1._prec)
ulrich_y's avatar
ulrich_y committed
290 291 292
    call check(res,ref)

    !ref = cmplx(0.01592795952537145)
293
    ref = (0.0485035531168477720320998551314479928648_prec,-2.82326885118931135859594797186107474877_prec)
ulrich_y's avatar
ulrich_y committed
294
    print*, '  ', 'testing GPL ', 'E.9', ' ...'
295
    res = G( (/ 3, 2 /), toinum((/ 1.3_prec, 1.1_prec /)), toinum(4._prec) )
ulrich_y's avatar
ulrich_y committed
296 297
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
298
    print*, '  ', 'testing GPL ', 'E.10', ' ...'
299
    res = G( (/ 3, 2 /), (/ 1.3_prec, 1.1_prec /), 4._prec )
ulrich_y's avatar
ulrich_y committed
300 301
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
302
    print*, '  ', 'testing GPL ', 'E.11', ' ...'
303 304 305
    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
306
    print*, '  ', 'testing GPL ', 'E.12', ' ...'
307 308
    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
309 310
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
311
    print*, '  ', 'testing GPL ', 'E.13', ' ...'
ulrich_y's avatar
ulrich_y committed
312 313 314 315
    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
316
    print*, '  ', 'testing GPL ', 'E.14', ' ...'
ulrich_y's avatar
ulrich_y committed
317 318 319 320
    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
321
    print*, '  ', 'testing GPL ', 'E.15', ' ...'
ulrich_y's avatar
ulrich_y committed
322 323 324 325
    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
326
    print*, '  ', 'testing GPL ', 'E.17', ' ...'
ulrich_y's avatar
ulrich_y committed
327 328 329 330
    ref = (0.190800137777535619036913153766083992418_prec, 0.)
    res = G((/ 6, 1, 1 /))
    call check(res,ref)

ulrich_y's avatar
ulrich_y committed
331
    print*, '  ', 'testing GPL ', 'E.18', ' ...'
ulrich_y's avatar
ulrich_y committed
332 333 334 335 336
    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
337
    print*, '  ', 'testing GPL ', 'E.19a', ' ...'
ulrich_y's avatar
ulrich_y committed
338 339
    res = G((/ inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
340
    print*, '  ', 'testing GPL ', 'E.19b', ' ...'
ulrich_y's avatar
ulrich_y committed
341 342
    res = G((/ inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
343
    print*, '  ', 'testing GPL ', 'E.19c', ' ...'
ulrich_y's avatar
ulrich_y committed
344 345
    res = G((/ inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
346
    print*, '  ', 'testing GPL ', 'E.19d', ' ...'
ulrich_y's avatar
ulrich_y committed
347 348 349
    res = G((/ inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
350
    print*, '  ', 'testing GPL ', 'E.19e', ' ...'
ulrich_y's avatar
ulrich_y committed
351 352
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
353
    print*, '  ', 'testing GPL ', 'E.19f', ' ...'
ulrich_y's avatar
ulrich_y committed
354 355
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
356
    print*, '  ', 'testing GPL ', 'E.19g', ' ...'
ulrich_y's avatar
ulrich_y committed
357 358
    res = G((/ inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
359
    print*, '  ', 'testing GPL ', 'E.19h', ' ...'
ulrich_y's avatar
ulrich_y committed
360 361
    res = G((/ inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
Luca Naterop's avatar
Luca Naterop committed
362

Luca Naterop's avatar
Luca Naterop committed
363

ulrich_y's avatar
ulrich_y committed
364
    ref = -(2.374395270272480200677499763071638424_prec, - 1.273806204919600530933131685580471698_prec)
ulrich_y's avatar
ulrich_y committed
365
    print*, '  ', 'testing GPL ', 'E.20a', ' ...'
ulrich_y's avatar
ulrich_y committed
366 367
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
368
    print*, '  ', 'testing GPL ', 'E.20b', ' ...'
ulrich_y's avatar
ulrich_y committed
369 370
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
371
    print*, '  ', 'testing GPL ', 'E.20c', ' ...'
ulrich_y's avatar
ulrich_y committed
372 373
    res = G((/ izero, inum(2._prec, +1) /), inum(3._prec, -1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
374
    print*, '  ', 'testing GPL ', 'E.20d', ' ...'
ulrich_y's avatar
ulrich_y committed
375 376 377
    res = G((/ izero, inum(2._prec, -1) /), inum(3._prec, -1))
    call check(res,conjg(ref))

ulrich_y's avatar
ulrich_y committed
378
    print*, '  ', 'testing GPL ', 'E.20e', ' ...'
ulrich_y's avatar
ulrich_y committed
379 380
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, +1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
381
    print*, '  ', 'testing GPL ', 'E.20f', ' ...'
ulrich_y's avatar
ulrich_y committed
382 383
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, +1))
    call check(res,ref)
ulrich_y's avatar
ulrich_y committed
384
    print*, '  ', 'testing GPL ', 'E.20g', ' ...'
ulrich_y's avatar
ulrich_y committed
385 386
    res = G((/ izero, inum(-2._prec, +1) /), inum(-3._prec, -1))
    call check(res,conjg(ref))
ulrich_y's avatar
ulrich_y committed
387
    print*, '  ', 'testing GPL ', 'E.20h', ' ...'
ulrich_y's avatar
ulrich_y committed
388 389
    res = G((/ izero, inum(-2._prec, -1) /), inum(-3._prec, -1))
    call check(res,ref)
390 391 392 393


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

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

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

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

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

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

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

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

    ref = (-1.11161035333074623447215317094270858897_prec,-0.875225967273157437459269290416981432294_prec)
ulrich_y's avatar
ulrich_y committed
418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443
    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')
444

ulrich_y's avatar
ulrich_y committed
445 446
    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')
447

ulrich_y's avatar
ulrich_y committed
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471

    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
472 473
  end subroutine do_GPL_tests
  
ulrich_y's avatar
ulrich_y committed
474 475 476



ulrich_y's avatar
ulrich_y committed
477
#ifdef HAVE_GINAC
ulrich_y's avatar
ulrich_y committed
478 479

  function evalt(arr, what)
ulrich_y's avatar
ulrich_y committed
480 481 482
#if KINDREAL==16
    use ieps, only: inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
483
    implicit none
ulrich_y's avatar
ulrich_y committed
484 485
    complex(kind=prec) :: arr(:), evalt
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
486 487 488 489 490
    integer what, i, l

    evalt =0.

    do i=1,size(arr)
491
      if (abs(arr(i)) .gt. 1.e10) then ! isnan?
ulrich_y's avatar
ulrich_y committed
492 493 494 495 496 497 498 499 500 501 502 503
        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
504 505 506
#if KINDREAL==16
      evalt = cmplx(geval(inum2inum(toinum(arr(1:l))),l),kind=prec)
#else
ulrich_y's avatar
ulrich_y committed
507
      evalt = geval(toinum(arr(1:l)),l)
ulrich_y's avatar
ulrich_y committed
508
#endif
ulrich_y's avatar
ulrich_y committed
509
    endif
ulrich_y's avatar
ulrich_y committed
510 511
  end function

ulrich_y's avatar
ulrich_y committed
512 513 514

#ifdef HAVE_MM

ulrich_y's avatar
ulrich_y committed
515
  subroutine perform_ginacv(n, args)
ulrich_y's avatar
ulrich_y committed
516
    use maths_functions, only:clearcache
ulrich_y's avatar
ulrich_y committed
517 518
    complex(kind=prec) :: args(:,:)
    integer i,n
ulrich_y's avatar
ulrich_y committed
519
    call clearcache
ulrich_y's avatar
ulrich_y committed
520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536

    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
537 538 539
    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
540

ulrich_y's avatar
ulrich_y committed
541
  end subroutine
ulrich_y's avatar
ulrich_y committed
542

ulrich_y's avatar
ulrich_y committed
543 544 545 546
  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
547 548
    integer(kind=8) cstart, cend, count_rate
    real(kind=prec) :: time(2), ttime(2)
ulrich_y's avatar
ulrich_y committed
549 550 551 552 553
    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
554
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
555 556 557
      do i=1,size(args,3)
        res=evalt(args(j,:,i),0)
      enddo
ulrich_y's avatar
ulrich_y committed
558
      call system_clock(cend, count_rate=count_rate)
559
      time(1) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
560
      if (time(1).lt.zero) print*,j, cend-cstart,count_rate
ulrich_y's avatar
ulrich_y committed
561 562
      ttime(1) = ttime(1) + time(1)

ulrich_y's avatar
ulrich_y committed
563
      call system_clock(cstart, count_rate=count_rate)
ulrich_y's avatar
ulrich_y committed
564 565 566
      do i=1,size(args,3)
        res=evalt(args(j,:,i),1)
      enddo
ulrich_y's avatar
ulrich_y committed
567
      call system_clock(cend, count_rate=count_rate)
568
      time(2) = real(cend-cstart)/real(count_rate,kind=prec)/size(args,3)
ulrich_y's avatar
ulrich_y committed
569 570 571 572 573 574 575 576
      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*,
577
    print*,ttime
ulrich_y's avatar
ulrich_y committed
578 579
    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
580 581 582 583 584 585
    ! 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)
586
    ttime(1) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
587 588 589 590 591 592

    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)
593
    ttime(2) = real(cend-cstart)/real(count_rate,kind=prec)
ulrich_y's avatar
ulrich_y committed
594 595 596

    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
597
900 FORMAT(a , 'Function ',i4,'/',i4,' for ',a)
598
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
599 600 601 602
  end subroutine

  subroutine do_timing_tests(n)
    use gtestchen   , only: inichen   =>args
ulrich_y's avatar
ulrich_y committed
603
    use gtestchenff , only: inichenff =>args
ulrich_y's avatar
ulrich_y committed
604 605 606 607 608 609
    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
610
    complex(kind=prec) :: fargs(  540,5,n)
ulrich_y's avatar
ulrich_y committed
611 612 613 614 615 616 617 618 619
    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
620 621
      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
622 623 624

      w = ran2(ranseed) ! 0<w<1
      z = ran2(ranseed) * (sqrt(1-w+w**2)-sqrt(w)) + sqrt(w)
625
      nargs(:,:,i) = inimuonenp(cmplx(w,kind=prec), cmplx(z,kind=prec))
ulrich_y's avatar
ulrich_y committed
626 627 628

      x = ran2(ranseed)
      y = ran2(ranseed)
629
      pargs(:,:,i) = inimuone(cmplx(x,kind=prec), cmplx(y,kind=prec))
ulrich_y's avatar
ulrich_y committed
630 631
    enddo

ulrich_y's avatar
ulrich_y committed
632 633
    cargs(1181,:,:)=1.e15
    fargs(367,:,:)=1.e15
ulrich_y's avatar
ulrich_y committed
634 635 636


    open(unit=9, file="stats.txt")
ulrich_y's avatar
ulrich_y committed
637 638
    write(9,*) "Chen form factor"
    call do_one_speed_test(fargs,9,"Chen FF")
ulrich_y's avatar
ulrich_y committed
639 640 641 642 643 644 645 646 647
    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
648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740

  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
741
#endif
ulrich_y's avatar
ulrich_y committed
742

ulrich_y's avatar
ulrich_y committed
743
  SUBROUTINE DO_LONG_TEST
ulrich_y's avatar
ulrich_y committed
744 745 746
#if KINDREAL==16
    use ieps,only:inum2inum
#endif
ulrich_y's avatar
ulrich_y committed
747 748 749 750 751
    implicit none
    integer,parameter :: nzero = 10
    integer,parameter :: nieps = 30
    integer,parameter :: ncmpl = 30

752
    integer,parameter :: perweight(4) = (/ 100, 500, 30000, 50000 /)
ulrich_y's avatar
ulrich_y committed
753 754 755 756
    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
757
    integer i, j, w, seed, oldseed
ulrich_y's avatar
ulrich_y committed
758
    integer(kind=1) i0
ulrich_y's avatar
ulrich_y committed
759 760 761
    real(kind=prec) :: v, maxd
    complex(kind=prec) :: ans(2)
    complex(kind=8) :: geval
ulrich_y's avatar
ulrich_y committed
762
    character,parameter :: cr = achar(13)
ulrich_y's avatar
ulrich_y committed
763
    maxd=0._prec
ulrich_y's avatar
ulrich_y committed
764 765

    seed = 112312
ulrich_y's avatar
ulrich_y committed
766
    tol = 0.01
ulrich_y's avatar
ulrich_y committed
767

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

ulrich_y's avatar
ulrich_y committed
812 813
    close(9)

ulrich_y's avatar
ulrich_y committed
814 815 816 817 818
    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
819 820 821
    close(9)

900 FORMAT(a , 'Testing ',i5,'/',i5,' GPLs with w=',i1)
ulrich_y's avatar
ulrich_y committed
822 823 824


  END SUBROUTINE
ulrich_y's avatar
ulrich_y committed
825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842
      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
843

844
#endif
ulrich_y's avatar
ulrich_y committed
845

846 847 848
  ! subroutine do_shuffle_tests() 
  !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
  !   complex(kind=prec) :: w(2) = cmplx((/3,4/))
849

850 851
  !   call print_matrix(shuffle_product(v,w))
  ! end subroutine do_shuffle_tests
852

Luca's avatar
Luca committed
853
END PROGRAM TEST
Luca's avatar
Luca committed
854 855

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