maths_functions.f90 19.9 KB
Newer Older
1 2 3

MODULE maths_functions
  use globals
ulrich_y's avatar
ulrich_y committed
4 5
  use ieps
  use utils, only: factorial
6
  implicit none
ulrich_y's avatar
ulrich_y committed
7 8 9
  interface polylog
    module procedure polylog1, polylog2
  end interface polylog
10

ulrich_y's avatar
ulrich_y committed
11 12 13 14 15
  real(kind=prec), parameter :: zeta(2:10) = (/ 1.6449340668482264364724151666460251892189499012068_prec,  &
      1.2020569031595942853997381615114499907649862923405_prec, 1.0823232337111381915160036965411679027747509519187_prec,  &
      1.0369277551433699263313654864570341680570809195019_prec, 1.0173430619844491397145179297909205279018174900329_prec,  &
      1.0083492773819228268397975498497967595998635605652_prec, 1.0040773561979443393786852385086524652589607906499_prec,  &
      1.0020083928260822144178527692324120604856058513949_prec, 1.0009945751278180853371459589003190170060195315645_prec /)
ulrich_y's avatar
ulrich_y committed
16 17 18 19 20 21 22
  type el
    type(inum) :: c
    complex(kind=prec) ans
  end type el

  type(el) :: cache(PolyLogCacheSize(1),PolyLogCacheSize(2))
  integer :: plcachesize(PolyLogCacheSize(1)) = 0
ulrich_y's avatar
ulrich_y committed
23
CONTAINS
24

ulrich_y's avatar
ulrich_y committed
25
  FUNCTION naive_polylog(m,x) result(res)
26 27 28
    ! Computes the classical polylogarithm Li_m(x) using series representation up to order n
    integer :: m
    complex(kind=prec) :: x, res
29
    integer(kind=ikin) :: i
ulrich_y's avatar
ulrich_y committed
30
    res=0._prec
31
    do i=1,PolylogInfinity
32
      if(i**m.lt.0) return ! roll over
ulrich_y's avatar
ulrich_y committed
33
      if(abs(x**i).lt.1.e-250_prec) return
34 35
      res = res+x**i/i**m
    enddo
36
  END FUNCTION naive_polylog
Luca's avatar
Luca committed
37

38 39 40 41
  FUNCTION Li2(x)

   !! Dilogarithm for arguments x < = 1.0

ulrich_y's avatar
ulrich_y committed
42
   real (kind=prec):: X,Y,T,S,A,ZERO,ONE,HALF,MALF,MONE,MTWO
ulrich_y's avatar
ulrich_y committed
43
   real (kind=prec):: C(0:42),H,ALFA,B0,B1,B2,LI2_OLD
44
   real (kind=prec):: Li2
ulrich_y's avatar
ulrich_y committed
45
   integer :: i, maxi
46 47 48 49 50

   DATA ZERO /0.0_prec/, ONE /1.0_prec/
   DATA HALF /0.5_prec/, MALF /-0.5_prec/ 
   DATA MONE /-1.0_prec/, MTWO /-2.0_prec/

ulrich_y's avatar
ulrich_y committed
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
   DATA C( 0) / 0.42996693560813697203703367869938799_prec/
   DATA C( 1) / 0.40975987533077105846826371092525528_prec/
   DATA C( 2) /-0.01858843665014591964764164914021227_prec/
   DATA C( 3) / 0.00145751084062267855367392841645949_prec/
   DATA C( 4) /-0.00014304184442340048774883012009088_prec/
   DATA C( 5) / 0.00001588415541879553236190550471677_prec/
   DATA C( 6) /-0.00000190784959386582722719632114209_prec/
   DATA C( 7) / 0.00000024195180854164749499461464343_prec/
   DATA C( 8) /-0.00000003193341274251783460496014143_prec/
   DATA C( 9) / 0.00000000434545062676912298795717848_prec/
   DATA C(10) /-0.00000000060578480118407444429705331_prec/
   DATA C(11) / 0.00000000008612097799359498244283685_prec/
   DATA C(12) /-0.00000000001244331659938867989642421_prec/
   DATA C(13) / 0.00000000000182255696235736330065548_prec/
   DATA C(14) /-0.00000000000027006766049114651808572_prec/
   DATA C(15) / 0.00000000000004042209263152664648833_prec/
   DATA C(16) /-0.00000000000000610325145269187950378_prec/
   DATA C(17) / 0.00000000000000092862975330195758613_prec/
   DATA C(18) /-0.00000000000000014226020855112446840_prec/
   DATA C(19) / 0.00000000000000002192631718153957354_prec/
   DATA C(20) /-0.00000000000000000339797324215897863_prec/
   DATA C(21) / 0.00000000000000000052919542448331471_prec/
   DATA C(22) /-0.00000000000000000008278580814278998_prec/
   DATA C(23) / 0.00000000000000000001300371734545560_prec/
   DATA C(24) /-0.00000000000000000000205022224255282_prec/
   DATA C(25) / 0.00000000000000000000032435785491489_prec/
   DATA C(26) /-0.00000000000000000000005147799903343_prec/
   DATA C(27) / 0.00000000000000000000000819387747717_prec/
   DATA C(28) /-0.00000000000000000000000130778354057_prec/
   DATA C(29) / 0.00000000000000000000000020925629306_prec/
   DATA C(30) /-0.00000000000000000000000003356166151_prec/
   DATA C(31) / 0.00000000000000000000000000539465777_prec/
   DATA C(32) /-0.00000000000000000000000000086891932_prec/
   DATA C(33) / 0.00000000000000000000000000014022817_prec/
   DATA C(34) /-0.00000000000000000000000000002267156_prec/
   DATA C(35) / 0.00000000000000000000000000000367174_prec/
   DATA C(36) /-0.00000000000000000000000000000059562_prec/
   DATA C(37) / 0.00000000000000000000000000000009677_prec/
   DATA C(38) /-0.00000000000000000000000000000001574_prec/
   DATA C(39) / 0.00000000000000000000000000000000257_prec/
   DATA C(40) /-0.00000000000000000000000000000000042_prec/
   DATA C(41) / 0.00000000000000000000000000000000007_prec/
   DATA C(42) /-0.00000000000000000000000000000000001_prec/
94 95 96 97 98 99 100 101

   if(X > 1.00000000001_prec) then
     print*, 'crashes because Li called with bad arguments'
   elseif(X > 1.0_prec) then
     X = 1._prec
   endif    

   IF(X > 0.999999_prec) THEN
ulrich_y's avatar
ulrich_y committed
102
    LI2_OLD=zeta(2)
103 104
    Li2 = Real(LI2_OLD,prec)
    RETURN
105
   ELSE IF(abs(x-MONE) < zero) THEN
ulrich_y's avatar
ulrich_y committed
106
    LI2_OLD=MALF*zeta(2)
107 108 109 110 111 112
    RETURN
   END IF
   T=-X
   IF(T .LE. MTWO) THEN
    Y=MONE/(ONE+T)
    S=ONE
ulrich_y's avatar
ulrich_y committed
113
    A=-2*zeta(2)+HALF*(LOG(-T)**2-LOG(ONE+ONE/T)**2)
114 115 116 117
   ELSE IF(T .LT. MONE) THEN
    Y=MONE-T
    S=MONE
    A=LOG(-T)
ulrich_y's avatar
ulrich_y committed
118
    A=-zeta(2)+A*(A+LOG(ONE+ONE/T))
119 120 121 122
   ELSE IF(T .LE. MALF) THEN
    Y=(MONE-T)/T
    S=ONE
    A=LOG(-T)
ulrich_y's avatar
ulrich_y committed
123
    A=-zeta(2)+A*(MALF*A+LOG(ONE+T))
124 125 126 127 128 129 130 131 132 133 134
   ELSE IF(T .LT. ZERO) THEN
    Y=-T/(ONE+T)
    S=MONE
    A=HALF*LOG(ONE+T)**2
   ELSE IF(T .LE. ONE) THEN
    Y=T
    S=ONE
    A=ZERO
   ELSE
    Y=ONE/T
    S=MONE
ulrich_y's avatar
ulrich_y committed
135
    A=zeta(2)+HALF*LOG(T)**2
136 137 138 139 140 141
   END IF

   H=Y+Y-ONE
   ALFA=H+H
   B1=ZERO
   B2=ZERO
ulrich_y's avatar
ulrich_y committed
142 143 144 145 146 147
   if (precision(1._prec) < 20) then
     maxi = 18
   else
     maxi = 42
   endif
   DO  I = maxi,0,-1
148 149 150 151 152 153 154 155
     B0=C(I)+ALFA*B1-B2
     B2=B1
     B1=B0
   ENDDO
   LI2_OLD=-(S*(B0-H*B2)+A)
         ! Artificial conversion           
   Li2 = Real(LI2_OLD,prec)
  END FUNCTION Li2
Luca's avatar
Luca committed
156

157
  RECURSIVE FUNCTION dilog(x) result(res)
ulrich_y's avatar
ulrich_y committed
158
    ! evaluates dilog for any argument |x|<1
159 160
    complex(kind=prec) :: res
    complex(kind=prec) :: x
Luca's avatar
Luca committed
161

ulrich_y's avatar
ulrich_y committed
162
    if(abs(aimag(x)) < zero ) then
ulrich_y's avatar
ulrich_y committed
163
      res = Li2(real(x,kind=prec))
164
    else
ulrich_y's avatar
ulrich_y committed
165 166
      res = naive_polylog(2,x)
    endif
167 168
  END FUNCTION dilog

Luca's avatar
Luca committed
169 170 171 172 173 174
  FUNCTION Li3(x)
    ! Trilogarithm for arguments x < = 1.0
    ! This was hacked from LI2 to also follow C332
    ! In theory this could also produce Re[Li [x]] for x>1

    real (kind=prec):: X,S,A
ulrich_y's avatar
ulrich_y committed
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271
    real (kind=prec):: CA(0:52),HA,ALFAA,BA0,BA1,BA2, YA
    real (kind=prec):: CB(0:52),HB,ALFAB,BB0,BB1,BB2, YB

    DATA CA( 0) / 0.46172939286012092817954516381760016_prec/
    DATA CA( 1) / 0.45017399588550288560580364647352070_prec/
    DATA CA( 2) /-0.01091284195229295374494914320402658_prec/
    DATA CA( 3) / 0.00059324547127243642952756961712713_prec/
    DATA CA( 4) /-0.00004479593219280756178998757870776_prec/
    DATA CA( 5) / 0.00000405154578580684540984293800468_prec/
    DATA CA( 6) /-0.00000041095398606214457668547736075_prec/
    DATA CA( 7) / 0.00000004513178770934181970313262557_prec/
    DATA CA( 8) /-0.00000000525466158515342604029419927_prec/
    DATA CA( 9) / 0.00000000063982547910449452549291936_prec/
    DATA CA(10) /-0.00000000008071938839872532971820424_prec/
    DATA CA(11) / 0.00000000001048073680087126094928657_prec/
    DATA CA(12) /-0.00000000000139365085138335067524094_prec/
    DATA CA(13) / 0.00000000000018907205037339730044704_prec/
    DATA CA(14) /-0.00000000000002609371657183250621931_prec/
    DATA CA(15) / 0.00000000000000365481859219879483309_prec/
    DATA CA(16) /-0.00000000000000051855842492271228151_prec/
    DATA CA(17) / 0.00000000000000007441491722173878908_prec/
    DATA CA(18) /-0.00000000000000001078686838424874221_prec/
    DATA CA(19) / 0.00000000000000000157774237809543778_prec/
    DATA CA(20) /-0.00000000000000000023264073800573828_prec/
    DATA CA(21) / 0.00000000000000000003455457587964154_prec/
    DATA CA(22) /-0.00000000000000000000516658458392580_prec/
    DATA CA(23) / 0.00000000000000000000077718849383139_prec/
    DATA CA(24) /-0.00000000000000000000011755815708807_prec/
    DATA CA(25) / 0.00000000000000000000001787262690583_prec/
    DATA CA(26) /-0.00000000000000000000000272999302683_prec/
    DATA CA(27) / 0.00000000000000000000000041881267359_prec/
    DATA CA(28) /-0.00000000000000000000000006451004176_prec/
    DATA CA(29) / 0.00000000000000000000000000997383916_prec/
    DATA CA(30) /-0.00000000000000000000000000154744603_prec/
    DATA CA(31) / 0.00000000000000000000000000024087296_prec/
    DATA CA(32) /-0.00000000000000000000000000003760889_prec/
    DATA CA(33) / 0.00000000000000000000000000000588900_prec/
    DATA CA(34) /-0.00000000000000000000000000000092463_prec/
    DATA CA(35) / 0.00000000000000000000000000000014555_prec/
    DATA CA(36) /-0.00000000000000000000000000000002297_prec/
    DATA CA(37) / 0.00000000000000000000000000000000363_prec/
    DATA CA(38) /-0.00000000000000000000000000000000058_prec/
    DATA CA(39) / 0.00000000000000000000000000000000009_prec/
    DATA CA(40) /-0.00000000000000000000000000000000001_prec/
    DATA CB( 0) /-0.01601618044919582873670691984756338_prec/
    DATA CB( 1) /-0.50364244007530129181209541016960792_prec/
    DATA CB( 2) /-0.01615099243050025888745446951929454_prec/
    DATA CB( 3) /-0.00124402421042449361265610524413112_prec/
    DATA CB( 4) /-0.00013757218124461673921971996409271_prec/
    DATA CB( 5) /-0.00001856381852603773316486795183129_prec/
    DATA CB( 6) /-0.00000284173534515440415934505790039_prec/
    DATA CB( 7) /-0.00000047459967905789937221638951390_prec/
    DATA CB( 8) /-0.00000008448038543781200676091819474_prec/
    DATA CB( 9) /-0.00000001578767124043400543246870475_prec/
    DATA CB(10) /-0.00000000306576207139903798128889004_prec/
    DATA CB(11) /-0.00000000061407921728125845808062189_prec/
    DATA CB(12) /-0.00000000012618830243156719690872484_prec/
    DATA CB(13) /-0.00000000002649314179819609957126783_prec/
    DATA CB(14) /-0.00000000000566470854636425926158812_prec/
    DATA CB(15) /-0.00000000000123041115779581117517467_prec/
    DATA CB(16) /-0.00000000000027093457836786768143960_prec/
    DATA CB(17) /-0.00000000000006038026463383701279197_prec/
    DATA CB(18) /-0.00000000000001360008993995749682352_prec/
    DATA CB(19) /-0.00000000000000309244740631856875855_prec/
    DATA CB(20) /-0.00000000000000070917249609207158220_prec/
    DATA CB(21) /-0.00000000000000016388083639226002471_prec/
    DATA CB(22) /-0.00000000000000003813464350168994613_prec/
    DATA CB(23) /-0.00000000000000000893010739611811656_prec/
    DATA CB(24) /-0.00000000000000000210331341599359416_prec/
    DATA CB(25) /-0.00000000000000000049802988416537866_prec/
    DATA CB(26) /-0.00000000000000000011850292695597351_prec/
    DATA CB(27) /-0.00000000000000000002832460494402074_prec/
    DATA CB(28) /-0.00000000000000000000679854955943073_prec/
    DATA CB(29) /-0.00000000000000000000163816629435900_prec/
    DATA CB(30) /-0.00000000000000000000039616291258646_prec/
    DATA CB(31) /-0.00000000000000000000009613022139972_prec/
    DATA CB(32) /-0.00000000000000000000002340035706102_prec/
    DATA CB(33) /-0.00000000000000000000000571315840877_prec/
    DATA CB(34) /-0.00000000000000000000000139876183805_prec/
    DATA CB(35) /-0.00000000000000000000000034336361321_prec/
    DATA CB(36) /-0.00000000000000000000000008449733573_prec/
    DATA CB(37) /-0.00000000000000000000000002084253881_prec/
    DATA CB(38) /-0.00000000000000000000000000515255292_prec/
    DATA CB(39) /-0.00000000000000000000000000127646290_prec/
    DATA CB(40) /-0.00000000000000000000000000031685555_prec/
    DATA CB(41) /-0.00000000000000000000000000007880228_prec/
    DATA CB(42) /-0.00000000000000000000000000001963363_prec/
    DATA CB(43) /-0.00000000000000000000000000000490016_prec/
    DATA CB(44) /-0.00000000000000000000000000000122499_prec/
    DATA CB(45) /-0.00000000000000000000000000000030671_prec/
    DATA CB(46) /-0.00000000000000000000000000000007691_prec/
    DATA CB(47) /-0.00000000000000000000000000000001931_prec/
    DATA CB(48) /-0.00000000000000000000000000000000486_prec/
    DATA CB(49) /-0.00000000000000000000000000000000122_prec/
    DATA CB(50) /-0.00000000000000000000000000000000031_prec/
    DATA CB(51) /-0.00000000000000000000000000000000008_prec/
    DATA CB(52) /-0.00000000000000000000000000000000002_prec/
Luca's avatar
Luca committed
272
    real (kind=prec):: Li3
ulrich_y's avatar
ulrich_y committed
273
    integer :: i, maxi
Luca's avatar
Luca committed
274 275 276 277 278 279 280 281 282


    if(x > 1.00000000001_prec) then
      print*, 'need to crash Li3, since not convergent'
    elseif(x > 1.0_prec) then
      x = 1._prec
    endif

    IF(X > 0.999999_prec) THEN
ulrich_y's avatar
ulrich_y committed
283
      LI3=zeta(3)
Luca's avatar
Luca committed
284
    RETURN
285
    ELSE IF( abs(x+1) < zero) THEN
ulrich_y's avatar
ulrich_y committed
286
      LI3=-0.75_prec*zeta(3)
Luca's avatar
Luca committed
287 288 289 290 291
    RETURN
    END IF
    IF(X .LE. -1._prec) THEN
      YA=1._prec/x ; YB=0._prec
      S=-1._prec
ulrich_y's avatar
ulrich_y committed
292
      A=-LOG(-X)*(zeta(2)+LOG(-x)**2/6._prec)
Luca's avatar
Luca committed
293 294 295 296 297 298 299 300 301 302 303
    ELSE IF(X .LE. 0._prec) THEN
      YA=x ; YB=0._prec
      S=-1._prec
      A=0._prec
    ELSE IF(X .LE. 0.5_prec) THEN
      YA=0._prec ; YB=x
      S=-1._prec
      A=0._prec
    ELSE IF(X .LE. 1._prec) THEN
      YA=(x-1._prec)/x ; YB=1._prec-x
      S=1._prec
ulrich_y's avatar
ulrich_y committed
304
      A=zeta(3) + zeta(2)*Log(x) - (Log(1._prec - X)*Log(X)**2)/2._prec + Log(X)**3/6._prec
Luca's avatar
Luca committed
305 306 307
    ELSE IF(X .LE. 2._prec) THEN
      YA=1._prec - X ; YB=(X-1._prec)/X
      S=1._prec
ulrich_y's avatar
ulrich_y committed
308
      A=zeta(3) + zeta(2)*Log(x) - (Log(X - 1._prec)*Log(X)**2)/2._prec + Log(X)**3/6._prec
Luca's avatar
Luca committed
309 310 311
    ELSE
      YA=0._prec ; YB=1._prec/X
      S=-1._prec
ulrich_y's avatar
ulrich_y committed
312
      A=2*zeta(2)*Log(x)-Log(x)**3/6._prec
Luca's avatar
Luca committed
313 314 315 316 317 318 319 320
    END IF


    HA=-2._prec*YA-1._prec ; HB= 2._prec*YB
    ALFAA=HA+HA ; ALFAB = HB+HB

    BA0 = 0. ; BA1=0. ; BA2=0.
    BB0 = 0. ; BB1=0. ; BB2=0.
ulrich_y's avatar
ulrich_y committed
321 322 323 324 325 326
    if (precision(1._prec) < 20) then
      maxi = 18
    else
      maxi = 42
    endif
    DO  I = maxi,0,-1
Luca's avatar
Luca committed
327 328 329 330 331
       BA0=CA(I)+ALFAA*BA1-BA2 ; BA2=BA1 ; BA1=BA0
       BB0=CB(I)+ALFAB*BB1-BB2 ; BB2=BB1 ; BB1=BB0
    ENDDO
    Li3 = A + S * (  (BA0 - HA*BA2) + (BB0 - HB*BB2) )
  END FUNCTION Li3
Luca's avatar
Luca committed
332 333

  FUNCTION trilog(x) result(res)
ulrich_y's avatar
ulrich_y committed
334
    ! evaluates trilog for any argument |x|<1
Luca's avatar
Luca committed
335 336
    complex(kind=prec) :: res
    complex(kind=prec) :: x
ulrich_y's avatar
ulrich_y committed
337
    if(abs(aimag(x)) < zero ) then
ulrich_y's avatar
ulrich_y committed
338
      res = Li3(real(x,kind=prec))
Luca's avatar
Luca committed
339
    else
ulrich_y's avatar
ulrich_y committed
340 341
      res = naive_polylog(3,x)
    endif
Luca's avatar
Luca committed
342 343
  END FUNCTION trilog

ulrich_y's avatar
ulrich_y committed
344
  FUNCTION BERNOULLI_POLYNOMIAL(n, x) result(res)
ulrich_y's avatar
ulrich_y committed
345
    integer, parameter :: maxn = 15
ulrich_y's avatar
ulrich_y committed
346
    integer n
ulrich_y's avatar
ulrich_y committed
347 348
    complex(kind=prec) :: x, res
    complex(kind=prec) :: xpow(maxn+1)
ulrich_y's avatar
ulrich_y committed
349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383
    integer, parameter :: coeffN(maxn+1, maxn) = reshape((/ &
        -   1, +   1,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0, &
        +   1, -   1, +   1,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0, &
            0, +   1, -   3, +   1,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0, &
        -   1,     0, +   1, -   2, +   1,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0, &
            0, -   1,     0, +   5, -   5, +   1,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0, &
        +   1,     0, -   1,     0, +   5, -   3, +   1,     0,     0,     0,     0,     0,     0,     0,     0,     0, &
            0, +   1,     0, -   7,     0, +   7, -   7, +   1,     0,     0,     0,     0,     0,     0,     0,     0, &
        -   1,     0, +   2,     0, -   7,     0, +  14, -   4, +   1,     0,     0,     0,     0,     0,     0,     0, &
            0, -   3,     0, +   2,     0, -  21,     0, +   6, -   9, +   1,     0,     0,     0,     0,     0,     0, &
        +   5,     0, -   3,     0, +   5,     0, -   7,     0, +  15, -   5, +   1,     0,     0,     0,     0,     0, &
            0, +   5,     0, -  11,     0, +  11,     0, -  11,     0, +  55, -  11, +   1,     0,     0,     0,     0, &
        - 691,     0, +   5,     0, -  33,     0, +  22,     0, -  33,     0, +  11, -   6, +   1,     0,     0,     0, &
            0, - 691,     0, +  65,     0, - 429,     0, + 286,     0, - 143,     0, +  13, -  13, +   1,     0,     0, &
        +   7,     0, - 691,     0, + 455,     0, -1001,     0, + 143,     0, -1001,     0, +  91, -   7, +   1,     0, &
            0, +  35,     0, - 691,     0, + 455,     0, - 429,     0, + 715,     0, -  91,     0, +  35, -  15, +   1 /), &
            (/maxn+1, maxn/))
    integer, parameter :: coeffD(maxn+1, maxn) = reshape((/ &
        +   2, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +   6, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +   1, +   2, +   2, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +  30, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +   1, +   6, +   1, +   3, +   2, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +  42, +   1, +   2, +   1, +   2, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +   1, +   6, +   1, +   6, +   1, +   2, +   2, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +  30, +   1, +   3, +   1, +   3, +   1, +   3, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +   1, +  10, +   1, +   1, +   1, +   5, +   1, +   1, +   2, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +  66, +   1, +   2, +   1, +   1, +   1, +   1, +   1, +   2, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +   1, +   6, +   1, +   2, +   1, +   1, +   1, +   1, +   1, +   6, +   2, +   1, +   1, +   1, +   1, +   1, &
        +2730, +   1, +   1, +   1, +   2, +   1, +   1, +   1, +   2, +   1, +   1, +   1, +   1, +   1, +   1, +   1, &
        +   1, + 210, +   1, +   3, +   1, +  10, +   1, +   7, +   1, +   6, +   1, +   1, +   2, +   1, +   1, +   1, &
        +   6, +   1, +  30, +   1, +   6, +   1, +  10, +   1, +   2, +   1, +  30, +   1, +   6, +   1, +   1, +   1, &
        +   1, +   2, +   1, +   6, +   1, +   2, +   1, +   2, +   1, +   6, +   1, +   2, +   1, +   2, +   2, +   1 /), &
            (/maxn+1, maxn/))

ulrich_y's avatar
ulrich_y committed
384
    real(kind=prec), parameter :: coeff(maxn+1,maxn) = coeffN/real(coeffD,kind=prec)
ulrich_y's avatar
ulrich_y committed
385 386 387 388 389 390
    integer i

    if (n>maxn) then
      print*,"Bernoulli beyond 15 is not implemented"
      stop
    endif
ulrich_y's avatar
ulrich_y committed
391

ulrich_y's avatar
ulrich_y committed
392 393
    xpow(1:n+1) = (/ ( x**i, i = 0, n ) /)
    res = sum( xpow(1:n+1) * coeff(1:n+1,n) )
ulrich_y's avatar
ulrich_y committed
394 395 396

  END FUNCTION

ulrich_y's avatar
ulrich_y committed
397
  RECURSIVE FUNCTION polylog1(m,x) result(res)
Luca Naterop's avatar
Luca Naterop committed
398 399
    ! computes the polylog
    
400
    integer :: m
ulrich_y's avatar
ulrich_y committed
401
    type(inum) :: x, inv
402
    complex(kind=prec) :: res
ulrich_y's avatar
ulrich_y committed
403 404
    integer i

405
    
ulrich_y's avatar
ulrich_y committed
406
#ifdef DEBUG
407
    if(verb >= 70) print*, 'called polylog(',m,',',x%c,x%i0,')'
ulrich_y's avatar
ulrich_y committed
408 409 410 411 412 413 414 415 416 417
#endif
#ifndef NOCACHE
    if (m.le.5) then
      do i=1,plcachesize(m)
        if( abs(cache(m,i)%c%c-x%c).lt.zero ) then
          res = cache(m,i)%ans
          return
        endif
      enddo
    endif
ulrich_y's avatar
ulrich_y committed
418
#endif
419
    if ((m.le.9).and.(abs(x%c-1.).lt.zero)) then
ulrich_y's avatar
ulrich_y committed
420
      res = zeta(m)
ulrich_y's avatar
ulrich_y committed
421 422
    else if ((m.le.9).and.(abs(x%c+1._prec).lt.zero)) then
      res = -(1._prec - 2._prec**(1-m))*zeta(m)
423
    else if (abs(x) .gt. 1) then
ulrich_y's avatar
ulrich_y committed
424
      inv = inum(1._prec/x%c, x%i0)
ulrich_y's avatar
ulrich_y committed
425
      res = (-1)**(m-1)*polylog(m,inv) &
ulrich_y's avatar
ulrich_y committed
426
          - (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*conjg(log(-x%c))/2._prec/pi) / factorial(m)
ulrich_y's avatar
ulrich_y committed
427
    else if(m == 2) then
428
      res = dilog(x%c)
Luca's avatar
Luca committed
429
    else if(m == 3) then
430
      res = trilog(x%c)
Luca's avatar
Luca committed
431
    else
432
      res = naive_polylog(m,x%c)
433
    end if
ulrich_y's avatar
ulrich_y committed
434 435 436 437 438 439 440 441 442

#ifndef NOCACHE
    if (m.le.PolyLogCacheSize(1)) then
      if (plcachesize(m).lt.PolyLogCacheSize(2)) then
        plcachesize(m) = plcachesize(m) + 1
        cache(m,plcachesize(m)) = el(x,res)
      endif
    endif
#endif
ulrich_y's avatar
ulrich_y committed
443 444 445 446 447 448 449 450 451
  END FUNCTION polylog1




  RECURSIVE FUNCTION polylog2(m,x,y) result(res)
    type(inum) :: x, y
    integer m
    complex(kind=prec) :: res
ulrich_y's avatar
ulrich_y committed
452
    !TODO!!
ulrich_y's avatar
ulrich_y committed
453 454 455 456 457 458 459 460 461
    res=polylog1(m,inum(x%c/y%c,di0))
  END FUNCTION POLYLOG2


  FUNCTION PLOG1(a,b)
  ! calculates log(1-a/b)
  implicit none
  type(inum) :: a,b
  complex(kind=prec) plog1
462 463 464 465 466 467 468 469 470 471 472 473 474

  if ( (abs(aimag(a)).lt.zero).and.(abs(aimag(b)).lt.zero) ) then
    ! Both arguments are real, only here does the ieps matter
    plog1 = log(abs(1.-a%c/b%c))
    ! this does not depend on the sign of a
    if (real(a).gt.real(b) .and. real(b).gt. 0) then
      plog1 = plog1 + b%i0*i_*pi
    elseif(real(a).lt.real(b) .and. real(b).lt. 0) then
      plog1 = plog1 - b%i0*i_*pi
    endif
  else
    plog1 = log(1.-a%c/b%c)
  endif
ulrich_y's avatar
ulrich_y committed
475
  END FUNCTION
Luca's avatar
Luca committed
476

ulrich_y's avatar
ulrich_y committed
477 478 479 480 481 482
#ifndef NOCACHE
  SUBROUTINE CLEARCACHE
  plcachesize=0
  END SUBROUTINE
#endif

483 484
END MODULE maths_functions

Luca's avatar
Luca committed
485 486 487 488 489 490 491
! PROGRAM test
!   use maths_functions
!   implicit none
!   complex(kind=prec) :: res
!   res = Li3(0.4d0)
!   print*, res
! END PROGRAM test
492