maths_functions.f90 21.6 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 /)
16 17 18 19 20 21

  real(kind=prec), parameter :: DirichletBeta(2:10) = (/ 0.91596559417721901505460351493238411077414937428167_prec, &
    0.96894614625936938048363484584691860006954026768391_prec, 0.98894455174110533610842263322837782131586088706273_prec, &
    0.99615782807708806400631936863097528151139552938826_prec, 0.99868522221843813544160078786020654967836454612651_prec, &
    0.99955450789053990949634654989905898300218848194998_prec, 0.99984999024682965633806705924046378147600743300743_prec, &
    0.99994968418722008982135887329384752737274799691796_prec, 0.99998316402619687740554072995833414145685781649717_prec /)
ulrich_y's avatar
ulrich_y committed
22
  type el
23
    complex(kind=prec) :: c
ulrich_y's avatar
ulrich_y committed
24 25 26 27 28
    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
29
CONTAINS
30

ulrich_y's avatar
ulrich_y committed
31
  FUNCTION naive_polylog(m,x) result(res)
32 33
    ! Computes the classical polylogarithm Li_m(x) using series representation up to order n
    integer :: m
34
    complex(kind=prec) :: x, res, del
35
    integer(kind=ikin) :: i
ulrich_y's avatar
ulrich_y committed
36
    res=0._prec
37 38 39 40

    i = 1
    del = 1._prec
    do while (abs(del) > zero)
41
      if(i**m.lt.0) return ! roll over
ulrich_y's avatar
ulrich_y committed
42
      if(abs(x**i).lt.1.e-250_prec) return
43
      del = x**i/i**m
ulrich_y's avatar
ulrich_y committed
44
      if (abs(del) < LiDelta) return
45 46 47
      res = res+del
      i = i+1
    end do
48
  END FUNCTION naive_polylog
Luca's avatar
Luca committed
49

50 51 52 53
  FUNCTION Li2(x)

   !! Dilogarithm for arguments x < = 1.0

ulrich_y's avatar
ulrich_y committed
54
   real (kind=prec):: X,Y,T,S,A,ZERO,ONE,HALF,MALF,MONE,MTWO
ulrich_y's avatar
ulrich_y committed
55
   real (kind=prec):: C(0:42),H,ALFA,B0,B1,B2,LI2_OLD
56
   real (kind=prec):: Li2
ulrich_y's avatar
ulrich_y committed
57
   integer :: i, maxi
58 59 60 61 62

   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
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 94 95 96 97 98 99 100 101 102 103 104 105
   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/
106 107 108 109 110 111 112 113

   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
114
    LI2_OLD=zeta(2)
115 116
    Li2 = Real(LI2_OLD,prec)
    RETURN
117
   ELSE IF(abs(x-MONE) < zero) THEN
ulrich_y's avatar
ulrich_y committed
118
    LI2_OLD=MALF*zeta(2)
119 120 121 122 123 124
    RETURN
   END IF
   T=-X
   IF(T .LE. MTWO) THEN
    Y=MONE/(ONE+T)
    S=ONE
ulrich_y's avatar
ulrich_y committed
125
    A=-2*zeta(2)+HALF*(LOG(-T)**2-LOG(ONE+ONE/T)**2)
126 127 128 129
   ELSE IF(T .LT. MONE) THEN
    Y=MONE-T
    S=MONE
    A=LOG(-T)
ulrich_y's avatar
ulrich_y committed
130
    A=-zeta(2)+A*(A+LOG(ONE+ONE/T))
131 132 133 134
   ELSE IF(T .LE. MALF) THEN
    Y=(MONE-T)/T
    S=ONE
    A=LOG(-T)
ulrich_y's avatar
ulrich_y committed
135
    A=-zeta(2)+A*(MALF*A+LOG(ONE+T))
136 137 138 139 140 141 142 143 144 145 146
   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
147
    A=zeta(2)+HALF*LOG(T)**2
148 149 150 151 152 153
   END IF

   H=Y+Y-ONE
   ALFA=H+H
   B1=ZERO
   B2=ZERO
ulrich_y's avatar
ulrich_y committed
154 155 156 157 158 159
   if (precision(1._prec) < 20) then
     maxi = 18
   else
     maxi = 42
   endif
   DO  I = maxi,0,-1
160 161 162 163 164 165 166 167
     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
168

169
  RECURSIVE FUNCTION dilog(x) result(res)
ulrich_y's avatar
ulrich_y committed
170
    ! evaluates dilog for any argument |x|<1
171 172
    complex(kind=prec) :: res
    complex(kind=prec) :: x
Luca's avatar
Luca committed
173

ulrich_y's avatar
ulrich_y committed
174
    if(abs(aimag(x)) < zero ) then
ulrich_y's avatar
ulrich_y committed
175
      res = Li2(real(x,kind=prec))
176
    else
ulrich_y's avatar
ulrich_y committed
177 178
      res = naive_polylog(2,x)
    endif
179 180
  END FUNCTION dilog

Luca's avatar
Luca committed
181 182 183 184 185 186
  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
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 272 273 274 275 276 277 278 279 280 281 282 283
    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
284
    real (kind=prec):: Li3
ulrich_y's avatar
ulrich_y committed
285
    integer :: i, maxi
Luca's avatar
Luca committed
286 287 288 289 290 291 292 293 294


    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
295
      LI3=zeta(3)
Luca's avatar
Luca committed
296
    RETURN
297
    ELSE IF( abs(x+1) < zero) THEN
ulrich_y's avatar
ulrich_y committed
298
      LI3=-0.75_prec*zeta(3)
Luca's avatar
Luca committed
299 300 301 302 303
    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
304
      A=-LOG(-X)*(zeta(2)+LOG(-x)**2/6._prec)
Luca's avatar
Luca committed
305 306 307 308 309 310 311 312 313 314 315
    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
316
      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
317 318 319
    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
320
      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
321 322 323
    ELSE
      YA=0._prec ; YB=1._prec/X
      S=-1._prec
ulrich_y's avatar
ulrich_y committed
324
      A=2*zeta(2)*Log(x)-Log(x)**3/6._prec
Luca's avatar
Luca committed
325 326 327 328 329 330 331 332
    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
333 334 335 336 337 338
    if (precision(1._prec) < 20) then
      maxi = 18
    else
      maxi = 42
    endif
    DO  I = maxi,0,-1
Luca's avatar
Luca committed
339 340 341 342 343
       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
344 345

  FUNCTION trilog(x) result(res)
ulrich_y's avatar
ulrich_y committed
346
    ! evaluates trilog for any argument |x|<1
Luca's avatar
Luca committed
347 348
    complex(kind=prec) :: res
    complex(kind=prec) :: x
ulrich_y's avatar
ulrich_y committed
349
    if(abs(aimag(x)) < zero ) then
ulrich_y's avatar
ulrich_y committed
350
      res = Li3(real(x,kind=prec))
Luca's avatar
Luca committed
351
    else
ulrich_y's avatar
ulrich_y committed
352 353
      res = naive_polylog(3,x)
    endif
Luca's avatar
Luca committed
354 355
  END FUNCTION trilog

ulrich_y's avatar
ulrich_y committed
356
  FUNCTION BERNOULLI_POLYNOMIAL(n, x) result(res)
ulrich_y's avatar
ulrich_y committed
357
    integer, parameter :: maxn = 15
ulrich_y's avatar
ulrich_y committed
358
    integer n
ulrich_y's avatar
ulrich_y committed
359 360
    complex(kind=prec) :: x, res
    complex(kind=prec) :: xpow(maxn+1)
ulrich_y's avatar
ulrich_y committed
361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395
    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
396
    real(kind=prec), parameter :: coeff(maxn+1,maxn) = coeffN/real(coeffD,kind=prec)
ulrich_y's avatar
ulrich_y committed
397 398 399 400 401 402
    integer i

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

ulrich_y's avatar
ulrich_y committed
404 405
    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
406 407 408

  END FUNCTION

ulrich_y's avatar
ulrich_y committed
409 410 411 412 413 414 415 416 417 418 419 420 421
  FUNCTION MYLOG(x)
    complex(kind=prec) :: x, mylog
    if (abs(aimag(x)) < zero) then
      if (real(x) > 0) then
        mylog = cmplx(log(real(+x,kind=prec)),     kind=prec)
      else
        mylog = cmplx(log(real(-x,kind=prec)), pi, kind=prec)
      endif
    else
      mylog = log(x)
    endif
  END FUNCTION

ulrich_y's avatar
ulrich_y committed
422
  RECURSIVE FUNCTION polylog1(m,x) result(res)
Luca Naterop's avatar
Luca Naterop committed
423 424
    ! computes the polylog
    
425
    integer :: m
426
    complex(kind=prec) :: x, inv
427
    complex(kind=prec) :: res
ulrich_y's avatar
ulrich_y committed
428 429
    integer i

430
    
ulrich_y's avatar
ulrich_y committed
431
#ifdef DEBUG
432
    if(verb >= 70) print*, 'called polylog(',m,',',x,')'
ulrich_y's avatar
ulrich_y committed
433 434 435 436
#endif
#ifndef NOCACHE
    if (m.le.5) then
      do i=1,plcachesize(m)
437
        if( abs(cache(m,i)%c-x).lt.zero ) then
ulrich_y's avatar
ulrich_y committed
438 439 440 441 442
          res = cache(m,i)%ans
          return
        endif
      enddo
    endif
ulrich_y's avatar
ulrich_y committed
443
#endif
444
    if ((m.le.9).and.(abs(x-1.).lt.zero)) then
ulrich_y's avatar
ulrich_y committed
445
      res = zeta(m)
446
    else if ((m.le.9).and.(abs(x+1._prec).lt.zero)) then
ulrich_y's avatar
ulrich_y committed
447
      res = -(1._prec - 2._prec**(1-m))*zeta(m)
448 449 450 451
    else if ((m.le.9).and.(abs(x-i_).lt.zero)) then
      res = -(0.5_prec**m - 0.5_prec**(2*m-1)) * zeta(m) + i_*dirichletbeta(m)
    else if ((m.le.9).and.(abs(x+i_).lt.zero)) then
      res = -(0.5_prec**m - 0.5_prec**(2*m-1)) * zeta(m) - i_*dirichletbeta(m)
452
    else if (abs(x) .gt. 1) then
453
      inv = 1._prec/x
ulrich_y's avatar
ulrich_y committed
454
      res = (-1)**(m-1)*polylog(m,inv) &
ulrich_y's avatar
ulrich_y committed
455
          - (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*mylog(-x)/2._prec/pi) / factorial(m)
ulrich_y's avatar
ulrich_y committed
456
    else if(m == 2) then
457
      res = dilog(x)
Luca's avatar
Luca committed
458
    else if(m == 3) then
459
      res = trilog(x)
Luca's avatar
Luca committed
460
    else
461
      res = naive_polylog(m,x)
462
    end if
ulrich_y's avatar
ulrich_y committed
463 464 465 466 467 468 469 470 471

#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
472 473 474 475 476 477 478 479 480
  END FUNCTION polylog1




  RECURSIVE FUNCTION polylog2(m,x,y) result(res)
    type(inum) :: x, y
    integer m
    complex(kind=prec) :: res
481 482 483 484 485 486 487 488 489 490 491
    res=polylog1(m,x%c/y%c)
    if ( (abs(aimag(x)).lt.zero).and.(abs(aimag(y)).lt.zero) ) then
      ! Both arguments are real, only here does the ieps matter
      ! FIXME this is rather ugly..
      if (real(x).gt.real(y) .and. real(y).gt. 0) then
        ! Force the sign to be -b%i0
        res = cmplx(real(res), sign(aimag(res), -y%i0*1._prec), kind=prec)
      elseif(real(x).lt.real(y) .and. real(y).lt. 0) then
        res = cmplx(real(res), sign(aimag(res), +y%i0*1._prec), kind=prec)
      endif
    endif
ulrich_y's avatar
ulrich_y committed
492 493 494 495 496 497 498 499
  END FUNCTION POLYLOG2


  FUNCTION PLOG1(a,b)
  ! calculates log(1-a/b)
  implicit none
  type(inum) :: a,b
  complex(kind=prec) plog1
500 501 502 503 504 505 506 507 508 509 510

  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
ulrich_y's avatar
ulrich_y committed
511
    plog1 = mylog(1.-a%c/b%c)
512
  endif
ulrich_y's avatar
ulrich_y committed
513
  END FUNCTION
Luca's avatar
Luca committed
514

ulrich_y's avatar
ulrich_y committed
515 516 517 518 519 520
#ifndef NOCACHE
  SUBROUTINE CLEARCACHE
  plcachesize=0
  END SUBROUTINE
#endif

521 522
END MODULE maths_functions

Luca's avatar
Luca committed
523 524 525 526 527 528 529
! PROGRAM test
!   use maths_functions
!   implicit none
!   complex(kind=prec) :: res
!   res = Li3(0.4d0)
!   print*, res
! END PROGRAM test
530