--- /dev/null
+// http://www.fractal-landscapes.co.uk/bigint.html\r
+\r
+using System;\r
+\r
+namespace BigNum\r
+{\r
+ /// <summary>\r
+ /// An arbitrary-precision floating-point class\r
+ /// \r
+ /// Format:\r
+ /// Each number is stored as an exponent (32-bit signed integer), and a mantissa\r
+ /// (n-bit) BigInteger. The sign of the number is stored in the BigInteger\r
+ /// \r
+ /// Applicability and Performance:\r
+ /// This class is designed to be used for small extended precisions. It may not be\r
+ /// safe (and certainly won't be fast) to use it with mixed-precision arguments.\r
+ /// It does support, but will not be efficient for, numbers over around 2048 bits.\r
+ /// \r
+ /// Notes:\r
+ /// All conversions to and from strings are slow.\r
+ /// \r
+ /// Conversions from simple integer types Int32, Int64, UInt32, UInt64 are performed\r
+ /// using the appropriate constructor, and are relatively fast.\r
+ /// \r
+ /// The class is written entirely in managed C# code, with not native or managed\r
+ /// assembler. The use of native assembler would speed up the multiplication operations\r
+ /// many times over, and therefore all higher-order operations too.\r
+ /// </summary>\r
+ public class BigFloat\r
+ {\r
+ /// <summary>\r
+ /// Floats can have 4 special value types:\r
+ /// \r
+ /// NaN: Not a number (cannot be changed using any operations)\r
+ /// Infinity: Positive infinity. Some operations e.g. Arctan() allow this input.\r
+ /// -Infinity: Negative infinity. Some operations allow this input.\r
+ /// Zero\r
+ /// </summary>\r
+ public enum SpecialValueType\r
+ {\r
+ /// <summary>\r
+ /// Not a special value\r
+ /// </summary>\r
+ NONE = 0,\r
+ /// <summary>\r
+ /// Zero\r
+ /// </summary>\r
+ ZERO,\r
+ /// <summary>\r
+ /// Positive infinity\r
+ /// </summary>\r
+ INF_PLUS,\r
+ /// <summary>\r
+ /// Negative infinity\r
+ /// </summary>\r
+ INF_MINUS,\r
+ /// <summary>\r
+ /// Not a number\r
+ /// </summary>\r
+ NAN\r
+ }\r
+\r
+ /// <summary>\r
+ /// This affects the ToString() method. \r
+ /// \r
+ /// With Trim rounding, all insignificant zero digits are drip\r
+ /// </summary>\r
+ public enum RoundingModeType\r
+ {\r
+ /// <summary>\r
+ /// Trim non-significant zeros from ToString output after rounding\r
+ /// </summary>\r
+ TRIM,\r
+ /// <summary>\r
+ /// Keep all non-significant zeroes in ToString output after rounding\r
+ /// </summary>\r
+ EXACT\r
+ }\r
+\r
+ /// <summary>\r
+ /// A wrapper for the signed exponent, avoiding overflow.\r
+ /// </summary>\r
+ protected struct ExponentAdaptor\r
+ {\r
+ /// <summary>\r
+ /// The 32-bit exponent\r
+ /// </summary>\r
+ public Int32 exponent\r
+ {\r
+ get { return expValue; }\r
+ set { expValue = value; }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Implicit cast to Int32\r
+ /// </summary>\r
+ public static implicit operator Int32(ExponentAdaptor adaptor)\r
+ {\r
+ return adaptor.expValue;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Implicit cast from Int32 to ExponentAdaptor\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <returns></returns>\r
+ public static implicit operator ExponentAdaptor(Int32 value)\r
+ {\r
+ ExponentAdaptor adaptor = new ExponentAdaptor();\r
+ adaptor.expValue = value;\r
+ return adaptor;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Overloaded increment operator\r
+ /// </summary>\r
+ public static ExponentAdaptor operator ++(ExponentAdaptor adaptor)\r
+ {\r
+ adaptor = adaptor + 1;\r
+ return adaptor;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Overloaded decrement operator\r
+ /// </summary>\r
+ public static ExponentAdaptor operator --(ExponentAdaptor adaptor)\r
+ {\r
+ adaptor = adaptor - 1;\r
+ return adaptor;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Overloaded addition operator\r
+ /// </summary>\r
+ public static ExponentAdaptor operator +(ExponentAdaptor a1, ExponentAdaptor a2)\r
+ {\r
+ if (a1.expValue == Int32.MaxValue) return a1;\r
+\r
+ Int64 temp = (Int64)a1.expValue;\r
+ temp += (Int64)(a2.expValue);\r
+\r
+ if (temp > (Int64)Int32.MaxValue)\r
+ {\r
+ a1.expValue = Int32.MaxValue;\r
+ }\r
+ else if (temp < (Int64)Int32.MinValue)\r
+ {\r
+ a1.expValue = Int32.MinValue;\r
+ }\r
+ else\r
+ {\r
+ a1.expValue = (Int32)temp;\r
+ }\r
+\r
+ return a1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Overloaded subtraction operator\r
+ /// </summary>\r
+ public static ExponentAdaptor operator -(ExponentAdaptor a1, ExponentAdaptor a2)\r
+ {\r
+ if (a1.expValue == Int32.MaxValue) return a1;\r
+\r
+ Int64 temp = (Int64)a1.expValue;\r
+ temp -= (Int64)(a2.expValue);\r
+\r
+ if (temp > (Int64)Int32.MaxValue)\r
+ {\r
+ a1.expValue = Int32.MaxValue;\r
+ }\r
+ else if (temp < (Int64)Int32.MinValue)\r
+ {\r
+ a1.expValue = Int32.MinValue;\r
+ }\r
+ else\r
+ {\r
+ a1.expValue = (Int32)temp;\r
+ }\r
+\r
+ return a1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Overloaded multiplication operator\r
+ /// </summary>\r
+ public static ExponentAdaptor operator *(ExponentAdaptor a1, ExponentAdaptor a2)\r
+ {\r
+ if (a1.expValue == Int32.MaxValue) return a1;\r
+\r
+ Int64 temp = (Int64)a1.expValue;\r
+ temp *= (Int64)a2.expValue;\r
+\r
+ if (temp > (Int64)Int32.MaxValue)\r
+ {\r
+ a1.expValue = Int32.MaxValue;\r
+ }\r
+ else if (temp < (Int64)Int32.MinValue)\r
+ {\r
+ a1.expValue = Int32.MinValue;\r
+ }\r
+ else\r
+ {\r
+ a1.expValue = (Int32)temp;\r
+ }\r
+\r
+ return a1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Overloaded division operator\r
+ /// </summary>\r
+ public static ExponentAdaptor operator /(ExponentAdaptor a1, ExponentAdaptor a2)\r
+ {\r
+ if (a1.expValue == Int32.MaxValue) return a1;\r
+\r
+ ExponentAdaptor res = new ExponentAdaptor();\r
+ res.expValue = a1.expValue / a2.expValue;\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Overloaded right-shift operator\r
+ /// </summary>\r
+ public static ExponentAdaptor operator >>(ExponentAdaptor a1, int shift)\r
+ {\r
+ if (a1.expValue == Int32.MaxValue) return a1;\r
+\r
+ ExponentAdaptor res = new ExponentAdaptor();\r
+ res.expValue = a1.expValue >> shift;\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Overloaded left-shift operator\r
+ /// </summary>\r
+ /// <param name="a1"></param>\r
+ /// <param name="shift"></param>\r
+ /// <returns></returns>\r
+ public static ExponentAdaptor operator <<(ExponentAdaptor a1, int shift)\r
+ {\r
+ if (a1.expValue == 0) return a1;\r
+\r
+ ExponentAdaptor res = new ExponentAdaptor();\r
+ res.expValue = a1.expValue;\r
+\r
+ if (shift > 31)\r
+ {\r
+ res.expValue = Int32.MaxValue;\r
+ }\r
+ else\r
+ {\r
+ Int64 temp = a1.expValue;\r
+ temp = temp << shift;\r
+\r
+ if (temp > (Int64)Int32.MaxValue)\r
+ {\r
+ res.expValue = Int32.MaxValue;\r
+ }\r
+ else if (temp < (Int64)Int32.MinValue)\r
+ {\r
+ res.expValue = Int32.MinValue;\r
+ }\r
+ else\r
+ {\r
+ res.expValue = (Int32)temp;\r
+ }\r
+ }\r
+\r
+ return res;\r
+ }\r
+\r
+ private Int32 expValue;\r
+ }\r
+\r
+ //************************ Constructors **************************\r
+\r
+ /// <summary>\r
+ /// Constructs a 128-bit BigFloat\r
+ /// \r
+ /// Sets the value to zero\r
+ /// </summary>\r
+ static BigFloat()\r
+ {\r
+ RoundingDigits = 3;\r
+ RoundingMode = RoundingModeType.TRIM;\r
+ scratch = new BigInt(new PrecisionSpec(128, PrecisionSpec.BaseType.BIN));\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigFloat of the required precision\r
+ /// \r
+ /// Sets the value to zero\r
+ /// </summary>\r
+ /// <param name="mantissaPrec"></param>\r
+ public BigFloat(PrecisionSpec mantissaPrec)\r
+ {\r
+ Init(mantissaPrec);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a big float from a UInt32 to the required precision\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <param name="mantissaPrec"></param>\r
+ public BigFloat(UInt32 value, PrecisionSpec mantissaPrec)\r
+ {\r
+ int mbWords = ((mantissaPrec.NumBits) >> 5);\r
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
+ int newManBits = mbWords << 5;\r
+\r
+ //For efficiency, we just use a 32-bit exponent\r
+ exponent = 0;\r
+\r
+ mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+ //scratch = new BigInt(mantissa.Precision);\r
+\r
+ int bit = BigInt.GetMSB(value);\r
+ if (bit == -1) return;\r
+\r
+ int shift = mantissa.Precision.NumBits - (bit + 1);\r
+ mantissa.LSH(shift);\r
+ exponent = bit;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigFloat from an Int32 to the required precision\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <param name="mantissaPrec"></param>\r
+ public BigFloat(Int32 value, PrecisionSpec mantissaPrec)\r
+ {\r
+ int mbWords = ((mantissaPrec.NumBits) >> 5);\r
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
+ int newManBits = mbWords << 5;\r
+\r
+ //For efficiency, we just use a 32-bit exponent\r
+ exponent = 0;\r
+ UInt32 uValue;\r
+ \r
+ if (value < 0)\r
+ {\r
+ if (value == Int32.MinValue)\r
+ {\r
+ uValue = 0x80000000;\r
+ }\r
+ else\r
+ {\r
+ uValue = (UInt32)(-value);\r
+ }\r
+ }\r
+ else\r
+ {\r
+ uValue = (UInt32)value;\r
+ }\r
+\r
+ mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+ //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+\r
+ int bit = BigInt.GetMSB(uValue);\r
+ if (bit == -1) return;\r
+\r
+ int shift = mantissa.Precision.NumBits - (bit + 1);\r
+ mantissa.LSH(shift);\r
+ exponent = bit;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigFloat from a 64-bit integer\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <param name="mantissaPrec"></param>\r
+ public BigFloat(Int64 value, PrecisionSpec mantissaPrec)\r
+ {\r
+ int mbWords = ((mantissaPrec.NumBits) >> 5);\r
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
+ int newManBits = mbWords << 5;\r
+\r
+ //For efficiency, we just use a 32-bit exponent\r
+ exponent = 0;\r
+ UInt64 uValue;\r
+\r
+ if (value < 0)\r
+ {\r
+ if (value == Int64.MinValue)\r
+ {\r
+ uValue = 0x80000000;\r
+ }\r
+ else\r
+ {\r
+ uValue = (UInt64)(-value);\r
+ }\r
+ }\r
+ else\r
+ {\r
+ uValue = (UInt64)value;\r
+ }\r
+\r
+ mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+ //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+\r
+ int bit = BigInt.GetMSB(uValue);\r
+ if (bit == -1) return;\r
+\r
+ int shift = mantissa.Precision.NumBits - (bit + 1);\r
+ if (shift > 0)\r
+ {\r
+ mantissa.LSH(shift);\r
+ }\r
+ else\r
+ {\r
+ mantissa.SetHighDigit((uint)(uValue >> (-shift)));\r
+ }\r
+ exponent = bit;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigFloat from a 64-bit unsigned integer\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <param name="mantissaPrec"></param>\r
+ public BigFloat(UInt64 value, PrecisionSpec mantissaPrec)\r
+ {\r
+ int mbWords = ((mantissaPrec.NumBits) >> 5);\r
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
+ int newManBits = mbWords << 5;\r
+\r
+ //For efficiency, we just use a 32-bit exponent\r
+ exponent = 0;\r
+\r
+ int bit = BigInt.GetMSB(value);\r
+\r
+ mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+ //scratch = new BigInt(mantissa.Precision);\r
+\r
+ int shift = mantissa.Precision.NumBits - (bit + 1);\r
+ if (shift > 0)\r
+ {\r
+ mantissa.LSH(shift);\r
+ }\r
+ else\r
+ {\r
+ mantissa.SetHighDigit((uint)(value >> (-shift)));\r
+ }\r
+ exponent = bit;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigFloat from a BigInt, using the specified precision\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <param name="mantissaPrec"></param>\r
+ public BigFloat(BigInt value, PrecisionSpec mantissaPrec)\r
+ {\r
+ if (value.IsZero())\r
+ {\r
+ Init(mantissaPrec);\r
+ SetZero();\r
+ return;\r
+ }\r
+\r
+ mantissa = new BigInt(value, mantissaPrec);\r
+ exponent = BigInt.GetMSB(value);\r
+ mantissa.Normalise();\r
+ }\r
+\r
+ /// <summary>\r
+ /// Construct a BigFloat from a double-precision floating point number\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <param name="mantissaPrec"></param>\r
+ public BigFloat(double value, PrecisionSpec mantissaPrec)\r
+ {\r
+ if (value == 0.0)\r
+ {\r
+ Init(mantissaPrec);\r
+ return;\r
+ }\r
+\r
+ bool sign = (value < 0) ? true : false;\r
+\r
+ long bits = BitConverter.DoubleToInt64Bits(value);\r
+ // Note that the shift is sign-extended, hence the test against -1 not 1\r
+ int valueExponent = (int)((bits >> 52) & 0x7ffL);\r
+ long valueMantissa = bits & 0xfffffffffffffL;\r
+\r
+ //The mantissa is stored with the top bit implied.\r
+ valueMantissa = valueMantissa | 0x10000000000000L;\r
+\r
+ //The exponent is biased by 1023.\r
+ exponent = valueExponent - 1023;\r
+\r
+ //Round the number of bits to the nearest word.\r
+ int mbWords = ((mantissaPrec.NumBits) >> 5);\r
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
+ int newManBits = mbWords << 5;\r
+\r
+ mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+ //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+\r
+ if (newManBits >= 64)\r
+ {\r
+ //The mantissa is 53 bits now, so add 11 to put it in the right place.\r
+ mantissa.SetHighDigits(valueMantissa << 11);\r
+ }\r
+ else\r
+ {\r
+ //To get the top word of the mantissa, shift up by 11 and down by 32 = down by 21\r
+ mantissa.SetHighDigit((uint)(valueMantissa >> 21));\r
+ }\r
+\r
+ mantissa.Sign = sign;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Copy constructor\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ public BigFloat(BigFloat value)\r
+ {\r
+ Init(value.mantissa.Precision);\r
+ exponent = value.exponent;\r
+ mantissa.Assign(value.mantissa);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Copy Constructor - constructs a new BigFloat with the specified precision, copying the old one.\r
+ /// \r
+ /// The value is rounded towards zero in the case where precision is decreased. The Round() function\r
+ /// should be used beforehand if a correctly rounded result is required.\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <param name="mantissaPrec"></param>\r
+ public BigFloat(BigFloat value, PrecisionSpec mantissaPrec)\r
+ {\r
+ Init(mantissaPrec);\r
+ exponent = value.exponent;\r
+ if (mantissa.AssignHigh(value.mantissa)) exponent++;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigFloat from a string\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <param name="mantissaPrec"></param>\r
+ public BigFloat(string value, PrecisionSpec mantissaPrec)\r
+ {\r
+ Init(mantissaPrec);\r
+\r
+ PrecisionSpec extendedPres = new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN);\r
+ BigFloat ten = new BigFloat(10, extendedPres);\r
+ BigFloat iPart = new BigFloat(extendedPres);\r
+ BigFloat fPart = new BigFloat(extendedPres);\r
+ BigFloat tenRCP = ten.Reciprocal();\r
+\r
+ if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol))\r
+ {\r
+ SetNaN();\r
+ return;\r
+ }\r
+ else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol))\r
+ {\r
+ SetInfPlus();\r
+ return;\r
+ }\r
+ else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol))\r
+ {\r
+ SetInfMinus();\r
+ return;\r
+ }\r
+\r
+ string decimalpoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;\r
+\r
+ char[] digitChars = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ',', '.' };\r
+\r
+ //Read in the integer part up the the decimal point.\r
+ bool sign = false;\r
+ value = value.Trim();\r
+\r
+ int i = 0;\r
+\r
+ if (value.Length > i && value[i] == '-')\r
+ {\r
+ sign = true;\r
+ i++;\r
+ }\r
+\r
+ if (value.Length > i && value[i] == '+')\r
+ {\r
+ i++;\r
+ }\r
+\r
+ for ( ; i < value.Length; i++)\r
+ {\r
+ //break on decimal point\r
+ if (value[i] == decimalpoint[0]) break;\r
+\r
+ int digit = Array.IndexOf(digitChars, value[i]);\r
+ if (digit < 0) break;\r
+\r
+ //Ignore place separators (assumed either , or .)\r
+ if (digit > 9) continue;\r
+\r
+ if (i > 0) iPart.Mul(ten);\r
+ iPart.Add(new BigFloat(digit, extendedPres));\r
+ }\r
+\r
+ //If we've run out of characters, assign everything and return\r
+ if (i == value.Length)\r
+ {\r
+ iPart.mantissa.Sign = sign;\r
+ exponent = iPart.exponent;\r
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
+ return;\r
+ }\r
+\r
+ //Assign the characters after the decimal point to fPart\r
+ if (value[i] == '.' && i < value.Length - 1)\r
+ {\r
+ BigFloat RecipToUse = new BigFloat(tenRCP);\r
+\r
+ for (i++; i < value.Length; i++)\r
+ {\r
+ int digit = Array.IndexOf(digitChars, value[i]);\r
+ if (digit < 0) break;\r
+ BigFloat temp = new BigFloat(digit, extendedPres);\r
+ temp.Mul(RecipToUse);\r
+ RecipToUse.Mul(tenRCP);\r
+ fPart.Add(temp);\r
+ }\r
+ }\r
+\r
+ //If we're run out of characters, add fPart and iPart and return\r
+ if (i == value.Length)\r
+ {\r
+ iPart.Add(fPart);\r
+ iPart.mantissa.Sign = sign;\r
+ exponent = iPart.exponent;\r
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
+ return;\r
+ }\r
+\r
+ if (value[i] == '+' || value[i] == '-') i++;\r
+\r
+ if (i == value.Length)\r
+ {\r
+ iPart.Add(fPart);\r
+ iPart.mantissa.Sign = sign;\r
+ exponent = iPart.exponent;\r
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
+ return;\r
+ }\r
+\r
+ //Look for exponential notation.\r
+ if ((value[i] == 'e' || value[i] == 'E') && i < value.Length - 1)\r
+ {\r
+ //Convert the exponent to an int.\r
+ int exp;\r
+\r
+ try\r
+ {\r
+ exp = System.Convert.ToInt32(new string(value.ToCharArray()));// i + 1, value.Length - (i + 1))));\r
+ }\r
+ catch (Exception)\r
+ {\r
+ iPart.Add(fPart);\r
+ iPart.mantissa.Sign = sign;\r
+ exponent = iPart.exponent;\r
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
+ return;\r
+ }\r
+\r
+ //Raise or lower 10 to the power of the exponent\r
+ BigFloat acc = new BigFloat(1, extendedPres);\r
+ BigFloat temp = new BigFloat(1, extendedPres);\r
+\r
+ int powerTemp = exp;\r
+\r
+ BigFloat multiplierToUse;\r
+\r
+ if (exp < 0)\r
+ {\r
+ multiplierToUse = new BigFloat(tenRCP);\r
+ powerTemp = -exp;\r
+ }\r
+ else\r
+ {\r
+ multiplierToUse = new BigFloat(ten);\r
+ }\r
+\r
+ //Fast power function\r
+ while (powerTemp != 0)\r
+ {\r
+ temp.Mul(multiplierToUse);\r
+ multiplierToUse.Assign(temp);\r
+\r
+ if ((powerTemp & 1) != 0)\r
+ {\r
+ acc.Mul(temp);\r
+ }\r
+\r
+ powerTemp >>= 1;\r
+ }\r
+\r
+ iPart.Add(fPart);\r
+ iPart.Mul(acc);\r
+ iPart.mantissa.Sign = sign;\r
+ exponent = iPart.exponent;\r
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
+\r
+ return;\r
+ }\r
+\r
+ iPart.Add(fPart);\r
+ iPart.mantissa.Sign = sign;\r
+ exponent = iPart.exponent;\r
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
+\r
+ }\r
+\r
+ private void Init(PrecisionSpec mantissaPrec)\r
+ {\r
+ int mbWords = ((mantissaPrec.NumBits) >> 5);\r
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
+ int newManBits = mbWords << 5;\r
+\r
+ //For efficiency, we just use a 32-bit exponent\r
+ exponent = 0;\r
+ mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+ //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
+ }\r
+\r
+ //************************** Properties *************************\r
+\r
+ /// <summary>\r
+ /// Read-only property. Returns the precision specification of the mantissa.\r
+ /// \r
+ /// Floating point numbers are represented as 2^exponent * mantissa, where the\r
+ /// mantissa and exponent are integers. Note that the exponent in this class is\r
+ /// always a 32-bit integer. The precision therefore specifies how many bits\r
+ /// the mantissa will have.\r
+ /// </summary>\r
+ public PrecisionSpec Precision\r
+ {\r
+ get { return mantissa.Precision; }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Writable property:\r
+ /// true iff the number is negative or in some cases zero (<0)\r
+ /// false iff the number if positive or in some cases zero (>0)\r
+ /// </summary>\r
+ public bool Sign \r
+ { \r
+ get { return mantissa.Sign; }\r
+ set { mantissa.Sign = value; }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Read-only property. \r
+ /// True if the number is NAN, INF_PLUS, INF_MINUS or ZERO\r
+ /// False if the number has any other value.\r
+ /// </summary>\r
+ public bool IsSpecialValue\r
+ {\r
+ get\r
+ {\r
+ return (exponent == Int32.MaxValue || mantissa.IsZero());\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Read-only property, returns the type of number this is. Special values include:\r
+ /// \r
+ /// NONE - a regular number\r
+ /// ZERO - zero\r
+ /// NAN - Not a Number (some operations will return this if their inputs are out of range)\r
+ /// INF_PLUS - Positive infinity, not really a number, but a valid input to and output of some functions.\r
+ /// INF_MINUS - Negative infinity, not really a number, but a valid input to and output of some functions.\r
+ /// </summary>\r
+ public SpecialValueType SpecialValue\r
+ {\r
+ get\r
+ {\r
+ if (exponent == Int32.MaxValue)\r
+ {\r
+ if (mantissa.IsZero())\r
+ {\r
+ if (mantissa.Sign) return SpecialValueType.INF_MINUS;\r
+ return SpecialValueType.INF_PLUS;\r
+ }\r
+\r
+ return SpecialValueType.NAN;\r
+ }\r
+ else\r
+ {\r
+ if (mantissa.IsZero()) return SpecialValueType.ZERO;\r
+ return SpecialValueType.NONE;\r
+ }\r
+ }\r
+ }\r
+\r
+ //******************** Mathematical Constants *******************\r
+\r
+ /// <summary>\r
+ /// Gets pi to the indicated precision\r
+ /// </summary>\r
+ /// <param name="precision">The precision to perform the calculation to</param>\r
+ /// <returns>pi (the ratio of the area of a circle to its diameter)</returns>\r
+ public static BigFloat GetPi(PrecisionSpec precision)\r
+ {\r
+ if (pi == null || precision.NumBits <= pi.mantissa.Precision.NumBits)\r
+ {\r
+ CalculatePi(precision.NumBits);\r
+ }\r
+\r
+ BigFloat ret = new BigFloat (precision);\r
+ ret.Assign(pi);\r
+\r
+ return ret;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Get e to the indicated precision\r
+ /// </summary>\r
+ /// <param name="precision">The preicision to perform the calculation to</param>\r
+ /// <returns>e (the number for which the d/dx(e^x) = e^x)</returns>\r
+ public static BigFloat GetE(PrecisionSpec precision)\r
+ {\r
+ if (eCache == null || eCache.mantissa.Precision.NumBits < precision.NumBits)\r
+ {\r
+ CalculateEOnly(precision.NumBits);\r
+ //CalculateFactorials(precision.NumBits);\r
+ }\r
+\r
+ BigFloat ret = new BigFloat(precision);\r
+ ret.Assign(eCache);\r
+\r
+ return ret;\r
+ }\r
+\r
+\r
+ //******************** Arithmetic Functions ********************\r
+\r
+ /// <summary>\r
+ /// Addition (this = this + n2)\r
+ /// </summary>\r
+ /// <param name="n2">The number to add</param>\r
+ public void Add(BigFloat n2)\r
+ {\r
+ if (SpecialValueAddTest(n2)) return;\r
+\r
+ if (scratch.Precision.NumBits != n2.mantissa.Precision.NumBits)\r
+ {\r
+ scratch = new BigInt(n2.mantissa.Precision);\r
+ }\r
+\r
+ if (exponent <= n2.exponent)\r
+ {\r
+ int diff = n2.exponent - exponent;\r
+ exponent = n2.exponent;\r
+\r
+ if (diff != 0)\r
+ {\r
+ mantissa.RSH(diff);\r
+ }\r
+\r
+ uint carry = mantissa.Add(n2.mantissa);\r
+\r
+ if (carry != 0)\r
+ {\r
+ mantissa.RSH(1);\r
+ mantissa.SetBit(mantissa.Precision.NumBits - 1);\r
+ exponent++;\r
+ }\r
+\r
+ exponent -= mantissa.Normalise();\r
+ }\r
+ else\r
+ {\r
+ int diff = exponent - n2.exponent;\r
+\r
+ scratch.Assign(n2.mantissa);\r
+ scratch.RSH(diff);\r
+\r
+ uint carry = scratch.Add(mantissa);\r
+\r
+ if (carry != 0)\r
+ {\r
+ scratch.RSH(1);\r
+ scratch.SetBit(mantissa.Precision.NumBits - 1);\r
+ exponent++;\r
+ }\r
+\r
+ mantissa.Assign(scratch);\r
+\r
+ exponent -= mantissa.Normalise();\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Subtraction (this = this - n2)\r
+ /// </summary>\r
+ /// <param name="n2">The number to subtract from this</param>\r
+ public void Sub(BigFloat n2)\r
+ {\r
+ n2.mantissa.Sign = !n2.mantissa.Sign;\r
+ Add(n2);\r
+ n2.mantissa.Sign = !n2.mantissa.Sign;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Multiplication (this = this * n2)\r
+ /// </summary>\r
+ /// <param name="n2">The number to multiply this by</param>\r
+ public void Mul(BigFloat n2)\r
+ {\r
+ if (SpecialValueMulTest(n2)) return;\r
+\r
+ //Anything times 0 = 0\r
+ if (n2.mantissa.IsZero())\r
+ {\r
+ mantissa.Assign(n2.mantissa);\r
+ exponent = 0;\r
+ return;\r
+ }\r
+\r
+ mantissa.MulHi(n2.mantissa);\r
+ int shift = mantissa.Normalise();\r
+ exponent = exponent + n2.exponent + 1 - shift;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Division (this = this / n2)\r
+ /// </summary>\r
+ /// <param name="n2">The number to divide this by</param>\r
+ public void Div(BigFloat n2)\r
+ {\r
+ if (SpecialValueDivTest(n2)) return;\r
+\r
+ if (mantissa.Precision.NumBits >= 8192)\r
+ {\r
+ BigFloat rcp = n2.Reciprocal();\r
+ Mul(rcp);\r
+ }\r
+ else\r
+ {\r
+ int shift = mantissa.DivAndShift(n2.mantissa);\r
+ exponent = exponent - (n2.exponent + shift);\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Multiply by a power of 2 (-ve implies division)\r
+ /// </summary>\r
+ /// <param name="pow2"></param>\r
+ public void MulPow2(int pow2)\r
+ {\r
+ exponent += pow2;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Division-based reciprocal, fastest for small precisions up to 15,000 bits.\r
+ /// </summary>\r
+ /// <returns>The reciprocal 1/this</returns>\r
+ public BigFloat Reciprocal()\r
+ {\r
+ if (mantissa.Precision.NumBits >= 8192) return ReciprocalNewton();\r
+\r
+ BigFloat reciprocal = new BigFloat(1u, mantissa.Precision);\r
+ reciprocal.Div(this);\r
+ return reciprocal;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.\r
+ /// </summary>\r
+ /// <returns>The reciprocal 1/this</returns>\r
+ public BigFloat ReciprocalNewton()\r
+ {\r
+ if (mantissa.IsZero())\r
+ {\r
+ exponent = Int32.MaxValue;\r
+ return null;\r
+ }\r
+\r
+ bool oldSign = mantissa.Sign;\r
+ int oldExponent = exponent;\r
+\r
+ //Kill exponent for now (will re-institute later)\r
+ exponent = 0;\r
+\r
+ bool topBit = mantissa.IsTopBitOnlyBit();\r
+\r
+ PrecisionSpec curPrec = new PrecisionSpec(32, PrecisionSpec.BaseType.BIN);\r
+\r
+ BigFloat reciprocal = new BigFloat(curPrec);\r
+ BigFloat constant2 = new BigFloat(curPrec);\r
+ BigFloat temp = new BigFloat(curPrec);\r
+ BigFloat thisPrec = new BigFloat(this, curPrec);\r
+\r
+ reciprocal.exponent = 1;\r
+ reciprocal.mantissa.SetHighDigit(3129112985u);\r
+\r
+ constant2.exponent = 1;\r
+ constant2.mantissa.SetHighDigit(0x80000000u);\r
+\r
+ //D is deliberately left negative for all the following operations.\r
+ thisPrec.mantissa.Sign = true;\r
+\r
+ //Initial estimate.\r
+ reciprocal.Add(thisPrec);\r
+\r
+ //mantissa.Sign = false;\r
+\r
+ //Shift down into 0.5 < this < 1 range\r
+ thisPrec.mantissa.RSH(1);\r
+\r
+ //Iteration.\r
+ int accuracyBits = 2;\r
+ int mantissaBits = mantissa.Precision.NumBits;\r
+\r
+ //Each iteration is a pass of newton's method for RCP.\r
+ //The is a substantial optimisation to be done here...\r
+ //You can double the number of bits for the calculations\r
+ //at each iteration, meaning that the whole process only\r
+ //takes some constant multiplier of the time for the\r
+ //full-scale multiplication.\r
+ while (accuracyBits < mantissaBits)\r
+ {\r
+ //Increase the precision as needed\r
+ if (accuracyBits >= curPrec.NumBits / 2)\r
+ {\r
+ int newBits = curPrec.NumBits * 2;\r
+ if (newBits > mantissaBits) newBits = mantissaBits;\r
+ curPrec = new PrecisionSpec(newBits, PrecisionSpec.BaseType.BIN);\r
+\r
+ reciprocal = new BigFloat(reciprocal, curPrec);\r
+\r
+ constant2 = new BigFloat(curPrec);\r
+ constant2.exponent = 1;\r
+ constant2.mantissa.SetHighDigit(0x80000000u);\r
+\r
+ temp = new BigFloat(temp, curPrec);\r
+\r
+ thisPrec = new BigFloat(this, curPrec);\r
+ thisPrec.mantissa.Sign = true;\r
+ thisPrec.mantissa.RSH(1);\r
+ }\r
+\r
+ //temp = Xn\r
+ temp.exponent = reciprocal.exponent;\r
+ temp.mantissa.Assign(reciprocal.mantissa);\r
+ //temp = -Xn * D\r
+ temp.Mul(thisPrec);\r
+ //temp = -Xn * D + 2 (= 2 - Xn * D)\r
+ temp.Add(constant2);\r
+ //reciprocal = X(n+1) = Xn * (2 - Xn * D)\r
+ reciprocal.Mul(temp);\r
+\r
+ accuracyBits *= 2;\r
+ }\r
+\r
+ //'reciprocal' is now the reciprocal of the shifted down, zero-exponent mantissa of 'this'\r
+ //Restore the mantissa.\r
+ //mantissa.LSH(1);\r
+ exponent = oldExponent;\r
+ //mantissa.Sign = oldSign;\r
+\r
+ if (topBit)\r
+ {\r
+ reciprocal.exponent = -(oldExponent);\r
+ }\r
+ else\r
+ {\r
+ reciprocal.exponent = -(oldExponent + 1);\r
+ }\r
+ reciprocal.mantissa.Sign = oldSign;\r
+\r
+ return reciprocal;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.\r
+ /// </summary>\r
+ /// <returns>The reciprocal 1/this</returns>\r
+ private BigFloat ReciprocalNewton2()\r
+ {\r
+ if (mantissa.IsZero())\r
+ {\r
+ exponent = Int32.MaxValue;\r
+ return null;\r
+ }\r
+\r
+ bool oldSign = mantissa.Sign;\r
+ int oldExponent = exponent;\r
+\r
+ //Kill exponent for now (will re-institute later)\r
+ exponent = 0;\r
+\r
+ BigFloat reciprocal = new BigFloat(mantissa.Precision);\r
+ BigFloat constant2 = new BigFloat(mantissa.Precision);\r
+ BigFloat temp = new BigFloat(mantissa.Precision);\r
+\r
+ reciprocal.exponent = 1;\r
+ reciprocal.mantissa.SetHighDigit(3129112985u);\r
+\r
+ constant2.exponent = 1;\r
+ constant2.mantissa.SetHighDigit(0x80000000u);\r
+\r
+ //D is deliberately left negative for all the following operations.\r
+ mantissa.Sign = true;\r
+\r
+ //Initial estimate.\r
+ reciprocal.Add(this);\r
+\r
+ //mantissa.Sign = false;\r
+\r
+ //Shift down into 0.5 < this < 1 range\r
+ mantissa.RSH(1);\r
+ \r
+ //Iteration.\r
+ int accuracyBits = 2;\r
+ int mantissaBits = mantissa.Precision.NumBits;\r
+\r
+ //Each iteration is a pass of newton's method for RCP.\r
+ //The is a substantial optimisation to be done here...\r
+ //You can double the number of bits for the calculations\r
+ //at each iteration, meaning that the whole process only\r
+ //takes some constant multiplier of the time for the\r
+ //full-scale multiplication.\r
+ while (accuracyBits < mantissaBits)\r
+ {\r
+ //temp = Xn\r
+ temp.exponent = reciprocal.exponent;\r
+ temp.mantissa.Assign(reciprocal.mantissa);\r
+ //temp = -Xn * D\r
+ temp.Mul(this);\r
+ //temp = -Xn * D + 2 (= 2 - Xn * D)\r
+ temp.Add(constant2);\r
+ //reciprocal = X(n+1) = Xn * (2 - Xn * D)\r
+ reciprocal.Mul(temp);\r
+\r
+ accuracyBits *= 2;\r
+ }\r
+\r
+ //'reciprocal' is now the reciprocal of the shifted down, zero-exponent mantissa of 'this'\r
+ //Restore the mantissa.\r
+ mantissa.LSH(1);\r
+ exponent = oldExponent;\r
+ mantissa.Sign = oldSign;\r
+\r
+ reciprocal.exponent = -(oldExponent + 1);\r
+ reciprocal.mantissa.Sign = oldSign;\r
+\r
+ return reciprocal;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Sets this equal to the input\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public void Assign(BigFloat n2)\r
+ {\r
+ exponent = n2.exponent;\r
+ if (mantissa.AssignHigh(n2.mantissa)) exponent++;\r
+ }\r
+\r
+\r
+ //********************* Comparison Functions *******************\r
+\r
+ /// <summary>\r
+ /// Greater than comparison\r
+ /// </summary>\r
+ /// <param name="n2">the number to compare this to</param>\r
+ /// <returns>true iff this is greater than n2 (this > n2)</returns>\r
+ public bool GreaterThan(BigFloat n2)\r
+ {\r
+ if (IsSpecialValue || n2.IsSpecialValue)\r
+ {\r
+ SpecialValueType s1 = SpecialValue;\r
+ SpecialValueType s2 = SpecialValue;\r
+\r
+ if (s1 == SpecialValueType.NAN || s2 == SpecialValueType.NAN) return false;\r
+ if (s1 == SpecialValueType.INF_MINUS) return false;\r
+ if (s2 == SpecialValueType.INF_PLUS) return false;\r
+ if (s1 == SpecialValueType.INF_PLUS) return true;\r
+ if (s2 == SpecialValueType.INF_MINUS) return true;\r
+\r
+ if (s1 == SpecialValueType.ZERO)\r
+ {\r
+ if (s2 != SpecialValueType.ZERO && n2.Sign)\r
+ {\r
+ return true;\r
+ }\r
+ else\r
+ {\r
+ return false;\r
+ }\r
+ }\r
+\r
+ if (s2 == SpecialValueType.ZERO)\r
+ {\r
+ return !Sign;\r
+ }\r
+ }\r
+\r
+ if (!mantissa.Sign && n2.mantissa.Sign) return true;\r
+ if (mantissa.Sign && !n2.mantissa.Sign) return false;\r
+ if (!mantissa.Sign)\r
+ {\r
+ if (exponent > n2.exponent) return true;\r
+ if (exponent < n2.exponent) return false;\r
+ }\r
+ if (mantissa.Sign)\r
+ {\r
+ if (exponent > n2.exponent) return false;\r
+ if (exponent < n2.exponent) return true;\r
+ }\r
+\r
+ return mantissa.GreaterThan(n2.mantissa);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Less than comparison\r
+ /// </summary>\r
+ /// <param name="n2">the number to compare this to</param>\r
+ /// <returns>true iff this is less than n2 (this < n2)</returns>\r
+ public bool LessThan(BigFloat n2)\r
+ {\r
+ if (IsSpecialValue || n2.IsSpecialValue)\r
+ {\r
+ SpecialValueType s1 = SpecialValue;\r
+ SpecialValueType s2 = SpecialValue;\r
+\r
+ if (s1 == SpecialValueType.NAN || s2 == SpecialValueType.NAN) return false;\r
+ if (s1 == SpecialValueType.INF_PLUS) return false;\r
+ if (s2 == SpecialValueType.INF_PLUS) return true;\r
+ if (s2 == SpecialValueType.INF_MINUS) return false;\r
+ if (s1 == SpecialValueType.INF_MINUS) return true;\r
+\r
+ if (s1 == SpecialValueType.ZERO)\r
+ {\r
+ if (s2 != SpecialValueType.ZERO && !n2.Sign)\r
+ {\r
+ return true;\r
+ }\r
+ else\r
+ {\r
+ return false;\r
+ }\r
+ }\r
+\r
+ if (s2 == SpecialValueType.ZERO)\r
+ {\r
+ return Sign;\r
+ }\r
+ }\r
+\r
+ if (!mantissa.Sign && n2.mantissa.Sign) return false;\r
+ if (mantissa.Sign && !n2.mantissa.Sign) return true;\r
+ if (!mantissa.Sign)\r
+ {\r
+ if (exponent > n2.exponent) return false;\r
+ if (exponent < n2.exponent) return true;\r
+ }\r
+ if (mantissa.Sign)\r
+ {\r
+ if (exponent > n2.exponent) return true;\r
+ if (exponent < n2.exponent) return false;\r
+ }\r
+\r
+ return mantissa.LessThan(n2.mantissa);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Greater than comparison\r
+ /// </summary>\r
+ /// <param name="i">the number to compare this to</param>\r
+ /// <returns>true iff this is greater than n2 (this > n2)</returns>\r
+ public bool GreaterThan(int i)\r
+ {\r
+ BigFloat integer = new BigFloat(i, mantissa.Precision);\r
+ return GreaterThan(integer);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Less than comparison\r
+ /// </summary>\r
+ /// <param name="i">the number to compare this to</param>\r
+ /// <returns>true iff this is less than n2 (this < n2)</returns>\r
+ public bool LessThan(int i)\r
+ {\r
+ BigFloat integer = new BigFloat(i, mantissa.Precision);\r
+ return LessThan(integer);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Compare to zero\r
+ /// </summary>\r
+ /// <returns>true if this is zero (this == 0)</returns>\r
+ public bool IsZero()\r
+ {\r
+ return (mantissa.IsZero());\r
+ }\r
+\r
+\r
+ //******************** Mathematical Functions ******************\r
+\r
+ /// <summary>\r
+ /// Sets the number to the biggest integer numerically closer to zero, if possible.\r
+ /// </summary>\r
+ public void Floor()\r
+ {\r
+ //Already an integer.\r
+ if (exponent >= mantissa.Precision.NumBits) return;\r
+\r
+ if (exponent < 0)\r
+ {\r
+ mantissa.ZeroBits(mantissa.Precision.NumBits);\r
+ exponent = 0;\r
+ return;\r
+ }\r
+\r
+ mantissa.ZeroBits(mantissa.Precision.NumBits - (exponent + 1));\r
+ }\r
+\r
+ /// <summary>\r
+ /// Sets the number to its fractional component (equivalent to 'this' - (int)'this')\r
+ /// </summary>\r
+ public void FPart()\r
+ {\r
+ //Already fractional\r
+ if (exponent < 0)\r
+ {\r
+ return;\r
+ }\r
+\r
+ //Has no fractional part\r
+ if (exponent >= mantissa.Precision.NumBits)\r
+ {\r
+ mantissa.Zero();\r
+ exponent = 0;\r
+ return;\r
+ }\r
+\r
+ mantissa.ZeroBitsHigh(exponent + 1);\r
+ exponent -= mantissa.Normalise();\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates tan(x)\r
+ /// </summary>\r
+ public void Tan()\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ //Tan(x) has no limit as x->inf\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
+ {\r
+ SetNaN();\r
+ }\r
+ else if (SpecialValue == SpecialValueType.ZERO)\r
+ {\r
+ SetZero();\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
+ {\r
+ CalculatePi(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ //Work out the sign change (involves replicating some rescaling).\r
+ bool sign = mantissa.Sign;\r
+ mantissa.Sign = false;\r
+\r
+ if (mantissa.IsZero())\r
+ {\r
+ return;\r
+ }\r
+\r
+ //Rescale into 0 <= x < pi\r
+ if (GreaterThan(pi))\r
+ {\r
+ //There will be an inherent loss of precision doing this.\r
+ BigFloat newAngle = new BigFloat(this);\r
+ newAngle.Mul(piRecip);\r
+ newAngle.FPart();\r
+ newAngle.Mul(pi);\r
+ Assign(newAngle);\r
+ }\r
+\r
+ //Rescale to -pi/2 <= x < pi/2\r
+ if (!LessThan(piBy2))\r
+ {\r
+ Sub(pi);\r
+ }\r
+\r
+ //Now the sign of the sin determines the sign of the tan.\r
+ //tan(x) = sin(x) / sqrt(1 - sin^2(x))\r
+ Sin();\r
+ BigFloat denom = new BigFloat(this);\r
+ denom.Mul(this);\r
+ denom.Sub(new BigFloat(1, mantissa.Precision));\r
+ denom.mantissa.Sign = !denom.mantissa.Sign;\r
+\r
+ if (denom.mantissa.Sign)\r
+ {\r
+ denom.SetZero();\r
+ }\r
+\r
+ denom.Sqrt();\r
+ Div(denom);\r
+ if (sign) mantissa.Sign = !mantissa.Sign;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates Cos(x)\r
+ /// </summary>\r
+ public void Cos()\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ //Cos(x) has no limit as x->inf\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
+ {\r
+ SetNaN();\r
+ }\r
+ else if (SpecialValue == SpecialValueType.ZERO)\r
+ {\r
+ Assign(new BigFloat(1, mantissa.Precision));\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
+ {\r
+ CalculatePi(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ Add(piBy2);\r
+ Sin();\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates Sin(x):\r
+ /// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].\r
+ /// </summary>\r
+ public void Sin()\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ //Sin(x) has no limit as x->inf\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
+ {\r
+ SetNaN();\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ //Convert to positive range (0 <= x < inf)\r
+ bool sign = mantissa.Sign;\r
+ mantissa.Sign = false;\r
+\r
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
+ {\r
+ CalculatePi(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ if (inverseFactorialCache == null || invFactorialCutoff != mantissa.Precision.NumBits)\r
+ {\r
+ CalculateFactorials(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ //Rescale into 0 <= x < 2*pi\r
+ if (GreaterThan(twoPi))\r
+ {\r
+ //There will be an inherent loss of precision doing this.\r
+ BigFloat newAngle = new BigFloat(this);\r
+ newAngle.Mul(twoPiRecip);\r
+ newAngle.FPart();\r
+ newAngle.Mul(twoPi);\r
+ Assign(newAngle);\r
+ }\r
+\r
+ //Rescale into range 0 <= x < pi\r
+ if (GreaterThan(pi))\r
+ {\r
+ //sin(pi + a) = sin(pi)cos(a) + sin(a)cos(pi) = 0 - sin(a) = -sin(a)\r
+ Sub(pi);\r
+ sign = !sign;\r
+ }\r
+\r
+ BigFloat temp = new BigFloat(mantissa.Precision);\r
+\r
+ //Rescale into range 0 <= x < pi/2\r
+ if (GreaterThan(piBy2))\r
+ {\r
+ temp.Assign(this);\r
+ Assign(pi);\r
+ Sub(temp);\r
+ }\r
+\r
+ //Rescale into range 0 <= x < pi/6 to accelerate convergence.\r
+ //This is done using sin(3x) = 3sin(x) - 4sin^3(x)\r
+ Mul(threeRecip);\r
+\r
+ if (mantissa.IsZero())\r
+ {\r
+ exponent = 0;\r
+ return;\r
+ }\r
+\r
+ BigFloat term = new BigFloat(this);\r
+\r
+ BigFloat square = new BigFloat(this);\r
+ square.Mul(term);\r
+\r
+ BigFloat sum = new BigFloat(this);\r
+\r
+ bool termSign = true;\r
+ int length = inverseFactorialCache.Length;\r
+ int numBits = mantissa.Precision.NumBits;\r
+\r
+ for (int i = 3; i < length; i += 2)\r
+ {\r
+ term.Mul(square);\r
+ temp.Assign(inverseFactorialCache[i]);\r
+ temp.Mul(term);\r
+ temp.mantissa.Sign = termSign;\r
+ termSign = !termSign;\r
+\r
+ if (temp.exponent < -numBits) break;\r
+\r
+ sum.Add(temp);\r
+ }\r
+\r
+ //Restore the triple-angle: sin(3x) = 3sin(x) - 4sin^3(x)\r
+ Assign(sum);\r
+ sum.Mul(this);\r
+ sum.Mul(this);\r
+ Mul(new BigFloat(3, mantissa.Precision));\r
+ sum.exponent += 2;\r
+ Sub(sum);\r
+\r
+ //Restore the sign\r
+ mantissa.Sign = sign;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Hyperbolic Sin (sinh) function\r
+ /// </summary>\r
+ public void Sinh()\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ return;\r
+ }\r
+\r
+ Exp();\r
+ Sub(Reciprocal());\r
+ exponent--;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Hyperbolic cosine (cosh) function\r
+ /// </summary>\r
+ public void Cosh()\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ if (SpecialValue == SpecialValueType.ZERO)\r
+ {\r
+ Assign(new BigFloat(1, mantissa.Precision));\r
+ }\r
+ else if (SpecialValue == SpecialValueType.INF_MINUS)\r
+ {\r
+ SetInfPlus();\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ Exp();\r
+ Add(Reciprocal());\r
+ exponent--;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Hyperbolic tangent function (tanh)\r
+ /// </summary>\r
+ public void Tanh()\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ if (SpecialValue == SpecialValueType.INF_MINUS)\r
+ {\r
+ Assign(new BigFloat(-1, mantissa.Precision));\r
+ }\r
+ else if (SpecialValue == SpecialValueType.INF_PLUS)\r
+ {\r
+ Assign(new BigFloat(1, mantissa.Precision));\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ exponent++;\r
+ Exp();\r
+ BigFloat temp = new BigFloat(this);\r
+ BigFloat one = new BigFloat(1, mantissa.Precision);\r
+ temp.Add(one);\r
+ Sub(one);\r
+ Div(temp);\r
+ }\r
+\r
+ /// <summary>\r
+ /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)\r
+ /// </summary>\r
+ public void Arcsin()\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)\r
+ {\r
+ SetNaN();\r
+ return;\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ BigFloat one = new BigFloat(1, mantissa.Precision);\r
+ BigFloat plusABit = new BigFloat(1, mantissa.Precision);\r
+ plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));\r
+ BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);\r
+ onePlusABit.Add(plusABit);\r
+\r
+ bool sign = mantissa.Sign;\r
+ mantissa.Sign = false;\r
+\r
+ if (GreaterThan(onePlusABit))\r
+ {\r
+ SetNaN();\r
+ }\r
+ else if (LessThan(one))\r
+ {\r
+ BigFloat temp = new BigFloat(this);\r
+ temp.Mul(this);\r
+ temp.Sub(one);\r
+ temp.mantissa.Sign = !temp.mantissa.Sign;\r
+ temp.Sqrt();\r
+ temp.Add(one);\r
+ Div(temp);\r
+ Arctan();\r
+ exponent++;\r
+ mantissa.Sign = sign;\r
+ }\r
+ else\r
+ {\r
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
+ {\r
+ CalculatePi(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ Assign(piBy2);\r
+ if (sign) mantissa.Sign = true;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// arccos(): the inverse function of cos(), range (0..pi)\r
+ /// </summary>\r
+ public void Arccos()\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)\r
+ {\r
+ SetNaN();\r
+ }\r
+ else if (SpecialValue == SpecialValueType.ZERO)\r
+ {\r
+ Assign(new BigFloat(1, mantissa.Precision));\r
+ exponent = 0;\r
+ Sign = false;\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ BigFloat one = new BigFloat(1, mantissa.Precision);\r
+ BigFloat plusABit = new BigFloat(1, mantissa.Precision);\r
+ plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));\r
+ BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);\r
+ onePlusABit.Add(plusABit);\r
+\r
+ bool sign = mantissa.Sign;\r
+ mantissa.Sign = false;\r
+\r
+ if (GreaterThan(onePlusABit))\r
+ {\r
+ SetNaN();\r
+ }\r
+ else if (LessThan(one))\r
+ {\r
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
+ {\r
+ CalculatePi(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ mantissa.Sign = sign;\r
+ BigFloat temp = new BigFloat(this);\r
+ Mul(temp);\r
+ Sub(one);\r
+ mantissa.Sign = !mantissa.Sign;\r
+ Sqrt();\r
+ temp.Add(one);\r
+ Div(temp);\r
+ Arctan();\r
+ exponent++;\r
+ }\r
+ else\r
+ {\r
+ if (sign)\r
+ {\r
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
+ {\r
+ CalculatePi(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ Assign(pi);\r
+ }\r
+ else\r
+ {\r
+ mantissa.Zero();\r
+ exponent = 0;\r
+ }\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// arctan(): the inverse function of sin(), range of (-pi/2..pi/2)\r
+ /// </summary>\r
+ public void Arctan()\r
+ {\r
+ //With 2 argument reductions, we increase precision by a minimum of 4 bits per term.\r
+ int numBits = mantissa.Precision.NumBits;\r
+ int maxTerms = numBits >> 2;\r
+\r
+ if (pi == null || pi.mantissa.Precision.NumBits != numBits)\r
+ {\r
+ CalculatePi(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ //Make domain positive\r
+ bool sign = mantissa.Sign;\r
+ mantissa.Sign = false;\r
+\r
+ if (IsSpecialValue)\r
+ {\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
+ {\r
+ Assign(piBy2);\r
+ mantissa.Sign = sign;\r
+ return;\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ if (reciprocals == null || reciprocals[0].mantissa.Precision.NumBits != numBits || reciprocals.Length < maxTerms)\r
+ {\r
+ CalculateReciprocals(numBits, maxTerms);\r
+ }\r
+\r
+ bool invert = false;\r
+ BigFloat one = new BigFloat(1, mantissa.Precision);\r
+\r
+ //Invert if outside of convergence\r
+ if (GreaterThan(one))\r
+ {\r
+ invert = true;\r
+ Assign(Reciprocal());\r
+ }\r
+\r
+ //Reduce using half-angle formula:\r
+ //arctan(2x) = 2 arctan (x / (1 + sqrt(1 + x)))\r
+\r
+ //First reduction (guarantees 2 bits per iteration)\r
+ BigFloat temp = new BigFloat(this);\r
+ temp.Mul(this);\r
+ temp.Add(one);\r
+ temp.Sqrt();\r
+ temp.Add(one);\r
+ this.Div(temp);\r
+\r
+ //Second reduction (guarantees 4 bits per iteration)\r
+ temp.Assign(this);\r
+ temp.Mul(this);\r
+ temp.Add(one);\r
+ temp.Sqrt();\r
+ temp.Add(one);\r
+ this.Div(temp);\r
+\r
+ //Actual series calculation\r
+ int length = reciprocals.Length;\r
+ BigFloat term = new BigFloat(this);\r
+\r
+ //pow = x^2\r
+ BigFloat pow = new BigFloat(this);\r
+ pow.Mul(this);\r
+\r
+ BigFloat sum = new BigFloat(this);\r
+\r
+ for (int i = 1; i < length; i++)\r
+ {\r
+ //u(n) = u(n-1) * x^2\r
+ //t(n) = u(n) / (2n+1)\r
+ term.Mul(pow);\r
+ term.Sign = !term.Sign;\r
+ temp.Assign(term);\r
+ temp.Mul(reciprocals[i]);\r
+\r
+ if (temp.exponent < -numBits) break;\r
+\r
+ sum.Add(temp);\r
+ }\r
+\r
+ //Undo the reductions.\r
+ Assign(sum);\r
+ exponent += 2;\r
+\r
+ if (invert)\r
+ {\r
+ //Assign(Reciprocal());\r
+ mantissa.Sign = true;\r
+ Add(piBy2);\r
+ }\r
+\r
+ if (sign)\r
+ {\r
+ mantissa.Sign = sign;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Arcsinh(): the inverse sinh function\r
+ /// </summary>\r
+ public void Arcsinh()\r
+ {\r
+ //Just let all special values fall through\r
+ if (IsSpecialValue)\r
+ {\r
+ return;\r
+ }\r
+\r
+ BigFloat temp = new BigFloat(this);\r
+ temp.Mul(this);\r
+ temp.Add(new BigFloat(1, mantissa.Precision));\r
+ temp.Sqrt();\r
+ Add(temp);\r
+ Log();\r
+ }\r
+\r
+ /// <summary>\r
+ /// Arccosh(): the inverse cosh() function\r
+ /// </summary>\r
+ public void Arccosh()\r
+ {\r
+ //acosh isn't defined for x < 1\r
+ if (IsSpecialValue)\r
+ {\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.ZERO)\r
+ {\r
+ SetNaN();\r
+ return;\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ BigFloat one = new BigFloat(1, mantissa.Precision);\r
+ if (LessThan(one))\r
+ {\r
+ SetNaN();\r
+ return;\r
+ }\r
+\r
+ BigFloat temp = new BigFloat(this);\r
+ temp.Mul(this);\r
+ temp.Sub(one);\r
+ temp.Sqrt();\r
+ Add(temp);\r
+ Log();\r
+ }\r
+\r
+ /// <summary>\r
+ /// Arctanh(): the inverse tanh function\r
+ /// </summary>\r
+ public void Arctanh()\r
+ {\r
+ //|x| <= 1 for a non-NaN output\r
+ if (IsSpecialValue)\r
+ {\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
+ {\r
+ SetNaN();\r
+ return;\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ BigFloat one = new BigFloat(1, mantissa.Precision);\r
+ BigFloat plusABit = new BigFloat(1, mantissa.Precision);\r
+ plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));\r
+ BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);\r
+ onePlusABit.Add(plusABit);\r
+\r
+ bool sign = mantissa.Sign;\r
+ mantissa.Sign = false;\r
+\r
+ if (GreaterThan(onePlusABit))\r
+ {\r
+ SetNaN();\r
+ }\r
+ else if (LessThan(one))\r
+ {\r
+ BigFloat temp = new BigFloat(this);\r
+ Add(one);\r
+ one.Sub(temp);\r
+ Div(one);\r
+ Log();\r
+ exponent--;\r
+ mantissa.Sign = sign;\r
+ }\r
+ else\r
+ {\r
+ if (sign)\r
+ {\r
+ SetInfMinus();\r
+ }\r
+ else\r
+ {\r
+ SetInfPlus();\r
+ }\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Two-variable iterative square root, taken from\r
+ /// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method\r
+ /// </summary>\r
+ public void Sqrt()\r
+ {\r
+ if (mantissa.Sign || IsSpecialValue)\r
+ {\r
+ if (SpecialValue == SpecialValueType.ZERO)\r
+ {\r
+ return;\r
+ }\r
+\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)\r
+ {\r
+ SetNaN();\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ BigFloat temp2;\r
+ BigFloat temp3 = new BigFloat(mantissa.Precision);\r
+ BigFloat three = new BigFloat(3, mantissa.Precision);\r
+\r
+ int exponentScale = 0;\r
+\r
+ //Rescale to 0.5 <= x < 2\r
+ if (exponent < -1)\r
+ {\r
+ int diff = -exponent;\r
+ if ((diff & 1) != 0)\r
+ {\r
+ diff--;\r
+ }\r
+\r
+ exponentScale = -diff;\r
+ exponent += diff;\r
+ }\r
+ else if (exponent > 0)\r
+ {\r
+ if ((exponent & 1) != 0)\r
+ {\r
+ exponentScale = exponent + 1;\r
+ exponent = -1;\r
+ }\r
+ else\r
+ {\r
+ exponentScale = exponent;\r
+ exponent = 0;\r
+ }\r
+ }\r
+\r
+ temp2 = new BigFloat(this);\r
+ temp2.Sub(new BigFloat(1, mantissa.Precision));\r
+\r
+ //if (temp2.mantissa.IsZero())\r
+ //{\r
+ // exponent += exponentScale;\r
+ // return;\r
+ //}\r
+\r
+ int numBits = mantissa.Precision.NumBits;\r
+\r
+ while ((exponent - temp2.exponent) < numBits && temp2.SpecialValue != SpecialValueType.ZERO)\r
+ {\r
+ //a(n+1) = an - an*cn / 2\r
+ temp3.Assign(this);\r
+ temp3.Mul(temp2);\r
+ temp3.MulPow2(-1);\r
+ this.Sub(temp3);\r
+\r
+ //c(n+1) = cn^2 * (cn - 3) / 4\r
+ temp3.Assign(temp2);\r
+ temp2.Sub(three);\r
+ temp2.Mul(temp3);\r
+ temp2.Mul(temp3);\r
+ temp2.MulPow2(-2);\r
+ }\r
+\r
+ exponent += (exponentScale >> 1);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The natural logarithm, ln(x)\r
+ /// </summary>\r
+ public void Log()\r
+ {\r
+ if (IsSpecialValue || mantissa.Sign)\r
+ {\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)\r
+ {\r
+ SetNaN();\r
+ }\r
+ else if (SpecialValue == SpecialValueType.ZERO)\r
+ {\r
+ SetInfMinus();\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ if (mantissa.Precision.NumBits >= 512)\r
+ {\r
+ LogAGM1();\r
+ return;\r
+ }\r
+\r
+ //Compute ln2.\r
+ if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)\r
+ {\r
+ CalculateLog2(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ Log2();\r
+ Mul(ln2cache);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Log to the base 10\r
+ /// </summary>\r
+ public void Log10()\r
+ {\r
+ if (IsSpecialValue || mantissa.Sign)\r
+ {\r
+ if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)\r
+ {\r
+ SetNaN();\r
+ }\r
+ else if (SpecialValue == SpecialValueType.ZERO)\r
+ {\r
+ SetInfMinus();\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ //Compute ln2.\r
+ if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)\r
+ {\r
+ CalculateLog2(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ Log();\r
+ Mul(log10recip);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The exponential function. Less accurate for high exponents, scales poorly with the number\r
+ /// of bits.\r
+ /// </summary>\r
+ public void Exp()\r
+ {\r
+ Exp(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Raises a number to an integer power (positive or negative). This is a very accurate and fast function,\r
+ /// comparable to or faster than division (although it is slightly slower for\r
+ /// negative powers, obviously)\r
+ /// \r
+ /// </summary>\r
+ /// <param name="power"></param>\r
+ public void Pow(int power)\r
+ {\r
+ BigFloat acc = new BigFloat(1, mantissa.Precision);\r
+ BigFloat temp = new BigFloat(1, mantissa.Precision);\r
+\r
+ int powerTemp = power;\r
+\r
+ if (power < 0)\r
+ {\r
+ Assign(Reciprocal());\r
+ powerTemp = -power;\r
+ }\r
+\r
+ //Fast power function\r
+ while (powerTemp != 0)\r
+ {\r
+ temp.Mul(this);\r
+ Assign(temp);\r
+\r
+ if ((powerTemp & 1) != 0)\r
+ {\r
+ acc.Mul(temp);\r
+ }\r
+\r
+ powerTemp >>= 1;\r
+ }\r
+\r
+ Assign(acc);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Raises to an aribitrary power. This is both slow (uses Log) and inaccurate. If you need to\r
+ /// raise e^x use exp(). If you need an integer power, use the integer power function Pow(int)\r
+ /// Accuracy Note:\r
+ /// The function is only ever accurate to a maximum of 4 decimal digits\r
+ /// For every 10x larger (or smaller) the power gets, you lose an additional decimal digit\r
+ /// If you really need a precise result, do the calculation with an extra 32-bits and round\r
+ /// Domain Note:\r
+ /// This only works for powers of positive real numbers. Negative numbers will fail.\r
+ /// </summary>\r
+ /// <param name="power"></param>\r
+ public void Pow(BigFloat power)\r
+ {\r
+ Log();\r
+ Mul(power);\r
+ Exp();\r
+ }\r
+\r
+\r
+ //******************** Static Math Functions *******************\r
+\r
+ /// <summary>\r
+ /// Returns the integer component of the input\r
+ /// </summary>\r
+ /// <param name="n1">The input number</param>\r
+ /// <remarks>The integer component returned will always be numerically closer to zero\r
+ /// than the input: an input of -3.49 for instance would produce a value of 3.</remarks>\r
+ public static BigFloat Floor(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Floor();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Returns the fractional (non-integer component of the input)\r
+ /// </summary>\r
+ /// <param name="n1">The input number</param>\r
+ public static BigFloat FPart(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.FPart();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates tan(x)\r
+ /// </summary>\r
+ /// <param name="n1">The angle (in radians) to find the tangent of</param>\r
+ public static BigFloat Tan(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Tan();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates Cos(x)\r
+ /// </summary>\r
+ /// <param name="n1">The angle (in radians) to find the cosine of</param>\r
+ /// <remarks>This is a reasonably fast function for smaller precisions, but\r
+ /// doesn't scale well for higher precision arguments</remarks>\r
+ public static BigFloat Cos(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Cos();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates Sin(x):\r
+ /// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].\r
+ /// </summary>\r
+ /// <param name="n1">The angle to find the sine of (in radians)</param>\r
+ /// <remarks>This is a resonably fast function, for smaller precision arguments, but doesn't\r
+ /// scale very well with the number of bits in the input.</remarks>\r
+ public static BigFloat Sin(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Sin();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Hyperbolic Sin (sinh) function\r
+ /// </summary>\r
+ /// <param name="n1">The number to find the hyperbolic sine of</param>\r
+ public static BigFloat Sinh(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Sinh();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Hyperbolic cosine (cosh) function\r
+ /// </summary>\r
+ /// <param name="n1">The number to find the hyperbolic cosine of</param>\r
+ public static BigFloat Cosh(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Cosh();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Hyperbolic tangent function (tanh)\r
+ /// </summary>\r
+ /// <param name="n1">The number to find the hyperbolic tangent of</param>\r
+ public static BigFloat Tanh(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Tanh();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)\r
+ /// </summary>\r
+ /// <param name="n1">The number to find the arcsine of (-pi/2..pi/2)</param>\r
+ /// <remarks>Note that inverse trig functions are only defined within a specific range.\r
+ /// Values outside this range will return NaN, although some margin for error is assumed.\r
+ /// </remarks>\r
+ public static BigFloat Arcsin(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Arcsin();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// arccos(): the inverse function of cos(), input range (0..pi)\r
+ /// </summary>\r
+ /// <param name="n1">The number to find the arccosine of (0..pi)</param>\r
+ /// <remarks>Note that inverse trig functions are only defined within a specific range.\r
+ /// Values outside this range will return NaN, although some margin for error is assumed.\r
+ /// </remarks>\r
+ public static BigFloat Arccos(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Arccos();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// arctan(): the inverse function of sin(), input range of (-pi/2..pi/2)\r
+ /// </summary>\r
+ /// <param name="n1">The number to find the arctangent of (-pi/2..pi/2)</param>\r
+ /// <remarks>Note that inverse trig functions are only defined within a specific range.\r
+ /// Values outside this range will return NaN, although some margin for error is assumed.\r
+ /// </remarks>\r
+ public static BigFloat Arctan(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Arctan();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Arcsinh(): the inverse sinh function\r
+ /// </summary>\r
+ /// <param name="n1">The number to find the inverse hyperbolic sine of</param>\r
+ public static BigFloat Arcsinh(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Arcsinh();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Arccosh(): the inverse cosh() function\r
+ /// </summary>\r
+ /// <param name="n1">The number to find the inverse hyperbolic cosine of</param>\r
+ public static BigFloat Arccosh(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Arccosh();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Arctanh(): the inverse tanh function\r
+ /// </summary>\r
+ /// <param name="n1">The number to fine the inverse hyperbolic tan of</param>\r
+ public static BigFloat Arctanh(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Arctanh();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Two-variable iterative square root, taken from\r
+ /// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method\r
+ /// </summary>\r
+ /// <remarks>This is quite a fast function, as elementary functions go. You can expect it to take\r
+ /// about twice as long as a floating-point division.\r
+ /// </remarks>\r
+ public static BigFloat Sqrt(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Sqrt();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// The natural logarithm, ln(x) (log base e)\r
+ /// </summary>\r
+ /// <remarks>This is a very slow function, despite repeated attempts at optimisation.\r
+ /// To make it any faster, different strategies would be needed for integer operations.\r
+ /// It does, however, scale well with the number of bits.\r
+ /// </remarks>\r
+ /// <param name="n1">The number to find the natural logarithm of</param>\r
+ public static BigFloat Log(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Log();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Base 10 logarithm of a number\r
+ /// </summary>\r
+ /// <remarks>This is a very slow function, despite repeated attempts at optimisation.\r
+ /// To make it any faster, different strategies would be needed for integer operations.\r
+ /// It does, however, scale well with the number of bits.\r
+ /// </remarks>\r
+ /// <param name="n1">The number to find the base 10 logarithm of</param>\r
+ public static BigFloat Log10(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Log10();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// The exponential function. Less accurate for high exponents, scales poorly with the number\r
+ /// of bits. This is quite fast for low-precision arguments.\r
+ /// </summary>\r
+ public static BigFloat Exp(BigFloat n1)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Exp();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Raises a number to an integer power (positive or negative). This is a very accurate and fast function,\r
+ /// comparable to or faster than division (although it is slightly slower for\r
+ /// negative powers, obviously).\r
+ /// </summary>\r
+ /// <param name="n1">The number to raise to the power</param>\r
+ /// <param name="power">The power to raise it to</param>\r
+ public static BigFloat Pow(BigFloat n1, int power)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Pow(power);\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Raises to an aribitrary power. This is both slow (uses Log) and inaccurate. If you need to\r
+ /// raise e^x use exp(). If you need an integer power, use the integer power function Pow(int)\r
+ /// </summary>\r
+ /// <remarks>\r
+ /// Accuracy Note:\r
+ /// The function is only ever accurate to a maximum of 4 decimal digits\r
+ /// For every 10x larger (or smaller) the power gets, you lose an additional decimal digit\r
+ /// If you really need a precise result, do the calculation with an extra 32-bits and round\r
+ /// \r
+ /// Domain Note:\r
+ /// This only works for powers of positive real numbers. Negative numbers will fail.\r
+ /// </remarks>\r
+ /// <param name="n1">The number to raise to a power</param>\r
+ /// <param name="power">The power to raise it to</param>\r
+ public static BigFloat Pow(BigFloat n1, BigFloat power)\r
+ {\r
+ BigFloat res = new BigFloat(n1);\r
+ n1.Pow(power);\r
+ return n1;\r
+ }\r
+\r
+ //********************** Static functions **********************\r
+\r
+ /// <summary>\r
+ /// Adds two numbers and returns the result\r
+ /// </summary>\r
+ public static BigFloat Add(BigFloat n1, BigFloat n2)\r
+ {\r
+ BigFloat ret = new BigFloat(n1);\r
+ ret.Add(n2);\r
+ return ret;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Subtracts two numbers and returns the result\r
+ /// </summary>\r
+ public static BigFloat Sub(BigFloat n1, BigFloat n2)\r
+ {\r
+ BigFloat ret = new BigFloat(n1);\r
+ ret.Sub(n2);\r
+ return ret;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Multiplies two numbers and returns the result\r
+ /// </summary>\r
+ public static BigFloat Mul(BigFloat n1, BigFloat n2)\r
+ {\r
+ BigFloat ret = new BigFloat(n1);\r
+ ret.Mul(n2);\r
+ return ret;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Divides two numbers and returns the result\r
+ /// </summary>\r
+ public static BigFloat Div(BigFloat n1, BigFloat n2)\r
+ {\r
+ BigFloat ret = new BigFloat(n1);\r
+ ret.Div(n2);\r
+ return ret;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Tests whether n1 is greater than n2\r
+ /// </summary>\r
+ public static bool GreaterThan(BigFloat n1, BigFloat n2)\r
+ {\r
+ return n1.GreaterThan(n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Tests whether n1 is less than n2\r
+ /// </summary>\r
+ public static bool LessThan(BigFloat n1, BigFloat n2)\r
+ {\r
+ return n1.LessThan(n2);\r
+ }\r
+\r
+\r
+ //******************* Fast static functions ********************\r
+\r
+ /// <summary>\r
+ /// Adds two numbers and assigns the result to res.\r
+ /// </summary>\r
+ /// <param name="res">a pre-existing BigFloat to take the result</param>\r
+ /// <param name="n1">the first number</param>\r
+ /// <param name="n2">the second number</param>\r
+ /// <returns>a handle to res</returns>\r
+ public static BigFloat Add(BigFloat res, BigFloat n1, BigFloat n2)\r
+ {\r
+ res.Assign(n1);\r
+ res.Add(n2);\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Subtracts two numbers and assigns the result to res.\r
+ /// </summary>\r
+ /// <param name="res">a pre-existing BigFloat to take the result</param>\r
+ /// <param name="n1">the first number</param>\r
+ /// <param name="n2">the second number</param>\r
+ /// <returns>a handle to res</returns>\r
+ public static BigFloat Sub(BigFloat res, BigFloat n1, BigFloat n2)\r
+ {\r
+ res.Assign(n1);\r
+ res.Sub(n2);\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Multiplies two numbers and assigns the result to res.\r
+ /// </summary>\r
+ /// <param name="res">a pre-existing BigFloat to take the result</param>\r
+ /// <param name="n1">the first number</param>\r
+ /// <param name="n2">the second number</param>\r
+ /// <returns>a handle to res</returns>\r
+ public static BigFloat Mul(BigFloat res, BigFloat n1, BigFloat n2)\r
+ {\r
+ res.Assign(n1);\r
+ res.Mul(n2);\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Divides two numbers and assigns the result to res.\r
+ /// </summary>\r
+ /// <param name="res">a pre-existing BigFloat to take the result</param>\r
+ /// <param name="n1">the first number</param>\r
+ /// <param name="n2">the second number</param>\r
+ /// <returns>a handle to res</returns>\r
+ public static BigFloat Div(BigFloat res, BigFloat n1, BigFloat n2)\r
+ {\r
+ res.Assign(n1);\r
+ res.Div(n2);\r
+ return res;\r
+ }\r
+\r
+\r
+ //************************* Operators **************************\r
+\r
+ /// <summary>\r
+ /// The addition operator\r
+ /// </summary>\r
+ public static BigFloat operator +(BigFloat n1, BigFloat n2)\r
+ {\r
+ return Add(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The subtraction operator\r
+ /// </summary>\r
+ public static BigFloat operator -(BigFloat n1, BigFloat n2)\r
+ {\r
+ return Sub(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The multiplication operator\r
+ /// </summary>\r
+ public static BigFloat operator *(BigFloat n1, BigFloat n2)\r
+ {\r
+ return Mul(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The division operator\r
+ /// </summary>\r
+ public static BigFloat operator /(BigFloat n1, BigFloat n2)\r
+ {\r
+ return Div(n1, n2);\r
+ }\r
+\r
+ //************************** Conversions *************************\r
+\r
+ /// <summary>\r
+ /// Converts a BigFloat to an BigInt with the specified precision\r
+ /// </summary>\r
+ /// <param name="n1">The number to convert</param>\r
+ /// <param name="precision">The precision to convert it with</param>\r
+ /// <param name="round">Do we round the number if we are truncating the mantissa?</param>\r
+ /// <returns></returns>\r
+ public static BigInt ConvertToInt(BigFloat n1, PrecisionSpec precision, bool round)\r
+ {\r
+ BigInt ret = new BigInt(precision);\r
+\r
+ int numBits = n1.mantissa.Precision.NumBits;\r
+ int shift = numBits - (n1.exponent + 1);\r
+\r
+ BigFloat copy = new BigFloat(n1);\r
+ bool inc = false;\r
+\r
+ //Rounding\r
+ if (copy.mantissa.Precision.NumBits > ret.Precision.NumBits)\r
+ {\r
+ inc = true;\r
+\r
+ for (int i = copy.exponent + 1; i <= ret.Precision.NumBits; i++)\r
+ {\r
+ if (copy.mantissa.GetBitFromTop(i) == 0)\r
+ {\r
+ inc = false;\r
+ break;\r
+ }\r
+ }\r
+ }\r
+\r
+ if (shift > 0)\r
+ {\r
+ copy.mantissa.RSH(shift);\r
+ }\r
+ else if (shift < 0)\r
+ {\r
+ copy.mantissa.LSH(-shift);\r
+ }\r
+\r
+ ret.Assign(copy.mantissa);\r
+\r
+ if (inc) ret.Increment();\r
+\r
+ return ret;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Returns a base-10 string representing the number.\r
+ /// \r
+ /// Note: This is inefficient and possibly inaccurate. Please use with enough\r
+ /// rounding digits (set using the RoundingDigits property) to ensure accuracy\r
+ /// </summary>\r
+ public override string ToString()\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ SpecialValueType s = SpecialValue;\r
+ if (s == SpecialValueType.ZERO)\r
+ {\r
+ return String.Format("0{0}0", System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator);\r
+ }\r
+ else if (s == SpecialValueType.INF_PLUS)\r
+ {\r
+ return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol;\r
+ }\r
+ else if (s == SpecialValueType.INF_MINUS)\r
+ {\r
+ return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol;\r
+ }\r
+ else if (s == SpecialValueType.NAN)\r
+ {\r
+ return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol;\r
+ }\r
+ else\r
+ {\r
+ return "Unrecognised special type";\r
+ }\r
+ }\r
+\r
+ if (scratch.Precision.NumBits != mantissa.Precision.NumBits)\r
+ {\r
+ scratch = new BigInt(mantissa.Precision);\r
+ }\r
+\r
+ //The mantissa expresses 1.xxxxxxxxxxx\r
+ //The highest possible value for the mantissa without the implicit 1. is 0.9999999...\r
+ scratch.Assign(mantissa);\r
+ //scratch.Round(3);\r
+ scratch.Sign = false;\r
+ BigInt denom = new BigInt("0", mantissa.Precision);\r
+ denom.SetBit(mantissa.Precision.NumBits - 1);\r
+\r
+ bool useExponentialNotation = false;\r
+ int halfBits = mantissa.Precision.NumBits / 2;\r
+ if (halfBits > 60) halfBits = 60;\r
+ int precDec = 10;\r
+\r
+ if (exponent > 0)\r
+ {\r
+ if (exponent < halfBits)\r
+ {\r
+ denom.RSH(exponent);\r
+ }\r
+ else\r
+ {\r
+ useExponentialNotation = true;\r
+ }\r
+ }\r
+ else if (exponent < 0)\r
+ {\r
+ int shift = -(exponent);\r
+ if (shift < precDec)\r
+ {\r
+ scratch.RSH(shift);\r
+ }\r
+ else\r
+ {\r
+ useExponentialNotation = true;\r
+ }\r
+ }\r
+\r
+ string output;\r
+\r
+ if (useExponentialNotation)\r
+ {\r
+ int absExponent = exponent;\r
+ if (absExponent < 0) absExponent = -absExponent;\r
+ int powerOf10 = (int)((double)absExponent * Math.Log10(2.0));\r
+\r
+ //Use 1 extra digit of precision (this is actually 32 bits more, nb)\r
+ BigFloat thisFloat = new BigFloat(this, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
+ thisFloat.mantissa.Sign = false;\r
+\r
+ //Multiplicative correction factor to bring number into range.\r
+ BigFloat one = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
+ BigFloat ten = new BigFloat(10, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
+ BigFloat tenRCP = ten.Reciprocal();\r
+\r
+ //Accumulator for the power of 10 calculation.\r
+ BigFloat acc = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
+\r
+ BigFloat tenToUse;\r
+\r
+ if (exponent > 0)\r
+ {\r
+ tenToUse = new BigFloat(tenRCP, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
+ }\r
+ else\r
+ {\r
+ tenToUse = new BigFloat(ten, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
+ }\r
+\r
+ BigFloat tenToPower = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
+\r
+ int powerTemp = powerOf10;\r
+\r
+ //Fast power function\r
+ while (powerTemp != 0)\r
+ {\r
+ tenToPower.Mul(tenToUse);\r
+ tenToUse.Assign(tenToPower);\r
+\r
+ if ((powerTemp & 1) != 0)\r
+ {\r
+ acc.Mul(tenToPower);\r
+ }\r
+\r
+ powerTemp >>= 1;\r
+ }\r
+\r
+ thisFloat.Mul(acc);\r
+\r
+ //If we are out of range, correct. \r
+ if (thisFloat.GreaterThan(ten))\r
+ {\r
+ thisFloat.Mul(tenRCP);\r
+ if (exponent > 0)\r
+ {\r
+ powerOf10++;\r
+ }\r
+ else\r
+ {\r
+ powerOf10--;\r
+ }\r
+ }\r
+ else if (thisFloat.LessThan(one))\r
+ {\r
+ thisFloat.Mul(ten);\r
+ if (exponent > 0)\r
+ {\r
+ powerOf10--;\r
+ }\r
+ else\r
+ {\r
+ powerOf10++;\r
+ }\r
+ }\r
+\r
+ //Restore the precision and the sign.\r
+ BigFloat printable = new BigFloat(thisFloat, mantissa.Precision);\r
+ printable.mantissa.Sign = mantissa.Sign;\r
+ output = printable.ToString();\r
+\r
+ if (exponent < 0) powerOf10 = -powerOf10;\r
+\r
+ output = String.Format("{0}E{1}", output, powerOf10);\r
+ }\r
+ else\r
+ {\r
+ BigInt bigDigit = BigInt.Div(scratch, denom);\r
+ bigDigit.Sign = false;\r
+ scratch.Sub(BigInt.Mul(denom, bigDigit));\r
+\r
+ if (mantissa.Sign)\r
+ {\r
+ output = String.Format("-{0}.", bigDigit);\r
+ }\r
+ else\r
+ {\r
+ output = String.Format("{0}.", bigDigit);\r
+ }\r
+\r
+ denom = BigInt.Div(denom, 10u);\r
+\r
+ while (!denom.IsZero())\r
+ {\r
+ uint digit = (uint)BigInt.Div(scratch, denom);\r
+ if (digit == 10) digit--;\r
+ scratch.Sub(BigInt.Mul(denom, digit));\r
+ output = String.Format("{0}{1}", output, digit);\r
+ denom = BigInt.Div(denom, 10u);\r
+ }\r
+\r
+ output = RoundString(output, RoundingDigits);\r
+ }\r
+\r
+ return output;\r
+ }\r
+\r
+ //**************** Special value handling for ops ***************\r
+\r
+ private void SetNaN()\r
+ {\r
+ exponent = Int32.MaxValue;\r
+ mantissa.SetBit(mantissa.Precision.NumBits - 1);\r
+ }\r
+\r
+ private void SetZero()\r
+ {\r
+ exponent = 0;\r
+ mantissa.Zero();\r
+ Sign = false;\r
+ }\r
+\r
+ private void SetInfPlus()\r
+ {\r
+ Sign = false;\r
+ exponent = Int32.MaxValue;\r
+ mantissa.Zero();\r
+ }\r
+\r
+ private void SetInfMinus()\r
+ {\r
+ Sign = true;\r
+ exponent = Int32.MaxValue;\r
+ mantissa.Zero();\r
+ }\r
+\r
+ private bool SpecialValueAddTest(BigFloat n2)\r
+ {\r
+ if (IsSpecialValue || n2.IsSpecialValue)\r
+ {\r
+ SpecialValueType s1 = SpecialValue;\r
+ SpecialValueType s2 = n2.SpecialValue;\r
+\r
+ if (s1 == SpecialValueType.NAN) return true;\r
+ if (s2 == SpecialValueType.NAN)\r
+ {\r
+ //Set NaN and return.\r
+ SetNaN();\r
+ return true;\r
+ }\r
+\r
+ if (s1 == SpecialValueType.INF_PLUS)\r
+ {\r
+ //INF+ + INF- = NAN\r
+ if (s2 == SpecialValueType.INF_MINUS)\r
+ {\r
+ SetNaN();\r
+ return true;\r
+ }\r
+\r
+ return true;\r
+ }\r
+\r
+ if (s1 == SpecialValueType.INF_MINUS)\r
+ {\r
+ //INF+ + INF- = NAN\r
+ if (s2 == SpecialValueType.INF_PLUS)\r
+ {\r
+ SetNaN();\r
+ return true;\r
+ }\r
+\r
+ return true;\r
+ }\r
+\r
+ if (s2 == SpecialValueType.ZERO)\r
+ {\r
+ return true;\r
+ }\r
+\r
+ if (s1 == SpecialValueType.ZERO)\r
+ {\r
+ Assign(n2);\r
+ return true;\r
+ }\r
+ }\r
+\r
+ return false;\r
+ }\r
+\r
+ private bool SpecialValueMulTest(BigFloat n2)\r
+ {\r
+ if (IsSpecialValue || n2.IsSpecialValue)\r
+ {\r
+ SpecialValueType s1 = SpecialValue;\r
+ SpecialValueType s2 = n2.SpecialValue;\r
+\r
+ if (s1 == SpecialValueType.NAN) return true;\r
+ if (s2 == SpecialValueType.NAN)\r
+ {\r
+ //Set NaN and return.\r
+ SetNaN();\r
+ return true;\r
+ }\r
+\r
+ if (s1 == SpecialValueType.INF_PLUS)\r
+ {\r
+ //Inf+ * Inf- = Inf-\r
+ if (s2 == SpecialValueType.INF_MINUS)\r
+ {\r
+ Assign(n2);\r
+ return true;\r
+ }\r
+\r
+ //Inf+ * 0 = NaN\r
+ if (s2 == SpecialValueType.ZERO)\r
+ {\r
+ //Set NaN and return.\r
+ SetNaN();\r
+ return true;\r
+ }\r
+\r
+ return true;\r
+ }\r
+\r
+ if (s1 == SpecialValueType.INF_MINUS)\r
+ {\r
+ //Inf- * Inf- = Inf+\r
+ if (s2 == SpecialValueType.INF_MINUS)\r
+ {\r
+ Sign = false;\r
+ return true;\r
+ }\r
+\r
+ //Inf- * 0 = NaN\r
+ if (s2 == SpecialValueType.ZERO)\r
+ {\r
+ //Set NaN and return.\r
+ SetNaN();\r
+ return true;\r
+ }\r
+\r
+ return true;\r
+ }\r
+\r
+ if (s2 == SpecialValueType.ZERO)\r
+ {\r
+ SetZero();\r
+ return true;\r
+ }\r
+\r
+ if (s1 == SpecialValueType.ZERO)\r
+ {\r
+ return true;\r
+ }\r
+ }\r
+\r
+ return false;\r
+ }\r
+\r
+ private bool SpecialValueDivTest(BigFloat n2)\r
+ {\r
+ if (IsSpecialValue || n2.IsSpecialValue)\r
+ {\r
+ SpecialValueType s1 = SpecialValue;\r
+ SpecialValueType s2 = n2.SpecialValue;\r
+\r
+ if (s1 == SpecialValueType.NAN) return true;\r
+ if (s2 == SpecialValueType.NAN)\r
+ {\r
+ //Set NaN and return.\r
+ SetNaN();\r
+ return true;\r
+ }\r
+\r
+ if ((s1 == SpecialValueType.INF_PLUS || s1 == SpecialValueType.INF_MINUS))\r
+ {\r
+ if (s2 == SpecialValueType.INF_PLUS || s2 == SpecialValueType.INF_MINUS)\r
+ {\r
+ //Set NaN and return.\r
+ SetNaN();\r
+ return true;\r
+ }\r
+\r
+ if (n2.Sign)\r
+ {\r
+ if (s1 == SpecialValueType.INF_PLUS)\r
+ {\r
+ SetInfMinus();\r
+ return true;\r
+ }\r
+\r
+ SetInfPlus();\r
+ return true;\r
+ }\r
+\r
+ //Keep inf\r
+ return true;\r
+ }\r
+\r
+ if (s2 == SpecialValueType.ZERO)\r
+ {\r
+ if (s1 == SpecialValueType.ZERO)\r
+ {\r
+ SetNaN();\r
+ return true;\r
+ }\r
+\r
+ if (Sign)\r
+ {\r
+ SetInfMinus();\r
+ return true;\r
+ }\r
+\r
+ SetInfPlus();\r
+ return true;\r
+ }\r
+ }\r
+\r
+ return false;\r
+ }\r
+\r
+ //****************** Internal helper functions *****************\r
+\r
+ /// <summary>\r
+ /// Used for fixed point speed-ups (where the extra precision is not required). Note that Denormalised\r
+ /// floats break the assumptions that underly Add() and Sub(), so they can only be used for multiplication\r
+ /// </summary>\r
+ /// <param name="targetExponent"></param>\r
+ private void Denormalise(int targetExponent)\r
+ {\r
+ int diff = targetExponent - exponent;\r
+ if (diff <= 0) return;\r
+\r
+ //This only works to reduce the precision, so if the difference implies an increase, we can't do anything.\r
+ mantissa.RSH(diff);\r
+ exponent += diff;\r
+ }\r
+\r
+ /// <summary>\r
+ /// The binary logarithm, log2(x) - for precisions above 1000 bits, use Log() and convert the base.\r
+ /// </summary>\r
+ private void Log2()\r
+ {\r
+ if (scratch.Precision.NumBits != mantissa.Precision.NumBits)\r
+ {\r
+ scratch = new BigInt(mantissa.Precision);\r
+ }\r
+\r
+ int bits = mantissa.Precision.NumBits;\r
+ BigFloat temp = new BigFloat(this);\r
+ BigFloat result = new BigFloat(exponent, mantissa.Precision);\r
+ BigFloat pow2 = new BigFloat(1, mantissa.Precision);\r
+ temp.exponent = 0;\r
+ int bitsCalculated = 0;\r
+\r
+ while (bitsCalculated < bits)\r
+ {\r
+ int i;\r
+ for (i = 0; (temp.exponent == 0); i++)\r
+ {\r
+ temp.mantissa.SquareHiFast(scratch);\r
+ int shift = temp.mantissa.Normalise();\r
+ temp.exponent += 1 - shift;\r
+ if (i + bitsCalculated >= bits) break;\r
+ }\r
+\r
+ pow2.MulPow2(-i);\r
+ result.Add(pow2);\r
+ temp.exponent = 0;\r
+ bitsCalculated += i;\r
+ }\r
+\r
+ this.Assign(result);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Tried the newton method for logs, but the exponential function is too slow to do it.\r
+ /// </summary>\r
+ private void LogNewton()\r
+ {\r
+ if (mantissa.IsZero() || mantissa.Sign)\r
+ {\r
+ return;\r
+ }\r
+\r
+ //Compute ln2.\r
+ if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)\r
+ {\r
+ CalculateLog2(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ int numBits = mantissa.Precision.NumBits;\r
+\r
+ //Use inverse exp function with Newton's method.\r
+ BigFloat xn = new BigFloat(this);\r
+ BigFloat oldExponent = new BigFloat(xn.exponent, mantissa.Precision);\r
+ xn.exponent = 0;\r
+ this.exponent = 0;\r
+ //Hack to subtract 1\r
+ xn.mantissa.ClearBit(numBits - 1);\r
+ //x0 = (x - 1) * log2 - this is a straight line fit between log(1) = 0 and log(2) = ln2\r
+ xn.Mul(ln2cache);\r
+ //x0 = (x - 1) * log2 + C - this corrects for minimum error over the range.\r
+ xn.Add(logNewtonConstant);\r
+ BigFloat term = new BigFloat(mantissa.Precision);\r
+ BigFloat one = new BigFloat(1, mantissa.Precision);\r
+\r
+ int precision = 32;\r
+ int normalPrecision = mantissa.Precision.NumBits;\r
+\r
+ int iterations = 0;\r
+\r
+ while (true)\r
+ {\r
+ term.Assign(xn);\r
+ term.mantissa.Sign = true;\r
+ term.Exp(precision);\r
+ term.Mul(this);\r
+ term.Sub(one);\r
+\r
+ iterations++;\r
+ if (term.exponent < -((precision >> 1) - 4))\r
+ {\r
+ if (precision == normalPrecision)\r
+ {\r
+ if (term.exponent < -(precision - 4)) break;\r
+ }\r
+ else\r
+ {\r
+ precision = precision << 1;\r
+ if (precision > normalPrecision) precision = normalPrecision;\r
+ }\r
+ }\r
+\r
+ xn.Add(term);\r
+ }\r
+\r
+ //log(2^n*s) = log(2^n) + log(s) = nlog(2) + log(s)\r
+ term.Assign(ln2cache);\r
+ term.Mul(oldExponent);\r
+\r
+ this.Assign(xn);\r
+ this.Add(term);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Log(x) implemented as an Arithmetic-Geometric Mean. Fast for high precisions.\r
+ /// </summary>\r
+ private void LogAGM1()\r
+ {\r
+ if (mantissa.IsZero() || mantissa.Sign)\r
+ {\r
+ return;\r
+ }\r
+\r
+ //Compute ln2.\r
+ if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)\r
+ {\r
+ CalculateLog2(mantissa.Precision.NumBits);\r
+ }\r
+\r
+ //Compute ln(x) using AGM formula\r
+\r
+ //1. Re-write the input as 2^n * (0.5 <= x < 1)\r
+ int power2 = exponent + 1;\r
+ exponent = -1;\r
+\r
+ //BigFloat res = new BigFloat(firstAGMcache);\r
+ BigFloat a0 = new BigFloat(1, mantissa.Precision);\r
+ BigFloat b0 = new BigFloat(pow10cache);\r
+ b0.Mul(this);\r
+\r
+ BigFloat r = R(a0, b0);\r
+\r
+ this.Assign(firstAGMcache);\r
+ this.Sub(r);\r
+\r
+ a0.Assign(ln2cache);\r
+ a0.Mul(new BigFloat(power2, mantissa.Precision));\r
+ this.Add(a0);\r
+ }\r
+\r
+ private void Exp(int numBits)\r
+ {\r
+ if (IsSpecialValue)\r
+ {\r
+ if (SpecialValue == SpecialValueType.ZERO)\r
+ {\r
+ //e^0 = 1\r
+ exponent = 0;\r
+ mantissa.SetHighDigit(0x80000000);\r
+ }\r
+ else if (SpecialValue == SpecialValueType.INF_MINUS)\r
+ {\r
+ //e^-inf = 0\r
+ SetZero();\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ PrecisionSpec prec = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);\r
+ numBits = prec.NumBits;\r
+\r
+ if (scratch.Precision.NumBits != prec.NumBits)\r
+ {\r
+ scratch = new BigInt(prec);\r
+ }\r
+\r
+ if (inverseFactorialCache == null || invFactorialCutoff < numBits)\r
+ {\r
+ CalculateFactorials(numBits);\r
+ }\r
+\r
+ //let x = 1 * 'this'.mantissa (i.e. 1 <= x < 2)\r
+ //exp(2^n * x) = e^(2^n * x) = (e^x)^2n = exp(x)^2n\r
+\r
+ int oldExponent = 0;\r
+\r
+ if (exponent > -4)\r
+ {\r
+ oldExponent = exponent + 4;\r
+ exponent = -4;\r
+ }\r
+\r
+ BigFloat thisSave = new BigFloat(this, prec);\r
+ BigFloat temp = new BigFloat(1, prec);\r
+ BigFloat temp2 = new BigFloat(this, prec);\r
+ BigFloat res = new BigFloat(1, prec);\r
+ int length = inverseFactorialCache.Length;\r
+\r
+ int iterations;\r
+ for (int i = 1; i < length; i++)\r
+ {\r
+ //temp = x^i\r
+ temp.Mul(thisSave);\r
+ temp2.Assign(inverseFactorialCache[i]);\r
+ temp2.Mul(temp);\r
+\r
+ if (temp2.exponent < -(numBits + 4)) { iterations = i; break; }\r
+\r
+ res.Add(temp2);\r
+ }\r
+\r
+ //res = exp(x)\r
+ //Now... x^(2^n) = (x^2)^(2^(n - 1))\r
+ for (int i = 0; i < oldExponent; i++)\r
+ {\r
+ res.mantissa.SquareHiFast(scratch);\r
+ int shift = res.mantissa.Normalise();\r
+ res.exponent = res.exponent << 1;\r
+ res.exponent += 1 - shift;\r
+ }\r
+\r
+ //Deal with +/- inf\r
+ if (res.exponent == Int32.MaxValue)\r
+ {\r
+ res.mantissa.Zero();\r
+ }\r
+\r
+ Assign(res);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates ln(2) and returns -10^(n/2 + a bit) for reuse, using the AGM method as described in\r
+ /// http://lacim.uqam.ca/~plouffe/articles/log2.pdf\r
+ /// </summary>\r
+ /// <param name="numBits"></param>\r
+ /// <returns></returns>\r
+ private static void CalculateLog2(int numBits)\r
+ {\r
+ //Use the AGM method formula to get log2 to N digits.\r
+ //R(a0, b0) = 1 / (1 - Sum(2^-n*(an^2 - bn^2)))\r
+ //log(1/2) = R(1, 10^-n) - R(1, 10^-n/2)\r
+ PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);\r
+ PrecisionSpec extendedPres = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);\r
+ BigFloat a0 = new BigFloat(1, extendedPres);\r
+ BigFloat b0 = TenPow(-(int)((double)((numBits >> 1) + 2) * 0.302), extendedPres);\r
+ BigFloat pow10saved = new BigFloat(b0);\r
+ BigFloat firstAGMcacheSaved = new BigFloat(extendedPres);\r
+\r
+ //save power of 10 (in normal precision)\r
+ pow10cache = new BigFloat(b0, normalPres);\r
+\r
+ ln2cache = R(a0, b0);\r
+\r
+ //save the first half of the log calculation\r
+ firstAGMcache = new BigFloat(ln2cache, normalPres);\r
+ firstAGMcacheSaved.Assign(ln2cache);\r
+\r
+ b0.MulPow2(-1);\r
+ ln2cache.Sub(R(a0, b0));\r
+\r
+ //Convert to log(2)\r
+ ln2cache.mantissa.Sign = false;\r
+\r
+ //Save magic constant for newton log\r
+ //First guess in range 1 <= x < 2 is x0 = ln2 * (x - 1) + C\r
+ logNewtonConstant = new BigFloat(ln2cache);\r
+ logNewtonConstant.Mul(new BigFloat(3, extendedPres));\r
+ logNewtonConstant.exponent--;\r
+ logNewtonConstant.Sub(new BigFloat(1, extendedPres));\r
+ logNewtonConstant = new BigFloat(logNewtonConstant, normalPres);\r
+\r
+ //Save the inverse.\r
+ log2ecache = new BigFloat(ln2cache);\r
+ log2ecache = new BigFloat(log2ecache.Reciprocal(), normalPres);\r
+\r
+ //Now cache log10\r
+ //Because the log functions call this function to the precision to which they\r
+ //are called, we cannot call them without causing an infinite loop, so we need\r
+ //to inline the code.\r
+ log10recip = new BigFloat(10, extendedPres);\r
+\r
+ {\r
+ int power2 = log10recip.exponent + 1;\r
+ log10recip.exponent = -1;\r
+\r
+ //BigFloat res = new BigFloat(firstAGMcache);\r
+ BigFloat ax = new BigFloat(1, extendedPres);\r
+ BigFloat bx = new BigFloat(pow10saved);\r
+ bx.Mul(log10recip);\r
+\r
+ BigFloat r = R(ax, bx);\r
+\r
+ log10recip.Assign(firstAGMcacheSaved);\r
+ log10recip.Sub(r);\r
+\r
+ ax.Assign(ln2cache);\r
+ ax.Mul(new BigFloat(power2, log10recip.mantissa.Precision));\r
+ log10recip.Add(ax);\r
+ }\r
+\r
+ log10recip = log10recip.Reciprocal();\r
+ log10recip = new BigFloat(log10recip, normalPres);\r
+\r
+\r
+ //Trim to n bits\r
+ ln2cache = new BigFloat(ln2cache, normalPres);\r
+ }\r
+\r
+ private static BigFloat TenPow(int power, PrecisionSpec precision)\r
+ {\r
+ BigFloat acc = new BigFloat(1, precision);\r
+ BigFloat temp = new BigFloat(1, precision);\r
+\r
+ int powerTemp = power;\r
+\r
+ BigFloat multiplierToUse = new BigFloat(10, precision);\r
+\r
+ if (power < 0)\r
+ {\r
+ multiplierToUse = multiplierToUse.Reciprocal();\r
+ powerTemp = -power;\r
+ }\r
+\r
+ //Fast power function\r
+ while (powerTemp != 0)\r
+ {\r
+ temp.Mul(multiplierToUse);\r
+ multiplierToUse.Assign(temp);\r
+\r
+ if ((powerTemp & 1) != 0)\r
+ {\r
+ acc.Mul(temp);\r
+ }\r
+\r
+ powerTemp >>= 1;\r
+ }\r
+\r
+ return acc;\r
+ }\r
+\r
+ private static BigFloat R(BigFloat a0, BigFloat b0)\r
+ {\r
+ //Precision extend taken out.\r
+ int bits = a0.mantissa.Precision.NumBits;\r
+ PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);\r
+ BigFloat an = new BigFloat(a0, extendedPres);\r
+ BigFloat bn = new BigFloat(b0, extendedPres);\r
+ BigFloat sum = new BigFloat(extendedPres);\r
+ BigFloat term = new BigFloat(extendedPres);\r
+ BigFloat temp1 = new BigFloat(extendedPres);\r
+ BigFloat one = new BigFloat(1, extendedPres);\r
+\r
+ int iteration = 0;\r
+\r
+ for (iteration = 0; ; iteration++)\r
+ {\r
+ //Get the sum term for this iteration.\r
+ term.Assign(an);\r
+ term.Mul(an);\r
+ temp1.Assign(bn);\r
+ temp1.Mul(bn);\r
+ //term = an^2 - bn^2\r
+ term.Sub(temp1);\r
+ //term = 2^(n-1) * (an^2 - bn^2)\r
+ term.exponent += iteration - 1;\r
+ sum.Add(term);\r
+\r
+ if (term.exponent < -(bits - 8)) break;\r
+\r
+ //Calculate the new AGM estimates.\r
+ temp1.Assign(an);\r
+ an.Add(bn);\r
+ //a(n+1) = (an + bn) / 2\r
+ an.MulPow2(-1);\r
+\r
+ //b(n+1) = sqrt(an*bn)\r
+ bn.Mul(temp1);\r
+ bn.Sqrt();\r
+ }\r
+\r
+ one.Sub(sum);\r
+ one = one.Reciprocal();\r
+ return new BigFloat(one, a0.mantissa.Precision);\r
+ }\r
+\r
+ private static void CalculateFactorials(int numBits)\r
+ {\r
+ System.Collections.Generic.List<BigFloat> list = new System.Collections.Generic.List<BigFloat>(64);\r
+ System.Collections.Generic.List<BigFloat> list2 = new System.Collections.Generic.List<BigFloat>(64);\r
+\r
+ PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);\r
+ PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);\r
+\r
+ BigFloat factorial = new BigFloat(1, extendedPrecision);\r
+ BigFloat reciprocal;\r
+\r
+ //Calculate e while we're at it\r
+ BigFloat e = new BigFloat(1, extendedPrecision);\r
+\r
+ list.Add(new BigFloat(factorial, normalPrecision));\r
+\r
+ for (int i = 1; i < Int32.MaxValue; i++)\r
+ {\r
+ BigFloat number = new BigFloat(i, extendedPrecision);\r
+ factorial.Mul(number);\r
+\r
+ if (factorial.exponent > numBits) break;\r
+\r
+ list2.Add(new BigFloat(factorial, normalPrecision));\r
+ reciprocal = factorial.Reciprocal();\r
+\r
+ e.Add(reciprocal);\r
+ list.Add(new BigFloat(reciprocal, normalPrecision));\r
+ }\r
+\r
+ //Set the cached static values.\r
+ inverseFactorialCache = list.ToArray();\r
+ factorialCache = list2.ToArray();\r
+ invFactorialCutoff = numBits;\r
+ eCache = new BigFloat(e, normalPrecision);\r
+ eRCPCache = new BigFloat(e.Reciprocal(), normalPrecision);\r
+ }\r
+\r
+ private static void CalculateEOnly(int numBits)\r
+ {\r
+ PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);\r
+ PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);\r
+\r
+ int iExponent = (int)(Math.Sqrt(numBits));\r
+\r
+ BigFloat factorial = new BigFloat(1, extendedPrecision);\r
+ BigFloat constant = new BigFloat(1, extendedPrecision);\r
+ constant.exponent -= iExponent;\r
+ BigFloat numerator = new BigFloat(constant);\r
+ BigFloat reciprocal;\r
+\r
+ //Calculate the 2^iExponent th root of e\r
+ BigFloat e = new BigFloat(1, extendedPrecision);\r
+\r
+ int i;\r
+ for (i = 1; i < Int32.MaxValue; i++)\r
+ {\r
+ BigFloat number = new BigFloat(i, extendedPrecision);\r
+ factorial.Mul(number);\r
+ reciprocal = factorial.Reciprocal();\r
+ reciprocal.Mul(numerator);\r
+\r
+ if (-reciprocal.exponent > numBits) break;\r
+\r
+ e.Add(reciprocal);\r
+ numerator.Mul(constant);\r
+ System.GC.Collect();\r
+ }\r
+\r
+ for (i = 0; i < iExponent; i++)\r
+ {\r
+ numerator.Assign(e);\r
+ e.Mul(numerator);\r
+ }\r
+\r
+ //Set the cached static values.\r
+ eCache = new BigFloat(e, normalPrecision);\r
+ eRCPCache = new BigFloat(e.Reciprocal(), normalPrecision);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Uses the Gauss-Legendre formula for pi\r
+ /// Taken from http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm\r
+ /// </summary>\r
+ /// <param name="numBits"></param>\r
+ private static void CalculatePi(int numBits)\r
+ {\r
+ int bits = numBits + 32;\r
+ //Precision extend taken out.\r
+ PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);\r
+ PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);\r
+\r
+ if (scratch.Precision.NumBits != bits)\r
+ {\r
+ scratch = new BigInt(extendedPres);\r
+ }\r
+\r
+ //a0 = 1\r
+ BigFloat an = new BigFloat(1, extendedPres);\r
+\r
+ //b0 = 1/sqrt(2)\r
+ BigFloat bn = new BigFloat(2, extendedPres);\r
+ bn.Sqrt();\r
+ bn.exponent--;\r
+\r
+ //to = 1/4\r
+ BigFloat tn = new BigFloat(1, extendedPres);\r
+ tn.exponent -= 2;\r
+\r
+ int pn = 0;\r
+\r
+ BigFloat anTemp = new BigFloat(extendedPres);\r
+\r
+ int iteration = 0;\r
+ int cutoffBits = numBits >> 5;\r
+\r
+ for (iteration = 0; ; iteration++)\r
+ {\r
+ //Save a(n)\r
+ anTemp.Assign(an);\r
+\r
+ //Calculate new an\r
+ an.Add(bn);\r
+ an.exponent--;\r
+\r
+ //Calculate new bn\r
+ bn.Mul(anTemp);\r
+ bn.Sqrt();\r
+\r
+ //Calculate new tn\r
+ anTemp.Sub(an);\r
+ anTemp.mantissa.SquareHiFast(scratch);\r
+ anTemp.exponent += anTemp.exponent + pn + 1 - anTemp.mantissa.Normalise();\r
+ tn.Sub(anTemp);\r
+\r
+ anTemp.Assign(an);\r
+ anTemp.Sub(bn);\r
+\r
+ if (anTemp.exponent < -(bits - cutoffBits)) break;\r
+\r
+ //New pn\r
+ pn++;\r
+ }\r
+\r
+ an.Add(bn);\r
+ an.mantissa.SquareHiFast(scratch);\r
+ an.exponent += an.exponent + 1 - an.mantissa.Normalise();\r
+ tn.exponent += 2;\r
+ an.Div(tn);\r
+\r
+ pi = new BigFloat(an, normalPres);\r
+ piBy2 = new BigFloat(pi);\r
+ piBy2.exponent--;\r
+ twoPi = new BigFloat(pi, normalPres);\r
+ twoPi.exponent++;\r
+ piRecip = new BigFloat(an.Reciprocal(), normalPres);\r
+ twoPiRecip = new BigFloat(piRecip);\r
+ twoPiRecip.exponent--;\r
+ //1/3 is going to be useful for sin.\r
+ threeRecip = new BigFloat((new BigFloat(3, extendedPres)).Reciprocal(), normalPres);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates the odd reciprocals of the natural numbers (for atan series)\r
+ /// </summary>\r
+ /// <param name="numBits"></param>\r
+ /// <param name="terms"></param>\r
+ private static void CalculateReciprocals(int numBits, int terms)\r
+ {\r
+ int bits = numBits + 32;\r
+ PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);\r
+ PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);\r
+\r
+ System.Collections.Generic.List<BigFloat> list = new System.Collections.Generic.List<BigFloat>(terms);\r
+\r
+ for (int i = 0; i < terms; i++)\r
+ {\r
+ BigFloat term = new BigFloat(i*2 + 1, extendedPres);\r
+ list.Add(new BigFloat(term.Reciprocal(), normalPres));\r
+ }\r
+\r
+ reciprocals = list.ToArray();\r
+ }\r
+\r
+ /// <summary>\r
+ /// Does decimal rounding, for numbers without E notation.\r
+ /// </summary>\r
+ /// <param name="input"></param>\r
+ /// <param name="places"></param>\r
+ /// <returns></returns>\r
+ private static string RoundString(string input, int places)\r
+ {\r
+ if (places <= 0) return input;\r
+ string trim = input.Trim();\r
+ char[] digits = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'};\r
+\r
+ /*\r
+ for (int i = 1; i <= places; i++)\r
+ {\r
+ //Skip decimal points.\r
+ if (trim[trim.Length - i] == '.')\r
+ {\r
+ places++;\r
+ continue;\r
+ }\r
+\r
+ int index = Array.IndexOf(digits, trim[trim.Length - i]);\r
+\r
+ if (index < 0) return input;\r
+\r
+ value += ten * index;\r
+ ten *= 10;\r
+ }\r
+ * */\r
+\r
+ //Look for a decimal point\r
+ string decimalPoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;\r
+\r
+ int indexPoint = trim.LastIndexOf(decimalPoint);\r
+ if (indexPoint < 0)\r
+ {\r
+ //We can't modify a string which doesn't have a decimal point.\r
+ return trim;\r
+ }\r
+\r
+ int trimPoint = trim.Length - places;\r
+ if (trimPoint < indexPoint) trimPoint = indexPoint;\r
+\r
+ bool roundDown = false;\r
+\r
+ if (trim[trimPoint] == '.')\r
+ {\r
+ if (trimPoint + 1 >= trim.Length)\r
+ {\r
+ roundDown = true;\r
+ }\r
+ else\r
+ {\r
+ int digit = Array.IndexOf(digits, trim[trimPoint + 1]);\r
+ if (digit < 5) roundDown = true;\r
+ }\r
+ }\r
+ else\r
+ {\r
+ int digit = Array.IndexOf(digits, trim[trimPoint]);\r
+ if (digit < 5) roundDown = true;\r
+ }\r
+\r
+ string output;\r
+\r
+ //Round down - just return a new string without the extra digits.\r
+ if (roundDown)\r
+ {\r
+ if (RoundingMode == RoundingModeType.EXACT)\r
+ {\r
+ return trim.Substring(0, trimPoint);\r
+ }\r
+ else\r
+ {\r
+ char[] trimChars = { '0' };\r
+ output = trim.Substring(0, trimPoint).TrimEnd(trimChars);\r
+ trimChars[0] = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator[0];\r
+ return output.TrimEnd(trimChars);\r
+ }\r
+ }\r
+ \r
+ //Round up - bit more complicated.\r
+ char [] arrayOutput = trim.ToCharArray();//0, trimPoint);\r
+\r
+ //Now, we round going from the back to the front.\r
+ int j;\r
+ for (j = trimPoint - 1; j >= 0; j--)\r
+ {\r
+ int index = Array.IndexOf(digits, arrayOutput[j]);\r
+\r
+ //Skip decimal points etc...\r
+ if (index < 0) continue;\r
+\r
+ if (index < 9)\r
+ {\r
+ arrayOutput[j] = digits[index + 1];\r
+ break;\r
+ }\r
+ else\r
+ {\r
+ arrayOutput[j] = digits[0];\r
+ }\r
+ }\r
+\r
+ output = new string(arrayOutput);\r
+\r
+ if (j < 0)\r
+ {\r
+ //Need to add a new digit.\r
+ output = String.Format("{0}{1}", "1", output);\r
+ }\r
+\r
+ if (RoundingMode == RoundingModeType.EXACT)\r
+ {\r
+ return output;\r
+ }\r
+ else\r
+ {\r
+ char[] trimChars = { '0' };\r
+ output = output.TrimEnd(trimChars);\r
+ trimChars[0] = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator[0];\r
+ return output.TrimEnd(trimChars);\r
+ }\r
+ }\r
+\r
+ //***************************** Data *****************************\r
+\r
+\r
+ //Side node - this way of doing things is far from optimal, both in terms of memory use and performance.\r
+ private ExponentAdaptor exponent;\r
+ private BigInt mantissa;\r
+\r
+ /// <summary>\r
+ /// Storage area for calculations.\r
+ /// </summary>\r
+ private static BigInt scratch;\r
+\r
+ private static BigFloat ln2cache; //Value of ln(2)\r
+ private static BigFloat log2ecache; //Value of log2(e) = 1/ln(2)\r
+ private static BigFloat pow10cache; //Cached power of 10 for AGM log calculation\r
+ private static BigFloat log10recip; //1/ln(10)\r
+ private static BigFloat firstAGMcache; //Cached half of AGM operation.\r
+ private static BigFloat[] factorialCache; //The values of n!\r
+ private static BigFloat[] inverseFactorialCache; //Values of 1/n! up to 2^-m where m = invFactorialCutoff (below)\r
+ private static int invFactorialCutoff; //The number of significant bits for the cutoff of the inverse factorials.\r
+ private static BigFloat eCache; //Value of e cached to invFactorialCutoff bits\r
+ private static BigFloat eRCPCache; //Reciprocal of e\r
+ private static BigFloat logNewtonConstant; //1.5*ln(2) - 1\r
+ private static BigFloat pi; //pi\r
+ private static BigFloat piBy2; //pi/2\r
+ private static BigFloat twoPi; //2*pi\r
+ private static BigFloat piRecip; //1/pi\r
+ private static BigFloat twoPiRecip; //1/2*pi\r
+ private static BigFloat threeRecip; //1/3\r
+ private static BigFloat[] reciprocals; //1/x\r
+ \r
+ /// <summary>\r
+ /// The number of decimal digits to round the output of ToString() by\r
+ /// </summary>\r
+ public static int RoundingDigits { get; set; }\r
+\r
+ /// <summary>\r
+ /// The way in which ToString() should deal with insignificant trailing zeroes\r
+ /// </summary>\r
+ public static RoundingModeType RoundingMode { get; set; }\r
+ }\r
+}
\ No newline at end of file
--- /dev/null
+// http://www.fractal-landscapes.co.uk/bigint.html\r
+\r
+using System;\r
+\r
+namespace BigNum\r
+{\r
+ /// <summary>\r
+ /// Specifies the desired precision for a BigInt or a BigFloat. \r
+ /// </summary>\r
+ public struct PrecisionSpec\r
+ {\r
+ /// <summary>\r
+ /// Precision can be specified in a choice of 8 bases.\r
+ /// Note that precision for decimals is approximate.\r
+ /// </summary>\r
+ public enum BaseType\r
+ {\r
+ /// <summary>\r
+ /// Binary base\r
+ /// </summary>\r
+ BIN,\r
+ /// <summary>\r
+ /// Octal base\r
+ /// </summary>\r
+ OCT,\r
+ /// <summary>\r
+ /// Decimal base\r
+ /// </summary>\r
+ DEC,\r
+ /// <summary>\r
+ /// Hexadecimal base\r
+ /// </summary>\r
+ HEX,\r
+ /// <summary>\r
+ /// 8-bits per digit\r
+ /// </summary>\r
+ BYTES,\r
+ /// <summary>\r
+ /// 16-bits per digit\r
+ /// </summary>\r
+ WORDS,\r
+ /// <summary>\r
+ /// 32-bits per digit\r
+ /// </summary>\r
+ DWORDS,\r
+ /// <summary>\r
+ /// 64-bits per digit\r
+ /// </summary>\r
+ QWORDS\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructor: Constructs a precision specification\r
+ /// </summary>\r
+ /// <param name="precision">The number of digits</param>\r
+ /// <param name="numberBase">The base of the digits</param>\r
+ public PrecisionSpec(int precision, BaseType numberBase)\r
+ {\r
+ this.prec = precision;\r
+ this.nB = numberBase;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Explicit cast from integer value.\r
+ /// </summary>\r
+ /// <param name="value">The value in bits for the new precision specification</param>\r
+ /// <returns>A new precision specification with the number of bits specified</returns>\r
+ public static explicit operator PrecisionSpec(int value)\r
+ {\r
+ return new PrecisionSpec(value, BaseType.BIN);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Equality test\r
+ /// </summary>\r
+ /// <param name="spec1">the first parameter</param>\r
+ /// <param name="spec2">the second parameter</param>\r
+ /// <returns>true iff both precisions have the same number of bits</returns>\r
+ public static bool operator ==(PrecisionSpec spec1, PrecisionSpec spec2)\r
+ {\r
+ return (spec1.NumBits == spec2.NumBits);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Inequality operator\r
+ /// </summary>\r
+ /// <param name="spec1">the first parameter</param>\r
+ /// <param name="spec2">the second parameter</param>\r
+ /// <returns>true iff the parameters do not have the same number of bits</returns>\r
+ public static bool operator !=(PrecisionSpec spec1, PrecisionSpec spec2)\r
+ {\r
+ return !(spec1 == spec2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Object equality override\r
+ /// </summary>\r
+ /// <param name="obj">the PrecisionSpec struct to compare</param>\r
+ /// <returns>true iff obj has the same number of bits as this</returns>\r
+ public override bool Equals(object obj)\r
+ {\r
+ return NumBits == ((PrecisionSpec)obj).NumBits; \r
+ }\r
+\r
+ /// <summary>\r
+ /// Override of the hash code\r
+ /// </summary>\r
+ /// <returns>A basic hash</returns>\r
+ public override int GetHashCode()\r
+ {\r
+ return NumBits * prec + NumBits;\r
+ }\r
+\r
+ /// <summary>\r
+ /// The precision in units specified by the base type (e.g. number of decimal digits)\r
+ /// </summary>\r
+ public int Precision\r
+ {\r
+ get { return prec; }\r
+ }\r
+\r
+ /// <summary>\r
+ /// The base type in which precision is specified\r
+ /// </summary>\r
+ public BaseType NumberBaseType\r
+ {\r
+ get { return nB; }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Converts the number-base to an integer\r
+ /// </summary>\r
+ public int NumberBase\r
+ {\r
+ get { return (int)nB; }\r
+ }\r
+\r
+ /// <summary>\r
+ /// The number of bits that this PrecisionSpec structure specifies.\r
+ /// </summary>\r
+ public int NumBits\r
+ {\r
+ get\r
+ {\r
+ if (nB == BaseType.BIN) return prec;\r
+ if (nB == BaseType.OCT) return prec * 3;\r
+ if (nB == BaseType.HEX) return prec * 4;\r
+ if (nB == BaseType.BYTES) return prec * 8;\r
+ if (nB == BaseType.WORDS) return prec * 16;\r
+ if (nB == BaseType.DWORDS) return prec * 32;\r
+ if (nB == BaseType.QWORDS) return prec * 64;\r
+\r
+ double factor = 3.322;\r
+ int bits = ((int)Math.Ceiling(factor * (double)prec));\r
+ return bits;\r
+ }\r
+ }\r
+\r
+ private int prec;\r
+ private BaseType nB;\r
+ }\r
+\r
+ \r
+ /// <summary>\r
+ /// An arbitrary-precision integer class\r
+ /// \r
+ /// Format:\r
+ /// Each number consists of an array of 32-bit unsigned integers, and a bool sign\r
+ /// value.\r
+ /// \r
+ /// Applicability and Performance:\r
+ /// This class is designed to be used for small extended precisions. It may not be\r
+ /// safe (and certainly won't be fast) to use it with mixed-precision arguments.\r
+ /// It does support, but will not be efficient for, numbers over around 2048 bits.\r
+ /// \r
+ /// Notes:\r
+ /// All conversions to and from strings are slow.\r
+ /// \r
+ /// Conversions from simple integer types Int32, Int64, UInt32, UInt64 are performed\r
+ /// using the appropriate constructor, and are relatively fast.\r
+ /// \r
+ /// The class is written entirely in managed C# code, with not native or managed\r
+ /// assembler. The use of native assembler would speed up the multiplication operations\r
+ /// many times over, and therefore all higher-order operations too.\r
+ /// </summary>\r
+ public class BigInt\r
+ {\r
+ //*************** Constructors ******************\r
+\r
+ /// <summary>\r
+ /// Constructs an empty BigInt to the desired precision.\r
+ /// </summary>\r
+ /// <param name="precision"></param>\r
+ public BigInt(PrecisionSpec precision)\r
+ {\r
+ Init(precision);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigInt from a string, using the string length to determine the precision\r
+ /// Note, operations on BigInts of non-matching precision are slow, so avoid using this constructor\r
+ /// </summary>\r
+ /// <param name="init"></param>\r
+ public BigInt(string init)\r
+ {\r
+ InitFromString(init, (PrecisionSpec)init.Length, 10);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructor for copying length and precision\r
+ /// </summary>\r
+ /// <param name="inputToCopy">The BigInt to copy</param>\r
+ /// <param name="precision">The precision of the new BigInt</param>\r
+ /// <param name="bCopyLengthOnly">decides whether to copy the actual input, or just its digit length</param>\r
+ /// <example><code>//Create an integer\r
+ /// BigInt four = new BigInt(4, new PrecisionSpec(128, PrecisionSpec.BaseType.BIN));\r
+ /// \r
+ /// //Pad four to double its usual number of digits (this does not affect the precision)\r
+ /// four.Pad();\r
+ /// \r
+ /// //Create a new, empty integer with matching precision, also padded to twice the usual length\r
+ /// BigInt newCopy = new BigInt(four, four.Precision, true);</code></example>\r
+ public BigInt(BigInt inputToCopy, PrecisionSpec precision, bool bCopyLengthOnly)\r
+ {\r
+ digitArray = new uint[inputToCopy.digitArray.Length];\r
+ workingSet = new uint[inputToCopy.digitArray.Length];\r
+ if (!bCopyLengthOnly) Array.Copy(inputToCopy.digitArray, digitArray, digitArray.Length);\r
+ sign = inputToCopy.sign;\r
+ pres = inputToCopy.pres;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a bigint from the string, with the desired precision, using base 10\r
+ /// </summary>\r
+ /// <param name="init"></param>\r
+ /// <param name="precision"></param>\r
+ public BigInt(string init, PrecisionSpec precision)\r
+ {\r
+ InitFromString(init, precision, 10);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigInt from a string, using the specified precision and base\r
+ /// </summary>\r
+ /// <param name="init"></param>\r
+ /// <param name="precision"></param>\r
+ /// <param name="numberBase"></param>\r
+ public BigInt(string init, PrecisionSpec precision, int numberBase)\r
+ {\r
+ InitFromString(init, precision, numberBase);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a bigint from the input.\r
+ /// </summary>\r
+ /// <param name="input"></param>\r
+ public BigInt(BigInt input)\r
+ {\r
+ digitArray = new uint[input.digitArray.Length];\r
+ workingSet = new uint[input.digitArray.Length];\r
+ Array.Copy(input.digitArray, digitArray, digitArray.Length);\r
+ sign = input.sign;\r
+ pres = input.pres;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a bigint from the input, matching the new precision provided\r
+ /// </summary>\r
+ public BigInt(BigInt input, PrecisionSpec precision)\r
+ {\r
+ //Casts the input to the new precision.\r
+ Init(precision);\r
+ int Min = (input.digitArray.Length < digitArray.Length) ? input.digitArray.Length : digitArray.Length;\r
+\r
+ for (int i = 0; i < Min; i++)\r
+ {\r
+ digitArray[i] = input.digitArray[i];\r
+ }\r
+\r
+ sign = input.sign;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigInt from a UInt32\r
+ /// </summary>\r
+ /// <param name="input"></param>\r
+ /// <param name="precision"></param>\r
+ public BigInt(UInt32 input, PrecisionSpec precision)\r
+ {\r
+ Init(precision);\r
+ digitArray[0] = input;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigInt from a UInt64\r
+ /// </summary>\r
+ /// <param name="input"></param>\r
+ /// <param name="precision"></param>\r
+ public BigInt(UInt64 input, PrecisionSpec precision)\r
+ {\r
+ Init(precision);\r
+ digitArray[0] = (UInt32)(input & 0xffffffff);\r
+ if (digitArray.Length > 1) digitArray[1] = (UInt32)(input >> 32);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigInt from an Int32\r
+ /// </summary>\r
+ /// <param name="input"></param>\r
+ /// <param name="precision"></param>\r
+ public BigInt(Int32 input, PrecisionSpec precision)\r
+ {\r
+ Init(precision);\r
+ if (input < 0)\r
+ {\r
+ sign = true;\r
+\r
+ if (input == Int32.MinValue)\r
+ {\r
+ digitArray[0] = 0x80000000;\r
+ }\r
+ else\r
+ {\r
+ digitArray[0] = (UInt32)(-input);\r
+ }\r
+ }\r
+ else\r
+ {\r
+ digitArray[0] = ((UInt32)input);\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Constructs a BigInt from a UInt32\r
+ /// </summary>\r
+ /// <param name="input"></param>\r
+ /// <param name="precision"></param>\r
+ public BigInt(Int64 input, PrecisionSpec precision)\r
+ {\r
+ Init(precision);\r
+ if (input < 0) sign = true;\r
+\r
+ digitArray[0] = (UInt32)(input & 0xffffffff);\r
+\r
+ if (digitArray.Length >= 2)\r
+ {\r
+ if (input == Int64.MinValue)\r
+ {\r
+ digitArray[1] = 0x80000000;\r
+ }\r
+ else\r
+ {\r
+ digitArray[1] = (UInt32)((input >> 32) & 0x7fffffff);\r
+ }\r
+ }\r
+ }\r
+\r
+ //***************** Properties *******************\r
+\r
+ /// <summary>\r
+ /// true iff the number is negative\r
+ /// </summary>\r
+ public bool Sign { get { return sign; } set { sign = value; } }\r
+\r
+ /// <summary>\r
+ /// The precision of the number.\r
+ /// </summary>\r
+ public PrecisionSpec Precision { get { return pres; } }\r
+\r
+ //*************** Utility Functions **************\r
+\r
+ /// <summary>\r
+ /// Casts a BigInt to the new precision provided.\r
+ /// Note: This will return the input if the precision already matches.\r
+ /// </summary>\r
+ /// <param name="input"></param>\r
+ /// <param name="precision"></param>\r
+ /// <returns></returns>\r
+ public static BigInt CastToPrecision(BigInt input, PrecisionSpec precision)\r
+ {\r
+ if (input.pres == precision) return input;\r
+ return new BigInt(input, precision);\r
+ }\r
+\r
+\r
+ //*************** Member Functions ***************\r
+\r
+ /// <summary>\r
+ /// Addition and assignment - without intermediate memory allocation.\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public uint Add(BigInt n2)\r
+ {\r
+ if (n2.digitArray.Length != digitArray.Length) MakeSafe(ref n2);\r
+\r
+ if (sign == n2.sign)\r
+ {\r
+ return AddInternalBits(n2.digitArray);\r
+ }\r
+ else\r
+ {\r
+ bool lessThan = LtInt(this, n2);\r
+\r
+ if (lessThan)\r
+ {\r
+ int Length = digitArray.Length;\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ workingSet[i] = digitArray[i];\r
+ digitArray[i] = n2.digitArray[i];\r
+ }\r
+\r
+ sign = !sign;\r
+ return SubInternalBits(workingSet);\r
+ }\r
+ else\r
+ {\r
+ return SubInternalBits(n2.digitArray);\r
+ }\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Subtraction and assignment - without intermediate memory allocation.\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public uint Sub(BigInt n2)\r
+ {\r
+ if (n2.digitArray.Length != digitArray.Length) MakeSafe(ref n2);\r
+\r
+ if (sign != n2.sign)\r
+ {\r
+ return AddInternalBits(n2.digitArray);\r
+ }\r
+ else\r
+ {\r
+ bool lessThan = LtInt(this, n2);\r
+\r
+ if (lessThan)\r
+ {\r
+ int Length = digitArray.Length;\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ workingSet[i] = digitArray[i];\r
+ digitArray[i] = n2.digitArray[i];\r
+ }\r
+\r
+ sign = !sign;\r
+ return SubInternalBits(workingSet);\r
+ }\r
+ else\r
+ {\r
+ return SubInternalBits(n2.digitArray);\r
+ }\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Multiplication and assignmnet - with minimal intermediate memory allocation.\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public void Mul(BigInt n2)\r
+ {\r
+ if (n2.digitArray.Length != digitArray.Length) MakeSafe(ref n2);\r
+\r
+ int Length = n2.digitArray.Length;\r
+\r
+ //Inner loop zero-mul avoidance\r
+ int maxDigit = 0;\r
+ for (int i = Length - 1; i >= 0; i--)\r
+ {\r
+ if (digitArray[i] != 0)\r
+ {\r
+ maxDigit = i + 1;\r
+ break;\r
+ }\r
+ }\r
+\r
+ //Result is zero, 'this' is unchanged\r
+ if (maxDigit == 0) return;\r
+\r
+ //Temp storage for source (both working sets are used by the calculation)\r
+ uint[] thisTemp = new uint [Length];\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ thisTemp[i] = digitArray[i];\r
+ digitArray[i] = 0;\r
+ }\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ //Clear the working set\r
+ for (int j = 0; j < i; j++)\r
+ {\r
+ workingSet[j] = 0;\r
+ n2.workingSet[j] = 0;\r
+ }\r
+\r
+ n2.workingSet[i] = 0;\r
+\r
+ ulong digit = n2.digitArray[i];\r
+ if (digit == 0) continue;\r
+\r
+ for (int j = 0; j + i < Length && j < maxDigit; j++)\r
+ {\r
+ //Multiply n1 by each of the integer digits of n2.\r
+ ulong temp = (ulong)thisTemp[j] * digit;\r
+ //n1.workingSet stores the low bits of each piecewise multiplication\r
+ workingSet[j + i] = (uint)(temp & 0xffffffff);\r
+ if (j + i + 1 < Length)\r
+ {\r
+ //n2.workingSet stores the high bits of each multiplication\r
+ n2.workingSet[j + i + 1] = (uint)(temp >> 32);\r
+ }\r
+ }\r
+\r
+ AddInternalBits(workingSet);\r
+ AddInternalBits(n2.workingSet);\r
+ }\r
+\r
+ sign = (sign != n2.sign);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Squares the number.\r
+ /// </summary>\r
+ public void Square()\r
+ {\r
+ int Length = digitArray.Length;\r
+\r
+ //Inner loop zero-mul avoidance\r
+ int maxDigit = 0;\r
+ for (int i = Length - 1; i >= 0; i--)\r
+ {\r
+ if (digitArray[i] != 0)\r
+ {\r
+ maxDigit = i + 1;\r
+ break;\r
+ }\r
+ }\r
+\r
+ //Result is zero, 'this' is unchanged\r
+ if (maxDigit == 0) return;\r
+\r
+ //Temp storage for source (both working sets are used by the calculation)\r
+ uint[] thisTemp = new uint[Length];\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ thisTemp[i] = digitArray[i];\r
+ digitArray[i] = 0;\r
+ }\r
+\r
+ UInt32 [] workingSet2 = new UInt32[Length];\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ //Clear the working set\r
+ for (int j = 0; j < i; j++)\r
+ {\r
+ workingSet[j] = 0;\r
+ workingSet2[j] = 0;\r
+ }\r
+\r
+ workingSet2[i] = 0;\r
+\r
+ ulong digit = thisTemp[i];\r
+ if (digit == 0) continue;\r
+\r
+ for (int j = 0; j + i < Length && j < maxDigit; j++)\r
+ {\r
+ //Multiply n1 by each of the integer digits of n2.\r
+ ulong temp = (ulong)thisTemp[j] * digit;\r
+ //n1.workingSet stores the low bits of each piecewise multiplication\r
+ workingSet[j + i] = (uint)(temp & 0xffffffff);\r
+ if (j + i + 1 < Length)\r
+ {\r
+ //n2.workingSet stores the high bits of each multiplication\r
+ workingSet2[j + i + 1] = (uint)(temp >> 32);\r
+ }\r
+ }\r
+\r
+ AddInternalBits(workingSet);\r
+ AddInternalBits(workingSet2);\r
+ }\r
+\r
+ sign = false;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Used for floating-point multiplication\r
+ /// Stores the high bits of the multiplication only (the carry bit from the\r
+ /// lower bits is missing, so the true answer might be 1 greater).\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public void MulHi(BigInt n2)\r
+ {\r
+ if (n2.digitArray.Length != digitArray.Length) MakeSafe(ref n2);\r
+\r
+ int Length = n2.digitArray.Length;\r
+\r
+ //Inner loop zero-mul avoidance\r
+ int maxDigit = 0;\r
+ for (int i = Length - 1; i >= 0; i--)\r
+ {\r
+ if (digitArray[i] != 0)\r
+ {\r
+ maxDigit = i + 1;\r
+ break;\r
+ }\r
+ }\r
+\r
+ //Result is zero, 'this' is unchanged\r
+ if (maxDigit == 0) return;\r
+\r
+ //Temp storage for source (both working sets are used by the calculation)\r
+ uint[] thisTemp = new uint[Length];\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ thisTemp[i] = digitArray[i];\r
+ digitArray[i] = 0;\r
+ }\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ //Clear the working set\r
+ for (int j = 0; j < Length; j++)\r
+ {\r
+ workingSet[j] = 0;\r
+ n2.workingSet[j] = 0;\r
+ }\r
+\r
+ n2.workingSet[i] = 0;\r
+\r
+ ulong digit = n2.digitArray[i];\r
+ if (digit == 0) continue;\r
+\r
+ //Only the high bits\r
+ if (maxDigit + i < Length - 1) continue;\r
+\r
+ for (int j = 0; j < maxDigit; j++)\r
+ {\r
+ if (j + i + 1 < Length) continue;\r
+ //Multiply n1 by each of the integer digits of n2.\r
+ ulong temp = (ulong)thisTemp[j] * digit;\r
+ //n1.workingSet stores the low bits of each piecewise multiplication\r
+ if (j + i >= Length)\r
+ {\r
+ workingSet[j + i - Length] = (uint)(temp & 0xffffffff);\r
+ }\r
+ \r
+ //n2.workingSet stores the high bits of each multiplication\r
+ n2.workingSet[j + i + 1 - Length] = (uint)(temp >> 32);\r
+ }\r
+\r
+ AddInternalBits(workingSet);\r
+ AddInternalBits(n2.workingSet);\r
+ }\r
+\r
+ sign = (sign != n2.sign);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Floating-point helper function.\r
+ /// Squares the number and keeps the high bits of the calculation.\r
+ /// Takes a temporary BigInt as a working set.\r
+ /// </summary>\r
+ public void SquareHiFast(BigInt scratch)\r
+ {\r
+ int Length = digitArray.Length;\r
+ uint[] tempDigits = scratch.digitArray;\r
+ uint[] workingSet2 = scratch.workingSet;\r
+\r
+ //Temp storage for source (both working sets are used by the calculation)\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ tempDigits[i] = digitArray[i];\r
+ digitArray[i] = 0;\r
+ }\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ //Clear the working set\r
+ for (int j = i; j < Length; j++)\r
+ {\r
+ workingSet[j] = 0;\r
+ workingSet2[j] = 0;\r
+ }\r
+\r
+ if (i - 1 >= 0) workingSet[i - 1] = 0;\r
+\r
+ ulong digit = tempDigits[i];\r
+ if (digit == 0) continue;\r
+\r
+ for (int j = 0; j < Length; j++)\r
+ {\r
+ if (j + i + 1 < Length) continue;\r
+ //Multiply n1 by each of the integer digits of n2.\r
+ ulong temp = (ulong)tempDigits[j] * digit;\r
+ //n1.workingSet stores the low bits of each piecewise multiplication\r
+ if (j + i >= Length)\r
+ {\r
+ workingSet[j + i - Length] = (uint)(temp & 0xffffffff);\r
+ }\r
+\r
+ //n2.workingSet stores the high bits of each multiplication\r
+ workingSet2[j + i + 1 - Length] = (uint)(temp >> 32);\r
+ }\r
+\r
+ AddInternalBits(workingSet);\r
+ AddInternalBits(workingSet2);\r
+ }\r
+\r
+ sign = false;\r
+ }\r
+\r
+ /// <summary>\r
+ /// This uses the schoolbook division algorithm, as decribed by http://www.treskal.com/kalle/exjobb/original-report.pdf\r
+ /// Algorithms 3.1 (implemented by Div_31) and 3.2 (implemented by Div_32)\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public void Div(BigInt n2)\r
+ {\r
+ if (n2.digitArray.Length != digitArray.Length) MakeSafe(ref n2);\r
+\r
+ int OldLength = digitArray.Length;\r
+\r
+ //First, we need to prepare the operands for division using Div_32, which requires\r
+ //That the most significant digit of n2 be set. To do this, we need to shift n2 (and therefore n1) up.\r
+ //This operation can potentially increase the precision of the operands.\r
+ int shift = MakeSafeDiv(this, n2);\r
+\r
+ BigInt Q = new BigInt(this, this.pres, true);\r
+ BigInt R = new BigInt(this, this.pres, true);\r
+\r
+ Div_32(this, n2, Q, R);\r
+\r
+ //Restore n2 to its pre-shift value\r
+ n2.RSH(shift);\r
+ AssignInt(Q);\r
+ sign = (sign != n2.sign);\r
+\r
+ //Reset the lengths of the operands\r
+ SetNumDigits(OldLength);\r
+ n2.SetNumDigits(OldLength);\r
+ }\r
+\r
+ /// <summary>\r
+ /// This function is used for floating-point division.\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ //Given two numbers:\r
+ // In floating point 1 <= a, b < 2, meaning that both numbers have their top bits set.\r
+ // To calculate a / b, maintaining precision, we:\r
+ // 1. Double the number of digits available to both numbers.\r
+ // 2. set a = a * 2^d (where d is the number of digits)\r
+ // 3. calculate the quotient a <- q: 2^(d-1) <= q < 2^(d+1)\r
+ // 4. if a >= 2^d, s = 1, else s = 0\r
+ // 6. shift a down by s, and undo the precision extension\r
+ // 7. return 1 - shift (change necessary to exponent)\r
+ public int DivAndShift(BigInt n2)\r
+ {\r
+ if (n2.IsZero()) return -1;\r
+ if (digitArray.Length != n2.digitArray.Length) MakeSafe(ref n2);\r
+\r
+ int oldLength = digitArray.Length;\r
+\r
+ //Double the number of digits, and shift a into the higher digits.\r
+ Pad();\r
+ n2.Extend();\r
+\r
+ //Do the divide (at double precision, ouch!)\r
+ Div(n2);\r
+\r
+ //Shift down if 'this' >= 2^d\r
+ int ret = 1;\r
+\r
+ if (digitArray[oldLength] != 0)\r
+ {\r
+ RSH(1);\r
+ ret--;\r
+ }\r
+\r
+ SetNumDigits(oldLength);\r
+ n2.SetNumDigits(oldLength);\r
+\r
+ return ret;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates 'this' mod n2 (using the schoolbook division algorithm as above)\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public void Mod(BigInt n2)\r
+ {\r
+ if (n2.digitArray.Length != digitArray.Length) MakeSafe(ref n2);\r
+\r
+ int OldLength = digitArray.Length;\r
+\r
+ //First, we need to prepare the operands for division using Div_32, which requires\r
+ //That the most significant digit of n2 be set. To do this, we need to shift n2 (and therefore n1) up.\r
+ //This operation can potentially increase the precision of the operands.\r
+ int shift = MakeSafeDiv(this, n2);\r
+\r
+ BigInt Q = new BigInt(this.pres);\r
+ BigInt R = new BigInt(this.pres);\r
+\r
+ Q.digitArray = new UInt32[this.digitArray.Length];\r
+ R.digitArray = new UInt32[this.digitArray.Length];\r
+\r
+ Div_32(this, n2, Q, R);\r
+\r
+ //Restore n2 to its pre-shift value\r
+ n2.RSH(shift);\r
+ R.RSH(shift);\r
+ R.sign = (sign != n2.sign);\r
+ AssignInt(R);\r
+\r
+ //Reset the lengths of the operands\r
+ SetNumDigits(OldLength);\r
+ n2.SetNumDigits(OldLength);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Logical left shift\r
+ /// </summary>\r
+ /// <param name="shift"></param>\r
+ public void LSH(int shift)\r
+ {\r
+ if (shift <= 0) return;\r
+ int length = digitArray.Length;\r
+ int digits = shift >> 5;\r
+ int rem = shift & 31;\r
+\r
+ for (int i = length - 1; i >= digits; i--)\r
+ {\r
+ digitArray[i] = digitArray[i - digits];\r
+ }\r
+\r
+ if (digits > 0)\r
+ {\r
+ for (int i = digits - 1; i >= 0; i--)\r
+ {\r
+ digitArray[i] = 0;\r
+ }\r
+ }\r
+\r
+ UInt64 lastShift = 0;\r
+\r
+ for (int i = 0; i < length; i++)\r
+ {\r
+ UInt64 temp = (((UInt64)digitArray[i]) << rem) | lastShift;\r
+ digitArray[i] = (UInt32)(temp & 0xffffffff);\r
+ lastShift = temp >> 32;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Logical right-shift\r
+ /// </summary>\r
+ /// <param name="shift"></param>\r
+ public void RSH(int shift)\r
+ {\r
+ if (shift < 0) return;\r
+ int length = digitArray.Length;\r
+ int digits = shift >> 5;\r
+ int rem = shift & 31;\r
+\r
+ for (int i = 0; i < length - digits; i++)\r
+ {\r
+ digitArray[i] = digitArray[i + digits];\r
+ }\r
+\r
+ int start = (length - digits);\r
+ if (start < 0) start = 0;\r
+\r
+ for (int i = start; i < length; i++)\r
+ {\r
+ digitArray[i] = 0;\r
+ }\r
+\r
+ UInt64 lastShift = 0;\r
+\r
+ for (int i = length - 1; i >= 0; i--)\r
+ {\r
+ UInt64 temp = ((((UInt64)digitArray[i]) << 32) >> rem) | lastShift;\r
+ digitArray[i] = (UInt32)(temp >> 32);\r
+ lastShift = (temp & 0xffffffff) << 32;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Changes the sign of the number\r
+ /// </summary>\r
+ public void Negate()\r
+ {\r
+ sign = !sign;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Increments the number.\r
+ /// </summary>\r
+ public void Increment()\r
+ {\r
+ if (sign)\r
+ {\r
+ DecrementInt();\r
+ }\r
+ else\r
+ {\r
+ IncrementInt();\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Decrements the number.\r
+ /// </summary>\r
+ public void Decrement()\r
+ {\r
+ if (sign)\r
+ {\r
+ IncrementInt();\r
+ }\r
+ else\r
+ {\r
+ DecrementInt();\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates the factorial 'this'!\r
+ /// </summary>\r
+ public void Factorial()\r
+ {\r
+ if (sign) return;\r
+\r
+ //Clamp to a reasonable range.\r
+ int factToUse = (int)(digitArray[0]);\r
+ if (factToUse > 65536) factToUse = 65536;\r
+\r
+ Zero();\r
+ digitArray[0] = 1;\r
+\r
+ for (uint i = 1; i <= factToUse; i++)\r
+ {\r
+ MulInternal(i);\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Calculates 'this'^power\r
+ /// </summary>\r
+ /// <param name="power"></param>\r
+ public void Power(BigInt power)\r
+ {\r
+ if (power.IsZero() || power.sign)\r
+ {\r
+ Zero();\r
+ digitArray[0] = 1;\r
+ return;\r
+ }\r
+\r
+ BigInt pow = new BigInt(power);\r
+ BigInt temp = new BigInt(this);\r
+ BigInt powTerm = new BigInt(this);\r
+\r
+ pow.Decrement();\r
+ for (; !pow.IsZero(); pow.RSH(1))\r
+ {\r
+ if ((pow.digitArray[0] & 1) == 1)\r
+ {\r
+ temp.Mul(powTerm);\r
+ }\r
+\r
+ powTerm.Square();\r
+ }\r
+\r
+ Assign(temp);\r
+ }\r
+\r
+ //***************** Comparison member functions *****************\r
+\r
+ /// <summary>\r
+ /// returns true if this bigint == 0\r
+ /// </summary>\r
+ /// <returns></returns>\r
+ public bool IsZero()\r
+ {\r
+ for (int i = 0; i < digitArray.Length; i++)\r
+ {\r
+ if (digitArray[i] != 0) return false;\r
+ }\r
+\r
+ return true;\r
+ }\r
+\r
+ /// <summary>\r
+ /// true iff n2 (precision adjusted to this) is less than 'this'\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public bool LessThan(BigInt n2)\r
+ {\r
+ if (digitArray.Length != n2.digitArray.Length) MakeSafe(ref n2);\r
+\r
+ if (sign)\r
+ {\r
+ if (!n2.sign) return true;\r
+ return GtInt(this, n2);\r
+ }\r
+ else\r
+ {\r
+ if (n2.sign) return false;\r
+ return LtInt(this, n2);\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// true iff n2 (precision adjusted to this) is greater than 'this'\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public bool GreaterThan(BigInt n2)\r
+ {\r
+ if (digitArray.Length != n2.digitArray.Length) MakeSafe(ref n2);\r
+\r
+ if (sign)\r
+ {\r
+ if (!n2.sign) return false;\r
+ return LtInt(this, n2);\r
+ }\r
+ else\r
+ {\r
+ if (n2.sign) return true;\r
+ return GtInt(this, n2);\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Override of base-class equals\r
+ /// </summary>\r
+ /// <param name="obj"></param>\r
+ /// <returns></returns>\r
+ public override bool Equals(object obj)\r
+ {\r
+ BigInt n2 = ((BigInt)obj);\r
+ return Equals(n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Get hash code\r
+ /// </summary>\r
+ /// <returns></returns>\r
+ public override int GetHashCode()\r
+ {\r
+ return (Int32)digitArray[0];\r
+ }\r
+\r
+ /// <summary>\r
+ /// True iff n2 (precision-adjusted to this) == n1\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public bool Equals(BigInt n2)\r
+ {\r
+ if (IsZero() && n2.IsZero()) return true;\r
+\r
+ if (sign != n2.sign) return false;\r
+\r
+ int Length = digitArray.Length;\r
+ if (n2.digitArray.Length != Length) MakeSafe(ref n2);\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ if (digitArray[i] != n2.digitArray[i]) return false;\r
+ }\r
+\r
+ return true;\r
+ }\r
+\r
+ //******************* Bitwise member functions ********************\r
+\r
+ /// <summary>\r
+ /// Takes the bitwise complement of the number\r
+ /// </summary>\r
+ public void Complement()\r
+ {\r
+ int Length = digitArray.Length;\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ digitArray[i] = ~digitArray[i];\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Bitwise And\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public void And(BigInt n2)\r
+ {\r
+ int Length = digitArray.Length;\r
+ if (n2.digitArray.Length != Length) MakeSafe(ref n2);\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ digitArray[i] &= n2.digitArray[i];\r
+ }\r
+\r
+ sign &= n2.sign;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Bitwise Or\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public void Or(BigInt n2)\r
+ {\r
+ int Length = digitArray.Length;\r
+ if (n2.digitArray.Length != Length) MakeSafe(ref n2);\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ digitArray[i] |= n2.digitArray[i];\r
+ }\r
+\r
+ sign |= n2.sign;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Bitwise Xor\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public void Xor(BigInt n2)\r
+ {\r
+ int Length = digitArray.Length;\r
+ if (n2.digitArray.Length != Length) MakeSafe(ref n2);\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ digitArray[i] ^= n2.digitArray[i];\r
+ }\r
+\r
+ sign ^= n2.sign;\r
+ }\r
+\r
+ //*************** Fast Static Arithmetic Functions *****************\r
+\r
+ /// <summary>\r
+ /// Adds n1 and n2 and puts result in dest, without intermediate memory allocation\r
+ /// (unsafe if n1 and n2 disagree in precision, safe even if dest is n1 or n2)\r
+ /// </summary>\r
+ /// <param name="dest"></param>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ public static void AddFast(BigInt dest, BigInt n1, BigInt n2)\r
+ {\r
+ //We cast to the highest input precision...\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+\r
+ //Then we up the output precision if less than the input precision.\r
+ if (dest.digitArray.Length < n1.digitArray.Length) n1.MakeSafe(ref dest);\r
+\r
+ int Length = n1.digitArray.Length;\r
+\r
+ if (n1.sign == n2.sign)\r
+ {\r
+ //Copies sources into digit array and working set for all cases, to avoid\r
+ //problems when dest is n1 or n2\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ dest.workingSet[i] = n2.digitArray[i];\r
+ dest.digitArray[i] = n1.digitArray[i];\r
+ }\r
+ dest.AddInternalBits(dest.workingSet);\r
+ dest.sign = n1.sign;\r
+ }\r
+ else\r
+ {\r
+ bool lessThan = LtInt(n1, n2);\r
+\r
+ if (lessThan)\r
+ {\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ dest.workingSet[i] = n1.digitArray[i];\r
+ dest.digitArray[i] = n2.digitArray[i];\r
+ }\r
+ dest.SubInternalBits(dest.workingSet);\r
+ dest.sign = !n1.sign;\r
+ }\r
+ else\r
+ {\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ dest.workingSet[i] = n2.digitArray[i];\r
+ dest.digitArray[i] = n1.digitArray[i];\r
+ }\r
+ dest.SubInternalBits(dest.workingSet);\r
+ dest.sign = n1.sign;\r
+ }\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Adds n1 and n2 and puts result in dest, without intermediate memory allocation\r
+ /// (unsafe if n1 and n2 disagree in precision, safe even if dest is n1 or n2)\r
+ /// </summary>\r
+ /// <param name="dest"></param>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ public static void SubFast(BigInt dest, BigInt n1, BigInt n2)\r
+ {\r
+ //We cast to the highest input precision...\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+\r
+ //Then we up the output precision if less than the input precision.\r
+ if (dest.digitArray.Length < n1.digitArray.Length) n1.MakeSafe(ref dest);\r
+\r
+ int Length = n1.digitArray.Length;\r
+\r
+ if (n1.sign != n2.sign)\r
+ {\r
+ //Copies sources into digit array and working set for all cases, to avoid\r
+ //problems when dest is n1 or n2\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ dest.workingSet[i] = n2.digitArray[i];\r
+ dest.digitArray[i] = n1.digitArray[i];\r
+ }\r
+ dest.AddInternalBits(dest.workingSet);\r
+ dest.sign = n1.sign;\r
+ }\r
+ else\r
+ {\r
+ bool lessThan = LtInt(n1, n2);\r
+\r
+ if (lessThan)\r
+ {\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ dest.workingSet[i] = n1.digitArray[i];\r
+ dest.digitArray[i] = n2.digitArray[i];\r
+ }\r
+ dest.SubInternalBits(dest.workingSet);\r
+ dest.sign = !n1.sign;\r
+ }\r
+ else\r
+ {\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ dest.workingSet[i] = n2.digitArray[i];\r
+ dest.digitArray[i] = n1.digitArray[i];\r
+ }\r
+ dest.SubInternalBits(dest.workingSet);\r
+ dest.sign = n1.sign;\r
+ }\r
+ }\r
+ }\r
+\r
+ //*************** Static Arithmetic Functions ***************\r
+\r
+ /// <summary>\r
+ /// signed addition of 2 numbers.\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt Add(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+ BigInt result; \r
+\r
+ if (n1.sign == n2.sign)\r
+ {\r
+ result = new BigInt(n1);\r
+ result.AddInternal(n2);\r
+ result.sign = n1.sign;\r
+ }\r
+ else\r
+ {\r
+ bool lessThan = LtInt(n1, n2);\r
+\r
+ if (lessThan)\r
+ {\r
+ result = new BigInt(n2);\r
+ result.SubInternal(n1);\r
+ result.sign = !n1.sign;\r
+ }\r
+ else\r
+ {\r
+ result = new BigInt(n1);\r
+ result.SubInternal(n2);\r
+ result.sign = n1.sign;\r
+ }\r
+ }\r
+\r
+ return result;\r
+ }\r
+\r
+ /// <summary>\r
+ /// signed subtraction of 2 numbers.\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt Sub(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+ BigInt result;\r
+\r
+ if ((n1.sign && !n2.sign) || (!n1.sign && n2.sign))\r
+ {\r
+ result = new BigInt(n1);\r
+ result.AddInternal(n2);\r
+ result.sign = n1.sign;\r
+ }\r
+ else\r
+ {\r
+ bool lessThan = LtInt(n1, n2);\r
+\r
+ if (lessThan)\r
+ {\r
+ result = new BigInt(n2);\r
+ result.SubInternal(n1);\r
+ result.sign = !n1.sign;\r
+ }\r
+ else\r
+ {\r
+ result = new BigInt(n1);\r
+ result.SubInternal(n2);\r
+ result.sign = n1.sign;\r
+ }\r
+ }\r
+\r
+ return result;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Multiplication of two BigInts\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt Mul(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+ \r
+ BigInt result = new BigInt(n1);\r
+ result.Mul(n2);\r
+ return result;\r
+ }\r
+\r
+ /// <summary>\r
+ /// True arbitrary precision divide.\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt Div(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+ BigInt res = new BigInt(n1);\r
+ res.Div(n2);\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// True arbitrary-precision mod operation\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt Mod(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+ BigInt res = new BigInt(n1);\r
+ res.Mod(n2);\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Unsigned multiplication of a BigInt by a small number\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt Mul(BigInt n1, uint n2)\r
+ {\r
+ BigInt result = new BigInt(n1);\r
+ result.MulInternal(n2);\r
+ return result;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Division of a BigInt by a small (unsigned) number\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt Div(BigInt n1, uint n2)\r
+ {\r
+ BigInt result = new BigInt(n1);\r
+ result.DivInternal(n2);\r
+ return result;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Division and remainder of a BigInt by a small (unsigned) number\r
+ /// n1 / n2 = div remainder mod\r
+ /// </summary>\r
+ /// <param name="n1">The number to divide (dividend)</param>\r
+ /// <param name="n2">The number to divide by (divisor)</param>\r
+ /// <param name="div">The quotient (output parameter)</param>\r
+ /// <param name="mod">The remainder (output parameter)</param>\r
+ public static void DivMod(BigInt n1, uint n2, out BigInt div, out BigInt mod)\r
+ {\r
+ div = Div(n1, n2);\r
+ mod = Mul(div, n2);\r
+ mod = Sub(n1, mod);\r
+ }\r
+\r
+ //**************** Static Comparison Functions ***************\r
+\r
+ /// <summary>\r
+ /// true iff n1 is less than n2\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static bool LessThan(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+\r
+ if (n1.sign)\r
+ {\r
+ if (!n2.sign) return true;\r
+ return GtInt(n1, n2);\r
+ }\r
+ else\r
+ {\r
+ if (n2.sign) return false;\r
+ return LtInt(n1, n2);\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// true iff n1 is greater than n2\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static bool GreaterThan(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+\r
+ if (n1.sign)\r
+ {\r
+ if (!n2.sign) return false;\r
+ return LtInt(n1, n2);\r
+ }\r
+ else\r
+ {\r
+ if (n2.sign) return true;\r
+ return GtInt(n1, n2);\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// true iff n1 == n2\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static bool Equals(BigInt n1, BigInt n2)\r
+ {\r
+ return n1.Equals(n2);\r
+ }\r
+\r
+ //***************** Static Bitwise Functions *****************\r
+\r
+ /// <summary>\r
+ /// Bitwise And\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt And(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+ BigInt res = new BigInt(n1);\r
+ res.And(n2);\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Bitwise Or\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt Or(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+ BigInt res = new BigInt(n1);\r
+ res.Or(n2);\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Bitwise Xor\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ public static BigInt Xor(BigInt n1, BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length != n2.digitArray.Length) MakeSafe(ref n1, ref n2);\r
+ BigInt res = new BigInt(n1);\r
+ res.And(n2);\r
+ return res;\r
+ }\r
+\r
+ //**************** Static Operator Overloads *****************\r
+\r
+ /// <summary>\r
+ /// Addition operator\r
+ /// </summary>\r
+ public static BigInt operator +(BigInt n1, BigInt n2)\r
+ {\r
+ return Add(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The subtraction operator\r
+ /// </summary>\r
+ public static BigInt operator -(BigInt n1, BigInt n2)\r
+ {\r
+ return Sub(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The multiplication operator\r
+ /// </summary>\r
+ public static BigInt operator *(BigInt n1, BigInt n2)\r
+ {\r
+ return Mul(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The division operator\r
+ /// </summary>\r
+ public static BigInt operator /(BigInt n1, BigInt n2)\r
+ {\r
+ return Div(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The remainder (mod) operator\r
+ /// </summary>\r
+ public static BigInt operator %(BigInt n1, BigInt n2)\r
+ {\r
+ return Mod(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The left-shift operator\r
+ /// </summary>\r
+ public static BigInt operator <<(BigInt n1, int n2)\r
+ {\r
+ BigInt res = new BigInt(n1);\r
+ res.LSH(n2);\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// The right-shift operator\r
+ /// </summary>\r
+ public static BigInt operator >>(BigInt n1, int n2)\r
+ {\r
+ BigInt res = new BigInt(n1);\r
+ res.RSH(n2);\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// The less than operator\r
+ /// </summary>\r
+ public static bool operator <(BigInt n1, BigInt n2)\r
+ {\r
+ return LessThan(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The greater than operator\r
+ /// </summary>\r
+ public static bool operator >(BigInt n1, BigInt n2)\r
+ {\r
+ return GreaterThan(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The less than or equal to operator\r
+ /// </summary>\r
+ public static bool operator <=(BigInt n1, BigInt n2)\r
+ {\r
+ return !GreaterThan(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The greater than or equal to operator\r
+ /// </summary>\r
+ public static bool operator >=(BigInt n1, BigInt n2)\r
+ {\r
+ return !LessThan(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The equality operator\r
+ /// </summary>\r
+ public static bool operator ==(BigInt n1, BigInt n2)\r
+ {\r
+ return Equals(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The inequality operator\r
+ /// </summary>\r
+ public static bool operator !=(BigInt n1, BigInt n2)\r
+ {\r
+ return !Equals(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The bitwise AND operator\r
+ /// </summary>\r
+ public static BigInt operator &(BigInt n1, BigInt n2)\r
+ {\r
+ return And(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The bitwise OR operator\r
+ /// </summary>\r
+ public static BigInt operator |(BigInt n1, BigInt n2)\r
+ {\r
+ return Or(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The bitwise eXclusive OR operator\r
+ /// </summary>\r
+ public static BigInt operator ^(BigInt n1, BigInt n2)\r
+ {\r
+ return Xor(n1, n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// The increment operator\r
+ /// </summary>\r
+ public static BigInt operator ++(BigInt n1)\r
+ {\r
+ n1.Increment();\r
+ return n1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// The decrement operator\r
+ /// </summary>\r
+ public static BigInt operator --(BigInt n1)\r
+ {\r
+ n1.Decrement();\r
+ return n1;\r
+ }\r
+\r
+ //**************** Private Member Functions *****************\r
+\r
+ /// <summary>\r
+ /// Unsigned multiplication and assignment by a small number\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ private void MulInternal(uint n2)\r
+ {\r
+ int Length = digitArray.Length;\r
+ ulong n2long = (ulong)n2;\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ workingSet[i] = 0;\r
+ }\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ if (digitArray[i] == 0) continue;\r
+ ulong temp = (ulong)digitArray[i] * n2long;\r
+ digitArray[i] = (uint)(temp & 0xffffffff);\r
+ if (i + 1 < Length) workingSet[i + 1] = (uint)(temp >> 32);\r
+ }\r
+\r
+ AddInternalBits(workingSet);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Unsigned division and assignment by a small number\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ private void DivInternal(uint n2)\r
+ {\r
+ int Length = digitArray.Length;\r
+ ulong carry = 0;\r
+\r
+ //Divide each digit by the small number.\r
+ for (int i = Length - 1; i >= 0; i--)\r
+ {\r
+ ulong temp = (ulong)digitArray[i] + (carry << 32);\r
+ digitArray[i] = (uint)(temp / (ulong)n2);\r
+ carry = temp % (ulong)n2;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Adds a signed integer to the number.\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ private void AddInternal(int n1)\r
+ {\r
+ if (n1 < 0)\r
+ {\r
+ SubInternal(-n1);\r
+ return;\r
+ }\r
+\r
+ uint carry = 0;\r
+ int length = digitArray.Length;\r
+\r
+ for (int i = 0; i < length && !(n1 == 0 && carry == 0); i++)\r
+ {\r
+ uint temp = digitArray[i];\r
+ digitArray[i] += (uint)n1 + carry;\r
+\r
+ carry = (digitArray[i] <= temp) ? 1u: 0u;\r
+ \r
+ n1 = 0;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Subtract a signed integer from the number.\r
+ /// This is internal because it will fail spectacularly if this number is negative or if n1 is bigger than this number.\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ private void SubInternal(int n1)\r
+ {\r
+ if (n1 < 0)\r
+ {\r
+ AddInternal(-n1);\r
+ return;\r
+ }\r
+\r
+ uint carry = 0;\r
+ int length = digitArray.Length;\r
+\r
+ for (int i = 0; i < length && !(n1 == 0 && carry == 0); i++)\r
+ {\r
+ uint temp = digitArray[i];\r
+ digitArray[i] -= ((uint)n1 + carry);\r
+\r
+ carry = (digitArray[i] >= temp) ? 1u: 0u;\r
+\r
+ n1 = 0;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Adds a signed integer to the number.\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ private bool AddInternal(uint n1)\r
+ {\r
+ uint carry = 0;\r
+ int length = digitArray.Length;\r
+\r
+ for (int i = 0; i < length && !(n1 == 0 && carry == 0); i++)\r
+ {\r
+ uint temp = digitArray[i];\r
+ digitArray[i] += n1 + carry;\r
+\r
+ carry = (digitArray[i] <= temp) ? 1u: 0u;\r
+\r
+ n1 = 0;\r
+ }\r
+\r
+ return (carry != 0);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Internally subtracts a uint from the number (sign insensitive)\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <returns></returns>\r
+ private bool SubInternal(uint n1)\r
+ {\r
+ uint carry = 0;\r
+ int length = digitArray.Length;\r
+\r
+ for (int i = 0; i < length && !(n1 == 0 && carry == 0); i++)\r
+ {\r
+ uint temp = digitArray[i];\r
+ digitArray[i] -= (n1 + carry);\r
+\r
+ carry = (digitArray[i] >= temp) ? 1u: 0u;\r
+\r
+ n1 = 0;\r
+ }\r
+\r
+ return (carry != 0);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Internal increment function (sign insensitive)\r
+ /// </summary>\r
+ private bool IncrementInt()\r
+ {\r
+ uint carry = 1;\r
+\r
+ int length = digitArray.Length;\r
+\r
+ for (int i = 0; i < length && carry != 0; i++)\r
+ {\r
+ uint temp = digitArray[i];\r
+ digitArray[i]++;\r
+\r
+ if (digitArray[i] > temp) carry = 0;\r
+ }\r
+\r
+ return (carry != 0);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Internal increment function (sign insensitive)\r
+ /// </summary>\r
+ private bool DecrementInt()\r
+ {\r
+ uint carry = 1;\r
+\r
+ int length = digitArray.Length;\r
+\r
+ for (int i = 0; i < length && carry != 0; i++)\r
+ {\r
+ uint temp = digitArray[i];\r
+ digitArray[i]--;\r
+\r
+ if (digitArray[i] < temp) carry = 0;\r
+ }\r
+\r
+ return (carry != 0);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Used to add a digit array to a big int.\r
+ /// </summary>\r
+ /// <param name="digitsToAdd"></param>\r
+ private uint AddInternalBits(uint[] digitsToAdd)\r
+ {\r
+ uint carry = 0;\r
+\r
+ int Length = digitArray.Length;\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ //Necessary because otherwise the carry calculation could go bad.\r
+ if (digitsToAdd[i] == 0 && carry == 0) continue;\r
+\r
+ uint temp = digitArray[i];\r
+ digitArray[i] += (digitsToAdd[i] + carry);\r
+\r
+ carry = 0;\r
+ if (digitArray[i] <= temp) carry = 1;\r
+ }\r
+\r
+ return carry;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Used to add with matching signs (true addition of the digit arrays)\r
+ /// This is internal because it will fail spectacularly if n1 is negative.\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ private uint AddInternal(BigInt n1)\r
+ {\r
+ return AddInternalBits(n1.digitArray);\r
+ }\r
+\r
+ private uint SubInternalBits(uint[] digitsToAdd)\r
+ {\r
+ uint carry = 0;\r
+\r
+ int Length = digitArray.Length;\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ //Necessary because otherwise the carry calculation could go bad.\r
+ if (digitsToAdd[i] == 0 && carry == 0) continue;\r
+\r
+ uint temp = digitArray[i];\r
+ digitArray[i] -= (digitsToAdd[i] + carry);\r
+\r
+ carry = 0;\r
+ if (digitArray[i] >= temp) carry = 1;\r
+ }\r
+\r
+ return carry;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Used to subtract n1 (true subtraction of digit arrays) - n1 must be less than or equal to this number\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ private uint SubInternal(BigInt n1)\r
+ {\r
+ return SubInternalBits(n1.digitArray);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Returns the length of the BigInt in 32-bit words for a given decimal precision\r
+ /// </summary>\r
+ /// <param name="precision"></param>\r
+ /// <returns></returns>\r
+ private static int GetRequiredDigitsForPrecision(PrecisionSpec precision)\r
+ {\r
+ int bits = precision.NumBits;\r
+ return ((bits - 1) >> 5) + 1;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Initialises the BigInt to a desired decimal precision\r
+ /// </summary>\r
+ /// <param name="precision"></param>\r
+ private void Init(PrecisionSpec precision)\r
+ {\r
+ int numDigits = GetRequiredDigitsForPrecision(precision);\r
+ digitArray = new uint[numDigits];\r
+ workingSet = new uint[numDigits];\r
+ pres = precision;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Initialises the BigInt from a string, given a base and precision\r
+ /// </summary>\r
+ /// <param name="init"></param>\r
+ /// <param name="precision"></param>\r
+ /// <param name="numberBase"></param>\r
+ private void InitFromString(string init, PrecisionSpec precision, int numberBase)\r
+ {\r
+ PrecisionSpec test;\r
+ if (numberBase == 2)\r
+ {\r
+ test = new PrecisionSpec(init.Length, PrecisionSpec.BaseType.BIN);\r
+ }\r
+ else if (numberBase == 8)\r
+ {\r
+ test = new PrecisionSpec(init.Length, PrecisionSpec.BaseType.OCT);\r
+ }\r
+ else if (numberBase == 10)\r
+ {\r
+ test = new PrecisionSpec(init.Length, PrecisionSpec.BaseType.DEC);\r
+ }\r
+ else if (numberBase == 16)\r
+ {\r
+ test = new PrecisionSpec(init.Length, PrecisionSpec.BaseType.HEX);\r
+ }\r
+ else\r
+ {\r
+ throw new ArgumentOutOfRangeException();\r
+ }\r
+\r
+ //if (test.NumBits > precision.NumBits) precision = test;\r
+ Init(precision);\r
+ FromStringInt(init, numberBase);\r
+ }\r
+\r
+ //************ Helper Functions for floating point *************\r
+\r
+ /// <summary>\r
+ /// Returns true if only the top bit is set: i.e. if the floating-point number is a power of 2\r
+ /// </summary>\r
+ /// <returns></returns>\r
+ public bool IsTopBitOnlyBit()\r
+ {\r
+ int length = digitArray.Length;\r
+\r
+ if (digitArray[length - 1] != 0x80000000u) return false;\r
+ length--;\r
+ for (int i = 0; i < length; i++)\r
+ {\r
+ if (digitArray[i] != 0) return false;\r
+ }\r
+\r
+ return true;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Zeroes the n most significant bits of the number\r
+ /// </summary>\r
+ /// <param name="bits"></param>\r
+ public void ZeroBitsHigh(int bits)\r
+ {\r
+ //Already done.\r
+ if (bits <= 0) return;\r
+\r
+ int length = digitArray.Length;\r
+\r
+ //The entire digit array.\r
+ if ((bits >> 5) > length)\r
+ {\r
+ bits = length << 5;\r
+ }\r
+\r
+ int remBits = (bits & 31);\r
+ int startDigit = length - ((bits >> 5) + 1);\r
+\r
+ if (remBits != 0)\r
+ {\r
+ digitArray[startDigit] = digitArray[startDigit] & (0xffffffffu >> remBits);\r
+ }\r
+\r
+ for (int i = startDigit + 1; i < length; i++)\r
+ {\r
+ digitArray[i] = 0;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Zeroes the least-significant n bits.\r
+ /// </summary>\r
+ /// <param name="bits"></param>\r
+ public void ZeroBits(int bits)\r
+ {\r
+ //Already done.\r
+ if (bits <= 0) return;\r
+\r
+ //The entire digit array.\r
+ if ((bits >> 5) > digitArray.Length)\r
+ {\r
+ bits = digitArray.Length << 5;\r
+ }\r
+\r
+ int remBits = (bits & 31);\r
+ int startDigit = bits >> 5;\r
+ \r
+ if (remBits != 0)\r
+ {\r
+ UInt32 startMask = 0xffffffffu & ~(UInt32)(((1 << remBits) - 1));\r
+ digitArray[startDigit] = digitArray[startDigit] & startMask;\r
+ }\r
+\r
+ for (int i = startDigit - 1; i >= 0; i--)\r
+ {\r
+ digitArray[i] = 0;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Sets the number to 0\r
+ /// </summary>\r
+ public void Zero()\r
+ {\r
+ int length = digitArray.Length;\r
+\r
+ for (int i = 0; i < length; i++)\r
+ {\r
+ digitArray[i] = 0;\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Rounds off the least significant bits of the number.\r
+ /// Can only round off up to 31 bits.\r
+ /// </summary>\r
+ /// <param name="bits">number of bits to round</param>\r
+ /// <returns></returns>\r
+ public bool Round(int bits)\r
+ {\r
+ //Always less than 32 bits, please!\r
+ if (bits < 32)\r
+ {\r
+ uint pow2 = 1u << bits;\r
+ uint test = digitArray[0] & (pow2 >> 1);\r
+\r
+ //Zero the lower bits\r
+ digitArray[0] = digitArray[0] & ~(pow2 - 1);\r
+\r
+ if (test != 0)\r
+ {\r
+ bool bRet = AddInternal(pow2);\r
+ digitArray[digitArray.Length - 1] = digitArray[digitArray.Length - 1] | 0x80000000;\r
+ return bRet;\r
+ }\r
+ }\r
+\r
+ return false;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Used for casting between BigFloats of different precisions - this assumes\r
+ /// that the number is a normalised mantissa.\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ /// <returns>true if a round-up caused the high bits to become zero</returns>\r
+ public bool AssignHigh(BigInt n2)\r
+ {\r
+ int length = digitArray.Length;\r
+ int length2 = n2.digitArray.Length;\r
+ int minWords = (length < length2) ? length: length2;\r
+ bool bRet = false;\r
+\r
+ for (int i = 1; i <= minWords; i++)\r
+ {\r
+ digitArray[length - i] = n2.digitArray[length2 - i];\r
+ }\r
+\r
+ if (length2 > length && n2.digitArray[length2 - (length + 1)] >= 0x80000000)\r
+ {\r
+ bRet = IncrementInt();\r
+\r
+ //Because we are assuming normalisation, we set the top bit (it will already be set if\r
+ //bRet is false.\r
+ digitArray[length - 1] = digitArray[length - 1] | 0x80000000;\r
+ }\r
+\r
+ sign = n2.sign;\r
+\r
+ return bRet;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Used for casting between long ints or doubles and floating-point numbers\r
+ /// </summary>\r
+ /// <param name="digits"></param>\r
+ public void SetHighDigits(Int64 digits)\r
+ {\r
+ digitArray[digitArray.Length - 1] = (uint)(digits >> 32);\r
+ if (digitArray.Length > 1) digitArray[digitArray.Length - 2] = (uint)(digits & 0xffffffff);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Used for casting between ints and doubles or floats.\r
+ /// </summary>\r
+ /// <param name="digit"></param>\r
+ public void SetHighDigit(UInt32 digit)\r
+ {\r
+ digitArray[digitArray.Length - 1] = digit;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Helper function for floating-point - extends the number to twice the precision\r
+ /// and shifts the digits into the upper bits.\r
+ /// </summary>\r
+ public void Pad()\r
+ {\r
+ int length = digitArray.Length;\r
+ int digits = length << 1;\r
+\r
+ UInt32[] newDigitArray = new UInt32[digits];\r
+ workingSet = new UInt32[digits];\r
+\r
+ for (int i = 0; i < length; i++)\r
+ {\r
+ newDigitArray[i + length] = digitArray[i];\r
+ }\r
+\r
+ digitArray = newDigitArray;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Helper function for floating-point - extends the number to twice the precision...\r
+ /// This is a necessary step in floating-point division.\r
+ /// </summary>\r
+ public void Extend()\r
+ {\r
+ SetNumDigits(digitArray.Length * 2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Gets the highest big of the integer (used for floating point stuff)\r
+ /// </summary>\r
+ /// <returns></returns>\r
+ public uint GetTopBit()\r
+ {\r
+ return (digitArray[digitArray.Length - 1] >> 31);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Used for floating point multiplication, this shifts the number so that\r
+ /// the highest bit is set, and returns the number of places shifted.\r
+ /// </summary>\r
+ /// <returns></returns>\r
+ public int Normalise()\r
+ {\r
+ if (IsZero()) return 0;\r
+\r
+ int MSD = GetMSD();\r
+ int digitShift = (digitArray.Length - (MSD + 1));\r
+ int bitShift = (31 - GetMSB(digitArray[MSD])) + (digitShift << 5);\r
+ LSH(bitShift);\r
+ return bitShift;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Gets the most significant bit\r
+ /// </summary>\r
+ /// <param name="value">the input to search for the MSB in</param>\r
+ /// <returns>-1 if the input was zero, the position of the MSB otherwise</returns>\r
+ public static int GetMSB(UInt32 value)\r
+ {\r
+ if (value == 0) return -1;\r
+\r
+ uint mask1 = 0xffff0000;\r
+ uint mask2 = 0xff00;\r
+ uint mask3 = 0xf0;\r
+ uint mask4 = 0xc; //1100 in binary\r
+ uint mask5 = 0x2; //10 in binary\r
+\r
+ int iPos = 0;\r
+\r
+ //Unrolled binary search for the most significant bit.\r
+ if ((value & mask1) != 0) iPos += 16;\r
+ mask2 <<= iPos;\r
+\r
+ if ((value & mask2) != 0) iPos += 8;\r
+ mask3 <<= iPos;\r
+\r
+ if ((value & mask3) != 0) iPos += 4;\r
+ mask4 <<= iPos;\r
+\r
+ if ((value & mask4) != 0) iPos += 2;\r
+ mask5 <<= iPos;\r
+\r
+ if ((value & mask5) != 0) iPos++;\r
+\r
+ return iPos;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Gets the most significant bit\r
+ /// </summary>\r
+ /// <param name="value">the input to search for the MSB in</param>\r
+ /// <returns>-1 if the input was zero, the position of the MSB otherwise</returns>\r
+ public static int GetMSB(UInt64 value)\r
+ {\r
+ if (value == 0) return -1;\r
+\r
+ UInt64 mask0 = 0xffffffff00000000ul;\r
+ UInt64 mask1 = 0xffff0000;\r
+ UInt64 mask2 = 0xff00;\r
+ UInt64 mask3 = 0xf0;\r
+ UInt64 mask4 = 0xc; //1100 in binary\r
+ UInt64 mask5 = 0x2; //10 in binary\r
+\r
+ int iPos = 0;\r
+\r
+ //Unrolled binary search for the most significant bit.\r
+ if ((value & mask0) != 0) iPos += 32;\r
+ mask1 <<= iPos;\r
+\r
+ if ((value & mask1) != 0) iPos += 16;\r
+ mask2 <<= iPos;\r
+\r
+ if ((value & mask2) != 0) iPos += 8;\r
+ mask3 <<= iPos;\r
+\r
+ if ((value & mask3) != 0) iPos += 4;\r
+ mask4 <<= iPos;\r
+\r
+ if ((value & mask4) != 0) iPos += 2;\r
+ mask5 <<= iPos;\r
+\r
+ if ((value & mask5) != 0) iPos++;\r
+\r
+ return iPos;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Gets the most significant bit\r
+ /// </summary>\r
+ /// <param name="value">the input to search for the MSB in</param>\r
+ /// <returns>-1 if the input was zero, the position of the MSB otherwise</returns>\r
+ public static int GetMSB(BigInt value)\r
+ {\r
+ int digit = value.GetMSD();\r
+ int bit = GetMSB(value.digitArray[digit]);\r
+ return (digit << 5) + bit;\r
+ }\r
+\r
+ //**************** Helper Functions for Div ********************\r
+\r
+ /// <summary>\r
+ /// Gets the index of the most significant digit\r
+ /// </summary>\r
+ /// <returns></returns>\r
+ private int GetMSD()\r
+ {\r
+ for (int i = digitArray.Length - 1; i >= 0; i--)\r
+ {\r
+ if (digitArray[i] != 0) return i;\r
+ }\r
+\r
+ return 0;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Gets the required bitshift for the Div_32 algorithm\r
+ /// </summary>\r
+ /// <returns></returns>\r
+ private int GetDivBitshift()\r
+ {\r
+ uint digit = digitArray[GetMSD()];\r
+ uint mask1 = 0xffff0000;\r
+ uint mask2 = 0xff00;\r
+ uint mask3 = 0xf0; \r
+ uint mask4 = 0xc; //1100 in binary\r
+ uint mask5 = 0x2; //10 in binary\r
+\r
+ int iPos = 0;\r
+\r
+ //Unrolled binary search for the most significant bit.\r
+ if ((digit & mask1) != 0) iPos += 16;\r
+ mask2 <<= iPos;\r
+\r
+ if ((digit & mask2) != 0) iPos += 8;\r
+ mask3 <<= iPos;\r
+ \r
+ if ((digit & mask3) != 0) iPos += 4;\r
+ mask4 <<= iPos;\r
+\r
+ if ((digit & mask4) != 0) iPos += 2;\r
+ mask5 <<= iPos;\r
+\r
+ if ((digit & mask5) != 0) return 30 - iPos;\r
+\r
+ return 31 - iPos;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Shifts and optionally precision-extends the arguments to prepare for Div_32\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ private static int MakeSafeDiv(BigInt n1, BigInt n2)\r
+ {\r
+ int shift = n2.GetDivBitshift();\r
+ int n1MSD = n1.GetMSD();\r
+\r
+ uint temp = n1.digitArray[n1MSD];\r
+ if (n1MSD == n1.digitArray.Length - 1 && ((temp << shift) >> shift) != n1.digitArray[n1MSD])\r
+ {\r
+ //Precision-extend n1 and n2 if necessary\r
+ int digits = n1.digitArray.Length;\r
+ n1.SetNumDigits(digits + 1);\r
+ n2.SetNumDigits(digits + 1);\r
+ }\r
+\r
+ //Logical left-shift n1 and n2\r
+ n1.LSH(shift);\r
+ n2.LSH(shift);\r
+\r
+ return shift;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Schoolbook division helper function.\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <param name="Q">Quotient output value</param>\r
+ /// <param name="R">Remainder output value</param>\r
+ private static void Div_31(BigInt n1, BigInt n2, BigInt Q, BigInt R)\r
+ {\r
+ int digitsN1 = n1.GetMSD() + 1;\r
+ int digitsN2 = n2.GetMSD() + 1; \r
+\r
+ if ((digitsN1 > digitsN2))\r
+ {\r
+ BigInt n1New = new BigInt(n2);\r
+ n1New.DigitShiftSelfLeft(1);\r
+\r
+ //If n1 >= n2 * 2^32\r
+ if (!LtInt(n1, n1New))\r
+ {\r
+ n1New.sign = n1.sign;\r
+ SubFast(n1New, n1, n1New);\r
+\r
+ Div_32(n1New, n2, Q, R);\r
+\r
+ //Q = (A - B*2^32)/B + 2^32\r
+ Q.Add2Pow32Self();\r
+ return;\r
+ }\r
+ }\r
+\r
+ UInt32 q = 0;\r
+\r
+ if (digitsN1 >= 2)\r
+ {\r
+ UInt64 q64 = ((((UInt64)n1.digitArray[digitsN1 - 1]) << 32) + n1.digitArray[digitsN1 - 2]) / (UInt64)n2.digitArray[digitsN2 - 1];\r
+\r
+ if (q64 > 0xfffffffful)\r
+ {\r
+ q = 0xffffffff;\r
+ }\r
+ else\r
+ {\r
+ q = (UInt32)q64;\r
+ }\r
+ }\r
+\r
+ BigInt temp = Mul(n2, q);\r
+\r
+ if (GtInt(temp, n1))\r
+ {\r
+ temp.SubInternalBits(n2.digitArray);\r
+ q--;\r
+\r
+ if (GtInt(temp, n1))\r
+ {\r
+ temp.SubInternalBits(n2.digitArray);\r
+ q--;\r
+ }\r
+ }\r
+\r
+ Q.Zero();\r
+ Q.digitArray[0] = q;\r
+ R.Assign(n1);\r
+ R.SubInternalBits(temp.digitArray);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Schoolbook division algorithm\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <param name="Q"></param>\r
+ /// <param name="R"></param>\r
+ private static void Div_32(BigInt n1, BigInt n2, BigInt Q, BigInt R)\r
+ {\r
+ int digitsN1 = n1.GetMSD() + 1;\r
+ int digitsN2 = n2.GetMSD() + 1;\r
+\r
+ //n2 is bigger than n1\r
+ if (digitsN1 < digitsN2)\r
+ {\r
+ R.AssignInt(n1);\r
+ Q.Zero();\r
+ return;\r
+ }\r
+\r
+ if (digitsN1 == digitsN2)\r
+ {\r
+ //n2 is bigger than n1\r
+ if (LtInt(n1, n2))\r
+ {\r
+ R.AssignInt(n1);\r
+ Q.Zero();\r
+ return;\r
+ }\r
+\r
+ //n2 >= n1, but less the 2x n1 (initial conditions make this certain)\r
+ Q.Zero();\r
+ Q.digitArray[0] = 1;\r
+ R.Assign(n1);\r
+ R.SubInternalBits(n2.digitArray);\r
+ return;\r
+ }\r
+\r
+ int digits = digitsN1 - (digitsN2 + 1);\r
+\r
+ //Algorithm Div_31 can be used to get the answer in O(n) time.\r
+ if (digits == 0)\r
+ {\r
+ Div_31(n1, n2, Q, R);\r
+ return;\r
+ }\r
+\r
+ BigInt n1New = DigitShiftRight(n1, digits);\r
+ BigInt s = DigitTruncate(n1, digits);\r
+\r
+ BigInt Q2 = new BigInt(n1, n1.pres, true);\r
+ BigInt R2 = new BigInt(n1, n1.pres, true);\r
+\r
+ Div_31(n1New, n2, Q2, R2);\r
+\r
+ R2.DigitShiftSelfLeft(digits);\r
+ R2.Add(s);\r
+\r
+ Div_32(R2, n2, Q, R);\r
+\r
+ Q2.DigitShiftSelfLeft(digits);\r
+ Q.Add(Q2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Sets the n-th bit of the number to 1\r
+ /// </summary>\r
+ /// <param name="bit">the index of the bit to set</param>\r
+ public void SetBit(int bit)\r
+ {\r
+ int digit = (bit >> 5);\r
+ if (digit >= digitArray.Length) return;\r
+ digitArray[digit] = digitArray[digit] | (1u << (bit - (digit << 5)));\r
+ }\r
+\r
+ /// <summary>\r
+ /// Sets the n-th bit of the number to 0\r
+ /// </summary>\r
+ /// <param name="bit">the index of the bit to set</param>\r
+ public void ClearBit(int bit)\r
+ {\r
+ int digit = (bit >> 5);\r
+ if (digit >= digitArray.Length) return;\r
+ digitArray[digit] = digitArray[digit] & (~(1u << (bit - (digit << 5))));\r
+ }\r
+\r
+ /// <summary>\r
+ /// Returns the n-th bit, counting from the MSB to the LSB\r
+ /// </summary>\r
+ /// <param name="bit">the index of the bit to return</param>\r
+ /// <returns>1 if the bit is 1, 0 otherwise</returns>\r
+ public uint GetBitFromTop(int bit)\r
+ {\r
+ if (bit < 0) return 0;\r
+ int wordCount = (bit >> 5);\r
+ int upBit = 31 - (bit & 31);\r
+ if (wordCount >= digitArray.Length) return 0;\r
+\r
+ return ((digitArray[digitArray.Length - (wordCount + 1)] & (1u << upBit)) >> upBit);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Assigns n2 to 'this'\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ public void Assign(BigInt n2)\r
+ {\r
+ if (digitArray.Length != n2.digitArray.Length) MakeSafe(ref n2);\r
+ sign = n2.sign;\r
+ AssignInt(n2);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Assign n2 to 'this', safe only if precision-matched\r
+ /// </summary>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ private void AssignInt(BigInt n2)\r
+ {\r
+ int Length = digitArray.Length;\r
+\r
+ for (int i = 0; i < Length; i++)\r
+ {\r
+ digitArray[i] = n2.digitArray[i];\r
+ }\r
+ }\r
+\r
+ private static BigInt DigitShiftRight(BigInt n1, int digits)\r
+ {\r
+ BigInt res = new BigInt(n1);\r
+\r
+ int Length = res.digitArray.Length;\r
+\r
+ for (int i = 0; i < Length - digits; i++)\r
+ {\r
+ res.digitArray[i] = res.digitArray[i + digits];\r
+ }\r
+\r
+ for (int i = Length - digits; i < Length; i++)\r
+ {\r
+ res.digitArray[i] = 0;\r
+ }\r
+\r
+ return res;\r
+ }\r
+\r
+ private void DigitShiftSelfRight(int digits)\r
+ {\r
+ for (int i = digits; i < digitArray.Length; i++)\r
+ {\r
+ digitArray[i - digits] = digitArray[i];\r
+ }\r
+\r
+ for (int i = digitArray.Length - digits; i < digitArray.Length; i++)\r
+ {\r
+ digitArray[i] = 0;\r
+ }\r
+ }\r
+\r
+ private void DigitShiftSelfLeft(int digits)\r
+ {\r
+ for (int i = digitArray.Length - 1; i >= digits; i--)\r
+ {\r
+ digitArray[i] = digitArray[i - digits];\r
+ }\r
+\r
+ for (int i = digits - 1; i >= 0; i--)\r
+ {\r
+ digitArray[i] = 0;\r
+ }\r
+ }\r
+\r
+ private static BigInt DigitTruncate(BigInt n1, int digits)\r
+ {\r
+ BigInt res = new BigInt(n1);\r
+\r
+ for (int i = res.digitArray.Length - 1; i >= digits; i--)\r
+ {\r
+ res.digitArray[i] = 0;\r
+ }\r
+\r
+ return res;\r
+ }\r
+\r
+ private void Add2Pow32Self()\r
+ {\r
+ int Length = digitArray.Length;\r
+\r
+ uint carry = 1;\r
+\r
+ for (int i = 1; i < Length; i++)\r
+ {\r
+ uint temp = digitArray[i];\r
+ digitArray[i] += carry;\r
+ if (digitArray[i] > temp) return;\r
+ }\r
+\r
+ return;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Sets the number of digits without changing the precision\r
+ /// This method is made public only to facilitate fixed-point operations\r
+ /// It should under no circumstances be used for anything else, because\r
+ /// it breaks the BigNum(PrecisionSpec precision) constructor in dangerous\r
+ /// and unpredictable ways.\r
+ /// </summary>\r
+ /// <param name="digits"></param>\r
+ public void SetNumDigits(int digits)\r
+ {\r
+ if (digits == digitArray.Length) return;\r
+\r
+ UInt32[] newDigitArray = new UInt32[digits];\r
+ workingSet = new UInt32[digits];\r
+\r
+ int numCopy = (digits < digitArray.Length) ? digits : digitArray.Length;\r
+\r
+ for (int i = 0; i < numCopy; i++)\r
+ {\r
+ newDigitArray[i] = digitArray[i];\r
+ }\r
+\r
+ digitArray = newDigitArray;\r
+ }\r
+\r
+ //********************** Explicit casts ***********************\r
+\r
+ /// <summary>\r
+ /// Cast to int\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <returns></returns>\r
+ public static explicit operator Int32(BigInt value)\r
+ {\r
+ if (value.digitArray[0] == 0x80000000 && value.sign) return Int32.MinValue;\r
+ int res = (int)(value.digitArray[0] & 0x7fffffff);\r
+ if (value.sign) res = -res;\r
+ return res;\r
+ }\r
+\r
+ /// <summary>\r
+ /// explicit cast to unsigned int.\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <returns></returns>\r
+ public static explicit operator UInt32(BigInt value)\r
+ {\r
+ if (value.sign) return (UInt32)((Int32)(value));\r
+ return (UInt32)value.digitArray[0];\r
+ }\r
+\r
+ /// <summary>\r
+ /// explicit cast to 64-bit signed integer.\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <returns></returns>\r
+ public static explicit operator Int64(BigInt value)\r
+ {\r
+ if (value.digitArray.Length < 2) return (value.sign ? -((Int64)value.digitArray[0]): ((Int64)value.digitArray[0]));\r
+ UInt64 ret = (((UInt64)value.digitArray[1]) << 32) + (UInt64)value.digitArray[0];\r
+ if (ret == 0x8000000000000000L && value.sign) return Int64.MinValue;\r
+ Int64 signedRet = (Int64)(ret & 0x7fffffffffffffffL);\r
+ if (value.sign) signedRet = -signedRet;\r
+ return signedRet;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Explicit cast to UInt64\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <returns></returns>\r
+ public static explicit operator UInt64(BigInt value)\r
+ {\r
+ if (value.sign) return (UInt64)((Int64)(value));\r
+ if (value.digitArray.Length < 2) return (UInt64)value.digitArray[0];\r
+ return ((((UInt64)value.digitArray[1]) << 32) + (UInt64)value.digitArray[0]);\r
+ }\r
+\r
+ /// <summary>\r
+ /// Cast to string\r
+ /// </summary>\r
+ /// <param name="value"></param>\r
+ /// <returns></returns>\r
+ public static explicit operator string(BigInt value)\r
+ {\r
+ return value.ToString();\r
+ }\r
+\r
+ /// <summary>\r
+ /// Cast from string - this is not wholly safe, because precision is not\r
+ /// specified. You should try to construct a BigInt with the appropriate\r
+ /// constructor instead.\r
+ /// </summary>\r
+ /// <param name="value">The decimal string to convert to a BigInt</param>\r
+ /// <returns>A BigInt of the precision required to encompass the string</returns>\r
+ public static explicit operator BigInt(string value)\r
+ {\r
+ return new BigInt(value);\r
+ }\r
+\r
+ //********************* ToString members **********************\r
+\r
+ /// <summary>\r
+ /// Converts this to a string, in the specified base\r
+ /// </summary>\r
+ /// <param name="numberBase">the base to use (min 2, max 16)</param>\r
+ /// <returns>a string representation of the number</returns>\r
+ public string ToString(int numberBase)\r
+ {\r
+ char[] digitChars = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B', 'C', 'D', 'E', 'F' };\r
+\r
+ string output = "";\r
+\r
+ BigInt clone = new BigInt(this);\r
+ clone.sign = false;\r
+\r
+ int numDigits = 0;\r
+ while (!clone.IsZero())\r
+ {\r
+ if (numberBase == 10 && (numDigits % 3) == 0 && numDigits != 0)\r
+ {\r
+ output = String.Format(",{0}", output);\r
+ }\r
+ else if (numberBase != 10 && (numDigits % 8) == 0 && numDigits != 0)\r
+ {\r
+ output = String.Format(" {0}", output);\r
+ }\r
+\r
+ BigInt div, mod;\r
+ DivMod(clone, (uint)numberBase, out div, out mod);\r
+ int iMod = (int)mod;\r
+ output = String.Format("{0}{1}", digitChars[(int)mod], output);\r
+\r
+ numDigits++;\r
+\r
+ clone = div;\r
+ }\r
+\r
+ if (output.Length == 0) output = String.Format("0");\r
+\r
+ if (sign) output = String.Format("-{0}", output);\r
+\r
+ return output;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Converts the number to a string, in base 10\r
+ /// </summary>\r
+ /// <returns>a string representation of the number in base 10</returns>\r
+ public override string ToString()\r
+ {\r
+ return ToString(10);\r
+ }\r
+\r
+ //***************** Internal helper functions *****************\r
+\r
+ private void FromStringInt(string init, int numberBase)\r
+ {\r
+ char [] digitChars = {'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'};\r
+\r
+ string formattedInput = init.Trim().ToUpper();\r
+\r
+ for (int i = 0; i < formattedInput.Length; i++)\r
+ {\r
+ int digitIndex = Array.IndexOf(digitChars, formattedInput[i]);\r
+\r
+ //Skip fractional part altogether\r
+ if (formattedInput[i] == '.') break;\r
+\r
+ //skip non-digit characters.\r
+ if (digitIndex < 0) continue;\r
+\r
+ //Multiply\r
+ MulInternal((uint)numberBase);\r
+ \r
+ //Add\r
+ AddInternal(digitIndex);\r
+ }\r
+\r
+ if (init.Length > 0 && init[0] == '-') sign = true;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Sign-insensitive less than comparison. \r
+ /// unsafe if n1 and n2 disagree in precision\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ private static bool LtInt(BigInt n1, BigInt n2)\r
+ {\r
+ //MakeSafe(ref n1, ref n2);\r
+\r
+ for (int i = n1.digitArray.Length - 1; i >= 0; i--)\r
+ {\r
+ if (n1.digitArray[i] < n2.digitArray[i]) return true;\r
+ if (n1.digitArray[i] > n2.digitArray[i]) return false;\r
+ }\r
+\r
+ //equal\r
+ return false;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Sign-insensitive greater than comparison. \r
+ /// unsafe if n1 and n2 disagree in precision\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ /// <returns></returns>\r
+ private static bool GtInt(BigInt n1, BigInt n2)\r
+ {\r
+ //MakeSafe(ref n1, ref n2);\r
+\r
+ for (int i = n1.digitArray.Length - 1; i >= 0; i--)\r
+ {\r
+ if (n1.digitArray[i] > n2.digitArray[i]) return true;\r
+ if (n1.digitArray[i] < n2.digitArray[i]) return false;\r
+ }\r
+\r
+ //equal\r
+ return false;\r
+ }\r
+\r
+ /// <summary>\r
+ /// Makes sure the numbers have matching precisions\r
+ /// </summary>\r
+ /// <param name="n1"></param>\r
+ /// <param name="n2"></param>\r
+ private static void MakeSafe(ref BigInt n1, ref BigInt n2)\r
+ {\r
+ if (n1.digitArray.Length == n2.digitArray.Length)\r
+ {\r
+ return;\r
+ }\r
+ else if (n1.digitArray.Length > n2.digitArray.Length)\r
+ {\r
+ n2 = new BigInt(n2, n1.pres);\r
+ }\r
+ else\r
+ {\r
+ n1 = new BigInt(n1, n2.pres);\r
+ }\r
+ }\r
+\r
+ /// <summary>\r
+ /// Makes sure the numbers have matching precisions\r
+ /// </summary>\r
+ /// <param name="n2">the number to match to this</param>\r
+ private void MakeSafe(ref BigInt n2)\r
+ {\r
+ n2 = new BigInt(n2, pres);\r
+ n2.SetNumDigits(digitArray.Length);\r
+ }\r
+\r
+\r
+ private PrecisionSpec pres;\r
+ private bool sign;\r
+ private uint[] digitArray;\r
+ private uint[] workingSet;\r
+ }\r
+\r
+}
\ No newline at end of file
OSDN Git Service
