/**********************************************************************

 *

 * GEOS - Geometry Engine Open Source

 * http://geos.osgeo.org

 *

 * Copyright (C) 2020 Crunchy Data

 *

 * This is free software; you can redistribute and/or modify it under

 * the terms of the GNU Lesser General Public Licence as published

 * by the Free Software Foundation.

 * See the COPYING file for more information.

 *

 **********************************************************************/

#include <cmath>

#include <geos/profiler.h>

#include <geos/math/DD.h>

namespace geos {

namespace math { // geos.util

/* private */

void DD::setNaN()

{

    hi = std::numeric_limits<double>::quiet_NaN();

}

/* public */

bool DD::isInf() const

{

    return std::isinf(hi);

}

/* public */

bool DD::isNaN() const

{

    return std::isnan(hi);

}

/* public */

bool DD::isNegative() const

{

    return hi < 0.0 || (hi == 0.0 && lo < 0.0);

}

/* public */

bool DD::isPositive() const

{

    return hi > 0.0 || (hi == 0.0 && lo > 0.0);

}

/* public */

bool DD::isZero() const

{

    return hi == 0.0 && lo == 0.0;

}

/* public */

double DD::doubleValue() const

{

    return hi + lo;

}

/* public */

int DD::intValue() const

{

    return (int) hi;

}

/* public */

void DD::selfAdd(const DD &y)

{

    return selfAdd(y.hi, y.lo);

}

/* public */

void DD::selfAdd(double yhi, double ylo)

{

    double H, h, T, t, S, s, e, f;

    S = hi + yhi;

    T = lo + ylo;

    e = S - hi;

    f = T - lo;

    s = S-e;

    t = T-f;

    s = (yhi-e)+(hi-s);

    t = (ylo-f)+(lo-t);

    e = s+T; H = S+e; h = e+(S-H); e = t+h;

    double zhi = H + e;

    double zlo = e + (H - zhi);

    hi = zhi;

    lo = zlo;

    return;

}

/* public */

void DD::selfAdd(double y)

{

    double H, h, S, s, e, f;

    S = hi + y;

    e = S - hi;

    s = S - e;

    s = (y - e) + (hi - s);

    f = s + lo;

    H = S + f;

    h = f + (S - H);

    hi = H + h;

    lo = h + (H - hi);

    return;

}

/* public */

DD operator+(const DD &lhs, const DD &rhs)

{

    DD rv(lhs.hi, lhs.lo);

    rv.selfAdd(rhs);

    return rv;

}

/* public */

DD operator+(const DD &lhs, double rhs)

{

    DD rv(lhs.hi, lhs.lo);

    rv.selfAdd(rhs);

    return rv;

}

/* public */

void DD::selfSubtract(const DD &d)

{

    return selfAdd(-1*d.hi, -1*d.lo);

}

/* public */

void DD::selfSubtract(double p_hi, double p_lo)

{

    return selfAdd(-1*p_hi, -1*p_lo);

}

/* public */

void DD::selfSubtract(double y)

{

    return selfAdd(-1*y, 0.0);

}

/* public */

DD operator-(const DD &lhs, const DD &rhs)

{

    DD rv(lhs.hi, lhs.lo);

    rv.selfSubtract(rhs);

    return rv;

}

/* public */

DD operator-(const DD &lhs, double rhs)

{

    DD rv(lhs.hi, lhs.lo);

    rv.selfSubtract(rhs);

    return rv;

}

/* public */

void DD::selfMultiply(double yhi, double ylo)

{

    double hx, tx, hy, ty, C, c;

    C = SPLIT * hi; hx = C-hi; c = SPLIT * yhi;

    hx = C-hx; tx = hi-hx; hy = c-yhi;

    C = hi*yhi; hy = c-hy; ty = yhi-hy;

    c = ((((hx*hy-C)+hx*ty)+tx*hy)+tx*ty)+(hi*ylo+lo*yhi);

    double zhi = C+c; hx = C-zhi;

    double zlo = c+hx;

    hi = zhi;

    lo = zlo;

    return;

}

/* public */

void DD::selfMultiply(DD const &d)

{

    return selfMultiply(d.hi, d.lo);

}

/* public */

void DD::selfMultiply(double y)

{

    return selfMultiply(y, 0.0);

}

/* public */

DD operator*(const DD &lhs, const DD &rhs)

{

    DD rv(lhs.hi, lhs.lo);

    rv.selfMultiply(rhs);

    return rv;

}

/* public */

DD operator*(const DD &lhs, double rhs)

{

    DD rv(lhs.hi, lhs.lo);

    rv.selfMultiply(rhs);

    return rv;

}

/* public */

void DD::selfDivide(double yhi, double ylo)

{

    double hc, tc, hy, ty, C, c, U, u;

    C = hi/yhi; c = SPLIT*C; hc =c-C;

    u = SPLIT*yhi; hc = c-hc;

    tc = C-hc; hy = u-yhi; U = C * yhi;

    hy = u-hy; ty = yhi-hy;

    u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;

    c = ((((hi-U)-u)+lo)-C*ylo)/yhi;

    u = C+c;

    hi = u;

    lo = (C-u)+c;

    return;

}

/* public */

void DD::selfDivide(const DD &d)

{

    return selfDivide(d.hi, d.lo);

}

/* public */

void DD::selfDivide(double y)

{

    return selfDivide(y, 0.0);

}

/* public */

DD operator/(const DD &lhs, const DD &rhs)

{

    DD rv(lhs.hi, lhs.lo);

    rv.selfDivide(rhs);

    return rv;

}

/* public */

DD operator/(const DD &lhs, double rhs)

{

    DD rv(lhs.hi, lhs.lo);

    rv.selfDivide(rhs);

    return rv;

}

DD operator- (const DD& x)

{

    return x.negate();

}

/* public */

DD DD::negate() const

{

    DD rv(hi, lo);

    if (rv.isNaN())

    {

        return rv;

    }

    rv.hi = -hi;

    rv.lo = -lo;

    return rv;

}

/* public static */

DD DD::reciprocal() const

{

    double  hc, tc, hy, ty, C, c, U, u;

    C = 1.0/hi;

    c = SPLIT*C;

    hc = c-C;

    u = SPLIT*hi;

    hc = c-hc; tc = C-hc; hy = u-hi; U = C*hi; hy = u-hy; ty = hi-hy;

    u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;

    c = ((((1.0-U)-u))-C*lo)/hi;

    double zhi = C+c;

    double zlo = (C-zhi)+c;

    return DD(zhi, zlo);

}

DD DD::floor() const

{

    DD rv(hi, lo);

    if (isNaN()) return rv;

    double fhi = std::floor(hi);

    double flo = 0.0;

    // Hi is already integral.  Floor the low word

    if (fhi == hi) {

      flo = std::floor(lo);

    }

      // do we need to renormalize here?

    rv.hi = fhi;

    rv.lo = flo;

    return rv;

}

DD DD::ceil() const

{

    DD rv(hi, lo);

    if (isNaN()) return rv;

    double fhi = std::ceil(hi);

    double flo = 0.0;

    // Hi is already integral.  Ceil the low word

    if (fhi == hi) {

      flo = std::ceil(lo);

      // do we need to renormalize here?

    }

    rv.hi = fhi;

    rv.lo = flo;

    return rv;

}

int DD::signum() const

{

    if (hi > 0) return 1;

    if (hi < 0) return -1;

    if (lo > 0) return 1;

    if (lo < 0) return -1;

    return 0;

}

DD DD::rint() const

{

    DD rv(hi, lo);

    if (isNaN()) return rv;

     return (rv + 0.5).floor();

}

/* public static */

DD DD::trunc(const DD &d)

{

    DD rv(d);

    if (rv.isNaN()) return rv;

    if (rv.isPositive())

        return rv.floor();

    return rv.ceil();

}

/* public static */

DD DD::abs(const DD &d)

{

    DD rv(d);

    if (rv.isNaN()) return rv;

    if (rv.isNegative())

        return rv.negate();

    return rv;

}

/* public static */

DD DD::determinant(const DD &x1, const DD &y1, const DD &x2, const DD &y2)

{

    return (x1 * y2) - (y1 * x2);

}

/* public static */

DD DD::determinant(double x1, double y1, double x2, double y2)

{

    return determinant(DD(x1), DD(y1), DD(x2), DD(y2) );

}

/**

* Computes the value of this number raised to an integral power.

* Follows semantics of Java Math.pow as closely as possible.

*/

/* public static */

DD DD::pow(const DD &d, int exp)

{

    if (exp == 0.0)

        return DD(1.0);

    DD r(d);

    DD s(1.0);

    int n = std::abs(exp);

    if (n > 1) {

        /* Use binary exponentiation */

        while (n > 0) {

        if (n % 2 == 1) {

            s.selfMultiply(r);

        }

        n /= 2;

        if (n > 0)

            r = r*r;

        }

    } else {

        s = r;

    }

    /* Compute the reciprocal if n is negative. */

    if (exp < 0)

        return s.reciprocal();

    return s;

}

}

}