spire.math.Real.scala Maven / Gradle / Ivy
package spire.math
import scala.annotation.tailrec
import scala.math.{ScalaNumber, ScalaNumericConversions}
import spire.algebra.{Order, Trig, Signed}
import spire.syntax.nroot._
sealed trait Real extends ScalaNumber with ScalaNumericConversions { x =>
import Real.{roundUp, Exact, Inexact}
def apply(p: Int): SafeLong
def toRational(p: Int): Rational = this match {
case Exact(n) => n
case _ => Rational(x(p), SafeLong.two.pow(p))
}
def toRational: Rational = toRational(Real.bits)
// ugh scala.math
def doubleValue(): Double = toRational.toDouble
def floatValue(): Float = toRational.toFloat
def intValue(): Int = toRational.toInt
def longValue(): Long = toRational.toLong
def underlying(): Object = this
override def isValidChar: Boolean = {
val r = toRational
r.isWhole && r.isValidChar
}
override def isValidByte: Boolean = {
val r = toRational
r.isWhole && r.isValidByte
}
override def isValidShort: Boolean = {
val r = toRational
r.isWhole && r.isValidShort
}
override def isValidInt: Boolean = {
val r = toRational
r.isWhole && r.isValidInt
}
def isValidLong: Boolean = {
val r = toRational
r.isWhole && r.isValidLong
}
override def hashCode(): Int = toRational.hashCode
override def equals(y: Any): Boolean = y match {
case y: Real => this === y
case y => toRational.equals(y)
}
def ===(y: Real): Boolean =
(x compare y) == 0
def =!=(y: Real): Boolean =
!(this === y)
def compare(y: Real): Int = (x, y) match {
case (Exact(nx), Exact(ny)) => nx compare ny
case _ => (x - y).signum
}
def min(y: Real): Real = (x, y) match {
case (Exact(nx), Exact(ny)) => Exact(nx min ny)
case _ => Real(p => x(p) min y(p))
}
def max(y: Real): Real = (x, y) match {
case (Exact(nx), Exact(ny)) => Exact(nx max ny)
case _ => Real(p => x(p) max y(p))
}
def abs(): Real = this match {
case Exact(n) => Exact(n.abs)
case _ => Real(p => x(p).abs)
}
def signum(): Int = this match {
case Exact(n) => n.signum
case _ => x(Real.bits).signum
}
def unary_-(): Real = this match {
case Exact(n) => Exact(-n)
case _ => Real(p => -x(p))
}
def reciprocal(): Real = {
def findNonzero(i: Int): Int =
if (SafeLong.three <= x(i).abs) i else findNonzero(i + 1)
this match {
case Exact(n) => Exact(n.reciprocal)
case _ => Real({p =>
val s = findNonzero(0)
roundUp(Rational(SafeLong.two.pow(2 * p + 2 * s + 2), x(p + 2 * s + 2)))
})
}
}
def +(y: Real): Real = (x, y) match {
case (Exact(nx), Exact(ny)) => Exact(nx + ny)
case (Exact(Rational.zero), _) => y
case (_, Exact(Rational.zero)) => x
case _ => Real(p => roundUp(Rational(x(p + 2) + y(p + 2), 4)))
}
def -(y: Real): Real = x + (-y)
def *(y: Real): Real = (x, y) match {
case (Exact(nx), Exact(ny)) => Exact(nx * ny)
case (Exact(Rational.zero), _) => Real.zero
case (_, Exact(Rational.zero)) => Real.zero
case (Exact(Rational.one), _) => y
case (_, Exact(Rational.one)) => x
case _ => Real({p =>
val x0 = x(0).abs + 2
val y0 = y(0).abs + 2
val sx = Real.sizeInBase(x0, 2) + 3
val sy = Real.sizeInBase(y0, 2) + 3
roundUp(Rational(x(p + sy) * y(p + sx), SafeLong.two.pow(p + sx + sy)))
})
}
def **(k: Int): Real = pow(k)
def pow(k: Int): Real = {
@tailrec
def loop(b: Real, k: Int, extra: Real): Real =
if (k == 1)
b * extra
else
loop(b * b, k >>> 1, if ((k & 1) == 1) b * extra else extra)
this match {
case Exact(n) =>
Exact(n.pow(k))
case _ =>
if (k < 0) {
reciprocal.pow(-k)
} else if (k == 0) {
Real.one
} else if (k == 1) {
this
} else {
loop(x, k - 1, x)
}
}
}
def /(y: Real): Real = x * y.reciprocal
def %(y: Real): Real = (x, y) match {
case (Exact(nx), Exact(ny)) => Exact(nx % ny)
case _ => Real({ p =>
val d = x / y
val s = d(2)
val d2 = if (s >= 0) d.floor else d.ceil
(x - d2 * y)(p)
})
}
def /~(y: Real): Real = (x, y) match {
case (Exact(nx), Exact(ny)) => Exact(nx /~ ny)
case _ => Real({ p =>
val d = x / y
val s = d(2)
val d2 = if (s >= 0) d.floor else d.ceil
d2(p)
})
}
def gcd(y: Real): Real = (x, y) match {
case (Exact(nx), Exact(ny)) => Exact(nx gcd ny)
case _ => Real({ p =>
val g = x.toRational(p) gcd y.toRational(p)
roundUp(g * SafeLong.two.pow(p))
})
}
def ceil(): Real = x match {
case Exact(n) => Exact(n.ceil)
case _ => Real({ p =>
val n = x(p)
val t = SafeLong.two.pow(p)
val m = n % t
if (m == 0) n
else if (n.signum >= 0) n + t - m
else n - m
})
}
def floor(): Real = x match {
case Exact(n) => Exact(n.floor)
case _ => Real({ p =>
val n = x(p)
val t = SafeLong.two.pow(p)
val m = n % t
if (n.signum >= 0) n - m else n - t - m
})
}
def round(): Real = x match {
case Exact(n) => Exact(n.round)
case _ => Real({ p =>
val n = x(p)
val t = SafeLong.two.pow(p)
val h = t / 2
val m = n % t
if (m < h) n - m else n - m + t
})
}
def isWhole(): Boolean = x match {
case Exact(n) =>
n.isWhole
case _ =>
val n = x(Real.bits)
val t = SafeLong.two.pow(Real.bits)
(n % t) == 0
}
def sqrt(): Real = Real(p => x(p * 2).sqrt)
def nroot(k: Int): Real = Real(p => x(p * k).nroot(k))
def fpow(r: Rational): Real =
Real({ p =>
val r2 = r.limitToInt
val n = r2.numerator
val d = r2.denominator
x.pow(n.toInt).nroot(d.toInt)(p)
})
// a bit hand-wavy
def fpow(y: Real): Real = y match {
case Exact(n) => x.fpow(n)
case _ => Real({ p =>
x.fpow(Rational(y(p), SafeLong.two.pow(p)))(p)
})
}
override def toString: String = x match {
case Exact(n) => n.toString
case _ => getString(Real.digits)
}
def repr: String = x match {
case Exact(n) => s"Exact(${n.toString})"
case _ => s"Inexact(${toRational})"
}
def getString(d: Int): String = {
val b = Real.digitsToBits(d)
val r = Rational(x(b) * SafeLong.ten.pow(d), SafeLong.two.pow(b))
val m = roundUp(r)
val (sign, str) = m.signum match {
case -1 => ("-", m.abs.toString)
case 0 => ("", "0")
case 1 => ("", m.toString)
}
val i = str.length - d
val s = if (i > 0) {
sign + str.substring(0, i) + "." + str.substring(i)
} else {
sign + "0." + ("0" * -i) + str
}
s.replaceAll("0+$", "").replaceAll("\\.$", "")
}
}
object Real extends RealInstances {
val zero: Real = Exact(Rational.zero)
val one: Real = Exact(Rational.one)
val two: Real = Exact(Rational(2))
val four: Real = Exact(Rational(4))
def apply(f: Int => SafeLong): Real = Inexact(f)
implicit def apply(n: Int): Real = Exact(Rational(n))
implicit def apply(n: Long): Real = Exact(Rational(n))
implicit def apply(n: BigInt): Real = Exact(Rational(n))
implicit def apply(n: SafeLong): Real = Exact(Rational(n))
implicit def apply(n: Rational): Real = Exact(n)
implicit def apply(n: Double): Real = Exact(Rational(n))
implicit def apply(n: BigDecimal): Real = Exact(Rational(n))
def apply(s: String): Real = Exact(Rational(s))
lazy val pi: Real =
Real(16) * atan(Real(Rational(1, 5))) - Real.four * atan(Real(Rational(1, 239)))
lazy val e: Real =
exp(Real.one)
lazy val phi: Real =
(Real.one + Real(5).sqrt) / Real.two
def log(x: Real): Real = {
val t = x(2)
val n = sizeInBase(t, 2) - 3
if (t < 0) throw new ArithmeticException("log of negative number")
else if (t < 4) -log(x.reciprocal)
else if (t < 8) logDr(x)
else logDr(div2n(x, n)) + Real(n) * log2
}
def exp(x: Real): Real = {
val u = x / log2
val n = u(0)
val s = x - Real(n) * log2
if (!n.isValidInt) throw new ArithmeticException("invalid power in exp")
else if (n < 0) div2n(expDr(s), -n.toInt)
else if (n > 0) mul2n(expDr(s), n.toInt)
else expDr(s)
}
def sin(x: Real): Real = {
val z = x / piBy4
val s = roundUp(Rational(z(2), 4))
val y = x - piBy4 * Real(s)
val m = (s % 8).toInt
val n = if (m < 0) m + 8 else m
n match {
case 0 => sinDr(y)
case 1 => sqrt1By2 * (cosDr(y) + sinDr(y))
case 2 => cosDr(y)
case 3 => sqrt1By2 * (cosDr(y) - sinDr(y))
case 4 => -sinDr(y)
case 5 => -sqrt1By2 * (cosDr(y) + sinDr(y))
case 6 => -cosDr(y)
case 7 => -sqrt1By2 * (cosDr(y) - sinDr(y))
}
}
def cos(x: Real): Real = {
val z = x / piBy4
val s = roundUp(Rational(z(2), 4))
val y = x - piBy4 * Real(s)
val m = (s % 8).toInt
val n = if (m < 0) m + 8 else m
n match {
case 0 => cosDr(y)
case 1 => sqrt1By2 * (cosDr(y) - sinDr(y))
case 2 => -sinDr(y)
case 3 => -sqrt1By2 * (cosDr(y) + sinDr(y))
case 4 => -cosDr(y)
case 5 => -sqrt1By2 * (cosDr(y) - sinDr(y))
case 6 => sinDr(y)
case 7 => sqrt1By2 * (cosDr(y) + sinDr(y))
}
}
def tan(x: Real): Real = sin(x) / cos(x)
def atan(x: Real): Real = {
val t = x(2)
val xp1 = x + Real.one
val xm1 = x - Real.one
if (t < -5) atanDr(-x.reciprocal) - piBy2
else if (t == -4) -piBy4 - atanDr(xp1 / xm1)
else if (t < 4) atanDr(x)
else if (t == 4) piBy4 + atanDr(xm1 / xp1)
else piBy2 - atanDr(x.reciprocal)
}
def atan2(y: Real, x: Real): Real = Real({ p =>
var pp = p
var sx = x(pp).signum
var sy = y(pp).signum
// val maxp = p * p
// while (sx == 0 && sy == 0 && pp < maxp) {
while (sx == 0 && sy == 0) {
sx = x(pp).signum
sy = y(pp).signum
pp += 1
}
if (sx > 0) {
atan(y / x)(p)
} else if (sy >= 0 && sx < 0) {
(atan(y / x) + Real.pi)(p)
} else if (sy < 0 && sx < 0) {
(atan(y / x) - Real.pi)(p)
} else if (sy > 0) {
(Real.pi / Real.two)(p)
} else if (sy < 0) {
(-Real.pi / Real.two)(p)
} else {
throw new IllegalArgumentException("atan2(0, 0) is undefined")
// // ugh
// Real.zero
// //sys.error("undefined sx=%s sy=%s" format (sx, sy))
}
})
def asin(x: Real): Real = {
val x0 = x(0)
val s = (Real.one - x * x).sqrt
x0.signum match {
case n if n > 0 => (Real.pi / Real.two) - atan(s / x)
case 0 => atan(x / s)
case _ => (-Real.pi / Real.two) - atan(s / x)
}
}
def acos(x: Real): Real = (Real.pi / Real.two) - asin(x)
def sinh(x: Real): Real = {
val y = exp(x)
(y - y.reciprocal) / Real.two
}
def cosh(x: Real): Real = {
val y = exp(x)
(y + y.reciprocal) / Real.two
}
def tanh(x: Real): Real = {
val y = exp(x);
val y2 = y.reciprocal
(y - y2) / (y + y2)
}
def asinh(x: Real): Real = log(x + (x * x + Real.one).sqrt)
def acosh(x: Real): Real = log(x + (x * x - Real.one).sqrt)
def atanh(x: Real): Real = log((Real.one + x) / (Real.one - x)) / Real.two
def digits: Int = 40
def bits: Int = digitsToBits(digits)
def digitsToBits(n: Int): Int =
spire.math.ceil(n * (spire.math.log(10.0) / spire.math.log(2.0))).toInt + 4
def sizeInBase(n: SafeLong, base: Int): Int = {
def loop(n: SafeLong, acc: Int): Int = if (n <= 1) acc + 1 else loop(n / base, acc + 1)
loop(n.abs, 0)
}
def roundUp(r: Rational): SafeLong = SafeLong(r.round.toBigInt)
def div2n(x: Real, n: Int): Real =
Real(p => if (p >= n) x(p - n) else roundUp(Rational(x(p), SafeLong.two.pow(n))))
def mul2n(x: Real, n: Int): Real =
Real(p => x(p + n))
lazy val piBy2 = div2n(pi, 1)
lazy val piBy4 = div2n(pi, 2)
lazy val log2 = div2n(logDrx(Real.two.reciprocal), 1)
lazy val sqrt1By2 = Real.two.reciprocal.sqrt
def accumulate(total: SafeLong, xs: Stream[SafeLong], cs: Stream[Rational]): SafeLong = {
(xs, cs) match {
case (_, Stream.Empty) => total
case (Stream.Empty, _) => sys.error("nooooo")
case (x #:: xs, c #:: cs) =>
val t = roundUp(c * Rational(x))
if (t == 0) total else accumulate(total + t, xs, cs)
}
}
private[spire] def powerSeries(ps: Stream[Rational], terms: Int => Int, x: Real): Real = {
Real({p =>
val t = terms(p)
val l2t = 2 * sizeInBase(SafeLong(t) + 1, 2) + 6
val p2 = p + l2t
val xr = x(p2)
val xn = SafeLong.two.pow(p2)
if (xn == 0) sys.error("oh no")
def g(yn: SafeLong): SafeLong = roundUp(Rational(yn * xr, xn))
val num = accumulate(SafeLong.zero, Stream.iterate(xn)(g), ps.take(t))
val denom = SafeLong.two.pow(l2t)
roundUp(Rational(num, denom))
})
}
private[spire] def accSeq(f: (Rational, SafeLong) => Rational): Stream[Rational] = {
def loop(r: Rational, n: SafeLong): Stream[Rational] =
r #:: loop(f(r, n), n + 1)
loop(Rational.one, SafeLong.one)
}
def expDr(x: Real): Real =
powerSeries(accSeq((r, n) => r / n), n => n, x)
def logDr(x: Real): Real = {
val y = (x - Real.one) / x
y * logDrx(y)
}
def logDrx(x: Real): Real = {
powerSeries(Stream.from(1).map(n => Rational(1, n)), _ + 1, x)
}
def sinDr(x: Real): Real =
x * powerSeries(accSeq((r, n) => -r * Rational(1, 2*n*(2*n+1))), n => n, x * x)
def cosDr(x: Real): Real =
powerSeries(accSeq((r, n) => -r * Rational(1, 2*n*(2*n-1))), n => n, x * x)
def atanDr(x: Real): Real = {
val y = x * x + Real(1)
(x / y) * atanDrx((x * x) / y)
}
def atanDrx(x: Real): Real =
//powerSeries(accSeq((r, n) => r * (Rational(2*n, 2*n + 1))), _ + 1, x)
powerSeries(accSeq((r, n) => r * (Rational(2*n, 2*n + 1))), _ * 2, x)
case class Exact(n: Rational) extends Real {
def apply(p: Int): SafeLong = Real.roundUp(Rational(2).pow(p) * n)
}
case class Inexact(f: Int => SafeLong) extends Real {
@volatile private[spire] var memo: Option[(Int, SafeLong)] = None
def apply(p: Int): SafeLong = memo match {
case Some((bits, value)) if bits >= p =>
Real.roundUp(Rational(value, SafeLong(2).pow(bits - p)))
case _ =>
val result = f(p)
memo = Some((p, result))
result
}
}
}
trait RealInstances {
implicit final val algebra = new RealAlgebra
import NumberTag._
implicit final val RealTag = new LargeTag[Real](Exact, Real.zero)
}
@SerialVersionUID(0L)
class RealAlgebra extends RealIsFractional {}
trait RealIsFractional extends Fractional[Real] with Order[Real] with Signed[Real] with Trig[Real] {
def abs(x: Real): Real = x.abs
def signum(x: Real): Int = x.signum
override def eqv(x: Real, y: Real): Boolean = x === y
def compare(x: Real, y: Real): Int = x compare y
def zero: Real = Real.zero
def one: Real = Real.one
def negate(x: Real): Real = -x
def plus(x: Real, y: Real): Real = x + y
override def minus(x: Real, y: Real): Real = x - y
def times(x: Real, y: Real): Real = x * y
def gcd(x: Real, y: Real): Real = x gcd y
def quot(x: Real, y: Real): Real = x /~ y
def mod(x: Real, y: Real): Real = x % y
override def reciprocal(x: Real): Real = x.reciprocal
def div(x: Real, y: Real): Real = x / y
override def sqrt(x: Real): Real = x.sqrt
def nroot(x: Real, k: Int): Real = x.nroot(k)
def fpow(x: Real, y: Real): Real = x fpow y
def acos(a: Real): Real = Real.acos(a)
def asin(a: Real): Real = Real.asin(a)
def atan(a: Real): Real = Real.atan(a)
def atan2(y: Real, x: Real): Real = Real.atan2(y, x)
def cos(a: Real): Real = Real.cos(a)
def cosh(x: Real): Real = Real.cosh(x)
def e: Real = Real.e
def exp(x: Real): Real = Real.exp(x)
def expm1(x: Real): Real = Real.exp(x) - Real.one
def log(x: Real): Real = Real.log(x)
def log1p(x: Real): Real = Real.log(Real.one + x)
def pi: Real = Real.pi
def sin(x: Real): Real = Real.sin(x)
def sinh(x: Real): Real = Real.sinh(x)
def tan(x: Real): Real = Real.tan(x)
def tanh(x: Real): Real = Real.tanh(x)
def toDegrees(a: Real): Real = a / (Real.two * Real.pi) * Real(360)
def toRadians(a: Real): Real = a / Real(360) * (Real.two * Real.pi)
def ceil(x: Real): Real = x.ceil
def floor(x: Real): Real = x.floor
def isWhole(x: Real): Boolean = x.isWhole
def round(x: Real): Real = x.round
def toByte(x: Real): Byte = x.toRational.toByte
def toInt(x: Real): Int = x.toRational.toInt
def toShort(x: Real): Short = x.toRational.toShort
def toLong(x: Real): Long = x.toRational.toLong
def toFloat(x: Real): Float = x.toRational.toFloat
def toDouble(x: Real): Double = x.toRational.toDouble
def toBigInt(x: Real): BigInt = x.toRational.toBigInt
def toBigDecimal(x: Real): BigDecimal = x.toRational.toBigDecimal(java.math.MathContext.DECIMAL64)
def toRational(x: Real): Rational = x.toRational
def toAlgebraic(x: Real): Algebraic = Algebraic(x.toRational) //FIXME
def toReal(x: Real): Real = x
def toNumber(x: Real): Number = Number(x.toRational)
def toString(x: Real): String = x.toString
def toType[B](x: Real)(implicit ev: ConvertableTo[B]): B =
ev.fromReal(x)
def fromByte(n: Byte): Real = Real(n)
def fromShort(n: Short): Real = Real(n)
def fromFloat(n: Float): Real = Real(n)
def fromLong(n: Long): Real = Real(n)
def fromBigInt(n: BigInt): Real = Real(n)
def fromBigDecimal(n: BigDecimal): Real = Real(n)
def fromRational(n: Rational): Real = Real(n)
def fromAlgebraic(n: Algebraic): Real = n.evaluateWith[Real]
def fromReal(n: Real): Real = n
def fromType[B](b: B)(implicit ev: ConvertableFrom[B]): Real =
ev.toReal(b)
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy