From 72418d802c006bc2fbe8123642eef1dddf0e4239 Mon Sep 17 00:00:00 2001 From: Milan Hanajik Date: Mon, 12 Dec 2016 12:51:00 +0100 Subject: [PATCH] UInt128 type added - Dirichlet.Numerics project/namespace --- Dirichlet.Numerics/Dirichlet.Numerics.csproj | 51 + Dirichlet.Numerics/Int128.cs | 1753 +++++++++++ Dirichlet.Numerics/Properties/AssemblyInfo.cs | 36 + Dirichlet.Numerics/UInt128.cs | 2638 +++++++++++++++++ Dirichlet.Numerics/license.txt | 13 + TestBenchFramework.sln | 12 + 6 files changed, 4503 insertions(+) create mode 100644 Dirichlet.Numerics/Dirichlet.Numerics.csproj create mode 100644 Dirichlet.Numerics/Int128.cs create mode 100644 Dirichlet.Numerics/Properties/AssemblyInfo.cs create mode 100644 Dirichlet.Numerics/UInt128.cs create mode 100644 Dirichlet.Numerics/license.txt diff --git a/Dirichlet.Numerics/Dirichlet.Numerics.csproj b/Dirichlet.Numerics/Dirichlet.Numerics.csproj new file mode 100644 index 000000000..6e1ff9f27 --- /dev/null +++ b/Dirichlet.Numerics/Dirichlet.Numerics.csproj @@ -0,0 +1,51 @@ + + + + Debug + AnyCPU + 8.0.30703 + 2.0 + {439D0878-C76E-452B-B17D-209A89E91D36} + Library + Properties + Dirichlet.Numerics + Dirichlet.Numerics + v4.0 + 512 + Client + + + true + full + false + bin\Debug\ + DEBUG;TRACE + prompt + 4 + + + pdbonly + true + bin\Release\ + TRACE + prompt + 4 + + + + + + + + + + + + + diff --git a/Dirichlet.Numerics/Int128.cs b/Dirichlet.Numerics/Int128.cs new file mode 100644 index 000000000..8473c9cc0 --- /dev/null +++ b/Dirichlet.Numerics/Int128.cs @@ -0,0 +1,1753 @@ +using System; +using System.Diagnostics; +using System.Globalization; +using System.Linq; +using System.Numerics; + +namespace Dirichlet.Numerics +{ + public struct Int128 : IFormattable, IComparable, IComparable, IEquatable + { + private UInt128 v; + + private static readonly Int128 minValue = (Int128)((UInt128)1 << 127); + private static readonly Int128 maxValue = (Int128)(((UInt128)1 << 127) - 1); + private static readonly Int128 zero = (Int128)0; + private static readonly Int128 one = (Int128)1; + private static readonly Int128 minusOne = (Int128)(-1); + + public static Int128 MinValue { get { return minValue; } } + public static Int128 MaxValue { get { return maxValue; } } + public static Int128 Zero { get { return zero; } } + public static Int128 One { get { return one; } } + public static Int128 MinusOne { get { return minusOne; } } + + public static Int128 Parse(string value) + { + Int128 c; + if (!TryParse(value, out c)) + throw new FormatException(); + return c; + } + + public static bool TryParse(string value, out Int128 result) + { + return TryParse(value, NumberStyles.Integer, NumberFormatInfo.CurrentInfo, out result); + } + + public static bool TryParse(string value, NumberStyles style, IFormatProvider format, out Int128 result) + { + BigInteger a; + if (!BigInteger.TryParse(value, style, format, out a)) + { + result = Int128.Zero; + return false; + } + UInt128.Create(out result.v, a); + return true; + } + + public Int128(long value) + { + UInt128.Create(out v, value); + } + + public Int128(ulong value) + { + UInt128.Create(out v, value); + } + + public Int128(double value) + { + UInt128.Create(out v, value); + } + + public Int128(decimal value) + { + UInt128.Create(out v, value); + } + + public Int128(BigInteger value) + { + UInt128.Create(out v, value); + } + + public ulong S0 { get { return v.S0; } } + public ulong S1 { get { return v.S1; } } + + public bool IsZero { get { return v.IsZero; } } + public bool IsOne { get { return v.IsOne; } } + public bool IsPowerOfTwo { get { return v.IsPowerOfTwo; } } + public bool IsEven { get { return v.IsEven; } } + public bool IsNegative { get { return v.S1 > long.MaxValue; } } + public int Sign { get { return IsNegative ? -1 : v.Sign; } } + + public override string ToString() + { + return ((BigInteger)this).ToString(); + } + + public string ToString(string format) + { + return ((BigInteger)this).ToString(format); + } + + public string ToString(IFormatProvider provider) + { + return ToString(null, provider); + } + + public string ToString(string format, IFormatProvider provider) + { + return ((BigInteger)this).ToString(format, provider); + } + + public static explicit operator Int128(double a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static implicit operator Int128(sbyte a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static implicit operator Int128(byte a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static implicit operator Int128(short a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static implicit operator Int128(ushort a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static implicit operator Int128(int a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static implicit operator Int128(uint a) + { + Int128 c; + UInt128.Create(out c.v, (ulong)a); + return c; + } + + public static implicit operator Int128(long a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static implicit operator Int128(ulong a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static explicit operator Int128(decimal a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static explicit operator Int128(UInt128 a) + { + Int128 c; + c.v = a; + return c; + } + + public static explicit operator UInt128(Int128 a) + { + return a.v; + } + + public static explicit operator Int128(BigInteger a) + { + Int128 c; + UInt128.Create(out c.v, a); + return c; + } + + public static explicit operator sbyte(Int128 a) + { + return (sbyte)a.v.S0; + } + + public static explicit operator byte(Int128 a) + { + return (byte)a.v.S0; + } + + public static explicit operator short(Int128 a) + { + return (short)a.v.S0; + } + + public static explicit operator ushort(Int128 a) + { + return (ushort)a.v.S0; + } + + public static explicit operator int(Int128 a) + { + return (int)a.v.S0; + } + + public static explicit operator uint(Int128 a) + { + return (uint)a.v.S0; + } + + public static explicit operator long(Int128 a) + { + return (long)a.v.S0; + } + + public static explicit operator ulong(Int128 a) + { + return a.v.S0; + } + + public static explicit operator decimal(Int128 a) + { + if (a.IsNegative) + { + UInt128 c; + UInt128.Negate(out c, ref a.v); + return -(decimal)c; + } + return (decimal)a.v; + } + + public static implicit operator BigInteger(Int128 a) + { + if (a.IsNegative) + { + UInt128 c; + UInt128.Negate(out c, ref a.v); + return -(BigInteger)c; + } + return (BigInteger)a.v; + } + + public static explicit operator float(Int128 a) + { + if (a.IsNegative) + { + UInt128 c; + UInt128.Negate(out c, ref a.v); + return -UInt128.ConvertToFloat(ref c); + } + return UInt128.ConvertToFloat(ref a.v); + } + + public static explicit operator double(Int128 a) + { + if (a.IsNegative) + { + UInt128 c; + UInt128.Negate(out c, ref a.v); + return -UInt128.ConvertToDouble(ref c); + } + return UInt128.ConvertToDouble(ref a.v); + } + + public static Int128 operator <<(Int128 a, int b) + { + Int128 c; + UInt128.LeftShift(out c.v, ref a.v, b); + return c; + } + + public static Int128 operator >>(Int128 a, int b) + { + Int128 c; + UInt128.ArithmeticRightShift(out c.v, ref a.v, b); + return c; + } + + public static Int128 operator &(Int128 a, Int128 b) + { + Int128 c; + UInt128.And(out c.v, ref a.v, ref b.v); + return c; + } + + public static int operator &(Int128 a, int b) + { + return (int)(a.v & (uint)b); + } + + public static int operator &(int a, Int128 b) + { + return (int)(b.v & (uint)a); + } + + public static long operator &(Int128 a, long b) + { + return (long)(a.v & (ulong)b); + } + + public static long operator &(long a, Int128 b) + { + return (long)(b.v & (ulong)a); + } + + public static Int128 operator |(Int128 a, Int128 b) + { + Int128 c; + UInt128.Or(out c.v, ref a.v, ref b.v); + return c; + } + + public static Int128 operator ^(Int128 a, Int128 b) + { + Int128 c; + UInt128.ExclusiveOr(out c.v, ref a.v, ref b.v); + return c; + } + + public static Int128 operator ~(Int128 a) + { + Int128 c; + UInt128.Not(out c.v, ref a.v); + return c; + } + + public static Int128 operator +(Int128 a, long b) + { + Int128 c; + if (b < 0) + UInt128.Subtract(out c.v, ref a.v, (ulong)(-b)); + else + UInt128.Add(out c.v, ref a.v, (ulong)b); + return c; + } + + public static Int128 operator +(long a, Int128 b) + { + Int128 c; + if (a < 0) + UInt128.Subtract(out c.v, ref b.v, (ulong)(-a)); + else + UInt128.Add(out c.v, ref b.v, (ulong)a); + return c; + } + + public static Int128 operator +(Int128 a, Int128 b) + { + Int128 c; + UInt128.Add(out c.v, ref a.v, ref b.v); + return c; + } + + public static Int128 operator ++(Int128 a) + { + Int128 c; + UInt128.Add(out c.v, ref a.v, 1); + return c; + } + + public static Int128 operator -(Int128 a, long b) + { + Int128 c; + if (b < 0) + UInt128.Add(out c.v, ref a.v, (ulong)(-b)); + else + UInt128.Subtract(out c.v, ref a.v, (ulong)b); + return c; + } + + public static Int128 operator -(Int128 a, Int128 b) + { + Int128 c; + UInt128.Subtract(out c.v, ref a.v, ref b.v); + return c; + } + + public static Int128 operator +(Int128 a) + { + return a; + } + + public static Int128 operator -(Int128 a) + { + Int128 c; + UInt128.Negate(out c.v, ref a.v); + return c; + } + + public static Int128 operator --(Int128 a) + { + Int128 c; + UInt128.Subtract(out c.v, ref a.v, 1); + return c; + } + + public static Int128 operator *(Int128 a, int b) + { + Int128 c; + Multiply(out c, ref a, b); + return c; + } + + public static Int128 operator *(int a, Int128 b) + { + Int128 c; + Multiply(out c, ref b, a); + return c; + } + + public static Int128 operator *(Int128 a, uint b) + { + Int128 c; + Multiply(out c, ref a, b); + return c; + } + + public static Int128 operator *(uint a, Int128 b) + { + Int128 c; + Multiply(out c, ref b, a); + return c; + } + + public static Int128 operator *(Int128 a, long b) + { + Int128 c; + Multiply(out c, ref a, b); + return c; + } + + public static Int128 operator *(long a, Int128 b) + { + Int128 c; + Multiply(out c, ref b, a); + return c; + } + + public static Int128 operator *(Int128 a, ulong b) + { + Int128 c; + Multiply(out c, ref a, b); + return c; + } + + public static Int128 operator *(ulong a, Int128 b) + { + Int128 c; + Multiply(out c, ref b, a); + return c; + } + + public static Int128 operator *(Int128 a, Int128 b) + { + Int128 c; + Multiply(out c, ref a, ref b); + return c; + } + + public static Int128 operator /(Int128 a, int b) + { + Int128 c; + Divide(out c, ref a, b); + return c; + } + + public static Int128 operator /(Int128 a, uint b) + { + Int128 c; + Divide(out c, ref a, b); + return c; + } + + public static Int128 operator /(Int128 a, long b) + { + Int128 c; + Divide(out c, ref a, b); + return c; + } + + public static Int128 operator /(Int128 a, ulong b) + { + Int128 c; + Divide(out c, ref a, b); + return c; + } + + public static Int128 operator /(Int128 a, Int128 b) + { + Int128 c; + Divide(out c, ref a, ref b); + return c; + } + + public static int operator %(Int128 a, int b) + { + return Remainder(ref a, b); + } + + public static int operator %(Int128 a, uint b) + { + return Remainder(ref a, b); + } + + public static long operator %(Int128 a, long b) + { + return Remainder(ref a, b); + } + + public static long operator %(Int128 a, ulong b) + { + return Remainder(ref a, b); + } + + public static Int128 operator %(Int128 a, Int128 b) + { + Int128 c; + Remainder(out c, ref a, ref b); + return c; + } + + public static bool operator <(Int128 a, UInt128 b) + { + return a.CompareTo(b) < 0; + } + + public static bool operator <(UInt128 a, Int128 b) + { + return b.CompareTo(a) > 0; + } + + public static bool operator <(Int128 a, Int128 b) + { + return LessThan(ref a.v, ref b.v); + } + + public static bool operator <(Int128 a, int b) + { + return LessThan(ref a.v, b); + } + + public static bool operator <(int a, Int128 b) + { + return LessThan(a, ref b.v); + } + + public static bool operator <(Int128 a, uint b) + { + return LessThan(ref a.v, b); + } + + public static bool operator <(uint a, Int128 b) + { + return LessThan(a, ref b.v); + } + + public static bool operator <(Int128 a, long b) + { + return LessThan(ref a.v, b); + } + + public static bool operator <(long a, Int128 b) + { + return LessThan(a, ref b.v); + } + + public static bool operator <(Int128 a, ulong b) + { + return LessThan(ref a.v, b); + } + + public static bool operator <(ulong a, Int128 b) + { + return LessThan(a, ref b.v); + } + + public static bool operator <=(Int128 a, UInt128 b) + { + return a.CompareTo(b) <= 0; + } + + public static bool operator <=(UInt128 a, Int128 b) + { + return b.CompareTo(a) >= 0; + } + + public static bool operator <=(Int128 a, Int128 b) + { + return !LessThan(ref b.v, ref a.v); + } + + public static bool operator <=(Int128 a, int b) + { + return !LessThan(b, ref a.v); + } + + public static bool operator <=(int a, Int128 b) + { + return !LessThan(ref b.v, a); + } + + public static bool operator <=(Int128 a, uint b) + { + return !LessThan(b, ref a.v); + } + + public static bool operator <=(uint a, Int128 b) + { + return !LessThan(ref b.v, a); + } + + public static bool operator <=(Int128 a, long b) + { + return !LessThan(b, ref a.v); + } + + public static bool operator <=(long a, Int128 b) + { + return !LessThan(ref b.v, a); + } + + public static bool operator <=(Int128 a, ulong b) + { + return !LessThan(b, ref a.v); + } + + public static bool operator <=(ulong a, Int128 b) + { + return !LessThan(ref b.v, a); + } + + public static bool operator >(Int128 a, UInt128 b) + { + return a.CompareTo(b) > 0; + } + + public static bool operator >(UInt128 a, Int128 b) + { + return b.CompareTo(a) < 0; + } + + public static bool operator >(Int128 a, Int128 b) + { + return LessThan(ref b.v, ref a.v); + } + + public static bool operator >(Int128 a, int b) + { + return LessThan(b, ref a.v); + } + + public static bool operator >(int a, Int128 b) + { + return LessThan(ref b.v, a); + } + + public static bool operator >(Int128 a, uint b) + { + return LessThan(b, ref a.v); + } + + public static bool operator >(uint a, Int128 b) + { + return LessThan(ref b.v, a); + } + + public static bool operator >(Int128 a, long b) + { + return LessThan(b, ref a.v); + } + + public static bool operator >(long a, Int128 b) + { + return LessThan(ref b.v, a); + } + + public static bool operator >(Int128 a, ulong b) + { + return LessThan(b, ref a.v); + } + + public static bool operator >(ulong a, Int128 b) + { + return LessThan(ref b.v, a); + } + + public static bool operator >=(Int128 a, UInt128 b) + { + return a.CompareTo(b) >= 0; + } + + public static bool operator >=(UInt128 a, Int128 b) + { + return b.CompareTo(a) <= 0; + } + + public static bool operator >=(Int128 a, Int128 b) + { + return !LessThan(ref a.v, ref b.v); + } + + public static bool operator >=(Int128 a, int b) + { + return !LessThan(ref a.v, b); + } + + public static bool operator >=(int a, Int128 b) + { + return !LessThan(a, ref b.v); + } + + public static bool operator >=(Int128 a, uint b) + { + return !LessThan(ref a.v, b); + } + + public static bool operator >=(uint a, Int128 b) + { + return !LessThan(a, ref b.v); + } + + public static bool operator >=(Int128 a, long b) + { + return !LessThan(ref a.v, b); + } + + public static bool operator >=(long a, Int128 b) + { + return !LessThan(a, ref b.v); + } + + public static bool operator >=(Int128 a, ulong b) + { + return !LessThan(ref a.v, b); + } + + public static bool operator >=(ulong a, Int128 b) + { + return !LessThan(a, ref b.v); + } + + public static bool operator ==(UInt128 a, Int128 b) + { + return b.Equals(a); + } + + public static bool operator ==(Int128 a, UInt128 b) + { + return a.Equals(b); + } + + public static bool operator ==(Int128 a, Int128 b) + { + return a.v.Equals(b.v); + } + + public static bool operator ==(Int128 a, int b) + { + return a.Equals(b); + } + + public static bool operator ==(int a, Int128 b) + { + return b.Equals(a); + } + + public static bool operator ==(Int128 a, uint b) + { + return a.Equals(b); + } + + public static bool operator ==(uint a, Int128 b) + { + return b.Equals(a); + } + + public static bool operator ==(Int128 a, long b) + { + return a.Equals(b); + } + + public static bool operator ==(long a, Int128 b) + { + return b.Equals(a); + } + + public static bool operator ==(Int128 a, ulong b) + { + return a.Equals(b); + } + + public static bool operator ==(ulong a, Int128 b) + { + return b.Equals(a); + } + + public static bool operator !=(UInt128 a, Int128 b) + { + return !b.Equals(a); + } + + public static bool operator !=(Int128 a, UInt128 b) + { + return !a.Equals(b); + } + + public static bool operator !=(Int128 a, Int128 b) + { + return !a.v.Equals(b.v); + } + + public static bool operator !=(Int128 a, int b) + { + return !a.Equals(b); + } + + public static bool operator !=(int a, Int128 b) + { + return !b.Equals(a); + } + + public static bool operator !=(Int128 a, uint b) + { + return !a.Equals(b); + } + + public static bool operator !=(uint a, Int128 b) + { + return !b.Equals(a); + } + + public static bool operator !=(Int128 a, long b) + { + return !a.Equals(b); + } + + public static bool operator !=(long a, Int128 b) + { + return !b.Equals(a); + } + + public static bool operator !=(Int128 a, ulong b) + { + return !a.Equals(b); + } + + public static bool operator !=(ulong a, Int128 b) + { + return !b.Equals(a); + } + + public int CompareTo(UInt128 other) + { + if (IsNegative) + return -1; + return v.CompareTo(other); + } + + public int CompareTo(Int128 other) + { + return SignedCompare(ref v, other.S0, other.S1); + } + + public int CompareTo(int other) + { + return SignedCompare(ref v, (ulong)other, (ulong)(other >> 31)); + } + + public int CompareTo(uint other) + { + return SignedCompare(ref v, (ulong)other, 0); + } + + public int CompareTo(long other) + { + return SignedCompare(ref v, (ulong)other, (ulong)(other >> 63)); + } + + public int CompareTo(ulong other) + { + return SignedCompare(ref v, other, 0); + } + + public int CompareTo(object obj) + { + if (obj == null) + return 1; + if (!(obj is Int128)) + throw new ArgumentException(); + return CompareTo((Int128)obj); + } + + private static bool LessThan(ref UInt128 a, ref UInt128 b) + { + var as1 = (long)a.S1; + var bs1 = (long)b.S1; + if (as1 != bs1) + return as1 < bs1; + return a.S0 < b.S0; + } + + private static bool LessThan(ref UInt128 a, long b) + { + var as1 = (long)a.S1; + var bs1 = b >> 63; + if (as1 != bs1) + return as1 < bs1; + return a.S0 < (ulong)b; + } + + private static bool LessThan(long a, ref UInt128 b) + { + var as1 = a >> 63; + var bs1 = (long)b.S1; + if (as1 != bs1) + return as1 < bs1; + return (ulong)a < b.S0; + } + + private static bool LessThan(ref UInt128 a, ulong b) + { + var as1 = (long)a.S1; + if (as1 != 0) + return as1 < 0; + return a.S0 < b; + } + + private static bool LessThan(ulong a, ref UInt128 b) + { + var bs1 = (long)b.S1; + if (0 != bs1) + return 0 < bs1; + return a < b.S0; + } + + private static int SignedCompare(ref UInt128 a, ulong bs0, ulong bs1) + { + var as1 = a.S1; + if (as1 != bs1) + return ((long)as1).CompareTo((long)bs1); + return a.S0.CompareTo(bs0); + } + + public bool Equals(UInt128 other) + { + return !IsNegative && v.Equals(other); + } + + public bool Equals(Int128 other) + { + return v.Equals(other.v); + } + + public bool Equals(int other) + { + if (other < 0) + return v.S1 == ulong.MaxValue && v.S0 == (uint)other; + return v.S1 == 0 && v.S0 == (uint)other; + } + + public bool Equals(uint other) + { + return v.S1 == 0 && v.S0 == other; + } + + public bool Equals(long other) + { + if (other < 0) + return v.S1 == ulong.MaxValue && v.S0 == (ulong)other; + return v.S1 == 0 && v.S0 == (ulong)other; + } + + public bool Equals(ulong other) + { + return v.S1 == 0 && v.S0 == other; + } + + public override bool Equals(object obj) + { + if (!(obj is Int128)) + return false; + return Equals((Int128)obj); + } + + public override int GetHashCode() + { + return v.GetHashCode(); + } + + public static void Multiply(out Int128 c, ref Int128 a, int b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b < 0) + UInt128.Multiply(out c.v, ref aneg, (uint)(-b)); + else + { + UInt128.Multiply(out c.v, ref aneg, (uint)b); + UInt128.Negate(ref c.v); + } + } + else + { + if (b < 0) + { + UInt128.Multiply(out c.v, ref a.v, (uint)(-b)); + UInt128.Negate(ref c.v); + } + else + UInt128.Multiply(out c.v, ref a.v, (uint)b); + } + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b); + } + + public static void Multiply(out Int128 c, ref Int128 a, uint b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + UInt128.Multiply(out c.v, ref aneg, b); + UInt128.Negate(ref c.v); + } + else + UInt128.Multiply(out c.v, ref a.v, b); + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b); + } + + public static void Multiply(out Int128 c, ref Int128 a, long b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b < 0) + UInt128.Multiply(out c.v, ref aneg, (ulong)(-b)); + else + { + UInt128.Multiply(out c.v, ref aneg, (ulong)b); + UInt128.Negate(ref c.v); + } + } + else + { + if (b < 0) + { + UInt128.Multiply(out c.v, ref a.v, (ulong)(-b)); + UInt128.Negate(ref c.v); + } + else + UInt128.Multiply(out c.v, ref a.v, (ulong)b); + } + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b); + } + + public static void Multiply(out Int128 c, ref Int128 a, ulong b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + UInt128.Multiply(out c.v, ref aneg, b); + UInt128.Negate(ref c.v); + } + else + UInt128.Multiply(out c.v, ref a.v, b); + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b); + } + + public static void Multiply(out Int128 c, ref Int128 a, ref Int128 b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b.IsNegative) + { + UInt128 bneg; + UInt128.Negate(out bneg, ref b.v); + UInt128.Multiply(out c.v, ref aneg, ref bneg); + } + else + { + UInt128.Multiply(out c.v, ref aneg, ref b.v); + UInt128.Negate(ref c.v); + } + } + else + { + if (b.IsNegative) + { + UInt128 bneg; + UInt128.Negate(out bneg, ref b.v); + UInt128.Multiply(out c.v, ref a.v, ref bneg); + UInt128.Negate(ref c.v); + } + else + UInt128.Multiply(out c.v, ref a.v, ref b.v); + } + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b); + } + + public static void Divide(out Int128 c, ref Int128 a, int b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b < 0) + UInt128.Multiply(out c.v, ref aneg, (uint)(-b)); + else + { + UInt128.Multiply(out c.v, ref aneg, (uint)b); + UInt128.Negate(ref c.v); + } + } + else + { + if (b < 0) + { + UInt128.Multiply(out c.v, ref a.v, (uint)(-b)); + UInt128.Negate(ref c.v); + } + else + UInt128.Multiply(out c.v, ref a.v, (uint)b); + } + Debug.Assert((BigInteger)c == (BigInteger)a / (BigInteger)b); + } + + public static void Divide(out Int128 c, ref Int128 a, uint b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + UInt128.Divide(out c.v, ref aneg, b); + UInt128.Negate(ref c.v); + } + else + UInt128.Divide(out c.v, ref a.v, b); + Debug.Assert((BigInteger)c == (BigInteger)a / (BigInteger)b); + } + + public static void Divide(out Int128 c, ref Int128 a, long b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b < 0) + UInt128.Divide(out c.v, ref aneg, (ulong)(-b)); + else + { + UInt128.Divide(out c.v, ref aneg, (ulong)b); + UInt128.Negate(ref c.v); + } + } + else + { + if (b < 0) + { + UInt128.Divide(out c.v, ref a.v, (ulong)(-b)); + UInt128.Negate(ref c.v); + } + else + UInt128.Divide(out c.v, ref a.v, (ulong)b); + } + Debug.Assert((BigInteger)c == (BigInteger)a / (BigInteger)b); + } + + public static void Divide(out Int128 c, ref Int128 a, ulong b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + UInt128.Divide(out c.v, ref aneg, b); + UInt128.Negate(ref c.v); + } + else + UInt128.Divide(out c.v, ref a.v, b); + Debug.Assert((BigInteger)c == (BigInteger)a / (BigInteger)b); + } + + public static void Divide(out Int128 c, ref Int128 a, ref Int128 b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b.IsNegative) + { + UInt128 bneg; + UInt128.Negate(out bneg, ref b.v); + UInt128.Divide(out c.v, ref aneg, ref bneg); + } + else + { + UInt128.Divide(out c.v, ref aneg, ref b.v); + UInt128.Negate(ref c.v); + } + } + else + { + if (b.IsNegative) + { + UInt128 bneg; + UInt128.Negate(out bneg, ref b.v); + UInt128.Divide(out c.v, ref a.v, ref bneg); + UInt128.Negate(ref c.v); + } + else + UInt128.Divide(out c.v, ref a.v, ref b.v); + } + Debug.Assert((BigInteger)c == (BigInteger)a / (BigInteger)b); + } + + public static int Remainder(ref Int128 a, int b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b < 0) + return (int)UInt128.Remainder(ref aneg, (uint)(-b)); + else + return -(int)UInt128.Remainder(ref aneg, (uint)b); + } + else + { + if (b < 0) + return -(int)UInt128.Remainder(ref a.v, (uint)(-b)); + else + return (int)UInt128.Remainder(ref a.v, (uint)b); + } + } + + public static int Remainder(ref Int128 a, uint b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + return -(int)UInt128.Remainder(ref aneg, b); + } + else + return (int)UInt128.Remainder(ref a.v, b); + } + + public static long Remainder(ref Int128 a, long b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b < 0) + return (long)UInt128.Remainder(ref aneg, (ulong)(-b)); + else + return -(long)UInt128.Remainder(ref aneg, (ulong)b); + } + else + { + if (b < 0) + return -(long)UInt128.Remainder(ref a.v, (ulong)(-b)); + else + return (long)UInt128.Remainder(ref a.v, (ulong)b); + } + } + + public static long Remainder(ref Int128 a, ulong b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + return -(long)UInt128.Remainder(ref aneg, b); + } + else + return (long)UInt128.Remainder(ref a.v, b); + } + + public static void Remainder(out Int128 c, ref Int128 a, ref Int128 b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b.IsNegative) + { + UInt128 bneg; + UInt128.Negate(out bneg, ref b.v); + UInt128.Remainder(out c.v, ref aneg, ref bneg); + } + else + { + UInt128.Remainder(out c.v, ref aneg, ref b.v); + UInt128.Negate(ref c.v); + } + } + else + { + if (b.IsNegative) + { + UInt128 bneg; + UInt128.Negate(out bneg, ref b.v); + UInt128.Remainder(out c.v, ref a.v, ref bneg); + UInt128.Negate(ref c.v); + } + else + UInt128.Remainder(out c.v, ref a.v, ref b.v); + } + Debug.Assert((BigInteger)c == (BigInteger)a % (BigInteger)b); + } + + public static Int128 Abs(Int128 a) + { + if (!a.IsNegative) + return a; + Int128 c; + UInt128.Negate(out c.v, ref a.v); + return c; + } + + public static Int128 Square(long a) + { + if (a < 0) + a = -a; + Int128 c; + UInt128.Square(out c.v, (ulong)a); + return c; + } + + public static Int128 Square(Int128 a) + { + Int128 c; + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + UInt128.Square(out c.v, ref aneg); + } + else + UInt128.Square(out c.v, ref a.v); + return c; + } + + public static Int128 Cube(long a) + { + Int128 c; + if (a < 0) + { + UInt128.Cube(out c.v, (ulong)(-a)); + UInt128.Negate(ref c.v); + } + else + UInt128.Cube(out c.v, (ulong)a); + return c; + } + + public static Int128 Cube(Int128 a) + { + Int128 c; + if (a < 0) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + UInt128.Cube(out c.v, ref aneg); + UInt128.Negate(ref c.v); + } + else + UInt128.Cube(out c.v, ref a.v); + return c; + } + + public static void Add(ref Int128 a, long b) + { + if (b < 0) + UInt128.Subtract(ref a.v, (ulong)(-b)); + else + UInt128.Add(ref a.v, (ulong)b); + } + + public static void Add(ref Int128 a, ref Int128 b) + { + UInt128.Add(ref a.v, ref b.v); + } + + public static void Subtract(ref Int128 a, long b) + { + if (b < 0) + UInt128.Add(ref a.v, (ulong)(-b)); + else + UInt128.Subtract(ref a.v, (ulong)b); + } + + public static void Subtract(ref Int128 a, ref Int128 b) + { + UInt128.Subtract(ref a.v, ref b.v); + } + + public static void Add(ref Int128 a, Int128 b) + { + UInt128.Add(ref a.v, ref b.v); + } + + public static void Subtract(ref Int128 a, Int128 b) + { + UInt128.Subtract(ref a.v, ref b.v); + } + + public static void AddProduct(ref Int128 a, ref UInt128 b, int c) + { + UInt128 product; + if (c < 0) + { + UInt128.Multiply(out product, ref b, (uint)(-c)); + UInt128.Subtract(ref a.v, ref product); + } + else + { + UInt128.Multiply(out product, ref b, (uint)c); + UInt128.Add(ref a.v, ref product); + } + } + + public static void AddProduct(ref Int128 a, ref UInt128 b, long c) + { + UInt128 product; + if (c < 0) + { + UInt128.Multiply(out product, ref b, (ulong)(-c)); + UInt128.Subtract(ref a.v, ref product); + } + else + { + UInt128.Multiply(out product, ref b, (ulong)c); + UInt128.Add(ref a.v, ref product); + } + } + + public static void SubtractProduct(ref Int128 a, ref UInt128 b, int c) + { + UInt128 d; + if (c < 0) + { + UInt128.Multiply(out d, ref b, (uint)(-c)); + UInt128.Add(ref a.v, ref d); + } + else + { + UInt128.Multiply(out d, ref b, (uint)c); + UInt128.Subtract(ref a.v, ref d); + } + } + + public static void SubtractProduct(ref Int128 a, ref UInt128 b, long c) + { + UInt128 d; + if (c < 0) + { + UInt128.Multiply(out d, ref b, (ulong)(-c)); + UInt128.Add(ref a.v, ref d); + } + else + { + UInt128.Multiply(out d, ref b, (ulong)c); + UInt128.Subtract(ref a.v, ref d); + } + } + + public static void AddProduct(ref Int128 a, UInt128 b, int c) + { + AddProduct(ref a, ref b, c); + } + + public static void AddProduct(ref Int128 a, UInt128 b, long c) + { + AddProduct(ref a, ref b, c); + } + + public static void SubtractProduct(ref Int128 a, UInt128 b, int c) + { + SubtractProduct(ref a, ref b, c); + } + + public static void SubtractProduct(ref Int128 a, UInt128 b, long c) + { + SubtractProduct(ref a, ref b, c); + } + + public static void Pow(out Int128 result, ref Int128 value, int exponent) + { + if (exponent < 0) + throw new ArgumentException("exponent must not be negative"); + if (value.IsNegative) + { + UInt128 valueneg; + UInt128.Negate(out valueneg, ref value.v); + if ((exponent & 1) == 0) + UInt128.Pow(out result.v, ref valueneg, (uint)exponent); + else + { + UInt128.Pow(out result.v, ref valueneg, (uint)exponent); + UInt128.Negate(ref result.v); + } + } + else + UInt128.Pow(out result.v, ref value.v, (uint)exponent); + } + + public static Int128 Pow(Int128 value, int exponent) + { + Int128 result; + Pow(out result, ref value, exponent); + return result; + } + + public static ulong FloorSqrt(Int128 a) + { + if (a.IsNegative) + throw new ArgumentException("argument must not be negative"); + return UInt128.FloorSqrt(a.v); + } + + public static ulong CeilingSqrt(Int128 a) + { + if (a.IsNegative) + throw new ArgumentException("argument must not be negative"); + return UInt128.CeilingSqrt(a.v); + } + + public static long FloorCbrt(Int128 a) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + return -(long)UInt128.FloorCbrt(aneg); + } + return (long)UInt128.FloorCbrt(a.v); + } + + public static long CeilingCbrt(Int128 a) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + return -(long)UInt128.CeilingCbrt(aneg); + } + return (long)UInt128.CeilingCbrt(a.v); + } + + public static Int128 Min(Int128 a, Int128 b) + { + if (LessThan(ref a.v, ref b.v)) + return a; + return b; + } + + public static Int128 Max(Int128 a, Int128 b) + { + if (LessThan(ref b.v, ref a.v)) + return a; + return b; + } + + public static double Log(Int128 a) + { + return Log(a, Math.E); + } + + public static double Log10(Int128 a) + { + return Log(a, 10); + } + + public static double Log(Int128 a, double b) + { + if (a.IsNegative || a.IsZero) + throw new ArgumentException("argument must be positive"); + return Math.Log(UInt128.ConvertToDouble(ref a.v), b); + } + + public static Int128 Add(Int128 a, Int128 b) + { + Int128 c; + UInt128.Add(out c.v, ref a.v, ref b.v); + return c; + } + + public static Int128 Subtract(Int128 a, Int128 b) + { + Int128 c; + UInt128.Subtract(out c.v, ref a.v, ref b.v); + return c; + } + + public static Int128 Multiply(Int128 a, Int128 b) + { + Int128 c; + Multiply(out c, ref a, ref b); + return c; + } + + public static Int128 Divide(Int128 a, Int128 b) + { + Int128 c; + Divide(out c, ref a, ref b); + return c; + } + + public static Int128 Remainder(Int128 a, Int128 b) + { + Int128 c; + Remainder(out c, ref a, ref b); + return c; + } + + public static Int128 DivRem(Int128 a, Int128 b, out Int128 remainder) + { + Int128 c; + Divide(out c, ref a, ref b); + Remainder(out remainder, ref a, ref b); + return c; + } + + public static Int128 Negate(Int128 a) + { + Int128 c; + UInt128.Negate(out c.v, ref a.v); + return c; + } + + public static Int128 GreatestCommonDivisor(Int128 a, Int128 b) + { + Int128 c; + GreatestCommonDivisor(out c, ref a, ref b); + return c; + } + + public static void GreatestCommonDivisor(out Int128 c, ref Int128 a, ref Int128 b) + { + if (a.IsNegative) + { + UInt128 aneg; + UInt128.Negate(out aneg, ref a.v); + if (b.IsNegative) + { + UInt128 bneg; + UInt128.Negate(out bneg, ref b.v); + UInt128.GreatestCommonDivisor(out c.v, ref aneg, ref bneg); + } + else + UInt128.GreatestCommonDivisor(out c.v, ref aneg, ref b.v); + } + else + { + if (b.IsNegative) + { + UInt128 bneg; + UInt128.Negate(out bneg, ref b.v); + UInt128.GreatestCommonDivisor(out c.v, ref a.v, ref bneg); + } + else + UInt128.GreatestCommonDivisor(out c.v, ref a.v, ref b.v); + } + } + + public static void LeftShift(ref Int128 c, int d) + { + UInt128.LeftShift(ref c.v, d); + } + + public static void LeftShift(ref Int128 c) + { + UInt128.LeftShift(ref c.v); + } + + public static void RightShift(ref Int128 c, int d) + { + UInt128.ArithmeticRightShift(ref c.v, d); + } + + public static void RightShift(ref Int128 c) + { + UInt128.ArithmeticRightShift(ref c.v); + } + + public static void Swap(ref Int128 a, ref Int128 b) + { + UInt128.Swap(ref a.v, ref b.v); + } + + public static int Compare(Int128 a, Int128 b) + { + return a.CompareTo(b); + } + + public static void Shift(out Int128 c, ref Int128 a, int d) + { + UInt128.ArithmeticShift(out c.v, ref a.v, d); + } + + public static void Shift(ref Int128 c, int d) + { + UInt128.ArithmeticShift(ref c.v, d); + } + + public static Int128 ModAdd(Int128 a, Int128 b, Int128 modulus) + { + Int128 c; + UInt128.ModAdd(out c.v, ref a.v, ref b.v, ref modulus.v); + return c; + } + + public static Int128 ModSub(Int128 a, Int128 b, Int128 modulus) + { + Int128 c; + UInt128.ModSub(out c.v, ref a.v, ref b.v, ref modulus.v); + return c; + } + + public static Int128 ModMul(Int128 a, Int128 b, Int128 modulus) + { + Int128 c; + UInt128.ModMul(out c.v, ref a.v, ref b.v, ref modulus.v); + return c; + } + + public static Int128 ModPow(Int128 value, Int128 exponent, Int128 modulus) + { +Int128 result; + UInt128.ModPow(out result.v, ref value.v, ref exponent.v, ref modulus.v); + return result; + } + } +} diff --git a/Dirichlet.Numerics/Properties/AssemblyInfo.cs b/Dirichlet.Numerics/Properties/AssemblyInfo.cs new file mode 100644 index 000000000..0c435d92b --- /dev/null +++ b/Dirichlet.Numerics/Properties/AssemblyInfo.cs @@ -0,0 +1,36 @@ +using System.Reflection; +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; + +// General Information about an assembly is controlled through the following +// set of attributes. Change these attribute values to modify the information +// associated with an assembly. +[assembly: AssemblyTitle("Dirichlet.Numerics")] +[assembly: AssemblyDescription("")] +[assembly: AssemblyConfiguration("")] +[assembly: AssemblyCompany("")] +[assembly: AssemblyProduct("Dirichlet.Numerics")] +[assembly: AssemblyCopyright("Copyright © 2012")] +[assembly: AssemblyTrademark("")] +[assembly: AssemblyCulture("")] + +// Setting ComVisible to false makes the types in this assembly not visible +// to COM components. If you need to access a type in this assembly from +// COM, set the ComVisible attribute to true on that type. +[assembly: ComVisible(false)] + +// The following GUID is for the ID of the typelib if this project is exposed to COM +[assembly: Guid("a588abc0-c631-4475-9eab-4163d9f81fed")] + +// Version information for an assembly consists of the following four values: +// +// Major Version +// Minor Version +// Build Number +// Revision +// +// You can specify all the values or you can default the Build and Revision Numbers +// by using the '*' as shown below: +// [assembly: AssemblyVersion("1.0.*")] +[assembly: AssemblyVersion("1.0.0.0")] +[assembly: AssemblyFileVersion("1.0.0.0")] diff --git a/Dirichlet.Numerics/UInt128.cs b/Dirichlet.Numerics/UInt128.cs new file mode 100644 index 000000000..c38671769 --- /dev/null +++ b/Dirichlet.Numerics/UInt128.cs @@ -0,0 +1,2638 @@ +using System; +using System.Diagnostics; +using System.Globalization; +using System.Linq; +using System.Numerics; + +namespace Dirichlet.Numerics +{ + public struct UInt128 : IFormattable, IComparable, IComparable, IEquatable + { + private struct UInt256 + { + public ulong s0; + public ulong s1; + public ulong s2; + public ulong s3; + + public uint r0 { get { return (uint)s0; } } + public uint r1 { get { return (uint)(s0 >> 32); } } + public uint r2 { get { return (uint)s1; } } + public uint r3 { get { return (uint)(s1 >> 32); } } + public uint r4 { get { return (uint)s2; } } + public uint r5 { get { return (uint)(s2 >> 32); } } + public uint r6 { get { return (uint)s3; } } + public uint r7 { get { return (uint)(s3 >> 32); } } + + public UInt128 t0 { get { UInt128 result; UInt128.Create(out result, s0, s1); return result; } } + public UInt128 t1 { get { UInt128 result; UInt128.Create(out result, s2, s3); return result; } } + + public static implicit operator BigInteger(UInt256 a) + { + return (BigInteger)a.s3 << 192 | (BigInteger)a.s2 << 128 | (BigInteger)a.s1 << 64 | a.s0; + } + + public override string ToString() + { + return ((BigInteger)this).ToString(); + } + } + + private ulong s0; + private ulong s1; + + private static readonly UInt128 maxValue = ~(UInt128)0; + private static readonly UInt128 zero = (UInt128)0; + private static readonly UInt128 one = (UInt128)1; + + public static UInt128 MinValue { get { return zero; } } + public static UInt128 MaxValue { get { return maxValue; } } + public static UInt128 Zero { get { return zero; } } + public static UInt128 One { get { return one; } } + + public static UInt128 Parse(string value) + { + UInt128 c; + if (!TryParse(value, out c)) + throw new FormatException(); + return c; + } + + public static bool TryParse(string value, out UInt128 result) + { + return TryParse(value, NumberStyles.Integer, NumberFormatInfo.CurrentInfo, out result); + } + + public static bool TryParse(string value, NumberStyles style, IFormatProvider provider, out UInt128 result) + { + BigInteger a; + if (!BigInteger.TryParse(value, style, provider, out a)) + { + result = Zero; + return false; + } + Create(out result, a); + return true; + } + + public UInt128(long value) + { + Create(out this, value); + } + + public UInt128(ulong value) + { + Create(out this, value); + } + + public UInt128(decimal value) + { + Create(out this, value); + } + + public UInt128(double value) + { + Create(out this, value); + } + + public UInt128(BigInteger value) + { + Create(out this, value); + } + + public static void Create(out UInt128 c, uint r0, uint r1, uint r2, uint r3) + { + c.s0 = (ulong)r1 << 32 | r0; + c.s1 = (ulong)r3 << 32 | r2; + } + + public static void Create(out UInt128 c, ulong s0, ulong s1) + { + c.s0 = s0; + c.s1 = s1; + } + + public static void Create(out UInt128 c, long a) + { + c.s0 = (ulong)a; + c.s1 = a < 0 ? ulong.MaxValue : 0; + } + + public static void Create(out UInt128 c, ulong a) + { + c.s0 = a; + c.s1 = 0; + } + + public static void Create(out UInt128 c, decimal a) + { + var bits = decimal.GetBits(decimal.Truncate(a)); + Create(out c, (uint)bits[0], (uint)bits[1], (uint)bits[2], 0); + if (a < 0) + Negate(ref c); + } + + public static void Create(out UInt128 c, BigInteger a) + { + var sign = a.Sign; + if (sign == -1) + a = -a; + c.s0 = (ulong)(a & ulong.MaxValue); + c.s1 = (ulong)(a >> 64); + if (sign == -1) + Negate(ref c); + } + + public static void Create(out UInt128 c, double a) + { + var negate = false; + if (a < 0) + { + negate = true; + a = -a; + } + if (a <= ulong.MaxValue) + { + c.s0 = (ulong)a; + c.s1 = 0; + } + else + { + var shift = Math.Max((int)Math.Ceiling(Math.Log(a, 2)) - 63, 0); + c.s0 = (ulong)(a / Math.Pow(2, shift)); + c.s1 = 0; + LeftShift(ref c, shift); + } + if (negate) + Negate(ref c); + } + + private uint r0 { get { return (uint)s0; } } + private uint r1 { get { return (uint)(s0 >> 32); } } + private uint r2 { get { return (uint)s1; } } + private uint r3 { get { return (uint)(s1 >> 32); } } + + public ulong S0 { get { return s0; } } + public ulong S1 { get { return s1; } } + + public bool IsZero { get { return (s0 | s1) == 0; } } + public bool IsOne { get { return (s1 ^ s0) == 1; } } + public bool IsPowerOfTwo { get { return (this & (this - 1)).IsZero; } } + public bool IsEven { get { return (s0 & 1) == 0; } } + public int Sign { get { return IsZero ? 0 : 1; } } + + public override string ToString() + { + return ((BigInteger)this).ToString(); + } + + public string ToString(string format) + { + return ((BigInteger)this).ToString(format); + } + + public string ToString(IFormatProvider provider) + { + return ToString(null, provider); + } + + public string ToString(string format, IFormatProvider provider) + { + return ((BigInteger)this).ToString(format, provider); + } + + public static explicit operator UInt128(double a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static explicit operator UInt128(sbyte a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static implicit operator UInt128(byte a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static explicit operator UInt128(short a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static implicit operator UInt128(ushort a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static explicit operator UInt128(int a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static implicit operator UInt128(uint a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static explicit operator UInt128(long a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static implicit operator UInt128(ulong a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static explicit operator UInt128(decimal a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static explicit operator UInt128(BigInteger a) + { + UInt128 c; + Create(out c, a); + return c; + } + + public static explicit operator float(UInt128 a) + { + return ConvertToFloat(ref a); + } + + public static explicit operator double(UInt128 a) + { + return ConvertToDouble(ref a); + } + + public static float ConvertToFloat(ref UInt128 a) + { + if (a.s1 == 0) + return a.s0; + return a.s1 * (float)ulong.MaxValue + a.s0; + } + + public static double ConvertToDouble(ref UInt128 a) + { + if (a.s1 == 0) + return a.s0; + return a.s1 * (double)ulong.MaxValue + a.s0; + } + + public static explicit operator sbyte(UInt128 a) + { + return (sbyte)a.s0; + } + + public static explicit operator byte(UInt128 a) + { + return (byte)a.s0; + } + + public static explicit operator short(UInt128 a) + { + return (short)a.s0; + } + + public static explicit operator ushort(UInt128 a) + { + return (ushort)a.s0; + } + + public static explicit operator int(UInt128 a) + { + return (int)a.s0; + } + + public static explicit operator uint(UInt128 a) + { + return (uint)a.s0; + } + + public static explicit operator long(UInt128 a) + { + return (long)a.s0; + } + + public static explicit operator ulong(UInt128 a) + { + return a.s0; + } + + public static explicit operator decimal(UInt128 a) + { + if (a.s1 == 0) + return a.s0; + var shift = Math.Max(0, 32 - GetBitLength(a.s1)); + UInt128 ashift; + RightShift(out ashift, ref a, shift); + return new decimal((int)a.r0, (int)a.r1, (int)a.r2, false, (byte)shift); + } + + public static implicit operator BigInteger(UInt128 a) + { + if (a.s1 == 0) + return a.s0; + return (BigInteger)a.s1 << 64 | a.s0; + } + + public static UInt128 operator <<(UInt128 a, int b) + { + UInt128 c; + LeftShift(out c, ref a, b); + return c; + } + + public static UInt128 operator >>(UInt128 a, int b) + { + UInt128 c; + RightShift(out c, ref a, b); + return c; + } + + public static UInt128 operator &(UInt128 a, UInt128 b) + { + UInt128 c; + And(out c, ref a, ref b); + return c; + } + + public static uint operator &(UInt128 a, uint b) + { + return (uint)a.s0 & b; + } + + public static uint operator &(uint a, UInt128 b) + { + return a & (uint)b.s0; + } + + public static ulong operator &(UInt128 a, ulong b) + { + return a.s0 & b; + } + + public static ulong operator &(ulong a, UInt128 b) + { + return a & b.s0; + } + + public static UInt128 operator |(UInt128 a, UInt128 b) + { + UInt128 c; + Or(out c, ref a, ref b); + return c; + } + + public static UInt128 operator ^(UInt128 a, UInt128 b) + { + UInt128 c; + ExclusiveOr(out c, ref a, ref b); + return c; + } + + public static UInt128 operator ~(UInt128 a) + { + UInt128 c; + Not(out c, ref a); + return c; + } + + public static UInt128 operator +(UInt128 a, UInt128 b) + { + UInt128 c; + Add(out c, ref a, ref b); + return c; + } + + public static UInt128 operator +(UInt128 a, ulong b) + { + UInt128 c; + Add(out c, ref a, b); + return c; + } + + public static UInt128 operator +(ulong a, UInt128 b) + { + UInt128 c; + Add(out c, ref b, a); + return c; + } + + public static UInt128 operator ++(UInt128 a) + { + UInt128 c; + Add(out c, ref a, 1); + return c; + } + + public static UInt128 operator -(UInt128 a, UInt128 b) + { + UInt128 c; + Subtract(out c, ref a, ref b); + return c; + } + + public static UInt128 operator -(UInt128 a, ulong b) + { + UInt128 c; + Subtract(out c, ref a, b); + return c; + } + + public static UInt128 operator -(ulong a, UInt128 b) + { + UInt128 c; + Subtract(out c, a, ref b); + return c; + } + + public static UInt128 operator --(UInt128 a) + { + UInt128 c; + Subtract(out c, ref a, 1); + return c; + } + + public static UInt128 operator +(UInt128 a) + { + return a; + } + + public static UInt128 operator *(UInt128 a, uint b) + { + UInt128 c; + Multiply(out c, ref a, b); + return c; + } + + public static UInt128 operator *(uint a, UInt128 b) + { + UInt128 c; + Multiply(out c, ref b, a); + return c; + } + + public static UInt128 operator *(UInt128 a, ulong b) + { + UInt128 c; + Multiply(out c, ref a, b); + return c; + } + + public static UInt128 operator *(ulong a, UInt128 b) + { + UInt128 c; + Multiply(out c, ref b, a); + return c; + } + + public static UInt128 operator *(UInt128 a, UInt128 b) + { + UInt128 c; + Multiply(out c, ref a, ref b); + return c; + } + + public static UInt128 operator /(UInt128 a, ulong b) + { + UInt128 c; + Divide(out c, ref a, b); + return c; + } + + public static UInt128 operator /(UInt128 a, UInt128 b) + { + UInt128 c; + Divide(out c, ref a, ref b); + return c; + } + + public static ulong operator %(UInt128 a, uint b) + { + return Remainder(ref a, b); + } + + public static ulong operator %(UInt128 a, ulong b) + { + return Remainder(ref a, b); + } + + public static UInt128 operator %(UInt128 a, UInt128 b) + { + UInt128 c; + Remainder(out c, ref a, ref b); + return c; + } + + public static bool operator <(UInt128 a, UInt128 b) + { + return LessThan(ref a, ref b); + } + + public static bool operator <(UInt128 a, int b) + { + return LessThan(ref a, b); + } + + public static bool operator <(int a, UInt128 b) + { + return LessThan(a, ref b); + } + + public static bool operator <(UInt128 a, uint b) + { + return LessThan(ref a, b); + } + + public static bool operator <(uint a, UInt128 b) + { + return LessThan(a, ref b); + } + + public static bool operator <(UInt128 a, long b) + { + return LessThan(ref a, b); + } + + public static bool operator <(long a, UInt128 b) + { + return LessThan(a, ref b); + } + + public static bool operator <(UInt128 a, ulong b) + { + return LessThan(ref a, b); + } + + public static bool operator <(ulong a, UInt128 b) + { + return LessThan(a, ref b); + } + + public static bool operator <=(UInt128 a, UInt128 b) + { + return !LessThan(ref b, ref a); + } + + public static bool operator <=(UInt128 a, int b) + { + return !LessThan(b, ref a); + } + + public static bool operator <=(int a, UInt128 b) + { + return !LessThan(ref b, a); + } + + public static bool operator <=(UInt128 a, uint b) + { + return !LessThan(b, ref a); + } + + public static bool operator <=(uint a, UInt128 b) + { + return !LessThan(ref b, a); + } + + public static bool operator <=(UInt128 a, long b) + { + return !LessThan(b, ref a); + } + + public static bool operator <=(long a, UInt128 b) + { + return !LessThan(ref b, a); + } + + public static bool operator <=(UInt128 a, ulong b) + { + return !LessThan(b, ref a); + } + + public static bool operator <=(ulong a, UInt128 b) + { + return !LessThan(ref b, a); + } + + public static bool operator >(UInt128 a, UInt128 b) + { + return LessThan(ref b, ref a); + } + + public static bool operator >(UInt128 a, int b) + { + return LessThan(b, ref a); + } + + public static bool operator >(int a, UInt128 b) + { + return LessThan(ref b, a); + } + + public static bool operator >(UInt128 a, uint b) + { + return LessThan(b, ref a); + } + + public static bool operator >(uint a, UInt128 b) + { + return LessThan(ref b, a); + } + + public static bool operator >(UInt128 a, long b) + { + return LessThan(b, ref a); + } + + public static bool operator >(long a, UInt128 b) + { + return LessThan(ref b, a); + } + + public static bool operator >(UInt128 a, ulong b) + { + return LessThan(b, ref a); + } + + public static bool operator >(ulong a, UInt128 b) + { + return LessThan(ref b, a); + } + + public static bool operator >=(UInt128 a, UInt128 b) + { + return !LessThan(ref a, ref b); + } + + public static bool operator >=(UInt128 a, int b) + { + return !LessThan(ref a, b); + } + + public static bool operator >=(int a, UInt128 b) + { + return !LessThan(a, ref b); + } + + public static bool operator >=(UInt128 a, uint b) + { + return !LessThan(ref a, b); + } + + public static bool operator >=(uint a, UInt128 b) + { + return !LessThan(a, ref b); + } + + public static bool operator >=(UInt128 a, long b) + { + return !LessThan(ref a, b); + } + + public static bool operator >=(long a, UInt128 b) + { + return !LessThan(a, ref b); + } + + public static bool operator >=(UInt128 a, ulong b) + { + return !LessThan(ref a, b); + } + + public static bool operator >=(ulong a, UInt128 b) + { + return !LessThan(a, ref b); + } + + public static bool operator ==(UInt128 a, UInt128 b) + { + return a.Equals(b); + } + + public static bool operator ==(UInt128 a, int b) + { + return a.Equals(b); + } + + public static bool operator ==(int a, UInt128 b) + { + return b.Equals(a); + } + + public static bool operator ==(UInt128 a, uint b) + { + return a.Equals(b); + } + + public static bool operator ==(uint a, UInt128 b) + { + return b.Equals(a); + } + + public static bool operator ==(UInt128 a, long b) + { + return a.Equals(b); + } + + public static bool operator ==(long a, UInt128 b) + { + return b.Equals(a); + } + + public static bool operator ==(UInt128 a, ulong b) + { + return a.Equals(b); + } + + public static bool operator ==(ulong a, UInt128 b) + { + return b.Equals(a); + } + + public static bool operator !=(UInt128 a, UInt128 b) + { + return !a.Equals(b); + } + + public static bool operator !=(UInt128 a, int b) + { + return !a.Equals(b); + } + + public static bool operator !=(int a, UInt128 b) + { + return !b.Equals(a); + } + + public static bool operator !=(UInt128 a, uint b) + { + return !a.Equals(b); + } + + public static bool operator !=(uint a, UInt128 b) + { + return !b.Equals(a); + } + + public static bool operator !=(UInt128 a, long b) + { + return !a.Equals(b); + } + + public static bool operator !=(long a, UInt128 b) + { + return !b.Equals(a); + } + + public static bool operator !=(UInt128 a, ulong b) + { + return !a.Equals(b); + } + + public static bool operator !=(ulong a, UInt128 b) + { + return !b.Equals(a); + } + + public int CompareTo(UInt128 other) + { + if (s1 != other.s1) + return s1.CompareTo(other.s1); + return s0.CompareTo(other.s0); + } + + public int CompareTo(int other) + { + if (s1 != 0 || other < 0) + return 1; + return s0.CompareTo((ulong)other); + } + + public int CompareTo(uint other) + { + if (s1 != 0) + return 1; + return s0.CompareTo((ulong)other); + } + + public int CompareTo(long other) + { + if (s1 != 0 || other < 0) + return 1; + return s0.CompareTo((ulong)other); + } + + public int CompareTo(ulong other) + { + if (s1 != 0) + return 1; + return s0.CompareTo(other); + } + + public int CompareTo(object obj) + { + if (obj == null) + return 1; + if (!(obj is UInt128)) + throw new ArgumentException(); + return CompareTo((UInt128)obj); + } + + private static bool LessThan(ref UInt128 a, long b) + { + return b >= 0 && a.s1 == 0 && a.s0 < (ulong)b; + } + + private static bool LessThan(long a, ref UInt128 b) + { + return a < 0 || b.s1 != 0 || (ulong)a < b.s0; + } + + private static bool LessThan(ref UInt128 a, ulong b) + { + return a.s1 == 0 && a.s0 < b; + } + + private static bool LessThan(ulong a, ref UInt128 b) + { + return b.s1 != 0 || a < b.s0; + } + + private static bool LessThan(ref UInt128 a, ref UInt128 b) + { + if (a.s1 != b.s1) + return a.s1 < b.s1; + return a.s0 < b.s0; + } + + public static bool Equals(ref UInt128 a, ref UInt128 b) + { + return a.s0 == b.s0 && a.s1 == b.s1; + } + + public bool Equals(UInt128 other) + { + return s0 == other.s0 && s1 == other.s1; + } + + public bool Equals(int other) + { + return other >= 0 && s0 == (uint)other && s1 == 0; + } + + public bool Equals(uint other) + { + return s0 == other && s1 == 0; + } + + public bool Equals(long other) + { + return other >= 0 && s0 == (ulong)other && s1 == 0; + } + + public bool Equals(ulong other) + { + return s0 == other && s1 == 0; + } + + public override bool Equals(object obj) + { + if (!(obj is UInt128)) + return false; + return Equals((UInt128)obj); + } + + public override int GetHashCode() + { + return s0.GetHashCode() ^ s1.GetHashCode(); + } + + public static void Multiply(out UInt128 c, ulong a, ulong b) + { + Multiply64(out c, a, b); + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b); + } + + public static void Multiply(out UInt128 c, ref UInt128 a, uint b) + { + if (a.s1 == 0) + Multiply64(out c, a.s0, b); + else + Multiply128(out c, ref a, b); + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b % ((BigInteger)1 << 128)); + } + + public static void Multiply(out UInt128 c, ref UInt128 a, ulong b) + { + if (a.s1 == 0) + Multiply64(out c, a.s0, b); + else + Multiply128(out c, ref a, b); + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b % ((BigInteger)1 << 128)); + } + + public static void Multiply(out UInt128 c, ref UInt128 a, ref UInt128 b) + { + if ((a.s1 | b.s1) == 0) + Multiply64(out c, a.s0, b.s0); + else if (a.s1 == 0) + Multiply128(out c, ref b, a.s0); + else if (b.s1 == 0) + Multiply128(out c, ref a, b.s0); + else + Multiply128(out c, ref a, ref b); + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b % ((BigInteger)1 << 128)); + } + + private static void Multiply(out UInt256 c, ref UInt128 a, ref UInt128 b) + { +#if true + UInt128 c00, c01, c10, c11; + Multiply64(out c00, a.s0, b.s0); + Multiply64(out c01, a.s0, b.s1); + Multiply64(out c10, a.s1, b.s0); + Multiply64(out c11, a.s1, b.s1); + var carry1 = (uint)0; + var carry2 = (uint)0; + c.s0 = c00.S0; + c.s1 = Add(Add(c00.s1, c01.s0, ref carry1), c10.s0, ref carry1); + c.s2 = Add(Add(Add(c01.s1, c10.s1, ref carry2), c11.s0, ref carry2), carry1, ref carry2); + c.s3 = c11.s1 + carry2; +#else + // Karatsuba method. + // Warning: doesn't correctly handle overflow. + UInt128 z0, z1, z2; + Multiply64(out z0, a.s0, b.s0); + Multiply64(out z2, a.s1, b.s1); + Multiply64(out z1, a.s0 + a.s1, b.s0 + b.s1); + Subtract(ref z1, ref z2); + Subtract(ref z1, ref z0); + var carry1 = (uint)0; + var carry2 = (uint)0; + c.s0 = z0.S0; + c.s1 = Add(z0.s1, z1.s0, ref carry1); + c.s2 = Add(Add(z1.s1, z2.s0, ref carry2), carry1, ref carry2); + c.s3 = z2.s1 + carry2; +#endif + Debug.Assert((BigInteger)c == (BigInteger)a * (BigInteger)b); + } + + public static UInt128 Abs(UInt128 a) + { + return a; + } + + public static UInt128 Square(ulong a) + { + UInt128 c; + Square(out c, a); + return c; + } + + public static UInt128 Square(UInt128 a) + { + UInt128 c; + Square(out c, ref a); + return c; + } + + public static void Square(out UInt128 c, ulong a) + { + Square64(out c, a); + } + + public static void Square(out UInt128 c, ref UInt128 a) + { + if (a.s1 == 0) + Square64(out c, a.s0); + else + Multiply128(out c, ref a, ref a); + } + + public static UInt128 Cube(ulong a) + { + UInt128 c; + Cube(out c, a); + return c; + } + + public static UInt128 Cube(UInt128 a) + { + UInt128 c; + Cube(out c, ref a); + return c; + } + + public static void Cube(out UInt128 c, ulong a) + { + UInt128 square; + Square(out square, a); + Multiply(out c, ref square, a); + } + + public static void Cube(out UInt128 c, ref UInt128 a) + { + UInt128 square; + if (a.s1 == 0) + { + Square64(out square, a.s0); + Multiply(out c, ref square, a.s0); + } + else + { + Multiply128(out square, ref a, ref a); + Multiply128(out c, ref square, ref a); + } + } + + public static void Add(out UInt128 c, ulong a, ulong b) + { + c.s0 = a + b; + c.s1 = 0; + if (c.s0 < a && c.s0 < b) + ++c.s1; + Debug.Assert((BigInteger)c == ((BigInteger)a + (BigInteger)b)); + } + + public static void Add(out UInt128 c, ref UInt128 a, ulong b) + { + c.s0 = a.s0 + b; + c.s1 = a.s1; + if (c.s0 < a.s0 && c.s0 < b) + ++c.s1; + Debug.Assert((BigInteger)c == ((BigInteger)a + (BigInteger)b) % ((BigInteger)1 << 128)); + } + + public static void Add(out UInt128 c, ref UInt128 a, ref UInt128 b) + { + c.s0 = a.s0 + b.s0; + c.s1 = a.s1 + b.s1; + if (c.s0 < a.s0 && c.s0 < b.s0) + ++c.s1; + Debug.Assert((BigInteger)c == ((BigInteger)a + (BigInteger)b) % ((BigInteger)1 << 128)); + } + + private static ulong Add(ulong a, ulong b, ref uint carry) + { + var c = a + b; + if (c < a && c < b) + ++carry; + return c; + } + + public static void Add(ref UInt128 a, ulong b) + { + var sum = a.s0 + b; + if (sum < a.s0 && sum < b) + ++a.s1; + a.s0 = sum; + } + + public static void Add(ref UInt128 a, ref UInt128 b) + { + var sum = a.s0 + b.s0; + if (sum < a.s0 && sum < b.s0) + ++a.s1; + a.s0 = sum; + a.s1 += b.s1; + } + + public static void Add(ref UInt128 a, UInt128 b) + { + Add(ref a, ref b); + } + + public static void Subtract(out UInt128 c, ref UInt128 a, ulong b) + { + c.s0 = a.s0 - b; + c.s1 = a.s1; + if (a.s0 < b) + --c.s1; + Debug.Assert((BigInteger)c == ((BigInteger)a - (BigInteger)b + ((BigInteger)1 << 128)) % ((BigInteger)1 << 128)); + } + + public static void Subtract(out UInt128 c, ulong a, ref UInt128 b) + { + c.s0 = a - b.s0; + c.s1 = 0 - b.s1; + if (a < b.s0) + --c.s1; + Debug.Assert((BigInteger)c == ((BigInteger)a - (BigInteger)b + ((BigInteger)1 << 128)) % ((BigInteger)1 << 128)); + } + + public static void Subtract(out UInt128 c, ref UInt128 a, ref UInt128 b) + { + c.s0 = a.s0 - b.s0; + c.s1 = a.s1 - b.s1; + if (a.s0 < b.s0) + --c.s1; + Debug.Assert((BigInteger)c == ((BigInteger)a - (BigInteger)b + ((BigInteger)1 << 128)) % ((BigInteger)1 << 128)); + } + + public static void Subtract(ref UInt128 a, ulong b) + { + if (a.s0 < b) + --a.s1; + a.s0 -= b; + } + + public static void Subtract(ref UInt128 a, ref UInt128 b) + { + if (a.s0 < b.s0) + --a.s1; + a.s0 -= b.s0; + a.s1 -= b.s1; + } + + public static void Subtract(ref UInt128 a, UInt128 b) + { + Subtract(ref a, ref b); + } + + private static void Square64(out UInt128 w, ulong u) + { + var u0 = (ulong)(uint)u; + var u1 = u >> 32; + var carry = u0 * u0; + var r0 = (uint)carry; + var u0u1 = u0 * u1; + carry = (carry >> 32) + u0u1; + var r2 = carry >> 32; + carry = (uint)carry + u0u1; + w.s0 = carry << 32 | r0; + w.s1 = (carry >> 32) + r2 + u1 * u1; + Debug.Assert((BigInteger)w == (BigInteger)u * u); + } + + private static void Multiply64(out UInt128 w, uint u, uint v) + { + w.s0 = (ulong)u * v; + w.s1 = 0; + Debug.Assert((BigInteger)w == (BigInteger)u * v); + } + + private static void Multiply64(out UInt128 w, ulong u, uint v) + { + var u0 = (ulong)(uint)u; + var u1 = u >> 32; + var carry = u0 * v; + var r0 = (uint)carry; + carry = (carry >> 32) + u1 * v; + w.s0 = carry << 32 | r0; + w.s1 = carry >> 32; + Debug.Assert((BigInteger)w == (BigInteger)u * v); + } + + private static void Multiply64(out UInt128 w, ulong u, ulong v) + { + var u0 = (ulong)(uint)u; + var u1 = u >> 32; + var v0 = (ulong)(uint)v; + var v1 = v >> 32; + var carry = u0 * v0; + var r0 = (uint)carry; + carry = (carry >> 32) + u0 * v1; + var r2 = carry >> 32; + carry = (uint)carry + u1 * v0; + w.s0 = carry << 32 | r0; + w.s1 = (carry >> 32) + r2 + u1 * v1; + Debug.Assert((BigInteger)w == (BigInteger)u * v); + } + + private static void Multiply64(out UInt128 w, ulong u, ulong v, ulong c) + { + var u0 = (ulong)(uint)u; + var u1 = u >> 32; + var v0 = (ulong)(uint)v; + var v1 = v >> 32; + var carry = u0 * v0 + (uint)c; + var r0 = (uint)carry; + carry = (carry >> 32) + u0 * v1 + (c >> 32); + var r2 = carry >> 32; + carry = (uint)carry + u1 * v0; + w.s0 = carry << 32 | r0; + w.s1 = (carry >> 32) + r2 + u1 * v1; + Debug.Assert((BigInteger)w == (BigInteger)u * v + c); + } + + private static ulong MultiplyHigh64(ulong u, ulong v, ulong c) + { + var u0 = (ulong)(uint)u; + var u1 = u >> 32; + var v0 = (ulong)(uint)v; + var v1 = v >> 32; + var carry = ((u0 * v0 + (uint)c) >> 32) + u0 * v1 + (c >> 32); + var r2 = carry >> 32; + carry = (uint)carry + u1 * v0; + return (carry >> 32) + r2 + u1 * v1; + } + + private static void Multiply128(out UInt128 w, ref UInt128 u, uint v) + { + Multiply64(out w, u.s0, v); + w.s1 += u.s1 * v; + Debug.Assert((BigInteger)w == (BigInteger)u * v % ((BigInteger)1 << 128)); + } + + private static void Multiply128(out UInt128 w, ref UInt128 u, ulong v) + { + Multiply64(out w, u.s0, v); + w.s1 += u.s1 * v; + Debug.Assert((BigInteger)w == (BigInteger)u * v % ((BigInteger)1 << 128)); + } + + private static void Multiply128(out UInt128 w, ref UInt128 u, ref UInt128 v) + { + Multiply64(out w, u.s0, v.s0); + w.s1 += u.s1 * v.s0 + u.s0 * v.s1; + Debug.Assert((BigInteger)w == (BigInteger)u * v % ((BigInteger)1 << 128)); + } + + public static void Divide(out UInt128 w, ref UInt128 u, uint v) + { + if (u.s1 == 0) + Divide64(out w, u.s0, v); + else if (u.s1 <= uint.MaxValue) + Divide96(out w, ref u, v); + else + Divide128(out w, ref u, v); + } + + public static void Divide(out UInt128 w, ref UInt128 u, ulong v) + { + if (u.s1 == 0) + Divide64(out w, u.s0, v); + else + { + var v0 = (uint)v; + if (v == v0) + { + if (u.s1 <= uint.MaxValue) + Divide96(out w, ref u, v0); + else + Divide128(out w, ref u, v0); + } + else + { + if (u.s1 <= uint.MaxValue) + Divide96(out w, ref u, v); + else + Divide128(out w, ref u, v); + } + } + } + + public static void Divide(out UInt128 c, ref UInt128 a, ref UInt128 b) + { + if (LessThan(ref a, ref b)) + c = Zero; + else if (b.s1 == 0) + Divide(out c, ref a, b.s0); + else if (b.s1 <= uint.MaxValue) + { + UInt128 rem; + Create(out c, DivRem96(out rem, ref a, ref b)); + } + else + { + UInt128 rem; + Create(out c, DivRem128(out rem, ref a, ref b)); + } + } + + public static uint Remainder(ref UInt128 u, uint v) + { + if (u.s1 == 0) + return (uint)(u.s0 % v); + if (u.s1 <= uint.MaxValue) + return Remainder96(ref u, v); + return Remainder128(ref u, v); + } + + public static ulong Remainder(ref UInt128 u, ulong v) + { + if (u.s1 == 0) + return u.s0 % v; + var v0 = (uint)v; + if (v == v0) + { + if (u.s1 <= uint.MaxValue) + return Remainder96(ref u, v0); + return Remainder128(ref u, v0); + } + if (u.s1 <= uint.MaxValue) + return Remainder96(ref u, v); + return Remainder128(ref u, v); + } + + public static void Remainder(out UInt128 c, ref UInt128 a, ref UInt128 b) + { + if (LessThan(ref a, ref b)) + c = a; + else if (b.s1 == 0) + Create(out c, Remainder(ref a, b.s0)); + else if (b.s1 <= uint.MaxValue) + DivRem96(out c, ref a, ref b); + else + DivRem128(out c, ref a, ref b); + } + + public static void Remainder(ref UInt128 a, ref UInt128 b) + { + UInt128 a2 = a; + Remainder(out a, ref a2, ref b); + } + + private static void Remainder(out UInt128 c, ref UInt256 a, ref UInt128 b) + { + if (b.r3 == 0) + Remainder192(out c, ref a, ref b); + else + Remainder256(out c, ref a, ref b); + } + + private static void Divide64(out UInt128 w, ulong u, ulong v) + { + w.s1 = 0; + w.s0 = u / v; + Debug.Assert((BigInteger)w == (BigInteger)u / v); + } + + private static void Divide96(out UInt128 w, ref UInt128 u, uint v) + { + var r2 = u.r2; + var w2 = r2 / v; + var u0 = (ulong)(r2 - w2 * v); + var u0u1 = u0 << 32 | u.r1; + var w1 = (uint)(u0u1 / v); + u0 = u0u1 - w1 * v; + u0u1 = u0 << 32 | u.r0; + var w0 = (uint)(u0u1 / v); + w.s1 = w2; + w.s0 = (ulong)w1 << 32 | w0; + Debug.Assert((BigInteger)w == (BigInteger)u / v); + } + + private static void Divide128(out UInt128 w, ref UInt128 u, uint v) + { + var r3 = u.r3; + var w3 = r3 / v; + var u0 = (ulong)(r3 - w3 * v); + var u0u1 = u0 << 32 | u.r2; + var w2 = (uint)(u0u1 / v); + u0 = u0u1 - w2 * v; + u0u1 = u0 << 32 | u.r1; + var w1 = (uint)(u0u1 / v); + u0 = u0u1 - w1 * v; + u0u1 = u0 << 32 | u.r0; + var w0 = (uint)(u0u1 / v); + w.s1 = (ulong)w3 << 32 | w2; + w.s0 = (ulong)w1 << 32 | w0; + Debug.Assert((BigInteger)w == (BigInteger)u / v); + } + + private static void Divide96(out UInt128 w, ref UInt128 u, ulong v) + { + w.s0 = w.s1 = 0; + var dneg = GetBitLength((uint)(v >> 32)); + var d = 32 - dneg; + var vPrime = v << d; + var v1 = (uint)(vPrime >> 32); + var v2 = (uint)vPrime; + var r0 = u.r0; + var r1 = u.r1; + var r2 = u.r2; + var r3 = (uint)0; + if (d != 0) + { + r3 = r2 >> dneg; + r2 = r2 << d | r1 >> dneg; + r1 = r1 << d | r0 >> dneg; + r0 <<= d; + } + var q1 = DivRem(r3, ref r2, ref r1, v1, v2); + var q0 = DivRem(r2, ref r1, ref r0, v1, v2); + w.s0 = (ulong)q1 << 32 | q0; + w.s1 = 0; + Debug.Assert((BigInteger)w == (BigInteger)u / v); + } + + private static void Divide128(out UInt128 w, ref UInt128 u, ulong v) + { + w.s0 = w.s1 = 0; + var dneg = GetBitLength((uint)(v >> 32)); + var d = 32 - dneg; + var vPrime = v << d; + var v1 = (uint)(vPrime >> 32); + var v2 = (uint)vPrime; + var r0 = u.r0; + var r1 = u.r1; + var r2 = u.r2; + var r3 = u.r3; + var r4 = (uint)0; + if (d != 0) + { + r4 = r3 >> dneg; + r3 = r3 << d | r2 >> dneg; + r2 = r2 << d | r1 >> dneg; + r1 = r1 << d | r0 >> dneg; + r0 <<= d; + } + w.s1 = DivRem(r4, ref r3, ref r2, v1, v2); + var q1 = DivRem(r3, ref r2, ref r1, v1, v2); + var q0 = DivRem(r2, ref r1, ref r0, v1, v2); + w.s0 = (ulong)q1 << 32 | q0; + Debug.Assert((BigInteger)w == (BigInteger)u / v); + } + + private static uint Remainder96(ref UInt128 u, uint v) + { + var u0 = (ulong)(u.r2 % v); + var u0u1 = u0 << 32 | u.r1; + u0 = u0u1 % v; + u0u1 = u0 << 32 | u.r0; + return (uint)(u0u1 % v); + } + + private static uint Remainder128(ref UInt128 u, uint v) + { + var u0 = (ulong)(u.r3 % v); + var u0u1 = u0 << 32 | u.r2; + u0 = u0u1 % v; + u0u1 = u0 << 32 | u.r1; + u0 = u0u1 % v; + u0u1 = u0 << 32 | u.r0; + return (uint)(u0u1 % v); + } + + private static ulong Remainder96(ref UInt128 u, ulong v) + { + var dneg = GetBitLength((uint)(v >> 32)); + var d = 32 - dneg; + var vPrime = v << d; + var v1 = (uint)(vPrime >> 32); + var v2 = (uint)vPrime; + var r0 = u.r0; + var r1 = u.r1; + var r2 = u.r2; + var r3 = (uint)0; + if (d != 0) + { + r3 = r2 >> dneg; + r2 = r2 << d | r1 >> dneg; + r1 = r1 << d | r0 >> dneg; + r0 <<= d; + } + DivRem(r3, ref r2, ref r1, v1, v2); + DivRem(r2, ref r1, ref r0, v1, v2); + return ((ulong)r1 << 32 | r0) >> d; + } + + private static ulong Remainder128(ref UInt128 u, ulong v) + { + var dneg = GetBitLength((uint)(v >> 32)); + var d = 32 - dneg; + var vPrime = v << d; + var v1 = (uint)(vPrime >> 32); + var v2 = (uint)vPrime; + var r0 = u.r0; + var r1 = u.r1; + var r2 = u.r2; + var r3 = u.r3; + var r4 = (uint)0; + if (d != 0) + { + r4 = r3 >> dneg; + r3 = r3 << d | r2 >> dneg; + r2 = r2 << d | r1 >> dneg; + r1 = r1 << d | r0 >> dneg; + r0 <<= d; + } + DivRem(r4, ref r3, ref r2, v1, v2); + DivRem(r3, ref r2, ref r1, v1, v2); + DivRem(r2, ref r1, ref r0, v1, v2); + return ((ulong)r1 << 32 | r0) >> d; + } + + private static ulong DivRem96(out UInt128 rem, ref UInt128 a, ref UInt128 b) + { + var d = 32 - GetBitLength(b.r2); + UInt128 v; + LeftShift64(out v, ref b, d); + var r4 = (uint)LeftShift64(out rem, ref a, d); + var v1 = v.r2; + var v2 = v.r1; + var v3 = v.r0; + var r3 = rem.r3; + var r2 = rem.r2; + var r1 = rem.r1; + var r0 = rem.r0; + var q1 = DivRem(r4, ref r3, ref r2, ref r1, v1, v2, v3); + var q0 = DivRem(r3, ref r2, ref r1, ref r0, v1, v2, v3); + Create(out rem, r0, r1, r2, 0); + var div = (ulong)q1 << 32 | q0; + RightShift64(ref rem, d); + Debug.Assert((BigInteger)div == (BigInteger)a / (BigInteger)b); + Debug.Assert((BigInteger)rem == (BigInteger)a % (BigInteger)b); + return div; + } + + private static uint DivRem128(out UInt128 rem, ref UInt128 a, ref UInt128 b) + { + var d = 32 - GetBitLength(b.r3); + UInt128 v; + LeftShift64(out v, ref b, d); + var r4 = (uint)LeftShift64(out rem, ref a, d); + var r3 = rem.r3; + var r2 = rem.r2; + var r1 = rem.r1; + var r0 = rem.r0; + var div = DivRem(r4, ref r3, ref r2, ref r1, ref r0, v.r3, v.r2, v.r1, v.r0); + Create(out rem, r0, r1, r2, r3); + RightShift64(ref rem, d); + Debug.Assert((BigInteger)div == (BigInteger)a / (BigInteger)b); + Debug.Assert((BigInteger)rem == (BigInteger)a % (BigInteger)b); + return div; + } + + private static void Remainder192(out UInt128 c, ref UInt256 a, ref UInt128 b) + { + var d = 32 - GetBitLength(b.r2); + UInt128 v; + LeftShift64(out v, ref b, d); + var v1 = v.r2; + var v2 = v.r1; + var v3 = v.r0; + UInt256 rem; + LeftShift64(out rem, ref a, d); + var r6 = rem.r6; + var r5 = rem.r5; + var r4 = rem.r4; + var r3 = rem.r3; + var r2 = rem.r2; + var r1 = rem.r1; + var r0 = rem.r0; + DivRem(r6, ref r5, ref r4, ref r3, v1, v2, v3); + DivRem(r5, ref r4, ref r3, ref r2, v1, v2, v3); + DivRem(r4, ref r3, ref r2, ref r1, v1, v2, v3); + DivRem(r3, ref r2, ref r1, ref r0, v1, v2, v3); + Create(out c, r0, r1, r2, 0); + RightShift64(ref c, d); + Debug.Assert((BigInteger)c == (BigInteger)a % (BigInteger)b); + } + + private static void Remainder256(out UInt128 c, ref UInt256 a, ref UInt128 b) + { + var d = 32 - GetBitLength(b.r3); + UInt128 v; + LeftShift64(out v, ref b, d); + var v1 = v.r3; + var v2 = v.r2; + var v3 = v.r1; + var v4 = v.r0; + UInt256 rem; + var r8 = (uint)LeftShift64(out rem, ref a, d); + var r7 = rem.r7; + var r6 = rem.r6; + var r5 = rem.r5; + var r4 = rem.r4; + var r3 = rem.r3; + var r2 = rem.r2; + var r1 = rem.r1; + var r0 = rem.r0; + DivRem(r8, ref r7, ref r6, ref r5, ref r4, v1, v2, v3, v4); + DivRem(r7, ref r6, ref r5, ref r4, ref r3, v1, v2, v3, v4); + DivRem(r6, ref r5, ref r4, ref r3, ref r2, v1, v2, v3, v4); + DivRem(r5, ref r4, ref r3, ref r2, ref r1, v1, v2, v3, v4); + DivRem(r4, ref r3, ref r2, ref r1, ref r0, v1, v2, v3, v4); + Create(out c, r0, r1, r2, r3); + RightShift64(ref c, d); + Debug.Assert((BigInteger)c == (BigInteger)a % (BigInteger)b); + } + + private static ulong Q(uint u0, uint u1, uint u2, uint v1, uint v2) + { + var u0u1 = (ulong)u0 << 32 | u1; + var qhat = u0 == v1 ? uint.MaxValue : u0u1 / v1; + var r = u0u1 - qhat * v1; + if (r == (uint)r && v2 * qhat > (r << 32 | u2)) + { + --qhat; + r += v1; + if (r == (uint)r && v2 * qhat > (r << 32 | u2)) + { + --qhat; + r += v1; + } + } + return qhat; + } + + private static uint DivRem(uint u0, ref uint u1, ref uint u2, uint v1, uint v2) + { + var qhat = Q(u0, u1, u2, v1, v2); + var carry = qhat * v2; + var borrow = (long)u2 - (uint)carry; + carry >>= 32; + u2 = (uint)borrow; + borrow >>= 32; + carry += qhat * v1; + borrow += (long)u1 - (uint)carry; + carry >>= 32; + u1 = (uint)borrow; + borrow >>= 32; + borrow += (long)u0 - (uint)carry; + if (borrow != 0) + { + --qhat; + carry = (ulong)u2 + v2; + u2 = (uint)carry; + carry >>= 32; + carry += (ulong)u1 + v1; + u1 = (uint)carry; + } + return (uint)qhat; + } + + private static uint DivRem(uint u0, ref uint u1, ref uint u2, ref uint u3, uint v1, uint v2, uint v3) + { + var qhat = Q(u0, u1, u2, v1, v2); + var carry = qhat * v3; + var borrow = (long)u3 - (uint)carry; + carry >>= 32; + u3 = (uint)borrow; + borrow >>= 32; + carry += qhat * v2; + borrow += (long)u2 - (uint)carry; + carry >>= 32; + u2 = (uint)borrow; + borrow >>= 32; + carry += qhat * v1; + borrow += (long)u1 - (uint)carry; + carry >>= 32; + u1 = (uint)borrow; + borrow >>= 32; + borrow += (long)u0 - (uint)carry; + if (borrow != 0) + { + --qhat; + carry = (ulong)u3 + v3; + u3 = (uint)carry; + carry >>= 32; + carry += (ulong)u2 + v2; + u2 = (uint)carry; + carry >>= 32; + carry += (ulong)u1 + v1; + u1 = (uint)carry; + } + return (uint)qhat; + } + + private static uint DivRem(uint u0, ref uint u1, ref uint u2, ref uint u3, ref uint u4, uint v1, uint v2, uint v3, uint v4) + { + var qhat = Q(u0, u1, u2, v1, v2); + var carry = qhat * v4; + var borrow = (long)u4 - (uint)carry; + carry >>= 32; + u4 = (uint)borrow; + borrow >>= 32; + carry += qhat * v3; + borrow += (long)u3 - (uint)carry; + carry >>= 32; + u3 = (uint)borrow; + borrow >>= 32; + carry += qhat * v2; + borrow += (long)u2 - (uint)carry; + carry >>= 32; + u2 = (uint)borrow; + borrow >>= 32; + carry += qhat * v1; + borrow += (long)u1 - (uint)carry; + carry >>= 32; + u1 = (uint)borrow; + borrow >>= 32; + borrow += (long)u0 - (uint)carry; + if (borrow != 0) + { + --qhat; + carry = (ulong)u4 + v4; + u4 = (uint)carry; + carry >>= 32; + carry += (ulong)u3 + v3; + u3 = (uint)carry; + carry >>= 32; + carry += (ulong)u2 + v2; + u2 = (uint)carry; + carry >>= 32; + carry += (ulong)u1 + v1; + u1 = (uint)carry; + } + return (uint)qhat; + } + + public static void ModAdd(out UInt128 c, ref UInt128 a, ref UInt128 b, ref UInt128 modulus) + { + Add(out c, ref a, ref b); + if (!LessThan(ref c, ref modulus) || LessThan(ref c, ref a) && LessThan(ref c, ref b)) + Subtract(ref c, ref modulus); + } + + public static void ModSub(out UInt128 c, ref UInt128 a, ref UInt128 b, ref UInt128 modulus) + { + Subtract(out c, ref a, ref b); + if (LessThan(ref a, ref b)) + Add(ref c, ref modulus); + } + + public static void ModMul(out UInt128 c, ref UInt128 a, ref UInt128 b, ref UInt128 modulus) + { + if (modulus.s1 == 0) + { + UInt128 product; + Multiply64(out product, a.s0, b.s0); + Create(out c, UInt128.Remainder(ref product, modulus.s0)); + } + else + { + UInt256 product; + Multiply(out product, ref a, ref b); + Remainder(out c, ref product, ref modulus); + } + } + + public static void ModMul(ref UInt128 a, ref UInt128 b, ref UInt128 modulus) + { + if (modulus.s1 == 0) + { + UInt128 product; + Multiply64(out product, a.s0, b.s0); + Create(out a, UInt128.Remainder(ref product, modulus.s0)); + } + else + { + UInt256 product; + Multiply(out product, ref a, ref b); + Remainder(out a, ref product, ref modulus); + } + } + + public static void ModPow(out UInt128 result, ref UInt128 value, ref UInt128 exponent, ref UInt128 modulus) + { + result = one; + var v = value; + var e = exponent.s0; + if (exponent.s1 != 0) + { + for (var i = 0; i < 64; i++) + { + if ((e & 1) != 0) + ModMul(ref result, ref v, ref modulus); + ModMul(ref v, ref v, ref modulus); + e >>= 1; + } + e = exponent.s1; + } + while (e != 0) + { + if ((e & 1) != 0) + ModMul(ref result, ref v, ref modulus); + if (e != 1) + ModMul(ref v, ref v, ref modulus); + e >>= 1; + } + Debug.Assert(BigInteger.ModPow(value, exponent, modulus) == result); + } + + public static void Shift(out UInt128 c, ref UInt128 a, int d) + { + if (d < 0) + RightShift(out c, ref a, -d); + else + LeftShift(out c, ref a, d); + } + + public static void ArithmeticShift(out UInt128 c, ref UInt128 a, int d) + { + if (d < 0) + ArithmeticRightShift(out c, ref a, -d); + else + LeftShift(out c, ref a, d); + } + + public static ulong LeftShift64(out UInt128 c, ref UInt128 a, int d) + { + if (d == 0) + { + c = a; + return 0; + } + var dneg = 64 - d; + c.s1 = a.s1 << d | a.s0 >> dneg; + c.s0 = a.s0 << d; + return a.s1 >> dneg; + } + + private static ulong LeftShift64(out UInt256 c, ref UInt256 a, int d) + { + if (d == 0) + { + c = a; + return 0; + } + var dneg = 64 - d; + c.s3 = a.s3 << d | a.s2 >> dneg; + c.s2 = a.s2 << d | a.s1 >> dneg; + c.s1 = a.s1 << d | a.s0 >> dneg; + c.s0 = a.s0 << d; + return a.s3 >> dneg; + } + + public static void LeftShift(out UInt128 c, ref UInt128 a, int b) + { + if (b < 64) + LeftShift64(out c, ref a, b); + else if (b == 64) + { + c.s0 = 0; + c.s1 = a.s0; + return; + } + else + { + c.s0 = 0; + c.s1 = a.s0 << (b - 64); + } + } + + public static void RightShift64(out UInt128 c, ref UInt128 a, int b) + { + if (b == 0) + c = a; + else + { + c.s0 = a.s0 >> b | a.s1 << (64 - b); + c.s1 = a.s1 >> b; + } + } + + public static void RightShift(out UInt128 c, ref UInt128 a, int b) + { + if (b < 64) + RightShift64(out c, ref a, b); + else if (b == 64) + { + c.s0 = a.s1; + c.s1 = 0; + } + else + { + c.s0 = a.s1 >> (b - 64); + c.s1 = 0; + } + } + + public static void ArithmeticRightShift64(out UInt128 c, ref UInt128 a, int b) + { + if (b == 0) + c = a; + else + { + c.s0 = a.s0 >> b | a.s1 << (64 - b); + c.s1 = (ulong)((long)a.s1 >> b); + } + } + + public static void ArithmeticRightShift(out UInt128 c, ref UInt128 a, int b) + { + if (b < 64) + ArithmeticRightShift64(out c, ref a, b); + else if (b == 64) + { + c.s0 = a.s1; + c.s1 = (ulong)((long)a.s1 >> 63); + } + else + { + c.s0 = a.s1 >> (b - 64); + c.s1 = (ulong)((long)a.s1 >> 63); + } + } + + public static void And(out UInt128 c, ref UInt128 a, ref UInt128 b) + { + c.s0 = a.s0 & b.s0; + c.s1 = a.s1 & b.s1; + } + + public static void Or(out UInt128 c, ref UInt128 a, ref UInt128 b) + { + c.s0 = a.s0 | b.s0; + c.s1 = a.s1 | b.s1; + } + + public static void ExclusiveOr(out UInt128 c, ref UInt128 a, ref UInt128 b) + { + c.s0 = a.s0 ^ b.s0; + c.s1 = a.s1 ^ b.s1; + } + + public static void Not(out UInt128 c, ref UInt128 a) + { + c.s0 = ~a.s0; + c.s1 = ~a.s1; + } + + public static void Negate(ref UInt128 a) + { + var s0 = a.s0; + a.s0 = 0 - s0; + a.s1 = 0 - a.s1; + if (s0 > 0) + --a.s1; + } + + + public static void Negate(out UInt128 c, ref UInt128 a) + { + c.s0 = 0 - a.s0; + c.s1 = 0 - a.s1; + if (a.s0 > 0) + --c.s1; + Debug.Assert((BigInteger)c == (BigInteger)(~a + 1)); + } + + public static void Pow(out UInt128 result, ref UInt128 value, uint exponent) + { + result = one; + while (exponent != 0) + { + + if ((exponent & 1) != 0) + { + var previous = result; + Multiply(out result, ref previous, ref value); + } + if (exponent != 1) + { + var previous = value; + Square(out value, ref previous); + } + exponent >>= 1; + } + } + + public static UInt128 Pow(UInt128 value, uint exponent) + { + UInt128 result; + Pow(out result, ref value, exponent); + return result; + } + + private const int maxRepShift = 53; + private static readonly ulong maxRep = (ulong)1 << maxRepShift; + private static readonly UInt128 maxRepSquaredHigh = (ulong)1 << (2 * maxRepShift - 64); + + public static ulong FloorSqrt(UInt128 a) + { + if (a.s1 == 0 && a.s0 <= maxRep) + return (ulong)Math.Sqrt(a.s0); + var s = (ulong)Math.Sqrt(ConvertToDouble(ref a)); + if (a.s1 < maxRepSquaredHigh) + { + UInt128 s2; + Square(out s2, s); + var r = a.s0 - s2.s0; + if (r > long.MaxValue) + --s; + else if (r - (s << 1) <= long.MaxValue) + ++s; + Debug.Assert((BigInteger)s * s <= a && (BigInteger)(s + 1) * (s + 1) > a); + return s; + } + s = FloorSqrt(ref a, s); + Debug.Assert((BigInteger)s * s <= a && (BigInteger)(s + 1) * (s + 1) > a); + return s; + } + + public static ulong CeilingSqrt(UInt128 a) + { + if (a.s1 == 0 && a.s0 <= maxRep) + return (ulong)Math.Ceiling(Math.Sqrt(a.s0)); + var s = (ulong)Math.Ceiling(Math.Sqrt(ConvertToDouble(ref a))); + if (a.s1 < maxRepSquaredHigh) + { + UInt128 s2; + Square(out s2, s); + var r = s2.s0 - a.s0; + if (r > long.MaxValue) + ++s; + else if (r - (s << 1) <= long.MaxValue) + --s; + Debug.Assert((BigInteger)(s - 1) * (s - 1) < a && (BigInteger)s * s >= a); + return s; + } + s = FloorSqrt(ref a, s); + UInt128 square; + Square(out square, s); + if (square.S0 != a.S0 || square.S1 != a.S1) + ++s; + Debug.Assert((BigInteger)(s - 1) * (s - 1) < a && (BigInteger)s * s >= a); + return s; + } + + private static ulong FloorSqrt(ref UInt128 a, ulong s) + { + var sprev = (ulong)0; + UInt128 div; + UInt128 sum; + while (true) + { + // Equivalent to: + // snext = (a / s + s) / 2; + Divide(out div, ref a, s); + Add(out sum, ref div, s); + var snext = sum.S0 >> 1; + if (sum.S1 != 0) + snext |= (ulong)1 << 63; + if (snext == sprev) + { + if (snext < s) + s = snext; + break; + } + sprev = s; + s = snext; + } + return s; + } + + public static ulong FloorCbrt(UInt128 a) + { + var s = (ulong)Math.Pow(ConvertToDouble(ref a), (double)1 / 3); + UInt128 s3; + Cube(out s3, s); + if (a < s3) + --s; + else + { + UInt128 sum; + Multiply(out sum, 3 * s, s + 1); + UInt128 diff; + Subtract(out diff, ref a, ref s3); + if (LessThan(ref sum, ref diff)) + ++s; + } + Debug.Assert((BigInteger)s * s * s <= a && (BigInteger)(s + 1) * (s + 1) * (s + 1) > a); + return s; + } + + public static ulong CeilingCbrt(UInt128 a) + { + var s = (ulong)Math.Ceiling(Math.Pow(ConvertToDouble(ref a), (double)1 / 3)); + UInt128 s3; + Cube(out s3, s); + if (s3 < a) + ++s; + else + { + UInt128 sum; + Multiply(out sum, 3 * s, s + 1); + UInt128 diff; + Subtract(out diff, ref s3, ref a); + if (LessThan(ref sum, ref diff)) + --s; + } + Debug.Assert((BigInteger)(s - 1) * (s - 1) * (s - 1) < a && (BigInteger)s * s * s >= a); + return s; + } + + public static UInt128 Min(UInt128 a, UInt128 b) + { + if (LessThan(ref a, ref b)) + return a; + return b; + } + + public static UInt128 Max(UInt128 a, UInt128 b) + { + if (LessThan(ref b, ref a)) + return a; + return b; + } + + public static double Log(UInt128 a) + { + return Log(a, Math.E); + } + + public static double Log10(UInt128 a) + { + return Log(a, 10); + } + + public static double Log(UInt128 a, double b) + { + return Math.Log(ConvertToDouble(ref a), b); + } + + public static UInt128 Add(UInt128 a, UInt128 b) + { + UInt128 c; + Add(out c, ref a, ref b); + return c; + } + + public static UInt128 Subtract(UInt128 a, UInt128 b) + { + UInt128 c; + Subtract(out c, ref a, ref b); + return c; + } + + public static UInt128 Multiply(UInt128 a, UInt128 b) + { + UInt128 c; + Multiply(out c, ref a, ref b); + return c; + } + + public static UInt128 Divide(UInt128 a, UInt128 b) + { + UInt128 c; + Divide(out c, ref a, ref b); + return c; + } + + public static UInt128 Remainder(UInt128 a, UInt128 b) + { + UInt128 c; + Remainder(out c, ref a, ref b); + return c; + } + + public static UInt128 DivRem(UInt128 a, UInt128 b, out UInt128 remainder) + { + UInt128 c; + Divide(out c, ref a, ref b); + Remainder(out remainder, ref a, ref b); + return c; + } + + public static UInt128 ModAdd(UInt128 a, UInt128 b, UInt128 modulus) + { + UInt128 c; + ModAdd(out c, ref a, ref b, ref modulus); + return c; + } + + public static UInt128 ModSub(UInt128 a, UInt128 b, UInt128 modulus) + { + UInt128 c; + ModSub(out c, ref a, ref b, ref modulus); + return c; + } + + public static UInt128 ModMul(UInt128 a, UInt128 b, UInt128 modulus) + { + UInt128 c; + ModMul(out c, ref a, ref b, ref modulus); + return c; + } + + public static UInt128 ModPow(UInt128 value, UInt128 exponent, UInt128 modulus) + { + UInt128 result; + ModPow(out result, ref value, ref exponent, ref modulus); + return result; + } + + public static UInt128 Negate(UInt128 a) + { + UInt128 c; + Negate(out c, ref a); + return c; + } + + public static UInt128 GreatestCommonDivisor(UInt128 a, UInt128 b) + { + UInt128 c; + GreatestCommonDivisor(out c, ref a, ref b); + return c; + } + + private static void RightShift64(ref UInt128 c, int d) + { + if (d == 0) + return; + c.s0 = c.s1 << (64 - d) | c.s0 >> d; + c.s1 >>= d; + } + + public static void RightShift(ref UInt128 c, int d) + { + if (d < 64) + RightShift64(ref c, d); + else + { + c.s0 = c.s1 >> (d - 64); + c.s1 = 0; + } + } + + public static void Shift(ref UInt128 c, int d) + { + if (d < 0) + RightShift(ref c, -d); + else + LeftShift(ref c, d); + } + + public static void ArithmeticShift(ref UInt128 c, int d) + { + if (d < 0) + ArithmeticRightShift(ref c, -d); + else + LeftShift(ref c, d); + } + + public static void RightShift(ref UInt128 c) + { + c.s0 = c.s1 << 63 | c.s0 >> 1; + c.s1 >>= 1; + } + + private static void ArithmeticRightShift64(ref UInt128 c, int d) + { + if (d == 0) + return; + c.s0 = c.s1 << (64 - d) | c.s0 >> d; + c.s1 = (ulong)((long)c.s1 >> d); + } + + public static void ArithmeticRightShift(ref UInt128 c, int d) + { + if (d < 64) + ArithmeticRightShift64(ref c, d); + else + { + c.s0 = (ulong)((long)c.s1 >> (d - 64)); + c.s1 = 0; + } + } + + public static void ArithmeticRightShift(ref UInt128 c) + { + c.s0 = c.s1 << 63 | c.s0 >> 1; + c.s1 = (ulong)((long)c.s1 >> 1); + } + + private static ulong LeftShift64(ref UInt128 c, int d) + { + if (d == 0) + return 0; + var dneg = 64 - d; + var result = c.s1 >> dneg; + c.s1 = c.s1 << d | c.s0 >> dneg; + c.s0 <<= d; + return result; + } + + public static void LeftShift(ref UInt128 c, int d) + { + if (d < 64) + LeftShift64(ref c, d); + else + { + c.s1 = c.s0 << (d - 64); + c.s0 = 0; + } + } + + public static void LeftShift(ref UInt128 c) + { + c.s1 = c.s1 << 1 | c.s0 >> 63; + c.s0 <<= 1; + } + + public static void Swap(ref UInt128 a, ref UInt128 b) + { + var as0 = a.s0; + var as1 = a.s1; + a.s0 = b.s0; + a.s1 = b.s1; + b.s0 = as0; + b.s1 = as1; + } + + public static void GreatestCommonDivisor(out UInt128 c, ref UInt128 a, ref UInt128 b) + { + // Check whether one number is > 64 bits and the other is <= 64 bits and both are non-zero. + UInt128 a1, b1; + if ((a.s1 == 0) != (b.s1 == 0) && !a.IsZero && !b.IsZero) + { + // Perform a normal step so that both a and b are <= 64 bits. + if (LessThan(ref a, ref b)) + { + a1 = a; + Remainder(out b1, ref b, ref a); + } + else + { + b1 = b; + Remainder(out a1, ref a, ref b); + } + } + else + { + a1 = a; + b1 = b; + } + + // Make sure neither is zero. + if (a1.IsZero) + { + c = b1; + return; + } + if (b1.IsZero) + { + c = a1; + return; + } + + // Ensure a >= b. + if (LessThan(ref a1, ref b1)) + Swap(ref a1, ref b1); + + // Lehmer-Euclid algorithm. + // See: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.31.693 + while (a1.s1 != 0 && !b.IsZero) + { + // Extract the high 63 bits of a and b. + var norm = 63 - GetBitLength(a1.s1); + UInt128 ahat, bhat; + Shift(out ahat, ref a1, norm); + Shift(out bhat, ref b1, norm); + var uhat = (long)ahat.s1; + var vhat = (long)bhat.s1; + + // Check whether q exceeds single-precision. + if (vhat == 0) + { + // Perform a normal step and try again. + UInt128 rem; + Remainder(out rem, ref a1, ref b1); + a1 = b1; + b1 = rem; + continue; + } + + // Perform steps using signed single-precision arithmetic. + var x0 = (long)1; + var y0 = (long)0; + var x1 = (long)0; + var y1 = (long)1; + var even = true; + while (true) + { + // Calculate quotient, cosquence pair, and update uhat and vhat. + var q = uhat / vhat; + var x2 = x0 - q * x1; + var y2 = y0 - q * y1; + var t = uhat; + uhat = vhat; + vhat = t - q * vhat; + even = !even; + + // Apply Jebelean's termination condition + // to check whether q is valid. + if (even) + { + if (vhat < -x2 || uhat - vhat < y2 - y1) + break; + } + else + { + if (vhat < -y2 || uhat - vhat < x2 - x1) + break; + } + + // Adjust cosequence history. + x0 = x1; y0 = y1; x1 = x2; y1 = y2; + } + + // Check whether a normal step is necessary. + if (x0 == 1 && y0 == 0) + { + UInt128 rem; + Remainder(out rem, ref a1, ref b1); + a1 = b1; + b1 = rem; + continue; + } + + // Back calculate a and b from the last valid cosequence pairs. + UInt128 anew, bnew; + if (even) + { + AddProducts(out anew, y0, ref b1, x0, ref a1); + AddProducts(out bnew, x1, ref a1, y1, ref b1); + } + else + { + AddProducts(out anew, x0, ref a1, y0, ref b1); + AddProducts(out bnew, y1, ref b1, x1, ref a1); + } + a1 = anew; + b1 = bnew; + } + + // Check whether we have any 64 bit work left. + if (!b1.IsZero) + { + var a2 = a1.s0; + var b2 = b1.s0; + + // Perform 64 bit steps. + while (a2 > uint.MaxValue && b2 != 0) + { + var t = a2 % b2; + a2 = b2; + b2 = t; + } + + // Check whether we have any 32 bit work left. + if (b2 != 0) + { + var a3 = (uint)a2; + var b3 = (uint)b2; + + // Perform 32 bit steps. + while (b3 != 0) + { + var t = a3 % b3; + a3 = b3; + b3 = t; + } + + Create(out c, a3); + } + else + Create(out c, a2); + } + else + c = a1; + } + + private static void AddProducts(out UInt128 result, long x, ref UInt128 u, long y, ref UInt128 v) + { + // Compute x * u + y * v assuming y is negative and the result is positive and fits in 128 bits. + UInt128 product1; + Multiply(out product1, ref u, (ulong)x); + UInt128 product2; + Multiply(out product2, ref v, (ulong)(-y)); + Subtract(out result, ref product1, ref product2); + } + + public static int Compare(UInt128 a, UInt128 b) + { + return a.CompareTo(b); + } + + private static byte[] bitLength = Enumerable.Range(0, byte.MaxValue + 1) + .Select(value => + { + int count; + for (count = 0; value != 0; count++) + value >>= 1; + return (byte)count; + }).ToArray(); + + private static int GetBitLength(uint value) + { + var tt = value >> 16; + if (tt != 0) + { + var t = tt >> 8; + if (t != 0) + return bitLength[t] + 24; + return bitLength[tt] + 16; + } + else + { + var t = value >> 8; + if (t != 0) + return bitLength[t] + 8; + return bitLength[value]; + } + } + + private static int GetBitLength(ulong value) + { + var r1 = value >> 32; + if (r1 != 0) + return GetBitLength((uint)r1) + 32; + return GetBitLength((uint)value); + } + + public static void Reduce(out UInt128 w, ref UInt128 u, ref UInt128 v, ref UInt128 n, ulong k0) + { + UInt128 carry; + Multiply64(out carry, u.s0, v.s0); + var t0 = carry.s0; + Multiply64(out carry, u.s1, v.s0, carry.s1); + var t1 = carry.s0; + var t2 = carry.s1; + + var m = t0 * k0; + Multiply64(out carry, m, n.s1, MultiplyHigh64(m, n.s0, t0)); + Add(ref carry, t1); + t0 = carry.s0; + Add(out carry, carry.s1, t2); + t1 = carry.s0; + t2 = carry.s1; + + Multiply64(out carry, u.s0, v.s1, t0); + t0 = carry.s0; + Multiply64(out carry, u.s1, v.s1, carry.s1); + Add(ref carry, t1); + t1 = carry.s0; + Add(out carry, carry.s1, t2); + t2 = carry.s0; + var t3 = carry.s1; + + m = t0 * k0; + Multiply64(out carry, m, n.s1, MultiplyHigh64(m, n.s0, t0)); + Add(ref carry, t1); + t0 = carry.s0; + Add(out carry, carry.s1, t2); + t1 = carry.s0; + t2 = t3 + carry.s1; + + Create(out w, t0, t1); + if (t2 != 0 || !LessThan(ref w, ref n)) + Subtract(ref w, ref n); + } + + public static void Reduce(out UInt128 w, ref UInt128 t, ref UInt128 n, ulong k0) + { + UInt128 carry; + var t0 = t.s0; + var t1 = t.s1; + var t2 = (ulong)0; + + for (var i = 0; i < 2; i++) + { + var m = t0 * k0; + Multiply64(out carry, m, n.s1, MultiplyHigh64(m, n.s0, t0)); + Add(ref carry, t1); + t0 = carry.s0; + Add(out carry, carry.s1, t2); + t1 = carry.s0; + t2 = carry.s1; + } + + Create(out w, t0, t1); + if (t2 != 0 || !LessThan(ref w, ref n)) + Subtract(ref w, ref n); + } + + public static UInt128 Reduce(UInt128 u, UInt128 v, UInt128 n, ulong k0) + { + UInt128 w; + Reduce(out w, ref u, ref v, ref n, k0); + return w; + } + + public static UInt128 Reduce(UInt128 t, UInt128 n, ulong k0) + { + UInt128 w; + Reduce(out w, ref t, ref n, k0); + return w; + } + } +} diff --git a/Dirichlet.Numerics/license.txt b/Dirichlet.Numerics/license.txt new file mode 100644 index 000000000..5e6188c6d --- /dev/null +++ b/Dirichlet.Numerics/license.txt @@ -0,0 +1,13 @@ +Copyright 2014 Rick Sladkey + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. \ No newline at end of file diff --git a/TestBenchFramework.sln b/TestBenchFramework.sln index 3ccd58b5e..c41e3316d 100644 --- a/TestBenchFramework.sln +++ b/TestBenchFramework.sln @@ -17,6 +17,8 @@ Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "ResultsBrowser", "ResultsBr EndProject Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "DeviceTest", "DeviceTest\DeviceTest.csproj", "{6D3384DC-4638-4A92-91A8-F39D900377C6}" EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Dirichlet.Numerics", "Dirichlet.Numerics\Dirichlet.Numerics.csproj", "{439D0878-C76E-452B-B17D-209A89E91D36}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|Any CPU = Debug|Any CPU @@ -83,6 +85,16 @@ Global {6D3384DC-4638-4A92-91A8-F39D900377C6}.Release|Mixed Platforms.Build.0 = Release|x86 {6D3384DC-4638-4A92-91A8-F39D900377C6}.Release|x86.ActiveCfg = Release|x86 {6D3384DC-4638-4A92-91A8-F39D900377C6}.Release|x86.Build.0 = Release|x86 + {439D0878-C76E-452B-B17D-209A89E91D36}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {439D0878-C76E-452B-B17D-209A89E91D36}.Debug|Any CPU.Build.0 = Debug|Any CPU + {439D0878-C76E-452B-B17D-209A89E91D36}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU + {439D0878-C76E-452B-B17D-209A89E91D36}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU + {439D0878-C76E-452B-B17D-209A89E91D36}.Debug|x86.ActiveCfg = Debug|Any CPU + {439D0878-C76E-452B-B17D-209A89E91D36}.Release|Any CPU.ActiveCfg = Release|Any CPU + {439D0878-C76E-452B-B17D-209A89E91D36}.Release|Any CPU.Build.0 = Release|Any CPU + {439D0878-C76E-452B-B17D-209A89E91D36}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU + {439D0878-C76E-452B-B17D-209A89E91D36}.Release|Mixed Platforms.Build.0 = Release|Any CPU + {439D0878-C76E-452B-B17D-209A89E91D36}.Release|x86.ActiveCfg = Release|Any CPU EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE