using System; using System.Collections.Generic; using System.Numerics; using System.Text; namespace Sandbox { public class Rational { public Rational(BigInteger num, BigInteger denom) { if (denom.IsZero) throw new ArgumentOutOfRangeException("denom"); if (denom < 0) { num = -num; denom = -denom; } var g = _Gcd(num, denom); Num = num / g; Denom = denom / g; } public BigInteger Num { get; private set; } public BigInteger Denom { get; private set; } public IEnumerable Approximations { get { var continuedFrac = new List(); var sgn = Num < 0 ? -1 : 1; var n = sgn * Num; var d = Denom; while (d > 0) { continuedFrac.Add(n / d); // n/d := 1 / (n/d % 1) = d / (n % d) var tmp = d; d = n % d; n = tmp; // TODO This can be improved Rational approx = continuedFrac[continuedFrac.Count - 1]; for (int i = continuedFrac.Count - 2; i >= 0; i--) approx = continuedFrac[i] + 1 / approx; yield return sgn * approx; } } } public Interval Bound(BigInteger targetDenominator) { var approximations = Approximations.GetEnumerator(); if (!approximations.MoveNext()) return this; var a = approximations.Current; Rational b = a; while (approximations.MoveNext()) { a = b; b = approximations.Current; if (b.Denom > targetDenominator) { break; } } return a < b ? new Interval(a, b) : new Interval(b, a); } public override string ToString() { return string.Format(Denom == 1 ? "{0}" : "{0}/{1}", Num, Denom); } public override int GetHashCode() { return Num.GetHashCode() * 37 + Denom.GetHashCode(); } public override bool Equals(object obj) { return obj is Rational && (this == (Rational)obj); } public static implicit operator Rational(int i) { return new Rational(i, 1); } public static implicit operator Rational(long i) { return new Rational(i, 1); } public static implicit operator Rational(BigInteger i) { return new Rational(i, 1); } public static implicit operator double(Rational r) { return (double)r.Num / (double)r.Denom; } public static Rational operator -(Rational a) { return new Rational(-a.Num, a.Denom); } public static Rational operator +(Rational a, Rational b) { return new Rational(a.Num * b.Denom + a.Denom * b.Num, a.Denom * b.Denom)$public static Rational operator -(Rational a, Rational b) { return new Rational(a.Num * b.Denom - a.Denom * b.Num, a.Denom * b.Denom)$ public static Rational operator *(Rational a, Rational b) { return new Rational(a.Num * b.Num, a.Denom * b.Denom); } public static Rational operator /(Rational a, Rational b) { return new Rational(a.Num * b.Denom, a.Denom * b.Num); } public static bool operator <(Rational a, Rational b) { return a.Num * b.Denom < a.Denom * b.Num; } public static bool operator >(Rational a, Rational b) { return a.Num * b.Denom > a.Denom * b.Num; } public static bool operator <=(Rational a, Rational b) { return a.Num * b.Denom <= a.Denom * b.Num; } public static bool operator >=(Rational a, Rational b) { return a.Num * b.Denom >= a.Denom * b.Num; } public static bool operator ==(Rational a, Rational b) { return a.Num == b.Num && a.Denom == b.Denom; } public static bool operator !=(Rational a, Rational b) { return !(a == b); } private static BigInteger _Gcd(BigInteger a, BigInteger b) { if (a < 0) a = -a; if (b < 0) b = -b; while (!a.IsZero) { var tmp = b % a; b = a; a = tmp; } return b; } } // https://arxiv.org/pdf/0708.3721.pdf class BigDecimalInterval { static readonly BigInteger Precision = BigInteger.Pow(2, 256); static readonly Rational Epsilon = new Rational(1, Precision); public BigDecimalInterval(Rational r) : this(r, Precision) { } private BigDecimalInterval(Rational r, BigInteger precision) { _Precision = precision; if (r.Denom <= precision) { Lower = Upper = r; return; } // TODO Should probably prefer continuants over treating the precision as a desired denominator r = r * precision; if (r.Denom == 1) { Lower = Upper = r; } else if (r.Num > 0) { var floor = r.Num / r.Denom; Lower = new Rational(floor, precision); Upper = new Rational(floor + 1, precision); } else { var ceil = -((-r.Num) / r.Denom); Lower = new Rational(ceil - 1, precision); Upper = new Rational(ceil, precision); } } public BigDecimalInterval(Rational lower, Rational upper) : this(lower, upper, Precision) { } private BigDecimalInterval(Rational lower, Rational upper, BigInteger precision) { if (lower > upper) throw new ArgumentOutOfRangeException(); _Precision = precision; Lower = new BigDecimalInterval(lower, precision).Lower; Upper = new BigDecimalInterval(upper, precision).Upper; } private BigInteger _Precision { get; set; } public Rational Lower { get; private set; } public Rational Upper { get; private set; } public override string ToString() => string.Format("[{0} ~= {1}, {2} ~= {3} width {4}]", Lower, (double)Lower, Upper, (double)Upper, (double)(Upper - Lower)); #region Operators public static implicit operator BigDecimalInterval(int r) => new BigDecimalInterval(r); public static implicit operator BigDecimalInterval(long r) => new BigDecimalInterval(r); public static implicit operator BigDecimalInterval(BigInteger r) => new BigDecimalInterval(r); public static implicit operator BigDecimalInterval(Rational r) => new BigDecimalInterval(r); public static BigDecimalInterval operator -(BigDecimalInterval a) => new BigDecimalInterval(-a.Upper, -a.Lower, a._Precision); public static BigDecimalInterval operator +(BigDecimalInterval a, BigDecimalInterval b) => new BigDecimalInterval(a.Lower + b.Lower, a.Upper + b.Upper, Max(a._Precision, b._Precision)); public static BigDecimalInterval operator -(BigDecimalInterval a, BigDecimalInterval b) => new BigDecimalInterval(a.Lower - b.Upper, a.Upper - b.Lower, Max(a._Precision, b._Precision)); // Assumes a > 0 public static BigDecimalInterval operator *(int a, BigDecimalInterval b) => new BigDecimalInterval(a * b.Lower, a * b.Upper, b._Precision); public static BigDecimalInterval operator /(BigDecimalInterval a, int b) => new BigDecimalInterval(a.Lower / b, a.Upper / b, a._Precision); #endregion private static BigInteger Max(BigInteger a, BigInteger b) => a > b ? a : b; public static BigDecimalInterval Atan(Rational x) { // atan(x) = \sum_{i=0}^\infty (-1)^i x^{2i+1} / (2i+1) Rational lower = x; Rational upper = x; var intermediatePrecision = Precision * Precision; var epsilon = new Rational(1, intermediatePrecision); BigDecimalInterval sum = 0; BigDecimalInterval term = new BigDecimalInterval(x, intermediatePrecision); var x2 = x * x; for (int i = 0; true; i++) { sum += term / (2 * i + 1); term = new BigDecimalInterval(term.Upper * -x2, term.Lower * -x2, intermediatePrecision); if (i % 2 == 0) upper = sum.Upper; else lower = sum.Lower; if (term.Lower.Abs() <= epsilon && term.Upper.Abs() <= epsilon) break; } return new BigDecimalInterval(lower, upper, Precision); } private static readonly BigDecimalInterval Pi = 16 * Atan(new Rational(1, 5)) - 4 * Atan(new Rational(1, 239)); private static BigDecimalInterval Sin(Rational x) { if (x < 0) return -Sin(-x); // sin(x) = sum_{i=0}^\infty (-1)^{i} x^{2i+1} / (2i+1)! Rational lower = x; Rational upper = x; var intermediatePrecision = Precision * Precision; var epsilon = new Rational(1, intermediatePrecision); BigDecimalInterval term = new BigDecimalInterval(x, intermediatePrecision); BigDecimalInterval sum = term; var x2 = x * x; for (int i = 1; true; i++) { var m = x2 / -((2 * i) * (2 * i + 1)); term = new BigDecimalInterval(term.Upper * m, term.Lower * m, intermediatePrecision); sum += term; if (i % 2 == 0) upper = sum.Upper; else lower = sum.Lower; if (term.Lower.Abs() <= epsilon && term.Upper.Abs() <= epsilon) break; } return new BigDecimalInterval(lower, upper, Precision); } private static BigDecimalInterval Cos(Rational x) { if (x < 0) x = -x; // cos(x) = sum_{i=0}^\infty (-1)^i x^{2i} / (2i)! Rational lower = 1; Rational upper = 1; var intermediatePrecision = Precision * Precision; var epsilon = new Rational(1, intermediatePrecision); BigDecimalInterval term = new BigDecimalInterval(1, intermediatePrecision); BigDecimalInterval sum = term; var x2 = x * x; for (int i = 1; true; i++) { var m = x2 / -((2 * i - 1) * (2 * i)); term = new BigDecimalInterval(term.Upper * m, term.Lower * m, intermediatePrecision); sum += term; if (i % 2 == 0) upper = sum.Upper; else lower = sum.Lower; if (term.Lower.Abs() <= epsilon && term.Upper.Abs() <= epsilon) break; } return new BigDecimalInterval(lower, upper, Precision); } public static BigDecimalInterval Tan(BigDecimalInterval x) { if (x.Upper < 0) return -Tan(-x); if (x.Lower >= Pi.Lower / -2 && x.Upper <= Pi.Lower / 2) { var lb = Sin(x.Lower).Lower / Cos(x.Lower).Upper; var ub = Sin(x.Upper).Upper / Cos(x.Upper).Lower; return new BigDecimalInterval(lb, ub); } throw new ArgumentOutOfRangeException(nameof(x)); } internal static void Main() { BigDecimalInterval x = 1; for (int i = 1; i < 1001; i++) { // Tan is fine for -pi/2 to pi/2. while (x.Upper > Pi.Lower / 2) x -= Pi; while (x.Lower < Pi.Lower / -2) x += Pi; x = Tan(x); double lb = x.Lower; double ub = x.Upper; double avg = (lb + ub) / 2; if (lb == ub) Console.WriteLine("{0}: {1}", i, lb); else { var sl = lb.ToString("F16"); var su = ub.ToString("F16"); var sa = avg.ToString("F16"); var sb = new StringBuilder(); for (int j = 0; j < sa.Length; j++) { if (j < sl.Length && j < su.Length && sl[j] == sa[j] && su[j] == sa[j]) sb.Append(sa[j]); else { sb.Append('(').Append(sa[j]).Append(')'); break; } } Console.WriteLine("{0}: {1}", i, sb); } } } } }