maths_functions.f90 13.3 KB
Newer Older
1 2 3 4 5

MODULE maths_functions
  use globals
  use utils
  implicit none
ulrich_y's avatar
ulrich_y committed
6 7 8
  interface polylog
    module procedure polylog1, polylog2
  end interface polylog
9

ulrich_y's avatar
ulrich_y committed
10
#ifndef HAVE_MPFUN
ulrich_y's avatar
ulrich_y committed
11
  real(kind=prec), parameter :: zeta(2:10) = (/1.6449340668482262, 1.2020569031595942, 1.0823232337111381, &
ulrich_y's avatar
ulrich_y committed
12 13
               1.03692775514337, 1.0173430619844488, 1.008349277381923, &
               1.0040773561979441, 1.0020083928260821, 1.000994575127818/)
ulrich_y's avatar
ulrich_y committed
14
#endif
ulrich_y's avatar
ulrich_y committed
15

ulrich_y's avatar
ulrich_y committed
16 17
  type el
    type(inum) :: c
ulrich_y's avatar
ulrich_y committed
18
    restype ans
ulrich_y's avatar
ulrich_y committed
19 20 21 22
  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

ulrich_y's avatar
ulrich_y committed
26
  FUNCTION naive_polylog(m,x) result(res)
27 28
    ! Computes the classical polylogarithm Li_m(x) using series representation up to order n
    integer :: m
ulrich_y's avatar
ulrich_y committed
29
    restype :: x, res
30
    integer :: i
ulrich_y's avatar
ulrich_y committed
31 32 33 34 35 36 37 38 39 40
#ifdef HAVE_MPFUN
    restype :: ii
    ii = nil
    res=0.
    do i=1,PolylogInfinity
      ii = ii + 1.
      if(abs(x**ii).lt.1.e-250) return !TODO
      res = res+x**ii/ii**m
    enddo
#else
41
    res=0.
42
    do i=1,PolylogInfinity
43 44 45 46
      if(i**m.lt.0) return ! roll over
      if(abs(x**i).lt.1.e-250) return
      res = res+x**i/i**m
    enddo
ulrich_y's avatar
ulrich_y committed
47
#endif
48
  END FUNCTION naive_polylog
Luca's avatar
Luca committed
49

ulrich_y's avatar
ulrich_y committed
50 51

#ifndef HAVE_MPFUN
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 94 95
  FUNCTION Li2(x)

   !! Dilogarithm for arguments x < = 1.0

   real (kind=prec):: X,Y,T,S,A,PI3,PI6,ZERO,ONE,HALF,MALF,MONE,MTWO
   real (kind=prec):: C(0:18),H,ALFA,B0,B1,B2,LI2_OLD
   real (kind=prec):: Li2
   integer :: i

   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/
   DATA PI3 /3.289868133696453_prec/, PI6 /1.644934066848226_prec/

   DATA C( 0) / 0.4299669356081370_prec/
   DATA C( 1) / 0.4097598753307711_prec/
   DATA C( 2) /-0.0185884366501460_prec/
   DATA C( 3) / 0.0014575108406227_prec/
   DATA C( 4) /-0.0001430418444234_prec/
   DATA C( 5) / 0.0000158841554188_prec/
   DATA C( 6) /-0.0000019078495939_prec/
   DATA C( 7) / 0.0000002419518085_prec/
   DATA C( 8) /-0.0000000319334127_prec/
   DATA C( 9) / 0.0000000043454506_prec/
   DATA C(10) /-0.0000000006057848_prec/
   DATA C(11) / 0.0000000000861210_prec/
   DATA C(12) /-0.0000000000124433_prec/
   DATA C(13) / 0.0000000000018226_prec/
   DATA C(14) /-0.0000000000002701_prec/
   DATA C(15) / 0.0000000000000404_prec/
   DATA C(16) /-0.0000000000000061_prec/
   DATA C(17) / 0.0000000000000009_prec/
   DATA C(18) /-0.0000000000000001_prec/

   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
    LI2_OLD=PI6
    Li2 = Real(LI2_OLD,prec)
    RETURN
96
   ELSE IF(abs(x-MONE) < zero) THEN
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
    LI2_OLD=MALF*PI6
    RETURN
   END IF
   T=-X
   IF(T .LE. MTWO) THEN
    Y=MONE/(ONE+T)
    S=ONE
    A=-PI3+HALF*(LOG(-T)**2-LOG(ONE+ONE/T)**2)
   ELSE IF(T .LT. MONE) THEN
    Y=MONE-T
    S=MONE
    A=LOG(-T)
    A=-PI6+A*(A+LOG(ONE+ONE/T))
   ELSE IF(T .LE. MALF) THEN
    Y=(MONE-T)/T
    S=ONE
    A=LOG(-T)
    A=-PI6+A*(MALF*A+LOG(ONE+T))
   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
    A=PI6+HALF*LOG(T)**2
   END IF

   H=Y+Y-ONE
   ALFA=H+H
   B1=ZERO
   B2=ZERO
   DO  I = 18,0,-1
     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
142

143
  RECURSIVE FUNCTION dilog(x) result(res)
ulrich_y's avatar
ulrich_y committed
144
    ! evaluates dilog for any argument |x|<1
145 146
    complex(kind=prec) :: res
    complex(kind=prec) :: x
Luca's avatar
Luca committed
147

ulrich_y's avatar
ulrich_y committed
148 149
    if(abs(aimag(x)) < zero ) then
      res = Li2(real(x))
150
    else
ulrich_y's avatar
ulrich_y committed
151 152
      res = naive_polylog(2,x)
    endif
153 154
  END FUNCTION dilog

Luca's avatar
Luca committed
155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 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
  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
    real (kind=prec):: CA(0:18),HA,ALFAA,BA0,BA1,BA2, YA
    real (kind=prec):: CB(0:18),HB,ALFAB,BB0,BB1,BB2, YB
    DATA CA(0) / 0.4617293928601208/
    DATA CA(1) / 0.4501739958855029/
    DATA CA(2) / -0.010912841952292843/
    DATA CA(3) / 0.0005932454712725702/
    DATA CA(4) / -0.00004479593219266303/
    DATA CA(5) / 4.051545785869334e-6/
    DATA CA(6) / -4.1095398602619446e-7/
    DATA CA(7) / 4.513178777974119e-8/
    DATA CA(8) / -5.254661564861129e-9/
    DATA CA(9) / 6.398255691618666e-10/
    DATA CA(10) / -8.071938105510391e-11/
    DATA CA(11) / 1.0480864927082917e-11/
    DATA CA(12) / -1.3936328400075057e-12/
    DATA CA(13) / 1.8919788723690422e-13/
    DATA CA(14) / -2.6097139622039465e-14/
    DATA CA(15) / 3.774985548158685e-15/
    DATA CA(16) / -5.671361978114946e-16/
    DATA CA(17) / 1.1023848202712794e-16/
    DATA CA(18) / -5.0940525990875006e-17/
    DATA CB(0) / -0.016016180449195803/
    DATA CB(1) / -0.5036424400753012/
    DATA CB(2) / -0.016150992430500253/
    DATA CB(3) / -0.0012440242104245127/
    DATA CB(4) / -0.00013757218124463538/
    DATA CB(5) / -0.000018563818526041144/
    DATA CB(6) / -2.841735345177361e-6/
    DATA CB(7) / -4.7459967908588557e-7/
    DATA CB(8) / -8.448038544563037e-8/
    DATA CB(9) / -1.5787671270014e-8/
    DATA CB(10) / -3.0657620579122164e-9/
    DATA CB(11) / -6.140791949281482e-10/
    DATA CB(12) / -1.2618831590198e-10/
    DATA CB(13) / -2.64931268635803e-11/
    DATA CB(14) / -5.664711482422879e-12/
    DATA CB(15) / -1.2303909436235178e-12/
    DATA CB(16) / -2.7089360852246495e-13/
    DATA CB(17) / -6.024075373994343e-14/
    DATA CB(18) / -1.2894320641440237e-14/
    real (kind=prec):: Li3
    integer :: i


    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
212
      LI3=zeta(3)
Luca's avatar
Luca committed
213
    RETURN
214
    ELSE IF( abs(x+1) < zero) THEN
ulrich_y's avatar
ulrich_y committed
215
      LI3=-0.75_prec*zeta(3)
Luca's avatar
Luca committed
216 217 218 219 220
    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
221
      A=-LOG(-X)*(zeta(2)+LOG(-x)**2/6._prec)
Luca's avatar
Luca committed
222 223 224 225 226 227 228 229 230 231 232
    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
233
      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
234 235 236
    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
237
      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
238 239 240
    ELSE
      YA=0._prec ; YB=1._prec/X
      S=-1._prec
ulrich_y's avatar
ulrich_y committed
241
      A=2*zeta(2)*Log(x)-Log(x)**3/6._prec
Luca's avatar
Luca committed
242 243 244 245 246 247 248 249 250 251 252 253 254 255
    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.
    DO  I = 18,0,-1
       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
256 257

  FUNCTION trilog(x) result(res)
ulrich_y's avatar
ulrich_y committed
258
    ! evaluates trilog for any argument |x|<1
Luca's avatar
Luca committed
259 260
    complex(kind=prec) :: res
    complex(kind=prec) :: x
ulrich_y's avatar
ulrich_y committed
261 262
    if(abs(aimag(x)) < zero ) then
      res = Li3(real(x))
Luca's avatar
Luca committed
263
    else
ulrich_y's avatar
ulrich_y committed
264 265
      res = naive_polylog(3,x)
    endif
Luca's avatar
Luca committed
266 267
  END FUNCTION trilog

ulrich_y's avatar
ulrich_y committed
268
  FUNCTION BERNOULLI_POLYNOMIAL(n, x) result(res)
ulrich_y's avatar
ulrich_y committed
269
    integer, parameter :: maxn = 15
ulrich_y's avatar
ulrich_y committed
270
    integer n
ulrich_y's avatar
ulrich_y committed
271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314
    restype :: x, res
    restype :: xpow(maxn+1)
    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/))

    real(kind=prec), parameter :: coeff(maxn+1,maxn) = coeffN/real(coeffD)
    integer i

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

ulrich_y's avatar
ulrich_y committed
316 317
    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
318 319 320

  END FUNCTION

ulrich_y's avatar
ulrich_y committed
321 322
#endif

ulrich_y's avatar
ulrich_y committed
323
  RECURSIVE FUNCTION polylog1(m,x) result(res)
Luca Naterop's avatar
Luca Naterop committed
324 325
    ! computes the polylog
    
326
    integer :: m
ulrich_y's avatar
ulrich_y committed
327
    type(inum) :: x, inv
ulrich_y's avatar
ulrich_y committed
328
    restype :: res
ulrich_y's avatar
ulrich_y committed
329 330
    integer i

331
    
ulrich_y's avatar
ulrich_y committed
332
#ifdef DEBUG
333
    if(verb >= 70) print*, 'called polylog(',m,',',x%c,x%i0,')'
ulrich_y's avatar
ulrich_y committed
334 335 336 337 338 339 340 341 342 343
#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
344
#endif
ulrich_y's avatar
ulrich_y committed
345 346

#ifndef HAVE_MPFUN
347
    if ((m.le.9).and.(abs(x%c-1.).lt.zero)) then
ulrich_y's avatar
ulrich_y committed
348
      res = zeta(m)
349
    else if ((m.le.9).and.(abs(x%c+1.).lt.zero)) then
ulrich_y's avatar
ulrich_y committed
350
      res = -(1. - 2.**(1-m))*zeta(m)
351
    else if (abs(x) .gt. 1) then
ulrich_y's avatar
ulrich_y committed
352 353 354
      inv = inum(1./x%c, x%i0)
      res = (-1)**(m-1)*polylog(m,inv) &
          - cmplx(0,2*pi)**m * bernoulli_polynomial(m, 0.5-cmplx(0.,1.)*conjg(log(-x%c))/2/pi) / factorial(m)
ulrich_y's avatar
ulrich_y committed
355
    else if(m == 2) then
356
      res = dilog(x%c)
Luca's avatar
Luca committed
357
    else if(m == 3) then
358
      res = trilog(x%c)
Luca's avatar
Luca committed
359
    else
360
      res = naive_polylog(m,x%c)
361
    end if
ulrich_y's avatar
ulrich_y committed
362
#endif
ulrich_y's avatar
ulrich_y committed
363 364 365 366 367 368 369 370 371

#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
372 373 374 375 376 377 378 379
  END FUNCTION polylog1




  RECURSIVE FUNCTION polylog2(m,x,y) result(res)
    type(inum) :: x, y
    integer m
ulrich_y's avatar
ulrich_y committed
380
    restype :: res
ulrich_y's avatar
ulrich_y committed
381
    !TODO!!
ulrich_y's avatar
ulrich_y committed
382 383 384 385 386 387 388 389
    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
ulrich_y's avatar
ulrich_y committed
390
  restype :: plog1
ulrich_y's avatar
ulrich_y committed
391
  !TODO!!
ulrich_y's avatar
ulrich_y committed
392 393 394
#ifdef HAVE_MPFUN
  plog1 = mpcmplx(log(1.-a%c/b%c))
#else
ulrich_y's avatar
ulrich_y committed
395
  plog1 = log(1.-a%c/b%c)
ulrich_y's avatar
ulrich_y committed
396
#endif
ulrich_y's avatar
ulrich_y committed
397
  END FUNCTION
Luca's avatar
Luca committed
398

ulrich_y's avatar
ulrich_y committed
399 400 401 402 403 404
#ifndef NOCACHE
  SUBROUTINE CLEARCACHE
  plcachesize=0
  END SUBROUTINE
#endif

405 406
END MODULE maths_functions

Luca's avatar
Luca committed
407 408 409 410 411 412 413
! PROGRAM test
!   use maths_functions
!   implicit none
!   complex(kind=prec) :: res
!   res = Li3(0.4d0)
!   print*, res
! END PROGRAM test
414