/*
 * expint.cpp: exponential integrals.
 */

#include <stdio.h>
#include <stdlib.h>

#include "spigot.h"
#include "funcs.h"
#include "cr.h"
#include "error.h"
#include "holefiller.h"

/*
 * Exponential integrals, logarithmic integrals, and the
 * Euler-Mascheroni constant (which we get for free as a by-product of
 * that lot).
 *
 * There's a slightly confusing family of exponential integral
 * functions. They're all based around the idea of wanting to
 * integrate e^x/x, but because that integrand has a pole at 0, things
 * get messy.
 *
 * If you're prepared to work in the complex numbers, then integrands
 * with poles are a very well understood problem, which you deal with
 * by introducing a branch cut. A typical place to put the branch cut
 * is along the negative real axis, which is less than helpful if you
 * want a function that restricts sensibly to R (or rather, R\{0}).
 * Also, even if you put the branch cut somewhere more out of the way,
 * you find that if you want your function to take real values on R^+
 * then it ends up taking values on R^- of the form (real + pi i) or
 * (real - pi i), depending which way you swerved round the origin.
 *
 * This is all a bit unsatisfactory for people who are only interested
 * in R. So the usual answer is to define a function called Ei, which
 * acts as an indefinite integral of e^x/x on each side of zero, so
 * that a _definite_ integral over any interval not containing zero
 * can be given as the difference of two values of Ei. (If the
 * interval does contain zero, of course, you go back to the person
 * who asked for it and demand more detail about what exactly they
 * were after.) To define this function fully, we must choose an
 * additive constant for R^- and another one for R^+ (since the
 * disconnected domain permits a separate choice in each component).
 *
 * The standard definition of Ei, as implemented by octave-specfun
 * under the name expint_Ei and by Mathematica under the name
 * ExpIntegralEi, chooses the negative domain constant so that the
 * limit of Ei at -infinity is zero, and then chooses the positive
 * domain constant by 'taking the Cauchy principal value', which
 * basically means looking at the integral you'd get if you zeroed out
 * a small interval [-epsilon,+epsilon] about the pole at zero, and
 * then took the limit as epsilon -> 0. An equivalent way of putting
 * that is to say that Ei's negative and positive constants are chosen
 * so that lim_{e->0} (Ei(-e) - Ei(+e)) = 0.
 *
 * The other basic idea of an exponential integral takes two
 * parameters, a non-negative integer n and a positive real x.
 * MathWorld calls this the 'E_n function':
 * http://mathworld.wolfram.com/En-Function.html
 * and defines it as
 *
 *            inf exp(-xt)
 *   E_n(x) = int -------- dt
 *             1    t^n
 *
 * That's not really an illuminating definition as far as I'm
 * concerned. I prefer the following inductive definition: E_0 is just
 * e^{-x}/x, and then E_{n+1} is the indefinite integral of (-E_n),
 * with the constant chosen every time so that every E_n tends to zero
 * as x->inf. (The flipping signs are slightly annoying, but the point
 * of them seems to be to make all the E_n positive and decreasing,
 * rather than having the sign of E_n depend on the parity of n.)
 *
 * Gnuplot implements E_n(x) under the name expint(n,x), and
 * Mathematica under the name ExpIntegralE[n,x].
 *
 * We implement Ei(x), En(n,x), and E1(x) = En(1,x). We also implement
 * a function that Wikipedia defines under the name 'Ein':
 * https://en.wikipedia.org/wiki/Exponential_integral which is defined
 * as the integral of (1-e^{-x})/x (with constant so that Ein(0)=0).
 * The point of that function is that it _doesn't_ have a pole at
 * zero, or indeed anywhere else: it's an actually sensible function
 * from R->R, and even from C->C, with an obvious power series that
 * converges everywhere (obtained from the power series of exp by
 * term-by-term integration). So it's easy to evaluate, and hence
 * makes a useful primitive to build the various other things out of.
 *
 * Also in this family are the logarithmic integrals li(x) and Li(x),
 * which are easy because they're basically just Ei(log x). Wikipedia
 * [https://en.wikipedia.org/wiki/Logarithmic_integral_function]
 * defines them to be translations of each other, such that li(0)=0
 * but Li(2)=0. (So Li lets you integrate up to large numbers without
 * the slight arbitrariness of having to invent a way to handle
 * crossing the pole at 1.)
 */

/* ----------------------------------------------------------------------
 * Continued-fraction expansion of E_n.
 */

class EnCfrac : public Source {
    /*
     * http://functions.wolfram.com/GammaBetaErf/ExpIntegralE/10/
     * says that
     *
     *            exp(-x)
     *   E_n(x) = -------------------------
     *            x + n
     *                ---------------------
     *                1 + 1
     *                    -----------------
     *                    x + n + 1
     *                        -------------
     *                        1 + 2
     *                            ---------
     *                            x + n + 2
     *                                -----
     *                                ...
     *
     * This class computes everything but the exp(-x) factor, which
     * the caller must multiply in afterwards.
     */

    bigint xn, xd, n, nim1, i;
    int crState;

  public:
    EnCfrac(const bigint &axn, const bigint &axd, const bigint &an)
        : xn(axn), xd(axd), n(an)
    {
        crState = -1;
    }

    virtual EnCfrac *clone() { return new EnCfrac(xn, xd, n); }

    bool gen_interval(bigint *low, bigint *high)
    {
        *low = 0;
        *high = 0;
        return false;
    }

    bool gen_matrix(bigint *matrix)
    {
        crBegin;

        /*
         * An initial matrix representing x |-> 1/x.
         */
        matrix[0] = 0;
        matrix[1] = 1;
        matrix[2] = 1;
        matrix[3] = 0;
        crReturn(false);

        /*
         * Then the regular series.
         */
        i = 1;
        nim1 = n;                      // n + i - 1

        while (1) {
            matrix[0] = xn;
            matrix[1] = nim1*xd;
            matrix[2] = xd;
            matrix[3] = 0;
            crReturn(false);

            matrix[0] = 1;
            matrix[1] = i;
            matrix[2] = 1;
            matrix[3] = 0;
            crReturn(false);

            ++i;
            ++nim1;
        }

        crEnd;
    }
};

class EinSeries : public Source {
    /*
     * Ein is computed by the following power series:
     *
     *            inf (-1)^{k+1} x^k
     *   Ein(x) = sum --------------
     *            k=1       k k!
     *
     * Let x = n/d. Then we have
     *
     *                 n       -n^2       +n^3
     *   Ein(n/d) = ------ + -------- + -------- + ...
     *              1 1! d   2 2! d^2   3 3! d^3
     *
     *              n       -n   1   -n   1
     *            = - ( 1 + -- ( - + -- ( - + ... )))
     *              d       2d   2   3d   3
     *
     * so our matrices go (with an anomalous non-negation of n in the
     * first one)
     *
     *   ( n n ) ( -2n -n ) ( -3n -n  ) ...
     *   ( 0 d ) (   0 4d ) (   0 9d )
     */

    bigint n, d, k;
    int crState;

  public:
    EinSeries(const bigint &an, const bigint &ad)
        : n(an), d(ad)
    {
        crState = -1;
    }

    virtual EinSeries *clone() { return new EinSeries(n, d); }

    bool gen_interval(bigint *low, bigint *high)
    {
        *low = -1;
        *high = 1;
        return true;
    }

    bool gen_matrix(bigint *matrix)
    {
        crBegin;

        /* Anomalous initial matrix without n negated */
        matrix[0] = n;
        matrix[1] = n;
        matrix[2] = 0;
        matrix[3] = d;
        crReturn(false);

        k = 2;
        while (1) {
            matrix[0] = -k*n;
            matrix[1] = -n;
            matrix[2] = 0;
            matrix[3] = k*k*d;
            ++k;
            crReturn(false);
        }

        crEnd;
    }
};

struct EnCfracConstructor : MonotoneConstructor {
    bigint n;
    EnCfracConstructor(const bigint &an) : n(an) {}
    MonotoneConstructor *clone() { return new EnCfracConstructor(n); }
    Spigot *construct(const bigint &xn, const bigint &xd) {
        return new EnCfrac(xn, xd, n);
    }
};
struct EinSeriesConstructor : MonotoneConstructor {
    MonotoneConstructor *clone() { return new EinSeriesConstructor(); }
    Spigot *construct(const bigint &n, const bigint &d) {
        return new EinSeries(n, d);
    }
};

static Spigot *spigot_En_internal(const bigint &n, Spigot *x)
{
    Spigot *expintcfrac = spigot_monotone(new EnCfracConstructor(n),
                                          x->clone());
    return spigot_mul(expintcfrac, spigot_exp(spigot_neg(x)));
}

Spigot *spigot_Ein(Spigot *x)
{
    return spigot_monotone(new EinSeriesConstructor, x->clone());
}

/*
 * Wikipedia shows the relation between E_1 and Ein as
 *
 *   E_1(z) = -gamma - log(z) + Ein(z)
 *
 * where gamma is the Euler-Mascheroni constant. This means that using
 * our continued fraction for E_1 and our series for Ein, we can
 * _compute_ the Euler-Mascheroni constant!
 */

Spigot *spigot_eulergamma()
{
    /*
     * This formula works for any z, so it's a matter of choosing one
     * which gives the best convergence.
     *
     * We care about generating a lot of digits of gamma, of course;
     * but we also care about not taking too long to _get started_,
     * because we'll also be using this as a component of some of the
     * user-facing exponential integral functions below. Benchmarking
     * suggests that as you turn up z, you get a gradual improvement
     * in the generation of lots of digits, but it gets more and more
     * gradual and at the same time the startup time begins to become
     * significant - by the time you're at z=4096, say, you're taking
     * (on my test machine) 12 seconds before _any_ digit is
     * generated, and not managing to overtake the z=512 run until
     * about 1000 digits have gone by, and even then, not by very
     * much.
     *
     * I suppose if you wanted a _really_ large number of digits, it
     * _might_ be worth rigging up a wrapper class which kept
     * restarting the computation from the beginning with bigger and
     * bigger z, akin to the strategies used by GammaResidualBase and
     * MonotoneHelper, but I think the gain would be marginal even so.
     * Based on my ad-hoc benchmarking, the sensible thing seems to be
     * to pick the largest value of z that we get to before the
     * startup time starts to become noticeable, and that's 256.
     */
    Spigot *z = spigot_integer(256);
    Spigot *Ein_minus_En = spigot_sub(spigot_Ein(z->clone()),
                                      spigot_En_internal(1, z->clone()));
    return spigot_sub(Ein_minus_En, spigot_log(z));
}

class EnHoleFiller : public HoleFiller {
    bigint n;

    virtual EnHoleFiller *clone() { return new EnHoleFiller(n, x->clone()); }
    virtual Spigot *xspecial(int i) {
        if (i == 0 && n > 1) return spigot_integer(0);
        return NULL;
    }
    virtual Spigot *yspecial(int i) {
        if (i == 0 && n > 1) return spigot_rational(1, n-1);
        return NULL;
    }
    virtual Spigot *ygeneral(Spigot *x) {
        return spigot_En_internal(n, x);
    }
  public:
    EnHoleFiller(const bigint &an, Spigot *ax) : HoleFiller(ax), n(an) {}
};

Spigot *spigot_En(Spigot *n, Spigot *x)
{
    bigint nn, nd;
    if (!n->is_rational(&nn, &nd) || nd != 1)
        throw spigot_error("first argument to En must be an integer");

    if (nn < 0)
        throw spigot_error("first argument to En must be non-negative");

    if (nn > 1)
        x = spigot_enforce(x, ENFORCE_GE, spigot_integer(0),
                           spigot_error("second argument to En must be "
                                        "non-negative"));
    else
        x = spigot_enforce(x, ENFORCE_GT, spigot_integer(0),
                           spigot_error("second argument to En must be "
                                        "positive"));

    return (new EnHoleFiller(nn, x))->replace();
}

Spigot *spigot_E1(Spigot *x)
{
    x = spigot_enforce(x, ENFORCE_GT, spigot_integer(0),
                       spigot_error("argument to E1 must be positive"));
    return spigot_En_internal(1, x);
}

Spigot *spigot_Ei(Spigot *x)
{
    Spigot *Ein_term = spigot_Ein(spigot_neg(x->clone()));
    return spigot_sub(spigot_add(spigot_eulergamma(),
                                 spigot_log(spigot_abs(x))),
                      Ein_term);
}

// Unusually lowercase class name, because li and Li really are
// different things :-)
class liHoleFiller : public HoleFiller {
    virtual liHoleFiller *clone() { return new liHoleFiller(x->clone()); }
    virtual Spigot *xspecial(int i) {
        if (i == 0) return spigot_integer(0);
        return NULL;
    }
    virtual Spigot *yspecial(int i) {
        if (i == 0) return spigot_integer(0);
        return NULL;
    }
    virtual Spigot *ygeneral(Spigot *x) {
        return spigot_Ei(spigot_log(x));
    }
  public:
    liHoleFiller(Spigot *ax) : HoleFiller(ax) {}
};

Spigot *spigot_li(Spigot *x)
{
    x = spigot_enforce(x, ENFORCE_GE, spigot_integer(0),
                       spigot_error("argument to li must be non-negative"));
    return (new liHoleFiller(x))->replace();
}

class LiHoleFiller : public HoleFiller {
    virtual LiHoleFiller *clone() { return new LiHoleFiller(x->clone()); }
    virtual Spigot *xspecial(int i) {
        if (i == 0) return spigot_integer(0);
        if (i == 1) return spigot_integer(2);
        return NULL;
    }
    virtual Spigot *yspecial(int i) {
        if (i == 0) return spigot_neg(spigot_li(spigot_integer(2)));
        if (i == 1) return spigot_integer(0);
        return NULL;
    }
    virtual Spigot *ygeneral(Spigot *x) {
        /*
         * Li is just li(x)-li(2), or equivalently, the integral from 2 to
         * x of t/log(t). But we can do better than actually subtracting
         * two values of li, because li is implemented in terms of Ei
         * which is in turn the sum of multiple terms, so we can simplify:
         *
         *   Li(x) = li(x) - li(2)
         *         = Ei(log x) - Ei(log 2)
         *         = gamma + log |log x| - Ein(-log x)
         *            - gamma - log |log 2| + Ein(-log 2)
         *         = log |log x| - log |log 2| - Ein(-log x) + Ein(-log 2)
         *         = log (|log_2 x|) - Ein(-log x) + Ein(-log 2)
         *
         * and now we don't have to compute gamma at all, plus the two
         * separate log terms have turned into a single log-base-2.
         */
        Spigot *log_log2_x = spigot_log(spigot_abs(spigot_log2(x->clone())));
        Spigot *Ein_log_x = spigot_Ein(spigot_neg(spigot_log(x)));
        Spigot *Ein_log_2 = spigot_Ein(spigot_log(spigot_rational(1,2)));
        return spigot_add(log_log2_x, spigot_sub(Ein_log_2, Ein_log_x));
    }
  public:
    LiHoleFiller(Spigot *ax) : HoleFiller(ax) {}
};

Spigot *spigot_Li(Spigot *x)
{
    x = spigot_enforce(x, ENFORCE_GE, spigot_integer(0),
                       spigot_error("argument to Li must be non-negative"));
    return (new LiHoleFiller(x))->replace();
}