use std;
use num::complex::Complex;

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

/// Returns the dilogarithm of a real number of type `f64`.
///
/// This function has been translated from the
/// [ROOT](https://root.cern.ch/) package.
///
/// # Example:
/// ```
/// use polylog::Li2;
///
/// let z = 1.0;
/// println!("Li2({}) = {}", z, z.li2());
/// ```
impl Li2<f64> for f64 {
    fn li2(&self) -> f64 {
        let pi  = 3.1415926535897932384626433832795;
        let pi2 = pi*pi;
        let pi3 = pi2/3.;
        let pi6 = pi2/6.;
        let coeffs = [0.42996693560813697, 0.40975987533077105,
           -0.01858843665014592, 0.00145751084062268,-0.00014304184442340,
            0.00001588415541880,-0.00000190784959387, 0.00000024195180854,
           -0.00000003193341274, 0.00000000434545063,-0.00000000060578480,
            0.00000000008612098,-0.00000000001244332, 0.00000000000182256,
           -0.00000000000027007, 0.00000000000004042,-0.00000000000000610,
            0.00000000000000093,-0.00000000000000014, 0.00000000000000002];

        if *self == 1.0 {
            pi6
        } else if *self == -1.0 {
            -pi2/12.
        } else {
            let t = -*self;
            let (y, s, a) = if t <= -2.0 {
                let b1 = (-t).ln();
                let b2 = (1.0 + 1.0/t).ln();
                (-1.0/(1.0 + t), 1.0, -pi3 + 0.5*(b1*b1 - b2*b2))
            } else if t < -1.0 {
                let a = (-t).ln();
                (-1.0 - t, -1.0, -pi6 + a*(a + (1.0 + 1.0/t).ln()))
            } else if t <= -0.5 {
                let a = (-t).ln();
                (-(1.0 + t)/t, 1.0, -pi6 + a*(-0.5*a + (1.0 + t).ln()))
            } else if t < 0.0 {
                let b1 = (1.0 + t).ln();
                (-t/(1.0 + t), -1.0, 0.5*b1*b1)
            } else if t <= 1.0 {
                (t, 1.0, 0.)
            } else {
                let b1 = t.ln();
                (1.0/t, -1.0, pi6 + 0.5*b1*b1)
            };

            let h      = y+y - 1.0;
            let alfa   = h+h;
            let mut b0 = 0.0;
            let mut b1 = 0.0;
            let mut b2 = 0.0;
            for c in coeffs.iter().rev() {
                b0 = c + alfa*b1 - b2;
                b2 = b1;
                b1 = b0;
            }
            -(s*(b0 - h*b2) + a)
        }
    }
}

/// 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());
/// }
/// ```
impl Li2<Complex<f64>> for Complex<f64> {
    fn li2(&self) -> Complex<f64> {
        let pi = 3.1415926535897932384626433832795;

        // 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./36.0e2,
              1./21168.0e1,
            - 1./108864.0e2,
              1./52690176.0e1,
            - 4.0647616451442255268059093862919666745470571274397078e-11,
              8.9216910204564525552179873167527488515142836130490451e-13,
            - 1.9939295860721075687236443477937897056306947496538801e-14,
              4.5189800296199181916504765528555932283968190144666184e-16,
            - 1.0356517612181247014483411542218656665960912381686505e-17,
              2.3952186210261867457402837430009803816789490019429743e-19,
            - 5.5817858743250093362830745056254199055670546676443981e-21,
              1.3091507554183212858123073991865923017498498387833038e-22,
            - 3.0874198024267402932422797648664624315955652561327457e-24,
              7.315975652702203420357905609252148591033401063690875e-26,
            - 1.7408456572340007409890551477597025453408414217542713e-27,
              4.1576356446138997196178996207752266734882541595115639e-29,
            - 9.9621484882846221031940067024558388498548600173944888e-31,
              2.3940344248961653005211679878937495629342791569329158e-32,
        ];

        let rz = self.re;
        let iz = self.im;
        let az = self.norm();

        // 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())
            }
        } else if az < std::f64::EPSILON {
            return *self;
        }

        let (cy, cz, jsgn, ipi12) = if rz <= 0.5 {
            if az > 1. {
                (-0.5 * sqr((-self).ln()), -(1. - 1. / self).ln(), -1., -2.)
            } else { // az <= 1.
                (Complex::new(0.,0.), -(1. - self).ln(), 1., 0.)
            }
        } else { // rz > 0.5
            if az <= (2.0*rz).sqrt() {
                let l = -(self).ln();
                (l * (1. - self).ln(), l, -1., 2.)
            } else { // az > sqrt(2*rz)
                (-0.5 * sqr((-self).ln()), -(1. - 1. / self).ln(), -1., -2.)
            }
        };

        // the dilogarithm
        let cz2 = sqr(cz);
        let sum =
            cz +
            cz2 * (bf[0] +
            cz  * (bf[1] +
            cz2 * (bf[2] +
            cz2 * (bf[3] +
            cz2 * (bf[4] +
            cz2 * (bf[5] +
            cz2 * (bf[6] +
            cz2 * (bf[7] +
            cz2 * (bf[8] +
            cz2 * (bf[9] +
            cz2 * (bf[10] +
            cz2 * (bf[11] +
            cz2 * (bf[12] +
            cz2 * (bf[13] +
            cz2 * (bf[14] +
            cz2 * (bf[15] +
            cz2 * (bf[16] +
            cz2 * (bf[17] +
            cz2 * (bf[18] +
            cz2 * (bf[19]))))))))))))))))))));

        jsgn * sum + cy + ipi12 * pi * pi / 12.
    }
}

fn sqr(x: Complex<f64>) -> Complex<f64> { x*x }