package dbldbl

import "math"

// NaN returns a “not-a-number” value.
func NaN() Number {
        return Number{y: math.NaN()}
}

// Inf returns positive infinity if sign >= 0, negative infinity if sign < 0.
func Inf(sign int) Number {
        return Number{y: math.Inf(sign)}
}

var (
        E   = Number{2.718281828459045, +0x1.4d57ee2b1013ap-53} // https://oeis.org/A001113
        Pi  = Number{3.141592653589793, +0x1.1a62633145c07p-53} // https://oeis.org/A000796
        Phi = Number{1.618033988749895, -0x1.f506319fcfd18p-55} // https://oeis.org/A001622

        Sqrt2   = Number{1.4142135623730951, -0x1.bdd3413b26456p-54} // https://oeis.org/A002193
        SqrtE   = Number{1.6487212707001282, -0x1.b4690082a4906p-55} // https://oeis.org/A019774
        SqrtPi  = Number{1.772453850905516, -0x1.618f13eb7ca89p-54}  // https://oeis.org/A002161
        SqrtPhi = Number{1.272019649514069, -0x1.6daabafeede14p-56}  // https://oeis.org/A139339

        Ln2  = Number{0.6931471805599453, +0x1.abc9e3b39804p-56} // https://oeis.org/A002162
        Ln10 = Number{2.302585092994046, -0x1.f48ad494ea3e9p-53} // https://oeis.org/A002392

        twoPi   = Number{Pi.y * 2, Pi.x * 2}
        halfPi  = Number{Pi.y / 2, Pi.x / 2}
        twoOfPi = Number{0.6366197723675814, -0x1.6b01ec5417056p-55} // https://oeis.org/A060294
)

package dbldbl

import "strconv"

// String implements [fmt.Stringer].
func (n Number) String() string {
        return strconv.FormatFloat(n.y, 'g', 15, 64)
}

// GoString implements [fmt.GoStringer].
func (n Number) GoString() string {
        y := strconv.FormatFloat(n.y, 'g', -1, 64)
        x := strconv.FormatFloat(n.x, 'x', -1, 64)
        sep := ", +"
        if x[0] == '-' {
                sep = sep[:2]
        }
        return "Number{" + y + sep + x + "}"
}

package dbldbl

import "math"

// Sinh returns the hyperbolic sine of n (approximate).
func Sinh(n Number) Number {
        switch {
        case n.y == 0:
                return n
        case n.y < 0:
                return Neg(Sinh(Neg(n)))
        }
        t := Expm1(n)
        if !isFinite(t.y) {
                return t
        }
        return shift(Add(t, Div(t, AddFloat(t, 1))), -1)
}

// Cosh returns the hyperbolic cosine of n (approximate).
func Cosh(n Number) Number {
        t := Exp(Abs(n))
        return shift(Add(t, Inv(t)), -1)
}

// Tanh returns the hyperbolic tangent of n (approximate).
func Tanh(n Number) Number {
        switch {
        case n.y == 0:
                return n
        case math.Abs(n.y) > 100:
                return Number{y: math.Copysign(1, n.y)}
        }
        t := Expm1(shift(n, 1))
        return Div(t, AddFloat(t, 2))
}

// Asinh returns the inverse hyperbolic sine of n (approximate).
func Asinh(n Number) Number {
        switch {
        case n.y == 0 || !isFinite(n.y):
                return n
        case n.y < 0:
                return Neg(Asinh(Neg(n)))
        case n.y > 100:
                return Add(Log(n), Ln2)
        }

        t := Sqr(n)
        return Log1p(Add(n, Div(t, AddFloat(Sqrt(AddFloat(t, 1)), 1))))
}

// Acosh returns the inverse hyperbolic cosine of n (approximate).
func Acosh(n Number) Number {
        switch {
        case n.y < 1:
                return NaN()
        case !isFinite(n.y):
                return n
        case n.y > 100:
                return Add(Log(n), Ln2)
        }

        t := AddFloat(n, -1)
        return Log1p(Add(t, Sqrt(Add(shift(t, 1), Sqr(t)))))
}

// Atanh returns the inverse hyperbolic tangent of n (approximate).
func Atanh(n Number) Number {
        switch {
        case n.y == 0:
                return n
        case n.y < 0:
                return Neg(Atanh(Neg(n)))
        case n.y > 1 || n.y == 1 && n.x > 0:
                return NaN()
        case n == Float(1):
                return Inf(0)
        }

        return shift(Log1p(Div(shift(n, 1), SubFloat(1, n))), -1)
}

package dbldbl

import "math"

func twoSumQuick(x, y float64) Number {
        // log₂|x| ≥ log₂|y|
        r := float64(x + y)
        e := y - float64(r-x)
        return Number{r, e}
}

func twoSum(x, y float64) Number {
        r := float64(x + y)
        t := float64(r - x)
        e := float64(x-float64(r-t)) + float64(y-t)
        return Number{r, e}
}

func twoDiff(x, y float64) Number {
        r := float64(x - y)
        t := float64(r - x)
        e := float64(x-float64(r-t)) - float64(y+t)
        return Number{r, e}
}

func twoProd(x, y float64) Number {
        r := float64(x * y)
        e := math.FMA(x, y, -r)
        return Number{r, e}
}

func twoFMA(x, y float64, n Number) Number {
        p := twoProd(x, y)
        s := twoSum(n.y, p.y)
        t := twoSum(n.x, p.x)
        s = twoSumQuick(s.y, s.x+t.y)
        s = twoSumQuick(s.y, s.x+t.x)
        return s
}

func twoFMAQuick(x, y float64, n Number) Number {
        p := twoProd(x, y)
        s := twoSum(n.y, p.y)
        return twoSumQuick(s.y, s.x+(n.x+p.x)) // approximation
}

package dbldbl

import "math"

// Log returns the natural logarithm of n (approximate).
func Log(n Number) Number {
        switch {
        case n.y < 0:
                return NaN()
        case n.y == 0:
                return Inf(-1)
        case n.y == 1:
                return Number{y: math.Log1p(n.x)}
        case !isFinite(n.y):
                return n
        case n.y > 1:
                return Neg(Log(Inv(n)))
        }

        // Reduce the argument until it is less than 2⁻⁵⁵.
        var halvings int8
        for n.y > 0x1p-55 {
                n = Sqr(n) // log(n) = log(n²)/2
                halvings++
        }

        // For n<2⁻⁵⁵ this is accurate to 107 bits:
        // log(1/n) ≈ π / 2⋅AGM(1, 4⋅n)
        n = Div(halfPi, agm(Float(1), shift(n, 2)))
        return Neg(shift(n, -halvings)) // log(n) = -log(1/n)
}

// Exp returns eⁿ, the base-e exponential of n (approximate).
func Exp(n Number) Number {
        y := math.Exp(n.y)
        if y == 0 || !isFinite(y) {
                return Number{y: y}
        }
        // Newton's method: y + y⋅(n-log(y))
        t := Sub(n, Log(Float(y)))
        return twoSumQuick(y, y*t.y)
}

// Log1p returns the natural logarithm of 1 plus n (approximate).
// It is more accurate than Log(AddFloat(n, 1)) when n is near zero.
func Log1p(n Number) Number {
        // https://www.johndcook.com/blog/2012/07/25/trick-for-computing-log1x/
        u := AddFloat(n, 1)
        switch {
        case u.y < 0:
                return NaN()
        case u == Float(1) || !isFinite(u.y):
                return n
        }
        // log(1+n) = n⋅log(u)/(u-1), u=1+n
        return Mul(n, Div(Log(u), AddFloat(u, -1)))
}

// Expm1 returns eⁿ-1, the base-e exponential of n minus 1 (approximate).
// It is more accurate than SubFloat(Exp(n), 1) when n is near zero.
func Expm1(n Number) Number {
        y := math.Expm1(n.y)
        switch {
        case y == 0 || !isFinite(y):
                return Number{y: y}
        case y == -1:
                return SubFloats(math.Exp(n.y), 1)
        }
        // Newton's method: y + (y+1)⋅(n-log1p(y))
        t := Sub(n, Log1p(Float(y)))
        return twoSumQuick(y, (y+1)*t.y)
}

func agm(a, g Number) Number {
        // https://en.wikipedia.org/wiki/Arithmetic–geometric_mean
        for {
                t := shift(Add(a, g), -1)
                if t == a {
                        return a
                }
                g = Sqrt(Mul(a, g))
                a = t
        }
}

// Package dbldbl provides double-double precision arithmetic.
package dbldbl

// Number is a double-double precision number.
type Number struct {
        y, x float64
}

// Float creates a Number from a float64.
func Float(a float64) Number {
        return Number{y: a}
}

// Int creates a Number from an int64.
func Int(a int64) Number {
        if y := a >> 48 << 48; y == 0 {
                return Number{y: float64(a)}
        } else {
                return twoSumQuick(float64(y), float64(a-y))
        }
}

// Uint creates a Number from a uint64.
func Uint(a uint64) Number {
        if y := a >> 48 << 48; y == 0 {
                return Number{y: float64(a)}
        } else {
                return twoSumQuick(float64(y), float64(a-y))
        }
}

// Float converts this Number to a float64.
func (n Number) Float() (_ float64, exact bool) {
        return n.y, n.x == 0
}

// Int converts this Number to an int64.
func (n Number) Int() int64 {
        return int64(n.y) + int64(n.x)
}

// Uint converts this Number to a uint64.
func (n Number) Uint() uint64 {
        return uint64(n.y) + uint64(n.x)
}

package dbldbl

import "math"

// Neg negates n (exact).
func Neg(n Number) Number {
        return Number{-n.y, 0 - n.x}
}

// Signbit reports whether x is negative or negative zero.
func Signbit(n Number) bool {
        return math.Signbit(n.y)
}

// Abs returns the absolute value of n (exact).
func Abs(n Number) Number {
        if Signbit(n) {
                return Neg(n)
        }
        return n
}

// IsNaN reports whether n is a “not-a-number” value.
func IsNaN(n Number) bool {
        return math.IsNaN(n.y)
}

// IsInf reports whether f is an infinity, according to sign.
// If sign > 0, IsInf reports whether f is positive infinity.
// If sign < 0, IsInf reports whether f is negative infinity.
// If sign == 0, IsInf reports whether f is either infinity.
func IsInf(n Number, sign int) bool {
        return math.IsInf(n.y, sign)
}

// Trunc returns the integer value of n (exact).
func Trunc(n Number) Number {
        y := math.Trunc(n.y)
        switch {
        case y != n.y:
                return Number{y: y}
        case y < 0:
                return Number{y: y, x: math.Ceil(n.x)}
        default:
                return Number{y: y, x: math.Floor(n.x)}
        }
}

// Floor returns the greatest integer value less than or equal to n (exact).
func Floor(n Number) Number {
        y := math.Floor(n.y)
        if y != n.y {
                return Number{y: y}
        }
        return Number{y, math.Floor(n.x)}
}

// Ceil returns the least integer value greater than or equal to n (exact).
func Ceil(n Number) Number {
        y := math.Ceil(n.y)
        if y != n.y {
                return Number{y: y}
        }
        return Number{y, math.Ceil(n.x)}
}

// Round returns the nearest integer, rounding half away from zero (exact).
func Round(n Number) Number {
        y := math.Round(n.y)

        switch d := y - n.y; {
        case math.IsNaN(d):
                return n

        case d != 0:
                switch {
                case n.x < 0 && d == +0.5:
                        return twoSumQuick(y, -1)
                case n.x > 0 && d == -0.5:
                        return twoSumQuick(y, +1)
                }
                return Number{y: y}
        }

        x := math.Round(n.x)
        switch {
        case n.y < 0 && x-n.x == +0.5:
                x--
        case n.y > 0 && x-n.x == -0.5:
                x++
        }
        return twoSumQuick(n.y, x)
}

// Ldexp returns the product of n by 2ⁱ (exact).
func Ldexp(n Number, i int) Number {
        y := math.Ldexp(n.y, i)
        if y == 0 || !isFinite(y) {
                return Number{y: y}
        }
        return Number{y, math.Ldexp(n.x, i)}
}

// AddFloats returns the sum of a and b (exact).
func AddFloats(a, b float64) Number {
        s := twoSum(a, b)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        return s
}

// AddFloat returns the sum of a and b (exactly rounded).
func AddFloat(a Number, b float64) Number {
        s := twoSum(a.y, b)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        return twoSumQuick(s.y, s.x+a.x)
}

// Add returns the sum of a and b (exactly rounded).
func Add(a, b Number) Number {
        s := twoSum(a.y, b.y)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        t := twoSum(a.x, b.x)
        s = twoSumQuick(s.y, s.x+t.y)
        s = twoSumQuick(s.y, s.x+t.x)
        return s
}

// SubFloats returns the difference of a and b (exact).
func SubFloats(a, b float64) Number {
        s := twoDiff(a, b)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        return s
}

// SubFloat returns the difference of a and b (exactly rounded).
func SubFloat(a float64, b Number) Number {
        s := twoDiff(a, b.y)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        return twoSumQuick(s.y, s.x-b.x)
}

// Sub returns the difference of a and b (exactly rounded).
func Sub(a, b Number) Number {
        s := twoDiff(a.y, b.y)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        t := twoDiff(a.x, b.x)
        s = twoSumQuick(s.y, s.x+t.y)
        s = twoSumQuick(s.y, s.x+t.x)
        return s
}

// MulFloats returns the product of a and b (exact).
func MulFloats(a, b float64) Number {
        s := twoProd(a, b)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        return s
}

// MulFloat returns the product of a and b (exactly rounded).
func MulFloat(a Number, b float64) Number {
        s := twoProd(a.y, b)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        s.x = math.FMA(a.x, b, s.x)
        return twoSumQuick(s.y, s.x)
}

// Mul returns the product of a and b (approximate).
func Mul(a, b Number) Number {
        s := twoProd(a.y, b.y)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        s.x = math.FMA(a.x, b.x, s.x)
        s.x = math.FMA(a.y, b.x, s.x)
        s.x = math.FMA(a.x, b.y, s.x)
        return twoSumQuick(s.y, s.x)
}

// Div returns the quotient of a and b (approximate).
func Div(a, b Number) Number {
        y := a.y / b.y
        if y == 0 || !isFinite(y) {
                return Number{y: y}
        }
        t := twoProd(y, b.y)
        x := (a.y - t.y - t.x + a.x - y*b.x) / b.y
        return twoSumQuick(y, x)
}

// Sqr returns the square of n (approximate).
func Sqr(n Number) Number {
        s := twoProd(n.y, n.y)
        if !isFinite(s.y) {
                return Number{y: s.y}
        }
        s.x = math.FMA(n.x, n.x, s.x)
        s.x = math.FMA(n.y, n.x+n.x, s.x)
        return twoSumQuick(s.y, s.x)
}

// Sqrt returns the square root of n (approximate).
func Sqrt(n Number) Number {
        y := math.Sqrt(n.y)
        if y == 0 || !isFinite(y) {
                return Number{y: y}
        }
        // Newton's method: y + (n-y²)/2y
        t := twoProd(y, y)
        x := (n.y - t.y - t.x + n.x) * 0.5 / y
        return twoSumQuick(y, x)
}

// Cbrt returns the cube root of n (approximate).
func Cbrt(n Number) Number {
        y := math.Cbrt(n.y)
        if y == 0 || !isFinite(y) {
                return Number{y: y}
        }
        // Newton's method: y + (n/y²-y)/3
        t := Div(n, twoProd(y, y))
        x := (t.y - y + t.x) / 3
        return twoSumQuick(y, x)
}

// Inv returns the reciprocal of n (approximate).
func Inv(n Number) Number {
        y := 1 / n.y
        if y == 0 || !isFinite(y) {
                return Number{y: y}
        }
        // Newton's method: y + (1-n⋅y)⋅y
        t := twoProd(y, n.y)
        x := (1 - t.y - t.x - y*n.x) * y
        return twoSumQuick(y, x)
}

// InvSqrt returns the reciprocal of the square root of n (approximate).
func InvSqrt(n Number) Number {
        y := 1 / math.Sqrt(n.y)
        if y == 0 || !isFinite(y) {
                return Number{y: y}
        }
        // Newton's method: y + (y-n⋅y³)/2
        t := Mul(MulFloat(n, y), twoProd(y, y))
        x := (y - t.y - t.x) * 0.5
        return twoSumQuick(y, x)
}

// FMAFloat returns a⋅b + c (exactly rounded).
func FMAFloat(a, b float64, c Number) Number {
        s := twoFMA(a, b, c)
        if !isFinite(s.y) {
                return Number{y: math.FMA(a, b, c.y)}
        }
        return s
}

// FMA returns a⋅b + c (approximate).
func FMA(a, b, c Number) Number {
        s := twoFMA(a.y, b.y, c)
        if !isFinite(s.y) {
                return Number{y: math.FMA(a.y, b.y, c.y)}
        }
        s = twoFMAQuick(a.x, b.x, s)
        s = twoFMAQuick(a.y, b.x, s)
        s = twoFMAQuick(a.x, b.y, s)
        return s
}

// Cmp compares x and y and returns:
//
//        -1 if x <  y (incl. NaN < !NaN)
//         0 if x == y (incl. -0 == 0, -Inf == -Inf, +Inf == +Inf, NaN == NaN)
//        +1 if x >  y (incl. !NaN > NaN)
func Cmp(a, b Number) int {
        switch {
        case a.y < b.y:
                return -1
        case a.y > b.y:
                return +1
        case a.x < b.x:
                return -1
        case a.x > b.x:
                return +1
        case a.y == b.y:
                return 0
        case !IsNaN(b):
                return -1
        case !IsNaN(a):
                return +1
        default:
                return 0
        }
}

func isFinite(x float64) bool {
        return math.Float64bits(x)>>52&0x7ff != 0x7ff
}

func shift(n Number, i int8) Number {
        // This is like ldexp, but nicer when i is a small constant,
        // because it gets inlined and skips some branching.
        e := math.Float64frombits((1023 + uint64(i)) << 52)
        return Number{e * n.y, e * n.x}
}

package dbldbl

import "math"

// Pow returns bⁿ, the base-b exponential of n (approximate).
func Pow(b Number, n Number) Number {
        switch {
        case n.y == 0 || b == Float(1):
                return Float(1)
        case n == Float(1):
                return b
        case IsNaN(b) || IsNaN(n):
                return NaN()
        case b.y == 0:
                if Signbit(n) {
                        if isOddInteger(n) && Signbit(b) {
                                return Inf(-1)
                        }
                        return Inf(0)
                } else {
                        if isOddInteger(n) {
                                return b
                        }
                        return Number{}
                }
        case IsInf(n, 0):
                switch {
                case b == Float(-1):
                        return Float(1)
                case (SubFloat(1, Abs(b)).y > 0) == IsInf(n, 1):
                        return Number{}
                default:
                        return Inf(0)
                }
        case IsInf(b, 0):
                switch {
                case IsInf(b, -1):
                        var zero float64
                        return Pow(Float(-zero), Neg(n))
                case n.y < 0:
                        return Number{}
                case n.y > 0:
                        return Inf(0)
                }
        case n == Float(+0.5):
                return Sqrt(b)
        case n == Float(-0.5):
                return InvSqrt(b)
        case n == Floor(n):
                return pow(b, n)
        case b.y < 0:
                return NaN()
        }

        return Exp(Mul(Log(b), n))
}

// Pow10 returns 10ⁱ, the base-10 exponential of i (approximate).
func Pow10(i int) Number {
        return pow(Float(10), Int(int64(i)))
}

func pow(b Number, n Number) Number {
        r := Float(1)
        i := Abs(n)
        for {
                if isOddInteger(i) {
                        i = twoSumQuick(i.y, i.x-1)
                        r = Mul(r, b)
                }
                if i.y == 0 {
                        break
                }
                i = shift(i, -1)
                b = Sqr(b)
        }
        if n.y < 0 {
                return Inv(r)
        }
        return r
}

func isOddInteger(n Number) bool {
        return isOddFloat(n.x) || isOddFloat(n.y)
}

func isOddFloat(x float64) bool {
        xi, xf := math.Modf(x)
        return xf == 0 && int64(xi)&1 == 1
}

package dbldbl

import "math"

// Sincos returns Sin(n), Cos(n) (approximate).
func Sincos(n Number) (sin, cos Number) {
        switch {
        case n.y == 0:
                return n, Float(1)
        case !isFinite(n.y):
                return NaN(), NaN()
        }

        // Range reduction modulo π/2.
        k := Round(Mul(n, twoOfPi))
        t := Sub(n, Mul(k, halfPi))

        // Halve the angle until it is less than 2⁻⁵³.
        var halvings int8
        if _, e := math.Frexp(t.y); t.y != 0 && e > -53 {
                halvings = int8(53 + e)
        }

        // For |θ|<2⁻⁵³ these are accurate to 107 bits.
        sin = shift(t, -halvings) // sin(θ) ≈ θ
        cos = Float(1)            // cos(θ) ≈ 1

        // Double-angle formulae.
        for range halvings {
                s, c := sin, cos
                sin = shift(Mul(s, c), 1)           // sin(2⋅t) = 2⋅sin(θ)⋅cos(θ)
                cos = SubFloat(1, shift(Sqr(s), 1)) // cos(2⋅t) = 1 - 2⋅sin²(θ)
        }

        yi, _ := math.Modf(k.y)
        xi, _ := math.Modf(k.x)
        switch (int64(xi) | int64(yi)) & 3 {
        default:
                return sin, cos
        case 1:
                return cos, Neg(sin)
        case 2:
                return Neg(sin), Neg(cos)
        case 3:
                return Neg(cos), sin
        }
}

// Sin returns the sine of the radian argument n (approximate).
func Sin(n Number) Number {
        sin, _ := Sincos(n)
        return sin
}

// Cos returns the cosine of the radian argument n (approximate).
func Cos(n Number) Number {
        _, cos := Sincos(n)
        return cos
}

// Tan returns the tangent of the radian argument n (approximate).
func Tan(n Number) Number {
        sin, cos := Sincos(n)
        return Div(sin, cos)
}

// Asin returns the arcsine, in radians, of n (approximate).
func Asin(n Number) Number {
        // asin(θ) = atan2(θ, √(1-θ²))
        return Atan2(n, Sqrt(SubFloat(1, Sqr(n))))
}

// Acos returns the arccosine, in radians, of n (approximate).
func Acos(n Number) Number {
        // acos(θ) = atan2(√(1-θ²), θ)
        return Atan2(Sqrt(SubFloat(1, Sqr(n))), n)
}

// Atan returns the arctangent, in radians, of n (approximate).
func Atan(n Number) Number {
        switch {
        case n.y == 0:
                return n
        case n.y < 0:
                // atan(θ) = -atan(-θ)
                return Neg(Atan(Neg(n)))
        case n.y > 1:
                // atan(θ) = π/2 - atan(1/θ), if θ>0
                return Sub(halfPi, Atan(Inv(n)))
        case IsNaN(n):
                return NaN()
        }

        // Reduce the argument until it is less than 2⁻⁵³.
        var doublings int8
        for math.Abs(n.y) > 0x1p-53 {
                // atan(θ) = 2·atan(θ/(1+√(1+θ²)))
                n = Div(n, AddFloat(Sqrt(AddFloat(Sqr(n), 1)), 1))
                doublings++
        }

        // For |θ|<2⁻⁵³ this is accurate to 107 bits
        return shift(n, doublings) // atan(θ) ≈ θ
}

// Atan2 returns the arc tangent of y/x, using the signs of the two
// to determine the quadrant of the return value (approximate).
func Atan2(y, x Number) Number {
        switch {
        case y.y == 0 && x.y == 0:
                switch {
                case !Signbit(x):
                        return y
                case !Signbit(y):
                        return Pi
                default:
                        return Neg(Pi)
                }
        case IsInf(y, 0) && IsInf(x, 0):
                y = Number{y: math.Copysign(1, y.y)}
                x = Number{y: math.Copysign(1, x.y)}
        }

        z := Atan(Div(y, x))

        switch {
        case !Signbit(x):
                return z
        case !Signbit(y):
                return Add(z, Pi)
        default:
                return Sub(z, Pi)
        }
}