/*
 * trig.cpp: Trigonometric functions (sin, cos, tan).
 */

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

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

class SinRational : public Source {
    /*
     * We compute the sine of a rational by translating the obvious
     * power series for sin(x) into a spigot description.
     *
     * Let x = n/d. Then we have
     *
     *              n     n^3      n^5
     *   sin(n/d) = - - ------ + ------ - ...
     *              d   3! d^5   5! d^5
     *
     *              n           n^2           n^2
     *            = - ( 1 - ------- ( 1 - ------- ( ... ) ) )
     *              d       2.3.d^2       4.5.d^2
     *
     * so our matrices go
     *
     *   ( n 0 ) ( -n^2 2.3.d^2 ) ( -n^2 4.5.d^2 ) ...
     *   ( 0 d ) (    0 2.3.d^2 ) (    0 4.5.d^2 )
     */

    bigint n, d, n2, d2, k, kk;
    bool negate;
    int crState;

  public:
    SinRational(const bigint &an, const bigint &ad, bool anegate)
        : n(an), d(ad), negate(anegate)
    {
        crState = -1;
    }

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

    bool gen_interval(bigint *low, bigint *high)
    {
        /* I totally made these numbers up, but they seem to work. Ahem. */
        *low = 0;
        *high = 2;
        return false;
    }

    bool gen_matrix(bigint *matrix)
    {
        crBegin;

        /*
         * The initial anomalous matrix.
         */
        matrix[1] = matrix[2] = 0;
        matrix[0] = (negate ? -n : n);
        matrix[3] = d;
        crReturn(false);

        /*
         * Then the regular series.
         */
        k = 1;
        n2 = n*n;
        d2 = d*d;
        while (1) {
            kk = ++k;
            kk *= ++k;
            matrix[0] = -n2;
            matrix[1] = matrix[3] = d2 * kk;
            matrix[2] = 0;
            crReturn(false);
        }

        crEnd;
    }
};

class CosSqrtRational : public Source {
    /*
     * The cosine series is closely related to the sine series.
     *
     * Let x = n/d. Then we have
     *
     *                    n^2      n^4
     *   cos(n/d) = 1 - ------ + ------ - ...
     *                  2! d^2   4! d^4
     *
     *                        n^2           n^2
     *            = ( 1 - ------- ( 1 - ------- ( ... ) ) )
     *                    1.2.d^2       3.4.d^2
     *
     * so our matrices go
     *
     *   ( -n^2 1.2.d^2 ) ( -n^2 3.4.d^2 ) ...
     *   (    0 1.2.d^2 ) (    0 3.4.d^2 )
     *
     * However, we're going to need this to be a monotonic function in
     * the interval it works in, and it's not, because of the turning
     * point at zero. To deal with this, we arrange that this class
     * receives n^2 and d^2 as parameters, i.e. it gets its input
     * pre-squared. And then we can do the squaring step in the input
     * value to spigot_monotone.
     */

    bigint n2, d2, k, kk;
    bool negate;
    int crState;

  public:
    CosSqrtRational(const bigint &an2, const bigint &ad2, bool anegate)
        : n2(an2), d2(ad2), negate(anegate)
    {
        crState = -1;
    }

    virtual CosSqrtRational *clone() {
        return new CosSqrtRational(n2, d2, negate);
    }

    bool gen_interval(bigint *low, bigint *high)
    {
        /* I totally made these numbers up, but they seem to work. Ahem. */
        *low = 0;
        *high = 2;
        return false;
    }

    bool gen_matrix(bigint *matrix)
    {
        crBegin;

        /*
         * Negate the entire thing if necessary.
         */
        if (negate) {
            matrix[1] = matrix[2] = 0;
            matrix[0] = -1;
            matrix[3] = 1;
            crReturn(false);
        }

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

        crEnd;
    }
};

class AtanRational : public Source {
    /*
     * We compute atan by the continued fraction method, which is
     * 
     *   atan(z) = z/(1+z^2/(3+4z^2/(5+...)))
     *
     * So let x = n/d. Then we have
     *
     *   atan(n/d) = (x |-> n/d(1+x)) o
     *      (x |-> n^2/(d^2(3+x) o (x |-> 4n^2/(d^2(5+x) o ...
     *               
     * so our matrices go
     *
     *   ( 0 n ) (  0   n^2 ) (  0  4n^2 ) (  0  9n^2 ) (  0  16n^2 ) ...
     *   ( d d ) ( d^2 3d^2 ) ( d^2 5d^2 ) ( d^2 7d^2 ) ( d^2  9d^2 )
     */

    bigint n, d, n2, d2, k, k2;
    int crState;

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

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

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

    bool gen_matrix(bigint *matrix)
    {
        crBegin;

        /*
         * The initial anomalous matrix.
         */
        matrix[0] = 0;
        matrix[1] = n;
        matrix[2] = matrix[3] = d;
        crReturn(false);

        /*
         * Then the regular series.
         */
        k = 3;
        k2 = 1;
        n2 = n*n;
        d2 = d*d;
        while (1) {
            matrix[0] = 0;
            matrix[1] = n2 * k2;
            matrix[2] = d2;
            matrix[3] = d2 * k;
            k2 += k;
            k += 2;
            crReturn(false);
        }

        crEnd;
    }
};

struct SinConstructor : MonotoneConstructor {
    MonotoneConstructor *clone() { return new SinConstructor(); }
    Spigot *construct(const bigint &n, const bigint &d) {
        return new SinRational(n, d, false);
    }
};
struct SinNegConstructor : MonotoneConstructor {
    MonotoneConstructor *clone() { return new SinNegConstructor(); }
    Spigot *construct(const bigint &n, const bigint &d) {
        return new SinRational(n, d, true);
    }
};
struct CosSqrtConstructor : MonotoneConstructor {
    MonotoneConstructor *clone() { return new CosSqrtConstructor(); }
    Spigot *construct(const bigint &n, const bigint &d) {
        return new CosSqrtRational(n, d, false);
    }
};
struct CosSqrtNegConstructor : MonotoneConstructor {
    MonotoneConstructor *clone() { return new CosSqrtNegConstructor(); }
    Spigot *construct(const bigint &n, const bigint &d) {
        return new CosSqrtRational(n, d, true);
    }
};
struct AtanConstructor : MonotoneConstructor {
    MonotoneConstructor *clone() { return new AtanConstructor(); }
    Spigot *construct(const bigint &n, const bigint &d) {
        return new AtanRational(n, d);
    }
};

Spigot *spigot_sincos(Spigot *a, int is_cos)
{
    bigint n;
    int quadrant;

    /*
     * Range-reduce by subtracting some multiple of pi/2.
     *
     * We absolutely don't want to do this _exactly_ right, by
     * rounding every number to the very nearest multiple of pi/2:
     * that would introduce an exactness hazard, so that sin(pi/4)
     * would hang without producing any output not because it's
     * fundamentally uncomputable but simply because the range
     * reducer couldn't decide which way to throw it.
     *
     * So instead, we find a near-enough rational approximation to the
     * answer, and range-reduce based on that.
     */
    {
        /*
         * Start by finding n such that n/64 is close (enough) to
         * a/pi.
         */
        StaticGenerator test(spigot_div(a->clone(), spigot_pi()));
        n = test.get_approximate_approximant(64);
    }

    /*
     * What we really want is the integer part of a / (pi/2), which is
     * close (enough) to (2*n/64) = n/32. Except that we want to round
     * to nearest, so add 16 first.
     */
    n = fdiv(n + 16, 32);
    if (n != 0) {
        /* Subtract off that multiple of pi/2. */
        a = spigot_sub(a, spigot_mul(spigot_rational(n, 2), spigot_pi()));
    }

    quadrant = (int)((n % 4U) + is_cos) & 3;

    if (quadrant == 0)
        return spigot_monotone(new SinConstructor, a);
    else if (quadrant == 1)
        return spigot_monotone(new CosSqrtConstructor, spigot_square(a));
    else if (quadrant == 2)
        return spigot_monotone(new SinNegConstructor, a);
    else /* (quadrant == 3) */
        return spigot_monotone(new CosSqrtNegConstructor, spigot_square(a));
}

Spigot *spigot_sin(Spigot *a)
{
    return spigot_sincos(a, 0);
}

Spigot *spigot_cos(Spigot *a)
{
    return spigot_sincos(a, 1);
}

Spigot *spigot_tan(Spigot *a)
{
    Spigot *ac = a->clone();
    return spigot_div(spigot_sin(a), spigot_cos(ac));
}

Spigot *spigot_atan(Spigot *a)
{
    /*
     * The only range reduction we need for atan is to spot numbers
     * with really big magnitude. Ideally, we'd spot anything with a
     * magnitude bigger than 1, and reduce it by taking the reciprocal
     * (using the identity atan(x) = pi/2 - atan(1/x), or the same
     * with -pi/2 for negative numbers).
     *
     * But, as usual, we don't need to get the right side of 1 in all
     * cases - it's fine to let borderline cases go in either
     * direction. So we avoid exactness hazards by doing only an
     * approximate check.
     */
    bool take_reciprocal = false;
    int sign = 0;

    {
        StaticGenerator test(a->clone());
        bigint n = test.get_approximate_approximant(32);
        if (n > 32) {
            sign = +1;
            take_reciprocal = true;
        } else if (n < -32) {
            sign = -1;
            take_reciprocal = true;
        }
    }

    if (take_reciprocal)
        a = spigot_reciprocal(a);

    a = spigot_monotone(new AtanConstructor, a);

    if (take_reciprocal)
        a = spigot_sub(spigot_rational_mul(spigot_pi(), sign, 2), a);

    return a;
}

Spigot *spigot_asin(Spigot *a)
{
    /*
     * The obvious way to compute asin(a) is as atan(a/sqrt(1-a^2)).
     * This is fine unless a is 1 or -1, or an expression which
     * non-obviously converges to that, in which case the argument to
     * atan becomes infinite and things go wrong. So we range-reduce
     * by first spotting things close to 1, and we generate those a
     * different way.
     */
    bool invert = false;
    int sign = 0;
    {
        StaticGenerator test(a->clone());
        bigint n = test.get_approximate_approximant(32);
        if (n > 16) {
            sign = +1;
            invert = true;
        } else if (n < -16) {
            sign = -1;
            invert = true;
        }
    }

    Spigot *sqrt_1_minus_a_squared =
        spigot_sqrt(spigot_quadratic(a->clone(), -1, 0, 1));

    if (!invert) {
        return spigot_atan(spigot_div(a, sqrt_1_minus_a_squared));
    } else {
        /*
         * For numbers close to 1 or -1, we invert the argument to
         * atan, i.e. we compute atan(sqrt(1-a^2)/a). Then we have to
         * subtract the result from pi/2 or -pi/2.
         */
        return spigot_sub(spigot_rational_mul(spigot_pi(), sign, 2),
                          spigot_atan(spigot_div(sqrt_1_minus_a_squared, a)));
    }
}

Spigot *spigot_acos(Spigot *a)
{
    /*
     * Simplest way to get acos right is just to subtract asin from
     * pi/2, having first special-cased for acos(1) = 0 exactly.
     */
    bigint n, d;
    if (a->is_rational(&n, &d) && n == d) {
        return spigot_integer(0);
    }
    return spigot_sub(spigot_rational_mul(spigot_pi(), 1, 2), spigot_asin(a));
}

static Spigot *spigot_atan2_inner(Spigot *y, Spigot *x,
                                  Spigot *(*pi)(),
                                  Spigot *(*scale)(Spigot *))
{
    bigint yn, yd, xn, xd;
    bool yrat, xrat, yzero, xzero;

    /*
     * Start by spotting exact zeroes among the inputs.
     */
    yrat = y->is_rational(&yn, &yd);
    yzero = yrat && yn == 0;
    xrat = x->is_rational(&xn, &xd);
    xzero = xrat && xn == 0;

    if (yzero && xzero) {
        throw spigot_error("atan2(0,0) has no answer");
    }
    if (yzero) {
        if (get_sign(x->clone()) < 0)
            return pi();
        else
            return spigot_integer(0);
    }
    if (xzero) {
        if (get_sign(y->clone()) < 0)
            return spigot_rational_mul(pi(), -1, 2);
        else
            return spigot_rational_mul(pi(), 1, 2);
    }

    if (yrat && xrat && yd == xd && bigint_abs(yn) == bigint_abs(xn)) {
        /*
         * For the case where we're returning an answer in degrees
         * (see below), we must give special handling to the case
         * where |y|=|x|, because it'll have to be a multiple of 45
         * degrees.
         */
        if (yn > 0 && xn > 0)
            return spigot_rational_mul(pi(), 1, 4);
        else if (yn > 0 && xn < 0)
            return spigot_rational_mul(pi(), 3, 4);
        else if (yn < 0 && xn > 0)
            return spigot_rational_mul(pi(), -1, 4);
        else if (yn < 0 && xn < 0)
            return spigot_rational_mul(pi(), -3, 4);
        else
            assert(!"Should have hit one of those cases");
    }

    /*
     * OK, those are the silly answers. Now for the sensible ones. We
     * begin by doing a _parallel_ sign test of y and x, in case one
     * of them is non-obviously zero. (If both are non-obviously zero,
     * we really do have to hang - there's nothing we can do.)
     */
    int s = parallel_sign_test(y->clone(), x->clone());
    if (s == +2) {
        /*
         * If x is positive, we're in the right half-plane, and can
         * just use normal atan.
         */
        return scale(spigot_atan(spigot_div(y, x)));
    } else if (s == +1) {
        /*
         * If y is positive, we're in the top half-plane, so we can
         * return pi/2 + atan(-x/y) = pi/2 - atan(x/y).
         */
        return spigot_sub(spigot_rational_mul(pi(), 1, 2),
                          scale(spigot_atan(spigot_div(x, y))));
    } else if (s == -1) {
        /*
         * If y is negative, we're in the bottom half-plane, so we can
         * do the same thing with -pi/2.
         */
        return spigot_sub(spigot_rational_mul(pi(), -1, 2),
                          scale(spigot_atan(spigot_div(x, y))));
    } else {
        /*
         * If x is negative, we really do need to know the exact sign
         * of y before we can decide whether to adjust atan(y/x) by
         * plus or minus pi.
         */
        if (get_sign(y->clone()) >= 0)  /* top left quadrant */
            return spigot_add(scale(spigot_atan(spigot_div(y, x))),
                              pi());
        else                            /* bottom left quadrant */
            return spigot_sub(scale(spigot_atan(spigot_div(y, x))),
                              pi());
    }
}

static Spigot *spigot_identity(Spigot *x)
{
    return x;
}

Spigot *spigot_atan2(Spigot *y, Spigot *x)
{
    return spigot_atan2_inner(y, x, spigot_pi, spigot_identity);
}

/* ----------------------------------------------------------------------
 * Trig functions taking arguments in degrees.
 *
 * These are basically equivalent to obvious things like sin(x*pi/180)
 * or asin(x)*180/pi, except that we take a little care to handle
 * rational input values which give rise to rational output values.
 *
 * The full set of such values is helpfully given for us by Niven's
 * theorem[1], which says (paraphrased) that the only rational values
 * of sin(x) arising from x being a rational multiple of pi are -1,
 * -1/2, 0, 1/2 and 1.
 *
 * This implies that the only input values we need to treat
 * interestingly for sind and cosd are multiples of 30 degrees.
 *
 * For tand, we can prove a corollary of Niven's theorem: the only
 * rational values of tan(x) arising from x being a rational multiple
 * of pi are -1, 0 and 1 (and infinity, which kind of counts in the
 * kind of sense that it's kind of a reciprocal of a rational, and
 * also counts in the more practical sense that we have to care about
 * it here).
 *
 * Proof: suppose tan(x) = p/q. Then the complex number q+ip has
 * argument x. If x is a rational multiple of pi, that implies that
 * q+ip has the same argument as a root of unity, i.e. q+ip is a real
 * multiple of a root of unity. So (q+ip)^2 is also a real multiple of
 * a root of unity. But (q+ip)^2 has modulus (p^2+q^2), which is an
 * integer, and hence it's a _rational_ multiple of a root of unity,
 * i.e. (q+ip)^2/(p^2+q^2) actually is a root of unity which has
 * _both_ real and imaginary parts rational. By Niven's theorem, those
 * real and imaginary parts must have absolute value 0, 1/2 or 1; the
 * 1/2 case is ruled out by _both_ parts having to be rational (any
 * angle which has sine 1/2 or -1/2 has an irrational cosine and vice
 * versa), so in fact both real and imaginary parts of that number
 * must be 0 or -1 or +1. Hence, (q+ip)^2 must be a real multiple of a
 * 4th root of unity; so (q+ip) must be a real multiple of an 8th root
 * of unity, i.e. its argument is a multiple of 45 degrees and so its
 * tangent is either 0, +1, -1, or infinite. []
 *
 * [1] http://en.wikipedia.org/wiki/Niven%27s_theorem
 */

Spigot *spigot_sind(Spigot *x)
{
    bigint n, d;

    if (x->is_rational(&n, &d) && d == 1 && n % 30U == 0) {
        // Reduce mod 360 using fdiv(), to get a non-negative result.
        switch ((int)(n - 360 * fdiv(n, 360))) {
          case 0: case 180:
            return spigot_integer(0);
          case 90:
            return spigot_integer(1);
          case 270:
            return spigot_integer(-1);
          case 30: case 150:
            return spigot_rational(1, 2);
          case 210: case 330:
            return spigot_rational(-1, 2);
        }
    }
    return spigot_sin(spigot_mul(spigot_pi(), spigot_rational_mul(x, 1, 180)));
}

Spigot *spigot_cosd(Spigot *x)
{
    bigint n, d;

    if (x->is_rational(&n, &d) && d == 1 && n % 30U == 0) {
        // Reduce mod 360 using fdiv(), to get a non-negative result.
        switch ((int)(n - 360 * fdiv(n, 360))) {
          case 90: case 270:
            return spigot_integer(0);
          case 0:
            return spigot_integer(1);
          case 180:
            return spigot_integer(-1);
          case 60: case 300:
            return spigot_rational(1, 2);
          case 120: case 240:
            return spigot_rational(-1, 2);
        }
    }
    return spigot_cos(spigot_mul(spigot_pi(), spigot_rational_mul(x, 1, 180)));
}

Spigot *spigot_tand(Spigot *x)
{
    bigint n, d;

    if (x->is_rational(&n, &d) && d == 1 && n % 45U == 0) {
        // Reduce mod 180 using fdiv(), to get a non-negative result.
        switch ((int)(n - 180 * fdiv(n, 180))) {
          case 0: case 180:
            return spigot_integer(0);
          case 45:
            return spigot_integer(1);
          case 135:
            return spigot_integer(-1);
          case 90:
            throw spigot_error("tand of an odd multiple of 90 degrees");
        }
    }
    return spigot_tan(spigot_mul(spigot_pi(), spigot_rational_mul(x, 1, 180)));
}

Spigot *spigot_asind(Spigot *x)
{
    bigint n, d;

    if (x->is_rational(&n, &d) && (bigint)2 % d == 0) {
        if (d == 1 && n == 0)
            return spigot_integer(0);
        else if (d == 1 && n == -1)
            return spigot_integer(-90);
        else if (d == 1 && n == +1)
            return spigot_integer(+90);
        else if (d == 2 && n == -1)
            return spigot_integer(-30);
        else if (d == 2 && n == +1)
            return spigot_integer(+30);
    }
    return spigot_div(spigot_rational_mul(spigot_asin(x), 180, 1),
                      spigot_pi());
}

Spigot *spigot_acosd(Spigot *x)
{
    bigint n, d;

    if (x->is_rational(&n, &d) && (bigint)2 % d == 0) {
        if (d == 1 && n == 0)
            return spigot_integer(90);
        else if (d == 1 && n == -1)
            return spigot_integer(180);
        else if (d == 1 && n == +1)
            return spigot_integer(+0);
        else if (d == 2 && n == -1)
            return spigot_integer(120);
        else if (d == 2 && n == +1)
            return spigot_integer(60);
    }
    return spigot_div(spigot_rational_mul(spigot_acos(x), 180, 1),
                      spigot_pi());
}

Spigot *spigot_atand(Spigot *x)
{
    bigint n, d;

    if (x->is_rational(&n, &d) && d == 1) {
        if (n == 0)
            return spigot_integer(0);
        else if (n == -1)
            return spigot_integer(-45);
        else if (n == +1)
            return spigot_integer(+45);
    }
    return spigot_div(spigot_rational_mul(spigot_atan(x), 180, 1),
                      spigot_pi());
}

static Spigot *one_hundred_and_eighty()
{
    return spigot_integer(180);
}
static Spigot *radians_to_degrees(Spigot *x)
{
    return spigot_div(spigot_rational_mul(x, 180, 1), spigot_pi());
}

Spigot *spigot_atan2d(Spigot *y, Spigot *x)
{
    /*
     * Most of the above functions are wrappers around the radian
     * versions with extra special-case handling on the front. This
     * one has enough complicated special cases already that I decided
     * it was easier to make the main version parametric in pi :-)
     */
    return spigot_atan2_inner(y, x, one_hundred_and_eighty,
                              radians_to_degrees);
}