Commit 3b2d170d authored by ulrich_y's avatar ulrich_y

Fixed bug related to branch cut in polylog

In bdc545ff we have introduced a conjg the
Bernoulli polynomial. It is not clear why.

Courtesy to Roman Zwicky and Ben Pullin for
pointing this out.
parent c624a50e
......@@ -423,7 +423,7 @@ CONTAINS
else if (abs(x) .gt. 1) then
inv = 1._prec/x
res = (-1)**(m-1)*polylog(m,inv) &
- (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*conjg(log(-x))/2._prec/pi) / factorial(m)
- (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*log(-x)/2._prec/pi) / factorial(m)
else if(m == 2) then
res = dilog(x)
else if(m == 3) then
......
  • It is now clear, why I did this. Fortran is sometimes inconsistent when it comes to certain logarithms

    complex(kind=prec) :: x
    real(kind=prec)
    x = cmplx(r, 0._prec)
    clog(x)

    However, IEEE754 supports signed zeroes (which in hind-sight would have been a fun idea instead of ieps). This is implemented assuming big-endian (>f)

    +---+------------------+----------------------------------------------------+
    | S |      Exponent    |                    Mantissa                        |
    +---+------------------+----------------------------------------------------+
    | 1 | 0 0 0 0  0 0 0 0 | 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 |  = -0.0
    +-+-----------------+--+---------------+-----------------+------------------+
    |        80         |        00        |        00       |        00        |
    +-------------------+------------------+-----------------+------------------+

    in contrast to +0.0 which is

    +---+------------------+----------------------------------------------------+
    | S |      Exponent    |                    Mantissa                        |
    +---+------------------+----------------------------------------------------+
    | 0 | 0 0 0 0  0 0 0 0 | 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 |  = +0.0
    +-+-----------------+--+---------------+-----------------+------------------+
    |        00         |        00        |        00       |        00        |
    +-------------------+------------------+-----------------+------------------+

    Fortran assigns meaning to the sign bit, even if the mantissa vanishes, i.e.

    clog( (-1._prec, -0._prec) ) = -pi
    clog( (-1._prec, +0._prec) ) = +pi

    When I implemented this for the first time, this is why I thought the conjg was necessary.

  • mentioned in commit 09592838

    Toggle commit list
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment