BigRational.cs

C# Arbitrary Precision Signed Big Rational Numbers

using System;
using System.Diagnostics;
using System.Globalization;
using System.Numerics;
using System.Runtime.InteropServices;
using System.Text;
[DebuggerDisplay("{" + nameof(DDisplay) + "}")]
[Serializable]
public struct BigRational : IComparable, IComparable<BigRational>, IEquatable<BigRational>
{
    [StructLayout(LayoutKind.Explicit)]
    internal struct DoubleUlong
    {
        [FieldOffset(0)] public double dbl;
        [FieldOffset(0)] public ulong  uu;
    }
    private const int DoubleMaxScale = 308;
    public static BigRational Pi = new BigRational(
        "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798162478513934506898440362801792706010179987216806726188740140466033567581311679376075335101609659171030644576233653027450257182803484658351860927270133809030914436823660262931162576284703194589395221866245992710817555393680237554917047871708932985106840785074833639247080859264327721882027979677397953754604196915619381410505600288856897761875941052867609089114345150157869223684881643245943313338421018485091403977277400743970527492816321894223953257584787737337170568053925027217102844351208765657302025589127695185039186644597240030541171074757870137431100579097277905612641495178817964173941740654985445918326928220945355416048444887050935562696866696019631573868714587428709669938320262709342763");
    public static           BigRational   E               = GetE(MaxFactorials);
    private static readonly BigInteger    DoublePrecision = BigInteger.Pow(10, DoubleMaxScale);
    private static readonly BigInteger    DoubleMaxValue  = (BigInteger) double.MaxValue;
    private static readonly BigInteger    DoubleMinValue  = (BigInteger) double.MinValue;
    private static          BigRational[] Factorials;
    static BigRational()
    {
    }
    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    private string DDisplay => AsDecimal(this);
    [StructLayout(LayoutKind.Explicit)]
    internal struct DecimalUInt32
    {
        [FieldOffset(0)] public decimal dec;
        [FieldOffset(0)] public int     flags;
    }
    private const int DecimalScaleMask = 0x00FF0000;
    private const int DecimalSignMask  = unchecked((int) 0x80000000);
    /// <summary>
    ///     Change here if more then 2048 bits are specified
    /// </summary>
    private const int DecimalMaxScale = 2048 / 64 * 20;
    private const           int        MaxFactorials    = 100;
    private static readonly BigInteger DecimalPrecision = BigInteger.Pow(10, DecimalMaxScale);
    private static readonly BigInteger DecimalMaxValue  = (BigInteger) decimal.MaxValue;
    private static readonly BigInteger DecimalMinValue  = (BigInteger) decimal.MinValue;
    private const           string     Solidus          = @"/";
    public static BigRational Zero
    {
        get;
    } = new BigRational(BigInteger.Zero);
    public static BigRational One
    {
        get;
    } = new BigRational(BigInteger.One);
    public static BigRational MinusOne
    {
        get;
    } = new BigRational(BigInteger.MinusOne);
    public int Sign => Numerator.Sign;
    public BigInteger Numerator
    {
        get;
        private set;
    }
    public BigInteger Denominator
    {
        get;
        private set;
    }
    public BigInteger GetWholePart => BigInteger.Divide(Numerator, Denominator);
    public bool IsFractionalPart
    {
        get
        {
            var fp = GetFractionPart;
            return fp.Numerator != 0 || fp.Denominator != 1;
        }
    }
    public BigInteger GetUnscaledAsDecimal => Numerator * DecimalPrecision / Denominator;
    public BigInteger Remainder            => Numerator                    % Denominator;
    public int DecimalPlaces
    {
        get
        {
            var a       = GetUnscaledAsDecimal;
            var dPlaces = 0;
            if (a.Sign == 0)
                return 1;
            if (a.Sign < 0)
                try
                {
                    a = -a;
                }
                catch (Exception ex)
                {
                    return 0;
                }
            var biRadix = new BigInteger(10);
            while (a > 0)
                try
                {
                    a /= biRadix;
                    dPlaces++;
                }
                catch (Exception ex)
                {
                    break;
                }
            return dPlaces;
        }
    }
    public static string AsDecimal(BigRational value)
    {
        var asd = new BigDecimal(value);
        return asd.ToString();
    }
    public static string CleanAsDecimal(BigRational value)
    {
        var fpas = AsDecimal(value);
        var rs   = fpas.Reverse();
        var fas  = "";
        foreach (var c in rs)
            if (c == '0')
                continue;
            else
                fas += c;
        return fas.Reverse();
    }
    public BigRational GetFractionPart
    {
        get
        {
            var rem = BigInteger.Remainder(Numerator, Denominator);
            return new BigRational(rem, Denominator);
        }
    }
    public override bool Equals(object obj)
    {
        if (obj == null)
            return false;
        if (!(obj is BigRational))
            return false;
        return Equals((BigRational) obj);
    }
    public override int GetHashCode()
    {
        return (Numerator / Denominator).GetHashCode();
    }
    int IComparable.CompareTo(object obj)
    {
        if (obj == null)
            return 1;
        if (!(obj is BigRational))
            throw new ArgumentException();
        return Compare(this, (BigRational) obj);
    }
    public int CompareTo(BigRational other)
    {
        return Compare(this, other);
    }
    public bool Equals(BigRational other)
    {
        if (Denominator == other.Denominator)
            return Numerator == other.Numerator;
        return Numerator * other.Denominator == Denominator * other.Numerator;
    }
    public BigRational(BigInteger numerator)
    {
        Numerator   = numerator;
        Denominator = BigInteger.One;
    }
    public BigRational(string n, string d)
    {
        Numerator   = new BigInteger().BigIntegerBase10(n);
        Denominator = new BigInteger().BigIntegerBase10(d);
    }
    public BigRational(string value)
    {
        if (!value.ContainsOnly("0123456789+-.eE"))
            throw new Exception(
                $"Input value must only contain these '0123456789+-.eE', value'{value}");
        var v1 = new BigDecimal(value);
        var (unscaledValue, scale) = v1.ToByteArrays();
        if (v1 == BigDecimal.Zero)
        {
            this = Zero;
            return;
        }
        Numerator   = new BigInteger(unscaledValue);
        Denominator = BigInteger.Pow(10, BitConverter.ToInt32(scale, 0));
        Simplify();
    }
    public static bool TryParse(string parse, out BigRational result)
    {
        result = default;
        if (!parse.ContainsOnly("0123456789+-.eE"))
            throw new Exception(
                $"Input value must only contain these '0123456789+-.eE', value'{parse}");
        try
        {
            result = new BigRational(parse);
        }
        catch
        {
            return false;
        }
        return true;
    }
    public BigRational(double value) : this((decimal) value)
    {
    }
    public BigRational(BigDecimal value)
    {
        var bits = value.ToByteArrays();
        if (value == BigDecimal.Zero)
        {
            this = Zero;
            return;
        }
        Numerator   = new BigInteger(bits.unscaledValue);
        Denominator = BigInteger.Pow(10, BitConverter.ToInt32(bits.scale, 0));
        Simplify();
    }
    public BigRational(decimal value)
    {
        var bits = decimal.GetBits(value);
        if (bits                                              == null || bits.Length != 4 ||
            (bits[3] & ~(DecimalSignMask | DecimalScaleMask)) != 0    ||
            (bits[3] & DecimalScaleMask)                      > 28 << 16)
            throw new ArgumentException();
        if (value == decimal.Zero)
        {
            this = Zero;
            return;
        }
        var ul = ((ulong) (uint) bits[2] << 32) | (uint) bits[1];
        Numerator = (new BigInteger(ul) << 32) | (uint) bits[0];
        var isNegative = (bits[3] & DecimalSignMask) != 0;
        if (isNegative)
            Numerator = BigInteger.Negate(Numerator);
        var scale = (bits[3] & DecimalScaleMask) >> 16;
        Denominator = BigInteger.Pow(10, scale);
        Simplify();
    }
    public BigRational(BigInteger numerator, BigInteger denominator)
    {
        if (denominator.Sign == 0)
            throw new DivideByZeroException();
        if (numerator.Sign == 0)
        {
            Numerator   = BigInteger.Zero;
            Denominator = BigInteger.One;
        }
        else if (denominator.Sign < 0)
        {
            Numerator   = BigInteger.Negate(numerator);
            Denominator = BigInteger.Negate(denominator);
        }
        else
        {
            Numerator   = numerator;
            Denominator = denominator;
        }
        Simplify();
    }
    public BigRational(BigInteger whole, BigInteger numerator, BigInteger denominator)
    {
        if (denominator.Sign == 0)
            throw new DivideByZeroException();
        if (numerator.Sign == 0 && whole.Sign == 0)
        {
            Numerator   = BigInteger.Zero;
            Denominator = BigInteger.One;
        }
        else if (denominator.Sign < 0)
        {
            Denominator = BigInteger.Negate(denominator);
            Numerator   = BigInteger.Negate(whole) * Denominator + BigInteger.Negate(numerator);
        }
        else
        {
            Denominator = denominator;
            Numerator   = whole * denominator + numerator;
        }
        Simplify();
    }
    public static BigRational Abs(BigRational r)
    {
        return r.Numerator.Sign < 0
            ? new BigRational(BigInteger.Abs(r.Numerator), r.Denominator)
            : r;
    }
    public static BigRational Negate(BigRational r)
    {
        return new BigRational(BigInteger.Negate(r.Numerator), r.Denominator);
    }
    public static BigRational Invert(BigRational r)
    {
        return new BigRational(r.Denominator, r.Numerator);
    }
    public static BigRational Add(BigRational x, BigRational y)
    {
        return x + y;
    }
    public static BigRational Subtract(BigRational x, BigRational y)
    {
        return x - y;
    }
    public static BigRational Multiply(BigRational x, BigRational y)
    {
        return x * y;
    }
    public static BigRational Divide(BigRational dividend, BigRational divisor)
    {
        return dividend / divisor;
    }
    public static BigRational DivRem(BigRational dividend,
        BigRational                              divisor,
        out BigRational                          remainder)
    {
        var ad = dividend.Numerator   * divisor.Denominator;
        var bc = dividend.Denominator * divisor.Numerator;
        var bd = dividend.Denominator * divisor.Denominator;
        remainder = new BigRational(ad % bc, bd);
        return new BigRational(ad, bc);
    }
    public static BigInteger LeastCommonDenominator(BigRational x, BigRational y)
    {
        return x.Denominator * y.Denominator /
               BigInteger.GreatestCommonDivisor(x.Denominator, y.Denominator);
    }
    public static int Compare(BigRational r1, BigRational r2)
    {
        return BigInteger.Compare(r1.Numerator * r2.Denominator, r2.Numerator * r1.Denominator);
    }
    public static bool operator ==(BigRational x, BigRational y)
    {
        return Compare(x, y) == 0;
    }
    public static bool operator !=(BigRational x, BigRational y)
    {
        return Compare(x, y) != 0;
    }
    public static bool operator <(BigRational x, BigRational y)
    {
        return Compare(x, y) < 0;
    }
    public static bool operator <=(BigRational x, BigRational y)
    {
        return Compare(x, y) <= 0;
    }
    public static bool operator >(BigRational x, BigRational y)
    {
        return Compare(x, y) > 0;
    }
    public static bool operator >=(BigRational x, BigRational y)
    {
        return Compare(x, y) >= 0;
    }
    public static BigRational operator +(BigRational r)
    {
        return r;
    }
    public static BigRational operator -(BigRational r)
    {
        return new BigRational(-r.Numerator, r.Denominator);
    }
    public static BigRational operator ++(BigRational r)
    {
        return r + One;
    }
    public static BigRational operator --(BigRational r)
    {
        return r - One;
    }
    public static BigRational operator +(BigRational r1, BigRational r2)
    {
        return new BigRational(r1.Numerator * r2.Denominator + r1.Denominator * r2.Numerator,
            r1.Denominator * r2.Denominator);
    }
    public static BigRational operator -(BigRational r1, BigRational r2)
    {
        return new BigRational(r1.Numerator * r2.Denominator - r1.Denominator * r2.Numerator,
            r1.Denominator * r2.Denominator);
    }
    public static BigRational operator *(BigRational r1, BigRational r2)
    {
        return new BigRational(r1.Numerator * r2.Numerator, r1.Denominator * r2.Denominator);
    }
    public static BigRational operator /(BigRational r1, BigRational r2)
    {
        return new BigRational(r1.Numerator * r2.Denominator, r1.Denominator * r2.Numerator);
    }
    public static BigRational operator %(BigRational r1, BigRational r2)
    {
        return new BigRational(r1.Numerator * r2.Denominator % (r1.Denominator * r2.Numerator),
            r1.Denominator                                   * r2.Denominator);
    }
    public static explicit operator sbyte(BigRational value)
    {
        return (sbyte) BigInteger.Divide(value.Numerator, value.Denominator);
    }
    public static explicit operator ushort(BigRational value)
    {
        return (ushort) BigInteger.Divide(value.Numerator, value.Denominator);
    }
    public static explicit operator uint(BigRational value)
    {
        return (uint) BigInteger.Divide(value.Numerator, value.Denominator);
    }
    public static explicit operator ulong(BigRational value)
    {
        return (ulong) BigInteger.Divide(value.Numerator, value.Denominator);
    }
    public static explicit operator byte(BigRational value)
    {
        return (byte) BigInteger.Divide(value.Numerator, value.Denominator);
    }
    public static explicit operator short(BigRational value)
    {
        return (short) BigInteger.Divide(value.Numerator, value.Denominator);
    }
    public static explicit operator int(BigRational value)
    {
        return (int) BigInteger.Divide(value.Numerator, value.Denominator);
    }
    public static explicit operator long(BigRational value)
    {
        return (long) BigInteger.Divide(value.Numerator, value.Denominator);
    }
    public static explicit operator BigInteger(BigRational value)
    {
        return BigInteger.Divide(value.Numerator, value.Denominator);
    }
    public static explicit operator float(BigRational value)
    {
        return (float) (double) value;
    }
    public static explicit operator double(BigRational value)
    {
        if (SafeCastToDouble(value.Numerator) && SafeCastToDouble(value.Denominator))
            return (double) value.Numerator / (double) value.Denominator;
        var denormalized = value.Numerator * DoublePrecision / value.Denominator;
        if (denormalized.IsZero)
            return value.Sign < 0
                ? BitConverter.Int64BitsToDouble(unchecked((long) 0x8000000000000000))
                : 0d;
        double result   = 0;
        var    isDouble = false;
        var    scale    = DoubleMaxScale;
        while (scale > 0)
        {
            if (!isDouble)
                if (SafeCastToDouble(denormalized))
                {
                    result   = (double) denormalized;
                    isDouble = true;
                }
                else
                {
                    denormalized = denormalized / 10;
                }
            result = result / 10;
            scale--;
        }
        if (!isDouble)
            return value.Sign < 0 ? double.NegativeInfinity : double.PositiveInfinity;
        return result;
    }
    public static explicit operator BigDecimal(BigRational value)
    {
        var denormalized = value.Numerator * DecimalPrecision / value.Denominator;
        return new BigDecimal(denormalized, DecimalMaxScale);
    }
    public static explicit operator decimal(BigRational value)
    {
        if (SafeCastToDecimal(value.Numerator) && SafeCastToDecimal(value.Denominator))
            return (decimal) value.Numerator / (decimal) value.Denominator;
        var denormalized = value.Numerator * DecimalPrecision / value.Denominator;
        if (denormalized.IsZero)
            return decimal.Zero;
        for (var scale = DecimalMaxScale; scale >= 0; scale--)
            if (!SafeCastToDecimal(denormalized))
            {
                denormalized /= 10;
            }
            else
            {
                var dec = new DecimalUInt32();
                dec.dec   = (decimal) denormalized;
                dec.flags = (dec.flags & ~DecimalScaleMask) | (scale << 16);
                return dec.dec;
            }
        throw new OverflowException();
    }
    public static implicit operator BigRational(sbyte value)
    {
        return new BigRational((BigInteger) value);
    }
    public static implicit operator BigRational(ushort value)
    {
        return new BigRational((BigInteger) value);
    }
    public static implicit operator BigRational(uint value)
    {
        return new BigRational((BigInteger) value);
    }
    public static implicit operator BigRational(ulong value)
    {
        return new BigRational((BigInteger) value);
    }
    public static implicit operator BigRational(byte value)
    {
        return new BigRational((BigInteger) value);
    }
    public static implicit operator BigRational(short value)
    {
        return new BigRational((BigInteger) value);
    }
    public static implicit operator BigRational(int value)
    {
        return new BigRational((BigInteger) value);
    }
    public static implicit operator BigRational(long value)
    {
        return new BigRational((BigInteger) value);
    }
    public static implicit operator BigRational(BigInteger value)
    {
        return new BigRational(value);
    }
    public static implicit operator BigRational(string value)
    {
        return new BigRational(value);
    }
    public static implicit operator BigRational(float value)
    {
        return new BigRational(value);
    }
    public static implicit operator BigRational(double value)
    {
        return new BigRational(value);
    }
    public static implicit operator BigRational(decimal value)
    {
        return new BigRational(value);
    }
    public static implicit operator BigRational(BigDecimal value)
    {
        return new BigRational(value);
    }
    private void Simplify()
    {
        if (Numerator == BigInteger.Zero)
            Denominator = BigInteger.One;
        var gcd = BigInteger.GreatestCommonDivisor(Numerator, Denominator);
        if (gcd > BigInteger.One)
        {
            Numerator   = Numerator   / gcd;
            Denominator = Denominator / gcd;
        }
    }
    private static bool SafeCastToDouble(BigInteger value)
    {
        return DoubleMinValue <= value && value <= DoubleMaxValue;
    }
    private static bool SafeCastToDecimal(BigInteger value)
    {
        return DecimalMinValue <= value && value <= DecimalMaxValue;
    }
    private static void SplitDoubleIntoParts(double dbl,
        out int                                     sign,
        out int                                     exp,
        out ulong                                   man,
        out bool                                    isFinite)
    {
        DoubleUlong du;
        du.uu  = 0;
        du.dbl = dbl;
        sign   = 1 - ((int) (du.uu >> 62) & 2);
        man    = du.uu               & 0x000FFFFFFFFFFFFF;
        exp    = (int) (du.uu >> 52) & 0x7FF;
        if (exp == 0)
        {
            isFinite = true;
            if (man != 0)
                exp = -1074;
        }
        else if (exp == 0x7FF)
        {
            isFinite = false;
            exp      = int.MaxValue;
        }
        else
        {
            isFinite =  true;
            man      |= 0x0010000000000000;
            exp      -= 1075;
        }
    }
    public static double GetDoubleFromParts(int sign, int exp, ulong man)
    {
        DoubleUlong du;
        du.dbl = 0;
        if (man == 0)
        {
            du.uu = 0;
        }
        else
        {
            var cbitShift = CbitHighZero(man) - 11;
            if (cbitShift < 0)
                man >>= -cbitShift;
            else
                man <<= cbitShift;
            exp += 1075;
            if (exp >= 0x7FF)
            {
                du.uu = 0x7FF0000000000000;
            }
            else if (exp <= 0)
            {
                exp--;
                if (exp < -52)
                    du.uu = 0;
                else
                    du.uu = man >> -exp;
            }
            else
            {
                du.uu = (man & 0x000FFFFFFFFFFFFF) | ((ulong) exp << 52);
            }
        }
        if (sign < 0)
            du.uu |= 0x8000000000000000;
        return du.dbl;
    }
    private static int CbitHighZero(ulong uu)
    {
        if ((uu & 0xFFFFFFFF00000000) == 0)
            return 32 + CbitHighZero((uint) uu);
        return CbitHighZero((uint) (uu >> 32));
    }
    private static int CbitHighZero(uint u)
    {
        if (u == 0)
            return 32;
        var cbit = 0;
        if ((u & 0xFFFF0000) == 0)
        {
            cbit +=  16;
            u    <<= 16;
        }
        if ((u & 0xFF000000) == 0)
        {
            cbit +=  8;
            u    <<= 8;
        }
        if ((u & 0xF0000000) == 0)
        {
            cbit +=  4;
            u    <<= 4;
        }
        if ((u & 0xC0000000) == 0)
        {
            cbit +=  2;
            u    <<= 2;
        }
        if ((u & 0x80000000) == 0)
            cbit += 1;
        return cbit;
    }
    private static (BigRational High, BigRational Low) SqrtLimits(BigInteger number)
    {
        if (number == BigInteger.Zero) return (0, 0);
        var high = number >> 1;
        var low  = BigInteger.Zero;
        while (high > low + 1)
        {
            var n = (high + low) >> 1;
            var p = n * n;
            if (number < p)
                high = n;
            else if (number > p)
                low = n;
            else
                break;
        }
        return (high, low);
    }
    public static BigRational Sqrt(BigRational value)
    {
        if (value == 0) return 0;
        var         hl = SqrtLimits(value.GetWholePart);
        BigRational n  = 0, p = 0;
        if (hl.High == 0 && hl.Low == 0)
            return 0;
        var high = hl.High;
        var low  = hl.Low;
        var d    = DecimalPrecision;
        var pp   = 1 / (BigRational) d;
        while (high > low + pp)
        {
            n = (high + low) / 2;
            p = n            * n;
            if (value < p)
                high = n;
            else if (value > p)
                low = n;
            else
                break;
        }
        var r = value == p ? n : low;
        return r;
    }
    public BigRational Sqrt()
    {
        return Sqrt(this);
    }
    public static BigRational ArcTangent(BigRational v, int n)
    {
        var retVal = v;
        for (var i = 1; i < n; i++)
        {
            var powRat = Pow(v, 2 * i + 1);
            retVal += new BigRational(powRat.Numerator * (BigInteger) Math.Pow(-1d, i),
                (2 * i + 1)                            * powRat.Denominator);
        }
        return retVal;
    }
    public static BigRational Reciprocal(BigRational v)
    {
        return new BigRational(v.Denominator, v.Numerator);
    }
    public static BigRational Round(BigRational number, int decimalPlaces)
    {
        BigRational power = BigInteger.Pow(10, decimalPlaces);
        number *= power;
        return number >= 0 ? (BigInteger) (number + 0.5) / power : (BigInteger) (number - 0.5) / power;
    }
    public void Round(int decimalPlaces)
    {
        var         number = this;
        BigRational power  = BigInteger.Pow(10, decimalPlaces);
        number *= power;
        var n = number >= 0 ? (BigInteger) (number + 0.5) / power : (BigInteger) (number - 0.5) / power;
        Numerator   = n.Numerator;
        Denominator = n.Denominator;
    }
    public static BigRational Pow(BigRational v, int e)
    {
        if (e < 1) throw new ArgumentException("Powers must be greater than or equal to one.");
        var retVal = new BigRational(v.Numerator, v.Denominator);
        for (var i = 1; i < e; i++)
        {
            retVal.Numerator   *= v.Numerator;
            retVal.Denominator *= v.Denominator;
        }
        return retVal;
    }
    public static BigRational Min(BigRational r, BigRational l)
    {
        return l < r ? l : r;
    }
    public static BigRational Max(BigRational r, BigRational l)
    {
        return l > r ? l : r;
    }
    /// <summary>
    ///     Set Pi before call
    /// </summary>
    public static BigRational ToRadians(BigRational degrees)
    {
        return degrees * Pi / 180;
    }
    /// <summary>
    ///     Set Pi before call
    /// </summary>
    public static BigRational ToDegrees(BigRational rads)
    {
        return rads * 180 / Pi;
    }
    private static BigRational Factorial(BigRational x)
    {
        BigRational r = 1;
        BigRational c = 1;
        while (c <= x)
        {
            r *= c;
            c++;
        }
        return r;
    }
    public static BigRational Exp(BigRational x)
    {
        BigRational r  = 0;
        BigRational r1 = 0;
        var         k  = 0;
        while (true)
        {
            r += Pow(x, k) / Factorial(k);
            if (r == r1)
                break;
            r1 = r;
            k++;
        }
        return r;
    }
    public static BigRational Sine(BigRational ar, int n)
    {
        if (Factorials == null)
        {
            Factorials = new BigRational[MaxFactorials];
            for (var i = 0; i < MaxFactorials; i++)
                Factorials[i] = new BigRational();
            for (var i = 1; i < MaxFactorials + 1; i++)
                Factorials[i - 1] = Factorial(i);
        }
        var sin = ar;
        for (var i = 1; i <= n; i++)
        {
            var trm = Pow(ar, i * 2 + 1);
            trm /= Factorials[i * 2];
            if ((i & 1) == 1)
                sin -= trm;
            else
                sin += trm;
        }
        return sin;
    }
    public static BigRational Atan(BigRational ar, int n)
    {
        var atan = ar;
        for (var i = 1; i <= n; i++)
        {
            var trm = Pow(ar, i * 2 + 1);
            trm /= i * 2;
            if ((i & 1) == 1)
                atan -= trm;
            else
                atan += trm;
        }
        return atan;
    }
    public static BigRational Cosine(BigRational ar, int n)
    {
        if (Factorials == null)
        {
            Factorials = new BigRational[MaxFactorials];
            for (var i = 0; i < MaxFactorials; i++)
                Factorials[i] = new BigRational();
            for (var i = 1; i < MaxFactorials + 1; i++)
                Factorials[i - 1] = Factorial(i);
        }
        BigRational cos = 1.0;
        for (var i = 1; i <= n; i++)
        {
            var trm = Pow(ar, i * 2);
            trm /= Factorials[i * 2 - 1];
            if ((i & 1) == 1)
                cos -= trm;
            else
                cos += trm;
        }
        return cos;
    }
    public static BigRational Tangent(BigRational ar, int n)
    {
        return Sine(ar, n) / Cosine(ar, n);
    }
    public static BigRational CoTangent(BigRational ar, int n)
    {
        return Cosine(ar, n) / Sine(ar, n);
    }
    public static BigRational Secant(BigRational ar, int n)
    {
        return 1.0 / Cosine(ar, n);
    }
    public static BigRational CoSecant(BigRational ar, int n)
    {
        return 1.0 / Sine(ar, n);
    }
    private static BigRational GetE(int n)
    {
        BigRational e = 1.0;
        var         c = n;
        while (c > 0)
        {
            BigRational f = 0;
            if (c == 1)
            {
                f = 1;
            }
            else
            {
                var i = c - 1;
                f = c;
                while (i > 0)
                {
                    f *= i;
                    i--;
                }
            }
            c--;
            e += 1.0 / f;
        }
        return e;
    }
    public static BigRational NthRoot(BigRational value, int nth)
    {
        BigRational lx;
        var         a = value;
        var         n = nth;
        BigRational s = 1.0;
        do
        {
            var t = s;
            lx = a / Pow(s, n - 1);
            var r = (n        - 1) * s;
            s = (lx + r) / n;
        } while (lx != s);
        return s;
    }
    public static BigRational LogN(BigRational value)
    {
        BigRational a;
        var         p = value;
        BigRational n = 0.0;
        while (p >= E)
        {
            p /= E;
            n++;
        }
        n += p / E;
        p =  value;
        do
        {
            a = n;
            var lx = p         / Exp(n - 1.0);
            var r  = (n - 1.0) * E;
            n = (lx + r) / E;
        } while (n != a);
        return n;
    }
    public static BigRational Log(BigRational n, int b)
    {
        return LogN(n) / LogN(b);
    }
    private static int ConversionIterations(BigRational v)
    {
        return (int) ((DecimalMaxScale + 1) / (2 * Math.Log10((double) Reciprocal(v))));
    }
    public static BigRational GetPI()
    {
        var oneFifth         = new BigRational(1, 5);
        var oneTwoThirtyNine = new BigRational(1, 239);
        var arcTanOneFifth   = ArcTangent(oneFifth, ConversionIterations(oneFifth));
        var arcTanOneTwoThirtyNine =
            ArcTangent(oneTwoThirtyNine, ConversionIterations(oneTwoThirtyNine));
        return arcTanOneFifth * 16 - arcTanOneTwoThirtyNine * 4;
    }
    public override string ToString()
    {
        var ret = new StringBuilder();
        ret.Append(Numerator.ToString("R", CultureInfo.InvariantCulture));
        ret.Append(Solidus);
        ret.Append(Denominator.ToString("R", CultureInfo.InvariantCulture));
        return ret.ToString();
    }
}

Random64csp.cs

64 Bit Random Number Generator Using RNGCryptoServiceProvider

using System;
using System.Security.Cryptography;
[Serializable]
public class Random64csp : RandomNumberGenerator
{
    private byte[] _buffer;
    private RNGCryptoServiceProvider _crng;
    private double _dBi;
    private ulong UpperLimit = ulong.MaxValue;
    public Random64csp()
    {
        SetDataUse = 8;
        _crng = new RNGCryptoServiceProvider();
    }
    
    public Random64csp(int dataSize)
    {
        SetDataUse = dataSize;
        _crng      = new RNGCryptoServiceProvider();
    }
    
    public int SetDataUse
    {
        get => _buffer.Length;
        set
        {
            var v = value;
            if (v < 1 || v > 8 || v == 3 || v == 5 || v == 6 || v == 7)
                throw new ArgumentException($"Value {v} must be either 1 or 2 or 4 or 8");
            switch (v)
            {
                case 1:
                    _dBi = 1.0D / byte.MaxValue;
                    UpperLimit = byte.MaxValue;
                    _buffer = new byte[1];
                    break;
                case 2:
                    _dBi = 1.0D / ushort.MaxValue;
                    UpperLimit = ushort.MaxValue;
                    _buffer = new byte[2];
                    break;
                case 4:
                    _dBi = 1.0D / uint.MaxValue;
                    UpperLimit = uint.MaxValue;
                    _buffer = new byte[4];
                    break;
                case 8:
                    _dBi = 1.0D / ulong.MaxValue;
                    UpperLimit = ulong.MaxValue;
                    _buffer = new byte[8];
                    break;
                default:
                    _dBi = 1.0D / ulong.MaxValue;
                    UpperLimit = ulong.MaxValue;
                    _buffer = new byte[8];
                    break;
            }
        }
    }
    /// <summary>
    ///     Will generate only odd values
    /// </summary>
    public bool OddsOnly
    {
        get;
        set;
    }
    private double Sample()
    {
        ulong Internal()
        {
            _crng.GetBytes(_buffer);
            return BufferToLong(_buffer);
        }
        return Internal() * _dBi;
    }
    public ulong Next(ulong minValue, ulong maxValue)
    {
        var sa = Sample();
        var fi = (double)(maxValue - minValue + minValue);
        var n = (ulong)(sa * fi);
        n = !OddsOnly ? n : n | 1;
        if (n < minValue)
            return Next(minValue, maxValue);
        return n;
    }
    public ulong Next(ulong maxValue)
    {
        return Next(0, maxValue);
    }
    public ulong Next()
    {
        return Next(0, UpperLimit);
    }
    public unsafe double NextDouble()
    {
        var buf = new byte[8];
        GetBytes(buf);
        fixed (byte* ptr = buf)
        {
            return *(ulong*)ptr * _dBi * ulong.MaxValue;
        }
    }
    public char[] GetNextCharArray(int size)
    {
        var xbc = new byte[1];
        var ca = new char[size];
        var ptr = 0;
        do
        {
            _crng.GetBytes(xbc);
            var c = xbc[0];
            if (c >= 32 && c <= 127)
                ca[ptr++] = (char)c;
        } while (ptr < size);
        return ca;
    }
    public char[] GetNextCharArrayX(int size)
    {
        var xbc = new byte[2];
        var ca = new char[size];
        var ptr = 0;
        do
        {
            _crng.GetBytes(xbc);
            ca[ptr++] = Convert.ToChar(xbc.ToShort());
        } while (ptr < size);
        return ca;
    }
    public byte[] GetNextByteArray(int size)
    {
        var ba = new byte[size];
        _crng.GetBytes(ba);
        return ba;
    }
    /// <summary>
    ///     Next(0,2)==0?false:true; has a distribution error of 1% weighted toward zero.
    ///     The distribution error here is 0.000000046% which is statistically insignificant in this context.
    /// </summary>
    /// <param name="size"></param>
    /// <returns></returns>
    public bool[] GetNextBoolArrayLimit(int size)
    {
        var ba = new bool[size];
        const uint ll = uint.MaxValue >> 1;
        for (var i = 0; i < size; ++i)
            ba[i] = Next(0, uint.MaxValue) > ll;
        return ba;
    }
    public byte[] GetNextByteArrayLimit(int size, ulong minValue, ulong maxValue)
    {
        var ba = new byte[size];
        for (var i = 0; i < size; ++i)
            ba[i] = (byte)Next(minValue, maxValue);
        return ba;
    }
    public ushort[] GetNextUShortArrayLimit(int size, ulong minValue, ulong maxValue)
    {
        var ba = new ushort[size];
        for (var i = 0; i < size; ++i)
            ba[i] = (ushort)Next(minValue, maxValue);
        return ba;
    }
    public uint[] GetNextUIntArrayLimit(int size, ulong minValue, ulong maxValue)
    {
        var ba = new uint[size];
        for (var i = 0; i < size; ++i)
            ba[i] = (uint)Next(minValue, maxValue);
        return ba;
    }
    public ulong[] GetNextULongArrayLimit(int size, ulong minValue, ulong maxValue)
    {
        var ba = new ulong[size];
        for (var i = 0; i < size; ++i)
            ba[i] = Next(minValue, maxValue);
        return ba;
    }
    public string GetRandomString(int minLen, int maxLen)
    {
        if (minLen == maxLen)
            return new string(GetNextCharArray(minLen));
        return new string(GetNextCharArray((int)Next((ulong)minLen, (ulong)maxLen)));
    }
    public override void GetBytes(byte[] data)
    {
        if (data == null)
            throw new ArgumentException("The buffer cannot be null.");
        _crng.GetBytes(data);
    }
    public void NextBytes(byte[] buffer)
    {
        if (buffer == null)
            throw new ArgumentNullException("The buffer cannot be null.");
        for (var index = 0; index < buffer.Length; ++index)
            buffer[index] = (byte)(Sample() * byte.MaxValue + 1);
    }
    public override void GetNonZeroBytes(byte[] buffer)
    {
        if (buffer == null)
            throw new ArgumentNullException("The buffer cannot be null.");
        var index = 0;
        do
        {
            var v = (byte)(Sample() * byte.MaxValue + 1);
            if (v > 0)
            {
                buffer[index] = v;
                index++;
            }
        } while (index < buffer.Length);
    }
    private unsafe ulong BufferToLong(byte[] buffer)
    {
        var len = buffer.Length;
        if (len < 1 || len > 8)
            throw new ArgumentException($"The array length {len} must be between 1 and 8");
        fixed (byte* Ptr = &buffer[0])
        {
            switch (len)
            {
                case 1:
                    return *Ptr;
                case 2:
                    return (uint)(*Ptr | (Ptr[1] << 8));
                case 3:
                    return (uint)(*Ptr | (Ptr[1] << 8) | (Ptr[2] << 16));
                case 4:
                    return (uint)(*Ptr | (Ptr[1] << 8) | (Ptr[2] << 16) | (Ptr[3] << 24));
                case 5:
                    return (uint)(*Ptr | (Ptr[1] << 8) | (Ptr[2] << 16) | (Ptr[3] << 24)) | ((ulong)Ptr[4] << 32);
                case 6:
                    return (uint)(*Ptr | (Ptr[1] << 8) | (Ptr[2] << 16) | (Ptr[3] << 24)) | ((ulong)(Ptr[4] | (Ptr[5] << 8)) << 32);
                case 7:
                    return (uint)(*Ptr | (Ptr[1] << 8) | (Ptr[2] << 16) | (Ptr[3] << 24)) | ((ulong)(Ptr[4] | (Ptr[5] << 8) | (Ptr[6] << 16)) << 32);
                case 8:
                    return (uint)(*Ptr | (Ptr[1] << 8) | (Ptr[2] << 16) | (Ptr[3] << 24)) | ((ulong)(Ptr[4] | (Ptr[5] << 8) | (Ptr[6] << 16) | (Ptr[7] << 24)) << 32);
                default:
                    return 0;
            }
        }
    }
}

Random64d.cs

64 Bit Version of Random

using System;
using System.Runtime.InteropServices;
using System.Security.Cryptography;
[Serializable]
public class Random64d : RandomNumberGenerator
{
    private const                    double                   DBi  = 1.0 / 9223372036854775807;
    private const                    double                   DBui = 1.0 / ulong.MaxValue;
    private readonly                 Random                   _r1;
    private readonly                 Random                   _r2;
    [NonSerialized] private readonly RNGCryptoServiceProvider _rng = new RNGCryptoServiceProvider();
    private                          RNumber                  _rNumber;
    public Random64d() : this(((long) (uint) Environment.TickCount << 32) | (uint) Environment.TickCount)
    {
    }
    public Random64d(long Seed)
    {
        var Seed0 = (int) (((Seed ^ 0x86B28CAB) << 16) & 0x7fffffff);
        var Seed1 = (int) (Seed                        & 0x7fffffff);
        _r1 = new Random(Seed0);
        _r2 = new Random(Seed1);
    }
    public long InternalSample()
    {
        _rNumber.n1 = _r1.Next();
        _rNumber.n2 = _r2.Next();
        return _rNumber.total;
    }
    private double GetSampleForLargeRange()
    {
        var value = InternalSample();
        if (InternalSample() % 2 == 0)
            value = -value;
        return (value + 4611686018427387903.0) / 9223372036854775807.0;
    }
    private double Sample()
    {
        return InternalSample() * DBi;
    }
    public long Next()
    {
        return InternalSample();
    }
    public long Next(long minValue, long maxValue)
    {
        if (minValue > maxValue)
            throw new ArgumentException("maxValue must be greater than or equal to minValue");
        var num = (ulong) maxValue - (ulong) minValue;
        if (num <= long.MaxValue)
            return (long) (Sample() * num) + minValue;
        return (long) ((ulong) (GetSampleForLargeRange() * num) + (ulong) minValue);
    }
    public ulong Next(ulong minValue, ulong maxValue)
    {
        if (minValue > maxValue)
            throw new ArgumentException("maxValue must be greater than or equal to minValue");
        var num = maxValue - minValue;
        return (ulong) (InternalSample() * DBui * num) + minValue;
    }
    public long Next(long maxValue)
    {
        if (maxValue < 0)
            throw new ArgumentException("maxValue must be greater than zero.");
        return (long) (Sample() * maxValue);
    }
    public ulong Next(ulong maxValue)
    {
        return (ulong) (InternalSample() * DBui * maxValue);
    }
    public double NextDouble()
    {
        return Sample();
    }
    public byte[] GetNextByteArray(int size)
    {
        var ba = new byte[size];
        _rng.GetBytes(ba);
        return ba;
    }
    public byte[] GetNextByteArrayFn(int size)
    {
        var ba = new byte[size];
        for (var i = 0; i < size; ++i)
            ba[i] = (byte) Next(0, 256);
        return ba;
    }
    public char[] GetNextCharArray(int size)
    {
        var ca = new char[size];
        for (var i = 0; i < size; ++i)
            ca[i] = (char)Next(32, 128);
        return ca;
    }
    public string GetRandomString(int minLen, int maxLen)
    {
        if (minLen == maxLen)
            return new string(GetNextCharArray(minLen));
        return new string(GetNextCharArray((int) Next(minLen, maxLen)));
    }
    public override void GetBytes(byte[] data)
    {
        if (data == null)
            throw new ArgumentException("The buffer cannot be null.");
        _rng.GetBytes(data);
    }
    [StructLayout(LayoutKind.Explicit, Size = 8, Pack = 1)]
    [Serializable]
    public struct RNumber
    {
        [FieldOffset(0)] public long n1;
        [FieldOffset(4)] public long n2;
        [FieldOffset(0)] public long total;
    }
}