1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 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
 94
 95
 96
 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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172

use num::complex::Complex;
use cln::CLn;

/// Provides the dilogarithm function `li2()` of a number of type `T`.
pub trait Li2<T> {
    fn li2(&self) -> T;
}

impl Li2<f64> for f64 {
    /// Returns the real dilogarithm of a real number of type `f64`.
    ///
    /// Implemented as an economized Pade approximation with a maximum
    /// error of 4.16e-18.
    ///
    /// # Example:
    /// ```
    /// use polylog::Li2;
    ///
    /// let z = 1.0;
    /// println!("Li2({}) = {}", z, z.li2());
    /// ```
    fn li2(&self) -> f64 {
        let pi = std::f64::consts::PI;
        let cp = [
            1.0706105563309304277e+0,
           -4.5353562730201404017e+0,
            7.4819657596286408905e+0,
           -6.0516124315132409155e+0,
            2.4733515209909815443e+0,
           -4.6937565143754629578e-1,
            3.1608910440687221695e-2,
           -2.4630612614645039828e-4
        ];
        let cq = [
            1.0000000000000000000e+0,
           -4.5355682121856044935e+0,
            8.1790029773247428573e+0,
           -7.4634190853767468810e+0,
            3.6245392503925290187e+0,
           -8.9936784740041174897e-1,
            9.8554565816757007266e-2,
           -3.2116618742475189569e-3
        ];

        let x = *self;

        // transform to [0, 1/2)
        let (y, r, s) = if x < -1. {
            let l = (1. - x).ln();
            (1./(1. - x), -pi*pi/6. + l*(0.5*l - (-x).ln()), 1.)
        } else if x == -1. {
            return -pi*pi/12.;
        } else if x < 0. {
            let l = (-x).ln_1p();
            (x/(x - 1.), -0.5*l*l, -1.)
        } else if x == 0. {
            return 0.;
        } else if x < 0.5 {
            (x, 0., 1.)
        } else if x < 1. {
            (1. - x, pi*pi/6. - x.ln()*(1. - x).ln(), -1.)
        } else if x == 1. {
            return pi*pi/6.;
        } else if x < 2. {
            let l = x.ln();
            (1. - 1./x, pi*pi/6. - l*((1. - 1./x).ln() + 0.5*l), 1.)
        } else {
            let l = x.ln();
            (1./x, pi*pi/3. - 0.5*l*l, -1.)
        };

        let z = y - 0.25;
        let p = horner(z, &cp);
        let q = horner(z, &cq);

        r + s*y*p/q
    }
}

impl Li2<Complex<f64>> for Complex<f64> {
    /// Returns the dilogarithm of a complex number of type
    /// `Complex<f64>`.
    ///
    /// This function has been translated from the
    /// [SPheno](https://spheno.hepforge.org/) package.
    ///
    /// # Example:
    /// ```
    /// extern crate num;
    /// extern crate polylog;
    /// use num::complex::Complex;
    /// use polylog::Li2;
    ///
    /// fn main() {
    ///     let z = Complex::new(1.0, 1.0);
    ///     println!("Li2({}) = {}", z, z.li2());
    /// }
    /// ```
    fn li2(&self) -> Complex<f64> {
        let pi = std::f64::consts::PI;

        // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
        // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 19}]
        let bf = [
            - 1./4.,
              1./36.,
            - 1./3600.,
              1./211680.,
            - 1./10886400.,
              1./526901760.,
            - 4.0647616451442255e-11,
              8.9216910204564526e-13,
            - 1.9939295860721076e-14,
              4.5189800296199182e-16,
        ];

        let rz = self.re;
        let iz = self.im;

        // special cases
        if iz == 0. {
            if rz <= 1. {
                return Complex::new(rz.li2(), 0.0)
            } else { // rz > 1.
                return Complex::new(rz.li2(), -pi*rz.ln())
            }
        }

        let nz = self.norm_sqr();

        if nz < std::f64::EPSILON {
            return *self;
        }

        let (u, rest, sgn) = if rz <= 0.5 {
            if nz > 1. {
                let l = (-self).cln();
                (-(1. - 1. / self).cln(), -0.5 * l * l - pi * pi / 6., -1.)
            } else { // nz <= 1.
                (-(1. - self).cln(), Complex::new(0.,0.), 1.)
            }
        } else { // rz > 0.5
            if nz <= 2.0*rz {
                let l = -(self).cln();
                (l, l * (1. - self).cln() + pi * pi / 6., -1.)
            } else { // nz > 2.0*rz
                let l = (-self).cln();
                (-(1. - 1. / self).cln(), -0.5 * l * l - pi * pi / 6., -1.)
            }
        };

        let u2 = u*u;
        let u4 = u2*u2;
        let sum =
            u +
            u2 * (bf[0] +
            u  * (bf[1] +
            u2 * (
                bf[2] +
                u2*bf[3] +
                u4*(bf[4] + u2*bf[5]) +
                u4*u4*(bf[6] + u2*bf[7] + u4*(bf[8] + u2*bf[9]))
            )));

        sgn * sum + rest
    }
}

/// evaluation of polynomial P(x) with coefficients `coeffs`
fn horner(x: f64, coeffs: &[f64]) -> f64 {
    coeffs.iter().rev().fold(0., |p, c| p*x + c)
}