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