X-Git-Url: http://git.osdn.jp/view?a=blobdiff_plain;f=dev5%2Fpsychlops%2Fextention%2Fmath%2FBigFloat.cs;fp=dev5%2Fpsychlops%2Fextention%2Fmath%2FBigFloat.cs;h=767f5b4a2d937bb87ef59da7a9311d30aaec70b2;hb=7fe25aa821826f09903fb14def74d6b0376e3b5a;hp=0000000000000000000000000000000000000000;hpb=af114e9e36fd9c2cdf79f95313e0b8712e253ed3;p=psychlops%2Fsilverlight.git
diff --git a/dev5/psychlops/extention/math/BigFloat.cs b/dev5/psychlops/extention/math/BigFloat.cs
new file mode 100644
index 0000000..767f5b4
--- /dev/null
+++ b/dev5/psychlops/extention/math/BigFloat.cs
@@ -0,0 +1,3789 @@
+// http://www.fractal-landscapes.co.uk/bigint.html
+
+using System;
+
+namespace BigNum
+{
+ ///
+ /// An arbitrary-precision floating-point class
+ ///
+ /// Format:
+ /// Each number is stored as an exponent (32-bit signed integer), and a mantissa
+ /// (n-bit) BigInteger. The sign of the number is stored in the BigInteger
+ ///
+ /// Applicability and Performance:
+ /// This class is designed to be used for small extended precisions. It may not be
+ /// safe (and certainly won't be fast) to use it with mixed-precision arguments.
+ /// It does support, but will not be efficient for, numbers over around 2048 bits.
+ ///
+ /// Notes:
+ /// All conversions to and from strings are slow.
+ ///
+ /// Conversions from simple integer types Int32, Int64, UInt32, UInt64 are performed
+ /// using the appropriate constructor, and are relatively fast.
+ ///
+ /// The class is written entirely in managed C# code, with not native or managed
+ /// assembler. The use of native assembler would speed up the multiplication operations
+ /// many times over, and therefore all higher-order operations too.
+ ///
+ public class BigFloat
+ {
+ ///
+ /// Floats can have 4 special value types:
+ ///
+ /// NaN: Not a number (cannot be changed using any operations)
+ /// Infinity: Positive infinity. Some operations e.g. Arctan() allow this input.
+ /// -Infinity: Negative infinity. Some operations allow this input.
+ /// Zero
+ ///
+ public enum SpecialValueType
+ {
+ ///
+ /// Not a special value
+ ///
+ NONE = 0,
+ ///
+ /// Zero
+ ///
+ ZERO,
+ ///
+ /// Positive infinity
+ ///
+ INF_PLUS,
+ ///
+ /// Negative infinity
+ ///
+ INF_MINUS,
+ ///
+ /// Not a number
+ ///
+ NAN
+ }
+
+ ///
+ /// This affects the ToString() method.
+ ///
+ /// With Trim rounding, all insignificant zero digits are drip
+ ///
+ public enum RoundingModeType
+ {
+ ///
+ /// Trim non-significant zeros from ToString output after rounding
+ ///
+ TRIM,
+ ///
+ /// Keep all non-significant zeroes in ToString output after rounding
+ ///
+ EXACT
+ }
+
+ ///
+ /// A wrapper for the signed exponent, avoiding overflow.
+ ///
+ protected struct ExponentAdaptor
+ {
+ ///
+ /// The 32-bit exponent
+ ///
+ public Int32 exponent
+ {
+ get { return expValue; }
+ set { expValue = value; }
+ }
+
+ ///
+ /// Implicit cast to Int32
+ ///
+ public static implicit operator Int32(ExponentAdaptor adaptor)
+ {
+ return adaptor.expValue;
+ }
+
+ ///
+ /// Implicit cast from Int32 to ExponentAdaptor
+ ///
+ ///
+ ///
+ public static implicit operator ExponentAdaptor(Int32 value)
+ {
+ ExponentAdaptor adaptor = new ExponentAdaptor();
+ adaptor.expValue = value;
+ return adaptor;
+ }
+
+ ///
+ /// Overloaded increment operator
+ ///
+ public static ExponentAdaptor operator ++(ExponentAdaptor adaptor)
+ {
+ adaptor = adaptor + 1;
+ return adaptor;
+ }
+
+ ///
+ /// Overloaded decrement operator
+ ///
+ public static ExponentAdaptor operator --(ExponentAdaptor adaptor)
+ {
+ adaptor = adaptor - 1;
+ return adaptor;
+ }
+
+ ///
+ /// Overloaded addition operator
+ ///
+ public static ExponentAdaptor operator +(ExponentAdaptor a1, ExponentAdaptor a2)
+ {
+ if (a1.expValue == Int32.MaxValue) return a1;
+
+ Int64 temp = (Int64)a1.expValue;
+ temp += (Int64)(a2.expValue);
+
+ if (temp > (Int64)Int32.MaxValue)
+ {
+ a1.expValue = Int32.MaxValue;
+ }
+ else if (temp < (Int64)Int32.MinValue)
+ {
+ a1.expValue = Int32.MinValue;
+ }
+ else
+ {
+ a1.expValue = (Int32)temp;
+ }
+
+ return a1;
+ }
+
+ ///
+ /// Overloaded subtraction operator
+ ///
+ public static ExponentAdaptor operator -(ExponentAdaptor a1, ExponentAdaptor a2)
+ {
+ if (a1.expValue == Int32.MaxValue) return a1;
+
+ Int64 temp = (Int64)a1.expValue;
+ temp -= (Int64)(a2.expValue);
+
+ if (temp > (Int64)Int32.MaxValue)
+ {
+ a1.expValue = Int32.MaxValue;
+ }
+ else if (temp < (Int64)Int32.MinValue)
+ {
+ a1.expValue = Int32.MinValue;
+ }
+ else
+ {
+ a1.expValue = (Int32)temp;
+ }
+
+ return a1;
+ }
+
+ ///
+ /// Overloaded multiplication operator
+ ///
+ public static ExponentAdaptor operator *(ExponentAdaptor a1, ExponentAdaptor a2)
+ {
+ if (a1.expValue == Int32.MaxValue) return a1;
+
+ Int64 temp = (Int64)a1.expValue;
+ temp *= (Int64)a2.expValue;
+
+ if (temp > (Int64)Int32.MaxValue)
+ {
+ a1.expValue = Int32.MaxValue;
+ }
+ else if (temp < (Int64)Int32.MinValue)
+ {
+ a1.expValue = Int32.MinValue;
+ }
+ else
+ {
+ a1.expValue = (Int32)temp;
+ }
+
+ return a1;
+ }
+
+ ///
+ /// Overloaded division operator
+ ///
+ public static ExponentAdaptor operator /(ExponentAdaptor a1, ExponentAdaptor a2)
+ {
+ if (a1.expValue == Int32.MaxValue) return a1;
+
+ ExponentAdaptor res = new ExponentAdaptor();
+ res.expValue = a1.expValue / a2.expValue;
+ return res;
+ }
+
+ ///
+ /// Overloaded right-shift operator
+ ///
+ public static ExponentAdaptor operator >>(ExponentAdaptor a1, int shift)
+ {
+ if (a1.expValue == Int32.MaxValue) return a1;
+
+ ExponentAdaptor res = new ExponentAdaptor();
+ res.expValue = a1.expValue >> shift;
+ return res;
+ }
+
+ ///
+ /// Overloaded left-shift operator
+ ///
+ ///
+ ///
+ ///
+ public static ExponentAdaptor operator <<(ExponentAdaptor a1, int shift)
+ {
+ if (a1.expValue == 0) return a1;
+
+ ExponentAdaptor res = new ExponentAdaptor();
+ res.expValue = a1.expValue;
+
+ if (shift > 31)
+ {
+ res.expValue = Int32.MaxValue;
+ }
+ else
+ {
+ Int64 temp = a1.expValue;
+ temp = temp << shift;
+
+ if (temp > (Int64)Int32.MaxValue)
+ {
+ res.expValue = Int32.MaxValue;
+ }
+ else if (temp < (Int64)Int32.MinValue)
+ {
+ res.expValue = Int32.MinValue;
+ }
+ else
+ {
+ res.expValue = (Int32)temp;
+ }
+ }
+
+ return res;
+ }
+
+ private Int32 expValue;
+ }
+
+ //************************ Constructors **************************
+
+ ///
+ /// Constructs a 128-bit BigFloat
+ ///
+ /// Sets the value to zero
+ ///
+ static BigFloat()
+ {
+ RoundingDigits = 3;
+ RoundingMode = RoundingModeType.TRIM;
+ scratch = new BigInt(new PrecisionSpec(128, PrecisionSpec.BaseType.BIN));
+ }
+
+ ///
+ /// Constructs a BigFloat of the required precision
+ ///
+ /// Sets the value to zero
+ ///
+ ///
+ public BigFloat(PrecisionSpec mantissaPrec)
+ {
+ Init(mantissaPrec);
+ }
+
+ ///
+ /// Constructs a big float from a UInt32 to the required precision
+ ///
+ ///
+ ///
+ public BigFloat(UInt32 value, PrecisionSpec mantissaPrec)
+ {
+ int mbWords = ((mantissaPrec.NumBits) >> 5);
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
+ int newManBits = mbWords << 5;
+
+ //For efficiency, we just use a 32-bit exponent
+ exponent = 0;
+
+ mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+ //scratch = new BigInt(mantissa.Precision);
+
+ int bit = BigInt.GetMSB(value);
+ if (bit == -1) return;
+
+ int shift = mantissa.Precision.NumBits - (bit + 1);
+ mantissa.LSH(shift);
+ exponent = bit;
+ }
+
+ ///
+ /// Constructs a BigFloat from an Int32 to the required precision
+ ///
+ ///
+ ///
+ public BigFloat(Int32 value, PrecisionSpec mantissaPrec)
+ {
+ int mbWords = ((mantissaPrec.NumBits) >> 5);
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
+ int newManBits = mbWords << 5;
+
+ //For efficiency, we just use a 32-bit exponent
+ exponent = 0;
+ UInt32 uValue;
+
+ if (value < 0)
+ {
+ if (value == Int32.MinValue)
+ {
+ uValue = 0x80000000;
+ }
+ else
+ {
+ uValue = (UInt32)(-value);
+ }
+ }
+ else
+ {
+ uValue = (UInt32)value;
+ }
+
+ mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+ //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+
+ int bit = BigInt.GetMSB(uValue);
+ if (bit == -1) return;
+
+ int shift = mantissa.Precision.NumBits - (bit + 1);
+ mantissa.LSH(shift);
+ exponent = bit;
+ }
+
+ ///
+ /// Constructs a BigFloat from a 64-bit integer
+ ///
+ ///
+ ///
+ public BigFloat(Int64 value, PrecisionSpec mantissaPrec)
+ {
+ int mbWords = ((mantissaPrec.NumBits) >> 5);
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
+ int newManBits = mbWords << 5;
+
+ //For efficiency, we just use a 32-bit exponent
+ exponent = 0;
+ UInt64 uValue;
+
+ if (value < 0)
+ {
+ if (value == Int64.MinValue)
+ {
+ uValue = 0x80000000;
+ }
+ else
+ {
+ uValue = (UInt64)(-value);
+ }
+ }
+ else
+ {
+ uValue = (UInt64)value;
+ }
+
+ mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+ //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+
+ int bit = BigInt.GetMSB(uValue);
+ if (bit == -1) return;
+
+ int shift = mantissa.Precision.NumBits - (bit + 1);
+ if (shift > 0)
+ {
+ mantissa.LSH(shift);
+ }
+ else
+ {
+ mantissa.SetHighDigit((uint)(uValue >> (-shift)));
+ }
+ exponent = bit;
+ }
+
+ ///
+ /// Constructs a BigFloat from a 64-bit unsigned integer
+ ///
+ ///
+ ///
+ public BigFloat(UInt64 value, PrecisionSpec mantissaPrec)
+ {
+ int mbWords = ((mantissaPrec.NumBits) >> 5);
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
+ int newManBits = mbWords << 5;
+
+ //For efficiency, we just use a 32-bit exponent
+ exponent = 0;
+
+ int bit = BigInt.GetMSB(value);
+
+ mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+ //scratch = new BigInt(mantissa.Precision);
+
+ int shift = mantissa.Precision.NumBits - (bit + 1);
+ if (shift > 0)
+ {
+ mantissa.LSH(shift);
+ }
+ else
+ {
+ mantissa.SetHighDigit((uint)(value >> (-shift)));
+ }
+ exponent = bit;
+ }
+
+ ///
+ /// Constructs a BigFloat from a BigInt, using the specified precision
+ ///
+ ///
+ ///
+ public BigFloat(BigInt value, PrecisionSpec mantissaPrec)
+ {
+ if (value.IsZero())
+ {
+ Init(mantissaPrec);
+ SetZero();
+ return;
+ }
+
+ mantissa = new BigInt(value, mantissaPrec);
+ exponent = BigInt.GetMSB(value);
+ mantissa.Normalise();
+ }
+
+ ///
+ /// Construct a BigFloat from a double-precision floating point number
+ ///
+ ///
+ ///
+ public BigFloat(double value, PrecisionSpec mantissaPrec)
+ {
+ if (value == 0.0)
+ {
+ Init(mantissaPrec);
+ return;
+ }
+
+ bool sign = (value < 0) ? true : false;
+
+ long bits = BitConverter.DoubleToInt64Bits(value);
+ // Note that the shift is sign-extended, hence the test against -1 not 1
+ int valueExponent = (int)((bits >> 52) & 0x7ffL);
+ long valueMantissa = bits & 0xfffffffffffffL;
+
+ //The mantissa is stored with the top bit implied.
+ valueMantissa = valueMantissa | 0x10000000000000L;
+
+ //The exponent is biased by 1023.
+ exponent = valueExponent - 1023;
+
+ //Round the number of bits to the nearest word.
+ int mbWords = ((mantissaPrec.NumBits) >> 5);
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
+ int newManBits = mbWords << 5;
+
+ mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+ //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+
+ if (newManBits >= 64)
+ {
+ //The mantissa is 53 bits now, so add 11 to put it in the right place.
+ mantissa.SetHighDigits(valueMantissa << 11);
+ }
+ else
+ {
+ //To get the top word of the mantissa, shift up by 11 and down by 32 = down by 21
+ mantissa.SetHighDigit((uint)(valueMantissa >> 21));
+ }
+
+ mantissa.Sign = sign;
+ }
+
+ ///
+ /// Copy constructor
+ ///
+ ///
+ public BigFloat(BigFloat value)
+ {
+ Init(value.mantissa.Precision);
+ exponent = value.exponent;
+ mantissa.Assign(value.mantissa);
+ }
+
+ ///
+ /// Copy Constructor - constructs a new BigFloat with the specified precision, copying the old one.
+ ///
+ /// The value is rounded towards zero in the case where precision is decreased. The Round() function
+ /// should be used beforehand if a correctly rounded result is required.
+ ///
+ ///
+ ///
+ public BigFloat(BigFloat value, PrecisionSpec mantissaPrec)
+ {
+ Init(mantissaPrec);
+ exponent = value.exponent;
+ if (mantissa.AssignHigh(value.mantissa)) exponent++;
+ }
+
+ ///
+ /// Constructs a BigFloat from a string
+ ///
+ ///
+ ///
+ public BigFloat(string value, PrecisionSpec mantissaPrec)
+ {
+ Init(mantissaPrec);
+
+ PrecisionSpec extendedPres = new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN);
+ BigFloat ten = new BigFloat(10, extendedPres);
+ BigFloat iPart = new BigFloat(extendedPres);
+ BigFloat fPart = new BigFloat(extendedPres);
+ BigFloat tenRCP = ten.Reciprocal();
+
+ if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol))
+ {
+ SetNaN();
+ return;
+ }
+ else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol))
+ {
+ SetInfPlus();
+ return;
+ }
+ else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol))
+ {
+ SetInfMinus();
+ return;
+ }
+
+ string decimalpoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;
+
+ char[] digitChars = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ',', '.' };
+
+ //Read in the integer part up the the decimal point.
+ bool sign = false;
+ value = value.Trim();
+
+ int i = 0;
+
+ if (value.Length > i && value[i] == '-')
+ {
+ sign = true;
+ i++;
+ }
+
+ if (value.Length > i && value[i] == '+')
+ {
+ i++;
+ }
+
+ for ( ; i < value.Length; i++)
+ {
+ //break on decimal point
+ if (value[i] == decimalpoint[0]) break;
+
+ int digit = Array.IndexOf(digitChars, value[i]);
+ if (digit < 0) break;
+
+ //Ignore place separators (assumed either , or .)
+ if (digit > 9) continue;
+
+ if (i > 0) iPart.Mul(ten);
+ iPart.Add(new BigFloat(digit, extendedPres));
+ }
+
+ //If we've run out of characters, assign everything and return
+ if (i == value.Length)
+ {
+ iPart.mantissa.Sign = sign;
+ exponent = iPart.exponent;
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
+ return;
+ }
+
+ //Assign the characters after the decimal point to fPart
+ if (value[i] == '.' && i < value.Length - 1)
+ {
+ BigFloat RecipToUse = new BigFloat(tenRCP);
+
+ for (i++; i < value.Length; i++)
+ {
+ int digit = Array.IndexOf(digitChars, value[i]);
+ if (digit < 0) break;
+ BigFloat temp = new BigFloat(digit, extendedPres);
+ temp.Mul(RecipToUse);
+ RecipToUse.Mul(tenRCP);
+ fPart.Add(temp);
+ }
+ }
+
+ //If we're run out of characters, add fPart and iPart and return
+ if (i == value.Length)
+ {
+ iPart.Add(fPart);
+ iPart.mantissa.Sign = sign;
+ exponent = iPart.exponent;
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
+ return;
+ }
+
+ if (value[i] == '+' || value[i] == '-') i++;
+
+ if (i == value.Length)
+ {
+ iPart.Add(fPart);
+ iPart.mantissa.Sign = sign;
+ exponent = iPart.exponent;
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
+ return;
+ }
+
+ //Look for exponential notation.
+ if ((value[i] == 'e' || value[i] == 'E') && i < value.Length - 1)
+ {
+ //Convert the exponent to an int.
+ int exp;
+
+ try
+ {
+ exp = System.Convert.ToInt32(new string(value.ToCharArray()));// i + 1, value.Length - (i + 1))));
+ }
+ catch (Exception)
+ {
+ iPart.Add(fPart);
+ iPart.mantissa.Sign = sign;
+ exponent = iPart.exponent;
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
+ return;
+ }
+
+ //Raise or lower 10 to the power of the exponent
+ BigFloat acc = new BigFloat(1, extendedPres);
+ BigFloat temp = new BigFloat(1, extendedPres);
+
+ int powerTemp = exp;
+
+ BigFloat multiplierToUse;
+
+ if (exp < 0)
+ {
+ multiplierToUse = new BigFloat(tenRCP);
+ powerTemp = -exp;
+ }
+ else
+ {
+ multiplierToUse = new BigFloat(ten);
+ }
+
+ //Fast power function
+ while (powerTemp != 0)
+ {
+ temp.Mul(multiplierToUse);
+ multiplierToUse.Assign(temp);
+
+ if ((powerTemp & 1) != 0)
+ {
+ acc.Mul(temp);
+ }
+
+ powerTemp >>= 1;
+ }
+
+ iPart.Add(fPart);
+ iPart.Mul(acc);
+ iPart.mantissa.Sign = sign;
+ exponent = iPart.exponent;
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
+
+ return;
+ }
+
+ iPart.Add(fPart);
+ iPart.mantissa.Sign = sign;
+ exponent = iPart.exponent;
+ if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
+
+ }
+
+ private void Init(PrecisionSpec mantissaPrec)
+ {
+ int mbWords = ((mantissaPrec.NumBits) >> 5);
+ if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
+ int newManBits = mbWords << 5;
+
+ //For efficiency, we just use a 32-bit exponent
+ exponent = 0;
+ mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+ //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
+ }
+
+ //************************** Properties *************************
+
+ ///
+ /// Read-only property. Returns the precision specification of the mantissa.
+ ///
+ /// Floating point numbers are represented as 2^exponent * mantissa, where the
+ /// mantissa and exponent are integers. Note that the exponent in this class is
+ /// always a 32-bit integer. The precision therefore specifies how many bits
+ /// the mantissa will have.
+ ///
+ public PrecisionSpec Precision
+ {
+ get { return mantissa.Precision; }
+ }
+
+ ///
+ /// Writable property:
+ /// true iff the number is negative or in some cases zero (<0)
+ /// false iff the number if positive or in some cases zero (>0)
+ ///
+ public bool Sign
+ {
+ get { return mantissa.Sign; }
+ set { mantissa.Sign = value; }
+ }
+
+ ///
+ /// Read-only property.
+ /// True if the number is NAN, INF_PLUS, INF_MINUS or ZERO
+ /// False if the number has any other value.
+ ///
+ public bool IsSpecialValue
+ {
+ get
+ {
+ return (exponent == Int32.MaxValue || mantissa.IsZero());
+ }
+ }
+
+ ///
+ /// Read-only property, returns the type of number this is. Special values include:
+ ///
+ /// NONE - a regular number
+ /// ZERO - zero
+ /// NAN - Not a Number (some operations will return this if their inputs are out of range)
+ /// INF_PLUS - Positive infinity, not really a number, but a valid input to and output of some functions.
+ /// INF_MINUS - Negative infinity, not really a number, but a valid input to and output of some functions.
+ ///
+ public SpecialValueType SpecialValue
+ {
+ get
+ {
+ if (exponent == Int32.MaxValue)
+ {
+ if (mantissa.IsZero())
+ {
+ if (mantissa.Sign) return SpecialValueType.INF_MINUS;
+ return SpecialValueType.INF_PLUS;
+ }
+
+ return SpecialValueType.NAN;
+ }
+ else
+ {
+ if (mantissa.IsZero()) return SpecialValueType.ZERO;
+ return SpecialValueType.NONE;
+ }
+ }
+ }
+
+ //******************** Mathematical Constants *******************
+
+ ///
+ /// Gets pi to the indicated precision
+ ///
+ /// The precision to perform the calculation to
+ /// pi (the ratio of the area of a circle to its diameter)
+ public static BigFloat GetPi(PrecisionSpec precision)
+ {
+ if (pi == null || precision.NumBits <= pi.mantissa.Precision.NumBits)
+ {
+ CalculatePi(precision.NumBits);
+ }
+
+ BigFloat ret = new BigFloat (precision);
+ ret.Assign(pi);
+
+ return ret;
+ }
+
+ ///
+ /// Get e to the indicated precision
+ ///
+ /// The preicision to perform the calculation to
+ /// e (the number for which the d/dx(e^x) = e^x)
+ public static BigFloat GetE(PrecisionSpec precision)
+ {
+ if (eCache == null || eCache.mantissa.Precision.NumBits < precision.NumBits)
+ {
+ CalculateEOnly(precision.NumBits);
+ //CalculateFactorials(precision.NumBits);
+ }
+
+ BigFloat ret = new BigFloat(precision);
+ ret.Assign(eCache);
+
+ return ret;
+ }
+
+
+ //******************** Arithmetic Functions ********************
+
+ ///
+ /// Addition (this = this + n2)
+ ///
+ /// The number to add
+ public void Add(BigFloat n2)
+ {
+ if (SpecialValueAddTest(n2)) return;
+
+ if (scratch.Precision.NumBits != n2.mantissa.Precision.NumBits)
+ {
+ scratch = new BigInt(n2.mantissa.Precision);
+ }
+
+ if (exponent <= n2.exponent)
+ {
+ int diff = n2.exponent - exponent;
+ exponent = n2.exponent;
+
+ if (diff != 0)
+ {
+ mantissa.RSH(diff);
+ }
+
+ uint carry = mantissa.Add(n2.mantissa);
+
+ if (carry != 0)
+ {
+ mantissa.RSH(1);
+ mantissa.SetBit(mantissa.Precision.NumBits - 1);
+ exponent++;
+ }
+
+ exponent -= mantissa.Normalise();
+ }
+ else
+ {
+ int diff = exponent - n2.exponent;
+
+ scratch.Assign(n2.mantissa);
+ scratch.RSH(diff);
+
+ uint carry = scratch.Add(mantissa);
+
+ if (carry != 0)
+ {
+ scratch.RSH(1);
+ scratch.SetBit(mantissa.Precision.NumBits - 1);
+ exponent++;
+ }
+
+ mantissa.Assign(scratch);
+
+ exponent -= mantissa.Normalise();
+ }
+ }
+
+ ///
+ /// Subtraction (this = this - n2)
+ ///
+ /// The number to subtract from this
+ public void Sub(BigFloat n2)
+ {
+ n2.mantissa.Sign = !n2.mantissa.Sign;
+ Add(n2);
+ n2.mantissa.Sign = !n2.mantissa.Sign;
+ }
+
+ ///
+ /// Multiplication (this = this * n2)
+ ///
+ /// The number to multiply this by
+ public void Mul(BigFloat n2)
+ {
+ if (SpecialValueMulTest(n2)) return;
+
+ //Anything times 0 = 0
+ if (n2.mantissa.IsZero())
+ {
+ mantissa.Assign(n2.mantissa);
+ exponent = 0;
+ return;
+ }
+
+ mantissa.MulHi(n2.mantissa);
+ int shift = mantissa.Normalise();
+ exponent = exponent + n2.exponent + 1 - shift;
+ }
+
+ ///
+ /// Division (this = this / n2)
+ ///
+ /// The number to divide this by
+ public void Div(BigFloat n2)
+ {
+ if (SpecialValueDivTest(n2)) return;
+
+ if (mantissa.Precision.NumBits >= 8192)
+ {
+ BigFloat rcp = n2.Reciprocal();
+ Mul(rcp);
+ }
+ else
+ {
+ int shift = mantissa.DivAndShift(n2.mantissa);
+ exponent = exponent - (n2.exponent + shift);
+ }
+ }
+
+ ///
+ /// Multiply by a power of 2 (-ve implies division)
+ ///
+ ///
+ public void MulPow2(int pow2)
+ {
+ exponent += pow2;
+ }
+
+ ///
+ /// Division-based reciprocal, fastest for small precisions up to 15,000 bits.
+ ///
+ /// The reciprocal 1/this
+ public BigFloat Reciprocal()
+ {
+ if (mantissa.Precision.NumBits >= 8192) return ReciprocalNewton();
+
+ BigFloat reciprocal = new BigFloat(1u, mantissa.Precision);
+ reciprocal.Div(this);
+ return reciprocal;
+ }
+
+ ///
+ /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.
+ ///
+ /// The reciprocal 1/this
+ public BigFloat ReciprocalNewton()
+ {
+ if (mantissa.IsZero())
+ {
+ exponent = Int32.MaxValue;
+ return null;
+ }
+
+ bool oldSign = mantissa.Sign;
+ int oldExponent = exponent;
+
+ //Kill exponent for now (will re-institute later)
+ exponent = 0;
+
+ bool topBit = mantissa.IsTopBitOnlyBit();
+
+ PrecisionSpec curPrec = new PrecisionSpec(32, PrecisionSpec.BaseType.BIN);
+
+ BigFloat reciprocal = new BigFloat(curPrec);
+ BigFloat constant2 = new BigFloat(curPrec);
+ BigFloat temp = new BigFloat(curPrec);
+ BigFloat thisPrec = new BigFloat(this, curPrec);
+
+ reciprocal.exponent = 1;
+ reciprocal.mantissa.SetHighDigit(3129112985u);
+
+ constant2.exponent = 1;
+ constant2.mantissa.SetHighDigit(0x80000000u);
+
+ //D is deliberately left negative for all the following operations.
+ thisPrec.mantissa.Sign = true;
+
+ //Initial estimate.
+ reciprocal.Add(thisPrec);
+
+ //mantissa.Sign = false;
+
+ //Shift down into 0.5 < this < 1 range
+ thisPrec.mantissa.RSH(1);
+
+ //Iteration.
+ int accuracyBits = 2;
+ int mantissaBits = mantissa.Precision.NumBits;
+
+ //Each iteration is a pass of newton's method for RCP.
+ //The is a substantial optimisation to be done here...
+ //You can double the number of bits for the calculations
+ //at each iteration, meaning that the whole process only
+ //takes some constant multiplier of the time for the
+ //full-scale multiplication.
+ while (accuracyBits < mantissaBits)
+ {
+ //Increase the precision as needed
+ if (accuracyBits >= curPrec.NumBits / 2)
+ {
+ int newBits = curPrec.NumBits * 2;
+ if (newBits > mantissaBits) newBits = mantissaBits;
+ curPrec = new PrecisionSpec(newBits, PrecisionSpec.BaseType.BIN);
+
+ reciprocal = new BigFloat(reciprocal, curPrec);
+
+ constant2 = new BigFloat(curPrec);
+ constant2.exponent = 1;
+ constant2.mantissa.SetHighDigit(0x80000000u);
+
+ temp = new BigFloat(temp, curPrec);
+
+ thisPrec = new BigFloat(this, curPrec);
+ thisPrec.mantissa.Sign = true;
+ thisPrec.mantissa.RSH(1);
+ }
+
+ //temp = Xn
+ temp.exponent = reciprocal.exponent;
+ temp.mantissa.Assign(reciprocal.mantissa);
+ //temp = -Xn * D
+ temp.Mul(thisPrec);
+ //temp = -Xn * D + 2 (= 2 - Xn * D)
+ temp.Add(constant2);
+ //reciprocal = X(n+1) = Xn * (2 - Xn * D)
+ reciprocal.Mul(temp);
+
+ accuracyBits *= 2;
+ }
+
+ //'reciprocal' is now the reciprocal of the shifted down, zero-exponent mantissa of 'this'
+ //Restore the mantissa.
+ //mantissa.LSH(1);
+ exponent = oldExponent;
+ //mantissa.Sign = oldSign;
+
+ if (topBit)
+ {
+ reciprocal.exponent = -(oldExponent);
+ }
+ else
+ {
+ reciprocal.exponent = -(oldExponent + 1);
+ }
+ reciprocal.mantissa.Sign = oldSign;
+
+ return reciprocal;
+ }
+
+ ///
+ /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.
+ ///
+ /// The reciprocal 1/this
+ private BigFloat ReciprocalNewton2()
+ {
+ if (mantissa.IsZero())
+ {
+ exponent = Int32.MaxValue;
+ return null;
+ }
+
+ bool oldSign = mantissa.Sign;
+ int oldExponent = exponent;
+
+ //Kill exponent for now (will re-institute later)
+ exponent = 0;
+
+ BigFloat reciprocal = new BigFloat(mantissa.Precision);
+ BigFloat constant2 = new BigFloat(mantissa.Precision);
+ BigFloat temp = new BigFloat(mantissa.Precision);
+
+ reciprocal.exponent = 1;
+ reciprocal.mantissa.SetHighDigit(3129112985u);
+
+ constant2.exponent = 1;
+ constant2.mantissa.SetHighDigit(0x80000000u);
+
+ //D is deliberately left negative for all the following operations.
+ mantissa.Sign = true;
+
+ //Initial estimate.
+ reciprocal.Add(this);
+
+ //mantissa.Sign = false;
+
+ //Shift down into 0.5 < this < 1 range
+ mantissa.RSH(1);
+
+ //Iteration.
+ int accuracyBits = 2;
+ int mantissaBits = mantissa.Precision.NumBits;
+
+ //Each iteration is a pass of newton's method for RCP.
+ //The is a substantial optimisation to be done here...
+ //You can double the number of bits for the calculations
+ //at each iteration, meaning that the whole process only
+ //takes some constant multiplier of the time for the
+ //full-scale multiplication.
+ while (accuracyBits < mantissaBits)
+ {
+ //temp = Xn
+ temp.exponent = reciprocal.exponent;
+ temp.mantissa.Assign(reciprocal.mantissa);
+ //temp = -Xn * D
+ temp.Mul(this);
+ //temp = -Xn * D + 2 (= 2 - Xn * D)
+ temp.Add(constant2);
+ //reciprocal = X(n+1) = Xn * (2 - Xn * D)
+ reciprocal.Mul(temp);
+
+ accuracyBits *= 2;
+ }
+
+ //'reciprocal' is now the reciprocal of the shifted down, zero-exponent mantissa of 'this'
+ //Restore the mantissa.
+ mantissa.LSH(1);
+ exponent = oldExponent;
+ mantissa.Sign = oldSign;
+
+ reciprocal.exponent = -(oldExponent + 1);
+ reciprocal.mantissa.Sign = oldSign;
+
+ return reciprocal;
+ }
+
+ ///
+ /// Sets this equal to the input
+ ///
+ ///
+ public void Assign(BigFloat n2)
+ {
+ exponent = n2.exponent;
+ if (mantissa.AssignHigh(n2.mantissa)) exponent++;
+ }
+
+
+ //********************* Comparison Functions *******************
+
+ ///
+ /// Greater than comparison
+ ///
+ /// the number to compare this to
+ /// true iff this is greater than n2 (this > n2)
+ public bool GreaterThan(BigFloat n2)
+ {
+ if (IsSpecialValue || n2.IsSpecialValue)
+ {
+ SpecialValueType s1 = SpecialValue;
+ SpecialValueType s2 = SpecialValue;
+
+ if (s1 == SpecialValueType.NAN || s2 == SpecialValueType.NAN) return false;
+ if (s1 == SpecialValueType.INF_MINUS) return false;
+ if (s2 == SpecialValueType.INF_PLUS) return false;
+ if (s1 == SpecialValueType.INF_PLUS) return true;
+ if (s2 == SpecialValueType.INF_MINUS) return true;
+
+ if (s1 == SpecialValueType.ZERO)
+ {
+ if (s2 != SpecialValueType.ZERO && n2.Sign)
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ }
+
+ if (s2 == SpecialValueType.ZERO)
+ {
+ return !Sign;
+ }
+ }
+
+ if (!mantissa.Sign && n2.mantissa.Sign) return true;
+ if (mantissa.Sign && !n2.mantissa.Sign) return false;
+ if (!mantissa.Sign)
+ {
+ if (exponent > n2.exponent) return true;
+ if (exponent < n2.exponent) return false;
+ }
+ if (mantissa.Sign)
+ {
+ if (exponent > n2.exponent) return false;
+ if (exponent < n2.exponent) return true;
+ }
+
+ return mantissa.GreaterThan(n2.mantissa);
+ }
+
+ ///
+ /// Less than comparison
+ ///
+ /// the number to compare this to
+ /// true iff this is less than n2 (this < n2)
+ public bool LessThan(BigFloat n2)
+ {
+ if (IsSpecialValue || n2.IsSpecialValue)
+ {
+ SpecialValueType s1 = SpecialValue;
+ SpecialValueType s2 = SpecialValue;
+
+ if (s1 == SpecialValueType.NAN || s2 == SpecialValueType.NAN) return false;
+ if (s1 == SpecialValueType.INF_PLUS) return false;
+ if (s2 == SpecialValueType.INF_PLUS) return true;
+ if (s2 == SpecialValueType.INF_MINUS) return false;
+ if (s1 == SpecialValueType.INF_MINUS) return true;
+
+ if (s1 == SpecialValueType.ZERO)
+ {
+ if (s2 != SpecialValueType.ZERO && !n2.Sign)
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ }
+
+ if (s2 == SpecialValueType.ZERO)
+ {
+ return Sign;
+ }
+ }
+
+ if (!mantissa.Sign && n2.mantissa.Sign) return false;
+ if (mantissa.Sign && !n2.mantissa.Sign) return true;
+ if (!mantissa.Sign)
+ {
+ if (exponent > n2.exponent) return false;
+ if (exponent < n2.exponent) return true;
+ }
+ if (mantissa.Sign)
+ {
+ if (exponent > n2.exponent) return true;
+ if (exponent < n2.exponent) return false;
+ }
+
+ return mantissa.LessThan(n2.mantissa);
+ }
+
+ ///
+ /// Greater than comparison
+ ///
+ /// the number to compare this to
+ /// true iff this is greater than n2 (this > n2)
+ public bool GreaterThan(int i)
+ {
+ BigFloat integer = new BigFloat(i, mantissa.Precision);
+ return GreaterThan(integer);
+ }
+
+ ///
+ /// Less than comparison
+ ///
+ /// the number to compare this to
+ /// true iff this is less than n2 (this < n2)
+ public bool LessThan(int i)
+ {
+ BigFloat integer = new BigFloat(i, mantissa.Precision);
+ return LessThan(integer);
+ }
+
+ ///
+ /// Compare to zero
+ ///
+ /// true if this is zero (this == 0)
+ public bool IsZero()
+ {
+ return (mantissa.IsZero());
+ }
+
+
+ //******************** Mathematical Functions ******************
+
+ ///
+ /// Sets the number to the biggest integer numerically closer to zero, if possible.
+ ///
+ public void Floor()
+ {
+ //Already an integer.
+ if (exponent >= mantissa.Precision.NumBits) return;
+
+ if (exponent < 0)
+ {
+ mantissa.ZeroBits(mantissa.Precision.NumBits);
+ exponent = 0;
+ return;
+ }
+
+ mantissa.ZeroBits(mantissa.Precision.NumBits - (exponent + 1));
+ }
+
+ ///
+ /// Sets the number to its fractional component (equivalent to 'this' - (int)'this')
+ ///
+ public void FPart()
+ {
+ //Already fractional
+ if (exponent < 0)
+ {
+ return;
+ }
+
+ //Has no fractional part
+ if (exponent >= mantissa.Precision.NumBits)
+ {
+ mantissa.Zero();
+ exponent = 0;
+ return;
+ }
+
+ mantissa.ZeroBitsHigh(exponent + 1);
+ exponent -= mantissa.Normalise();
+ }
+
+ ///
+ /// Calculates tan(x)
+ ///
+ public void Tan()
+ {
+ if (IsSpecialValue)
+ {
+ //Tan(x) has no limit as x->inf
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
+ {
+ SetNaN();
+ }
+ else if (SpecialValue == SpecialValueType.ZERO)
+ {
+ SetZero();
+ }
+
+ return;
+ }
+
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
+ {
+ CalculatePi(mantissa.Precision.NumBits);
+ }
+
+ //Work out the sign change (involves replicating some rescaling).
+ bool sign = mantissa.Sign;
+ mantissa.Sign = false;
+
+ if (mantissa.IsZero())
+ {
+ return;
+ }
+
+ //Rescale into 0 <= x < pi
+ if (GreaterThan(pi))
+ {
+ //There will be an inherent loss of precision doing this.
+ BigFloat newAngle = new BigFloat(this);
+ newAngle.Mul(piRecip);
+ newAngle.FPart();
+ newAngle.Mul(pi);
+ Assign(newAngle);
+ }
+
+ //Rescale to -pi/2 <= x < pi/2
+ if (!LessThan(piBy2))
+ {
+ Sub(pi);
+ }
+
+ //Now the sign of the sin determines the sign of the tan.
+ //tan(x) = sin(x) / sqrt(1 - sin^2(x))
+ Sin();
+ BigFloat denom = new BigFloat(this);
+ denom.Mul(this);
+ denom.Sub(new BigFloat(1, mantissa.Precision));
+ denom.mantissa.Sign = !denom.mantissa.Sign;
+
+ if (denom.mantissa.Sign)
+ {
+ denom.SetZero();
+ }
+
+ denom.Sqrt();
+ Div(denom);
+ if (sign) mantissa.Sign = !mantissa.Sign;
+ }
+
+ ///
+ /// Calculates Cos(x)
+ ///
+ public void Cos()
+ {
+ if (IsSpecialValue)
+ {
+ //Cos(x) has no limit as x->inf
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
+ {
+ SetNaN();
+ }
+ else if (SpecialValue == SpecialValueType.ZERO)
+ {
+ Assign(new BigFloat(1, mantissa.Precision));
+ }
+
+ return;
+ }
+
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
+ {
+ CalculatePi(mantissa.Precision.NumBits);
+ }
+
+ Add(piBy2);
+ Sin();
+ }
+
+ ///
+ /// Calculates Sin(x):
+ /// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].
+ ///
+ public void Sin()
+ {
+ if (IsSpecialValue)
+ {
+ //Sin(x) has no limit as x->inf
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
+ {
+ SetNaN();
+ }
+
+ return;
+ }
+
+ //Convert to positive range (0 <= x < inf)
+ bool sign = mantissa.Sign;
+ mantissa.Sign = false;
+
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
+ {
+ CalculatePi(mantissa.Precision.NumBits);
+ }
+
+ if (inverseFactorialCache == null || invFactorialCutoff != mantissa.Precision.NumBits)
+ {
+ CalculateFactorials(mantissa.Precision.NumBits);
+ }
+
+ //Rescale into 0 <= x < 2*pi
+ if (GreaterThan(twoPi))
+ {
+ //There will be an inherent loss of precision doing this.
+ BigFloat newAngle = new BigFloat(this);
+ newAngle.Mul(twoPiRecip);
+ newAngle.FPart();
+ newAngle.Mul(twoPi);
+ Assign(newAngle);
+ }
+
+ //Rescale into range 0 <= x < pi
+ if (GreaterThan(pi))
+ {
+ //sin(pi + a) = sin(pi)cos(a) + sin(a)cos(pi) = 0 - sin(a) = -sin(a)
+ Sub(pi);
+ sign = !sign;
+ }
+
+ BigFloat temp = new BigFloat(mantissa.Precision);
+
+ //Rescale into range 0 <= x < pi/2
+ if (GreaterThan(piBy2))
+ {
+ temp.Assign(this);
+ Assign(pi);
+ Sub(temp);
+ }
+
+ //Rescale into range 0 <= x < pi/6 to accelerate convergence.
+ //This is done using sin(3x) = 3sin(x) - 4sin^3(x)
+ Mul(threeRecip);
+
+ if (mantissa.IsZero())
+ {
+ exponent = 0;
+ return;
+ }
+
+ BigFloat term = new BigFloat(this);
+
+ BigFloat square = new BigFloat(this);
+ square.Mul(term);
+
+ BigFloat sum = new BigFloat(this);
+
+ bool termSign = true;
+ int length = inverseFactorialCache.Length;
+ int numBits = mantissa.Precision.NumBits;
+
+ for (int i = 3; i < length; i += 2)
+ {
+ term.Mul(square);
+ temp.Assign(inverseFactorialCache[i]);
+ temp.Mul(term);
+ temp.mantissa.Sign = termSign;
+ termSign = !termSign;
+
+ if (temp.exponent < -numBits) break;
+
+ sum.Add(temp);
+ }
+
+ //Restore the triple-angle: sin(3x) = 3sin(x) - 4sin^3(x)
+ Assign(sum);
+ sum.Mul(this);
+ sum.Mul(this);
+ Mul(new BigFloat(3, mantissa.Precision));
+ sum.exponent += 2;
+ Sub(sum);
+
+ //Restore the sign
+ mantissa.Sign = sign;
+ }
+
+ ///
+ /// Hyperbolic Sin (sinh) function
+ ///
+ public void Sinh()
+ {
+ if (IsSpecialValue)
+ {
+ return;
+ }
+
+ Exp();
+ Sub(Reciprocal());
+ exponent--;
+ }
+
+ ///
+ /// Hyperbolic cosine (cosh) function
+ ///
+ public void Cosh()
+ {
+ if (IsSpecialValue)
+ {
+ if (SpecialValue == SpecialValueType.ZERO)
+ {
+ Assign(new BigFloat(1, mantissa.Precision));
+ }
+ else if (SpecialValue == SpecialValueType.INF_MINUS)
+ {
+ SetInfPlus();
+ }
+
+ return;
+ }
+
+ Exp();
+ Add(Reciprocal());
+ exponent--;
+ }
+
+ ///
+ /// Hyperbolic tangent function (tanh)
+ ///
+ public void Tanh()
+ {
+ if (IsSpecialValue)
+ {
+ if (SpecialValue == SpecialValueType.INF_MINUS)
+ {
+ Assign(new BigFloat(-1, mantissa.Precision));
+ }
+ else if (SpecialValue == SpecialValueType.INF_PLUS)
+ {
+ Assign(new BigFloat(1, mantissa.Precision));
+ }
+
+ return;
+ }
+
+ exponent++;
+ Exp();
+ BigFloat temp = new BigFloat(this);
+ BigFloat one = new BigFloat(1, mantissa.Precision);
+ temp.Add(one);
+ Sub(one);
+ Div(temp);
+ }
+
+ ///
+ /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)
+ ///
+ public void Arcsin()
+ {
+ if (IsSpecialValue)
+ {
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)
+ {
+ SetNaN();
+ return;
+ }
+
+ return;
+ }
+
+ BigFloat one = new BigFloat(1, mantissa.Precision);
+ BigFloat plusABit = new BigFloat(1, mantissa.Precision);
+ plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));
+ BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);
+ onePlusABit.Add(plusABit);
+
+ bool sign = mantissa.Sign;
+ mantissa.Sign = false;
+
+ if (GreaterThan(onePlusABit))
+ {
+ SetNaN();
+ }
+ else if (LessThan(one))
+ {
+ BigFloat temp = new BigFloat(this);
+ temp.Mul(this);
+ temp.Sub(one);
+ temp.mantissa.Sign = !temp.mantissa.Sign;
+ temp.Sqrt();
+ temp.Add(one);
+ Div(temp);
+ Arctan();
+ exponent++;
+ mantissa.Sign = sign;
+ }
+ else
+ {
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
+ {
+ CalculatePi(mantissa.Precision.NumBits);
+ }
+
+ Assign(piBy2);
+ if (sign) mantissa.Sign = true;
+ }
+ }
+
+ ///
+ /// arccos(): the inverse function of cos(), range (0..pi)
+ ///
+ public void Arccos()
+ {
+ if (IsSpecialValue)
+ {
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)
+ {
+ SetNaN();
+ }
+ else if (SpecialValue == SpecialValueType.ZERO)
+ {
+ Assign(new BigFloat(1, mantissa.Precision));
+ exponent = 0;
+ Sign = false;
+ }
+
+ return;
+ }
+
+ BigFloat one = new BigFloat(1, mantissa.Precision);
+ BigFloat plusABit = new BigFloat(1, mantissa.Precision);
+ plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));
+ BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);
+ onePlusABit.Add(plusABit);
+
+ bool sign = mantissa.Sign;
+ mantissa.Sign = false;
+
+ if (GreaterThan(onePlusABit))
+ {
+ SetNaN();
+ }
+ else if (LessThan(one))
+ {
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
+ {
+ CalculatePi(mantissa.Precision.NumBits);
+ }
+
+ mantissa.Sign = sign;
+ BigFloat temp = new BigFloat(this);
+ Mul(temp);
+ Sub(one);
+ mantissa.Sign = !mantissa.Sign;
+ Sqrt();
+ temp.Add(one);
+ Div(temp);
+ Arctan();
+ exponent++;
+ }
+ else
+ {
+ if (sign)
+ {
+ if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
+ {
+ CalculatePi(mantissa.Precision.NumBits);
+ }
+
+ Assign(pi);
+ }
+ else
+ {
+ mantissa.Zero();
+ exponent = 0;
+ }
+ }
+ }
+
+ ///
+ /// arctan(): the inverse function of sin(), range of (-pi/2..pi/2)
+ ///
+ public void Arctan()
+ {
+ //With 2 argument reductions, we increase precision by a minimum of 4 bits per term.
+ int numBits = mantissa.Precision.NumBits;
+ int maxTerms = numBits >> 2;
+
+ if (pi == null || pi.mantissa.Precision.NumBits != numBits)
+ {
+ CalculatePi(mantissa.Precision.NumBits);
+ }
+
+ //Make domain positive
+ bool sign = mantissa.Sign;
+ mantissa.Sign = false;
+
+ if (IsSpecialValue)
+ {
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
+ {
+ Assign(piBy2);
+ mantissa.Sign = sign;
+ return;
+ }
+
+ return;
+ }
+
+ if (reciprocals == null || reciprocals[0].mantissa.Precision.NumBits != numBits || reciprocals.Length < maxTerms)
+ {
+ CalculateReciprocals(numBits, maxTerms);
+ }
+
+ bool invert = false;
+ BigFloat one = new BigFloat(1, mantissa.Precision);
+
+ //Invert if outside of convergence
+ if (GreaterThan(one))
+ {
+ invert = true;
+ Assign(Reciprocal());
+ }
+
+ //Reduce using half-angle formula:
+ //arctan(2x) = 2 arctan (x / (1 + sqrt(1 + x)))
+
+ //First reduction (guarantees 2 bits per iteration)
+ BigFloat temp = new BigFloat(this);
+ temp.Mul(this);
+ temp.Add(one);
+ temp.Sqrt();
+ temp.Add(one);
+ this.Div(temp);
+
+ //Second reduction (guarantees 4 bits per iteration)
+ temp.Assign(this);
+ temp.Mul(this);
+ temp.Add(one);
+ temp.Sqrt();
+ temp.Add(one);
+ this.Div(temp);
+
+ //Actual series calculation
+ int length = reciprocals.Length;
+ BigFloat term = new BigFloat(this);
+
+ //pow = x^2
+ BigFloat pow = new BigFloat(this);
+ pow.Mul(this);
+
+ BigFloat sum = new BigFloat(this);
+
+ for (int i = 1; i < length; i++)
+ {
+ //u(n) = u(n-1) * x^2
+ //t(n) = u(n) / (2n+1)
+ term.Mul(pow);
+ term.Sign = !term.Sign;
+ temp.Assign(term);
+ temp.Mul(reciprocals[i]);
+
+ if (temp.exponent < -numBits) break;
+
+ sum.Add(temp);
+ }
+
+ //Undo the reductions.
+ Assign(sum);
+ exponent += 2;
+
+ if (invert)
+ {
+ //Assign(Reciprocal());
+ mantissa.Sign = true;
+ Add(piBy2);
+ }
+
+ if (sign)
+ {
+ mantissa.Sign = sign;
+ }
+ }
+
+ ///
+ /// Arcsinh(): the inverse sinh function
+ ///
+ public void Arcsinh()
+ {
+ //Just let all special values fall through
+ if (IsSpecialValue)
+ {
+ return;
+ }
+
+ BigFloat temp = new BigFloat(this);
+ temp.Mul(this);
+ temp.Add(new BigFloat(1, mantissa.Precision));
+ temp.Sqrt();
+ Add(temp);
+ Log();
+ }
+
+ ///
+ /// Arccosh(): the inverse cosh() function
+ ///
+ public void Arccosh()
+ {
+ //acosh isn't defined for x < 1
+ if (IsSpecialValue)
+ {
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.ZERO)
+ {
+ SetNaN();
+ return;
+ }
+
+ return;
+ }
+
+ BigFloat one = new BigFloat(1, mantissa.Precision);
+ if (LessThan(one))
+ {
+ SetNaN();
+ return;
+ }
+
+ BigFloat temp = new BigFloat(this);
+ temp.Mul(this);
+ temp.Sub(one);
+ temp.Sqrt();
+ Add(temp);
+ Log();
+ }
+
+ ///
+ /// Arctanh(): the inverse tanh function
+ ///
+ public void Arctanh()
+ {
+ //|x| <= 1 for a non-NaN output
+ if (IsSpecialValue)
+ {
+ if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
+ {
+ SetNaN();
+ return;
+ }
+
+ return;
+ }
+
+ BigFloat one = new BigFloat(1, mantissa.Precision);
+ BigFloat plusABit = new BigFloat(1, mantissa.Precision);
+ plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));
+ BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);
+ onePlusABit.Add(plusABit);
+
+ bool sign = mantissa.Sign;
+ mantissa.Sign = false;
+
+ if (GreaterThan(onePlusABit))
+ {
+ SetNaN();
+ }
+ else if (LessThan(one))
+ {
+ BigFloat temp = new BigFloat(this);
+ Add(one);
+ one.Sub(temp);
+ Div(one);
+ Log();
+ exponent--;
+ mantissa.Sign = sign;
+ }
+ else
+ {
+ if (sign)
+ {
+ SetInfMinus();
+ }
+ else
+ {
+ SetInfPlus();
+ }
+ }
+ }
+
+ ///
+ /// Two-variable iterative square root, taken from
+ /// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method
+ ///
+ public void Sqrt()
+ {
+ if (mantissa.Sign || IsSpecialValue)
+ {
+ if (SpecialValue == SpecialValueType.ZERO)
+ {
+ return;
+ }
+
+ if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)
+ {
+ SetNaN();
+ }
+
+ return;
+ }
+
+ BigFloat temp2;
+ BigFloat temp3 = new BigFloat(mantissa.Precision);
+ BigFloat three = new BigFloat(3, mantissa.Precision);
+
+ int exponentScale = 0;
+
+ //Rescale to 0.5 <= x < 2
+ if (exponent < -1)
+ {
+ int diff = -exponent;
+ if ((diff & 1) != 0)
+ {
+ diff--;
+ }
+
+ exponentScale = -diff;
+ exponent += diff;
+ }
+ else if (exponent > 0)
+ {
+ if ((exponent & 1) != 0)
+ {
+ exponentScale = exponent + 1;
+ exponent = -1;
+ }
+ else
+ {
+ exponentScale = exponent;
+ exponent = 0;
+ }
+ }
+
+ temp2 = new BigFloat(this);
+ temp2.Sub(new BigFloat(1, mantissa.Precision));
+
+ //if (temp2.mantissa.IsZero())
+ //{
+ // exponent += exponentScale;
+ // return;
+ //}
+
+ int numBits = mantissa.Precision.NumBits;
+
+ while ((exponent - temp2.exponent) < numBits && temp2.SpecialValue != SpecialValueType.ZERO)
+ {
+ //a(n+1) = an - an*cn / 2
+ temp3.Assign(this);
+ temp3.Mul(temp2);
+ temp3.MulPow2(-1);
+ this.Sub(temp3);
+
+ //c(n+1) = cn^2 * (cn - 3) / 4
+ temp3.Assign(temp2);
+ temp2.Sub(three);
+ temp2.Mul(temp3);
+ temp2.Mul(temp3);
+ temp2.MulPow2(-2);
+ }
+
+ exponent += (exponentScale >> 1);
+ }
+
+ ///
+ /// The natural logarithm, ln(x)
+ ///
+ public void Log()
+ {
+ if (IsSpecialValue || mantissa.Sign)
+ {
+ if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)
+ {
+ SetNaN();
+ }
+ else if (SpecialValue == SpecialValueType.ZERO)
+ {
+ SetInfMinus();
+ }
+
+ return;
+ }
+
+ if (mantissa.Precision.NumBits >= 512)
+ {
+ LogAGM1();
+ return;
+ }
+
+ //Compute ln2.
+ if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
+ {
+ CalculateLog2(mantissa.Precision.NumBits);
+ }
+
+ Log2();
+ Mul(ln2cache);
+ }
+
+ ///
+ /// Log to the base 10
+ ///
+ public void Log10()
+ {
+ if (IsSpecialValue || mantissa.Sign)
+ {
+ if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)
+ {
+ SetNaN();
+ }
+ else if (SpecialValue == SpecialValueType.ZERO)
+ {
+ SetInfMinus();
+ }
+
+ return;
+ }
+
+ //Compute ln2.
+ if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
+ {
+ CalculateLog2(mantissa.Precision.NumBits);
+ }
+
+ Log();
+ Mul(log10recip);
+ }
+
+ ///
+ /// The exponential function. Less accurate for high exponents, scales poorly with the number
+ /// of bits.
+ ///
+ public void Exp()
+ {
+ Exp(mantissa.Precision.NumBits);
+ }
+
+ ///
+ /// Raises a number to an integer power (positive or negative). This is a very accurate and fast function,
+ /// comparable to or faster than division (although it is slightly slower for
+ /// negative powers, obviously)
+ ///
+ ///
+ ///
+ public void Pow(int power)
+ {
+ BigFloat acc = new BigFloat(1, mantissa.Precision);
+ BigFloat temp = new BigFloat(1, mantissa.Precision);
+
+ int powerTemp = power;
+
+ if (power < 0)
+ {
+ Assign(Reciprocal());
+ powerTemp = -power;
+ }
+
+ //Fast power function
+ while (powerTemp != 0)
+ {
+ temp.Mul(this);
+ Assign(temp);
+
+ if ((powerTemp & 1) != 0)
+ {
+ acc.Mul(temp);
+ }
+
+ powerTemp >>= 1;
+ }
+
+ Assign(acc);
+ }
+
+ ///
+ /// Raises to an aribitrary power. This is both slow (uses Log) and inaccurate. If you need to
+ /// raise e^x use exp(). If you need an integer power, use the integer power function Pow(int)
+ /// Accuracy Note:
+ /// The function is only ever accurate to a maximum of 4 decimal digits
+ /// For every 10x larger (or smaller) the power gets, you lose an additional decimal digit
+ /// If you really need a precise result, do the calculation with an extra 32-bits and round
+ /// Domain Note:
+ /// This only works for powers of positive real numbers. Negative numbers will fail.
+ ///
+ ///
+ public void Pow(BigFloat power)
+ {
+ Log();
+ Mul(power);
+ Exp();
+ }
+
+
+ //******************** Static Math Functions *******************
+
+ ///
+ /// Returns the integer component of the input
+ ///
+ /// The input number
+ /// The integer component returned will always be numerically closer to zero
+ /// than the input: an input of -3.49 for instance would produce a value of 3.
+ public static BigFloat Floor(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Floor();
+ return n1;
+ }
+
+ ///
+ /// Returns the fractional (non-integer component of the input)
+ ///
+ /// The input number
+ public static BigFloat FPart(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.FPart();
+ return n1;
+ }
+
+ ///
+ /// Calculates tan(x)
+ ///
+ /// The angle (in radians) to find the tangent of
+ public static BigFloat Tan(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Tan();
+ return n1;
+ }
+
+ ///
+ /// Calculates Cos(x)
+ ///
+ /// The angle (in radians) to find the cosine of
+ /// This is a reasonably fast function for smaller precisions, but
+ /// doesn't scale well for higher precision arguments
+ public static BigFloat Cos(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Cos();
+ return n1;
+ }
+
+ ///
+ /// Calculates Sin(x):
+ /// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].
+ ///
+ /// The angle to find the sine of (in radians)
+ /// This is a resonably fast function, for smaller precision arguments, but doesn't
+ /// scale very well with the number of bits in the input.
+ public static BigFloat Sin(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Sin();
+ return n1;
+ }
+
+ ///
+ /// Hyperbolic Sin (sinh) function
+ ///
+ /// The number to find the hyperbolic sine of
+ public static BigFloat Sinh(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Sinh();
+ return n1;
+ }
+
+ ///
+ /// Hyperbolic cosine (cosh) function
+ ///
+ /// The number to find the hyperbolic cosine of
+ public static BigFloat Cosh(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Cosh();
+ return n1;
+ }
+
+ ///
+ /// Hyperbolic tangent function (tanh)
+ ///
+ /// The number to find the hyperbolic tangent of
+ public static BigFloat Tanh(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Tanh();
+ return n1;
+ }
+
+ ///
+ /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)
+ ///
+ /// The number to find the arcsine of (-pi/2..pi/2)
+ /// Note that inverse trig functions are only defined within a specific range.
+ /// Values outside this range will return NaN, although some margin for error is assumed.
+ ///
+ public static BigFloat Arcsin(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Arcsin();
+ return n1;
+ }
+
+ ///
+ /// arccos(): the inverse function of cos(), input range (0..pi)
+ ///
+ /// The number to find the arccosine of (0..pi)
+ /// Note that inverse trig functions are only defined within a specific range.
+ /// Values outside this range will return NaN, although some margin for error is assumed.
+ ///
+ public static BigFloat Arccos(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Arccos();
+ return n1;
+ }
+
+ ///
+ /// arctan(): the inverse function of sin(), input range of (-pi/2..pi/2)
+ ///
+ /// The number to find the arctangent of (-pi/2..pi/2)
+ /// Note that inverse trig functions are only defined within a specific range.
+ /// Values outside this range will return NaN, although some margin for error is assumed.
+ ///
+ public static BigFloat Arctan(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Arctan();
+ return n1;
+ }
+
+ ///
+ /// Arcsinh(): the inverse sinh function
+ ///
+ /// The number to find the inverse hyperbolic sine of
+ public static BigFloat Arcsinh(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Arcsinh();
+ return n1;
+ }
+
+ ///
+ /// Arccosh(): the inverse cosh() function
+ ///
+ /// The number to find the inverse hyperbolic cosine of
+ public static BigFloat Arccosh(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Arccosh();
+ return n1;
+ }
+
+ ///
+ /// Arctanh(): the inverse tanh function
+ ///
+ /// The number to fine the inverse hyperbolic tan of
+ public static BigFloat Arctanh(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Arctanh();
+ return n1;
+ }
+
+ ///
+ /// Two-variable iterative square root, taken from
+ /// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method
+ ///
+ /// This is quite a fast function, as elementary functions go. You can expect it to take
+ /// about twice as long as a floating-point division.
+ ///
+ public static BigFloat Sqrt(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Sqrt();
+ return n1;
+ }
+
+ ///
+ /// The natural logarithm, ln(x) (log base e)
+ ///
+ /// This is a very slow function, despite repeated attempts at optimisation.
+ /// To make it any faster, different strategies would be needed for integer operations.
+ /// It does, however, scale well with the number of bits.
+ ///
+ /// The number to find the natural logarithm of
+ public static BigFloat Log(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Log();
+ return n1;
+ }
+
+ ///
+ /// Base 10 logarithm of a number
+ ///
+ /// This is a very slow function, despite repeated attempts at optimisation.
+ /// To make it any faster, different strategies would be needed for integer operations.
+ /// It does, however, scale well with the number of bits.
+ ///
+ /// The number to find the base 10 logarithm of
+ public static BigFloat Log10(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Log10();
+ return n1;
+ }
+
+ ///
+ /// The exponential function. Less accurate for high exponents, scales poorly with the number
+ /// of bits. This is quite fast for low-precision arguments.
+ ///
+ public static BigFloat Exp(BigFloat n1)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Exp();
+ return n1;
+ }
+
+ ///
+ /// Raises a number to an integer power (positive or negative). This is a very accurate and fast function,
+ /// comparable to or faster than division (although it is slightly slower for
+ /// negative powers, obviously).
+ ///
+ /// The number to raise to the power
+ /// The power to raise it to
+ public static BigFloat Pow(BigFloat n1, int power)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Pow(power);
+ return n1;
+ }
+
+ ///
+ /// Raises to an aribitrary power. This is both slow (uses Log) and inaccurate. If you need to
+ /// raise e^x use exp(). If you need an integer power, use the integer power function Pow(int)
+ ///
+ ///
+ /// Accuracy Note:
+ /// The function is only ever accurate to a maximum of 4 decimal digits
+ /// For every 10x larger (or smaller) the power gets, you lose an additional decimal digit
+ /// If you really need a precise result, do the calculation with an extra 32-bits and round
+ ///
+ /// Domain Note:
+ /// This only works for powers of positive real numbers. Negative numbers will fail.
+ ///
+ /// The number to raise to a power
+ /// The power to raise it to
+ public static BigFloat Pow(BigFloat n1, BigFloat power)
+ {
+ BigFloat res = new BigFloat(n1);
+ n1.Pow(power);
+ return n1;
+ }
+
+ //********************** Static functions **********************
+
+ ///
+ /// Adds two numbers and returns the result
+ ///
+ public static BigFloat Add(BigFloat n1, BigFloat n2)
+ {
+ BigFloat ret = new BigFloat(n1);
+ ret.Add(n2);
+ return ret;
+ }
+
+ ///
+ /// Subtracts two numbers and returns the result
+ ///
+ public static BigFloat Sub(BigFloat n1, BigFloat n2)
+ {
+ BigFloat ret = new BigFloat(n1);
+ ret.Sub(n2);
+ return ret;
+ }
+
+ ///
+ /// Multiplies two numbers and returns the result
+ ///
+ public static BigFloat Mul(BigFloat n1, BigFloat n2)
+ {
+ BigFloat ret = new BigFloat(n1);
+ ret.Mul(n2);
+ return ret;
+ }
+
+ ///
+ /// Divides two numbers and returns the result
+ ///
+ public static BigFloat Div(BigFloat n1, BigFloat n2)
+ {
+ BigFloat ret = new BigFloat(n1);
+ ret.Div(n2);
+ return ret;
+ }
+
+ ///
+ /// Tests whether n1 is greater than n2
+ ///
+ public static bool GreaterThan(BigFloat n1, BigFloat n2)
+ {
+ return n1.GreaterThan(n2);
+ }
+
+ ///
+ /// Tests whether n1 is less than n2
+ ///
+ public static bool LessThan(BigFloat n1, BigFloat n2)
+ {
+ return n1.LessThan(n2);
+ }
+
+
+ //******************* Fast static functions ********************
+
+ ///
+ /// Adds two numbers and assigns the result to res.
+ ///
+ /// a pre-existing BigFloat to take the result
+ /// the first number
+ /// the second number
+ /// a handle to res
+ public static BigFloat Add(BigFloat res, BigFloat n1, BigFloat n2)
+ {
+ res.Assign(n1);
+ res.Add(n2);
+ return res;
+ }
+
+ ///
+ /// Subtracts two numbers and assigns the result to res.
+ ///
+ /// a pre-existing BigFloat to take the result
+ /// the first number
+ /// the second number
+ /// a handle to res
+ public static BigFloat Sub(BigFloat res, BigFloat n1, BigFloat n2)
+ {
+ res.Assign(n1);
+ res.Sub(n2);
+ return res;
+ }
+
+ ///
+ /// Multiplies two numbers and assigns the result to res.
+ ///
+ /// a pre-existing BigFloat to take the result
+ /// the first number
+ /// the second number
+ /// a handle to res
+ public static BigFloat Mul(BigFloat res, BigFloat n1, BigFloat n2)
+ {
+ res.Assign(n1);
+ res.Mul(n2);
+ return res;
+ }
+
+ ///
+ /// Divides two numbers and assigns the result to res.
+ ///
+ /// a pre-existing BigFloat to take the result
+ /// the first number
+ /// the second number
+ /// a handle to res
+ public static BigFloat Div(BigFloat res, BigFloat n1, BigFloat n2)
+ {
+ res.Assign(n1);
+ res.Div(n2);
+ return res;
+ }
+
+
+ //************************* Operators **************************
+
+ ///
+ /// The addition operator
+ ///
+ public static BigFloat operator +(BigFloat n1, BigFloat n2)
+ {
+ return Add(n1, n2);
+ }
+
+ ///
+ /// The subtraction operator
+ ///
+ public static BigFloat operator -(BigFloat n1, BigFloat n2)
+ {
+ return Sub(n1, n2);
+ }
+
+ ///
+ /// The multiplication operator
+ ///
+ public static BigFloat operator *(BigFloat n1, BigFloat n2)
+ {
+ return Mul(n1, n2);
+ }
+
+ ///
+ /// The division operator
+ ///
+ public static BigFloat operator /(BigFloat n1, BigFloat n2)
+ {
+ return Div(n1, n2);
+ }
+
+ //************************** Conversions *************************
+
+ ///
+ /// Converts a BigFloat to an BigInt with the specified precision
+ ///
+ /// The number to convert
+ /// The precision to convert it with
+ /// Do we round the number if we are truncating the mantissa?
+ ///
+ public static BigInt ConvertToInt(BigFloat n1, PrecisionSpec precision, bool round)
+ {
+ BigInt ret = new BigInt(precision);
+
+ int numBits = n1.mantissa.Precision.NumBits;
+ int shift = numBits - (n1.exponent + 1);
+
+ BigFloat copy = new BigFloat(n1);
+ bool inc = false;
+
+ //Rounding
+ if (copy.mantissa.Precision.NumBits > ret.Precision.NumBits)
+ {
+ inc = true;
+
+ for (int i = copy.exponent + 1; i <= ret.Precision.NumBits; i++)
+ {
+ if (copy.mantissa.GetBitFromTop(i) == 0)
+ {
+ inc = false;
+ break;
+ }
+ }
+ }
+
+ if (shift > 0)
+ {
+ copy.mantissa.RSH(shift);
+ }
+ else if (shift < 0)
+ {
+ copy.mantissa.LSH(-shift);
+ }
+
+ ret.Assign(copy.mantissa);
+
+ if (inc) ret.Increment();
+
+ return ret;
+ }
+
+ ///
+ /// Returns a base-10 string representing the number.
+ ///
+ /// Note: This is inefficient and possibly inaccurate. Please use with enough
+ /// rounding digits (set using the RoundingDigits property) to ensure accuracy
+ ///
+ public override string ToString()
+ {
+ if (IsSpecialValue)
+ {
+ SpecialValueType s = SpecialValue;
+ if (s == SpecialValueType.ZERO)
+ {
+ return String.Format("0{0}0", System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator);
+ }
+ else if (s == SpecialValueType.INF_PLUS)
+ {
+ return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol;
+ }
+ else if (s == SpecialValueType.INF_MINUS)
+ {
+ return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol;
+ }
+ else if (s == SpecialValueType.NAN)
+ {
+ return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol;
+ }
+ else
+ {
+ return "Unrecognised special type";
+ }
+ }
+
+ if (scratch.Precision.NumBits != mantissa.Precision.NumBits)
+ {
+ scratch = new BigInt(mantissa.Precision);
+ }
+
+ //The mantissa expresses 1.xxxxxxxxxxx
+ //The highest possible value for the mantissa without the implicit 1. is 0.9999999...
+ scratch.Assign(mantissa);
+ //scratch.Round(3);
+ scratch.Sign = false;
+ BigInt denom = new BigInt("0", mantissa.Precision);
+ denom.SetBit(mantissa.Precision.NumBits - 1);
+
+ bool useExponentialNotation = false;
+ int halfBits = mantissa.Precision.NumBits / 2;
+ if (halfBits > 60) halfBits = 60;
+ int precDec = 10;
+
+ if (exponent > 0)
+ {
+ if (exponent < halfBits)
+ {
+ denom.RSH(exponent);
+ }
+ else
+ {
+ useExponentialNotation = true;
+ }
+ }
+ else if (exponent < 0)
+ {
+ int shift = -(exponent);
+ if (shift < precDec)
+ {
+ scratch.RSH(shift);
+ }
+ else
+ {
+ useExponentialNotation = true;
+ }
+ }
+
+ string output;
+
+ if (useExponentialNotation)
+ {
+ int absExponent = exponent;
+ if (absExponent < 0) absExponent = -absExponent;
+ int powerOf10 = (int)((double)absExponent * Math.Log10(2.0));
+
+ //Use 1 extra digit of precision (this is actually 32 bits more, nb)
+ BigFloat thisFloat = new BigFloat(this, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
+ thisFloat.mantissa.Sign = false;
+
+ //Multiplicative correction factor to bring number into range.
+ BigFloat one = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
+ BigFloat ten = new BigFloat(10, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
+ BigFloat tenRCP = ten.Reciprocal();
+
+ //Accumulator for the power of 10 calculation.
+ BigFloat acc = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
+
+ BigFloat tenToUse;
+
+ if (exponent > 0)
+ {
+ tenToUse = new BigFloat(tenRCP, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
+ }
+ else
+ {
+ tenToUse = new BigFloat(ten, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
+ }
+
+ BigFloat tenToPower = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
+
+ int powerTemp = powerOf10;
+
+ //Fast power function
+ while (powerTemp != 0)
+ {
+ tenToPower.Mul(tenToUse);
+ tenToUse.Assign(tenToPower);
+
+ if ((powerTemp & 1) != 0)
+ {
+ acc.Mul(tenToPower);
+ }
+
+ powerTemp >>= 1;
+ }
+
+ thisFloat.Mul(acc);
+
+ //If we are out of range, correct.
+ if (thisFloat.GreaterThan(ten))
+ {
+ thisFloat.Mul(tenRCP);
+ if (exponent > 0)
+ {
+ powerOf10++;
+ }
+ else
+ {
+ powerOf10--;
+ }
+ }
+ else if (thisFloat.LessThan(one))
+ {
+ thisFloat.Mul(ten);
+ if (exponent > 0)
+ {
+ powerOf10--;
+ }
+ else
+ {
+ powerOf10++;
+ }
+ }
+
+ //Restore the precision and the sign.
+ BigFloat printable = new BigFloat(thisFloat, mantissa.Precision);
+ printable.mantissa.Sign = mantissa.Sign;
+ output = printable.ToString();
+
+ if (exponent < 0) powerOf10 = -powerOf10;
+
+ output = String.Format("{0}E{1}", output, powerOf10);
+ }
+ else
+ {
+ BigInt bigDigit = BigInt.Div(scratch, denom);
+ bigDigit.Sign = false;
+ scratch.Sub(BigInt.Mul(denom, bigDigit));
+
+ if (mantissa.Sign)
+ {
+ output = String.Format("-{0}.", bigDigit);
+ }
+ else
+ {
+ output = String.Format("{0}.", bigDigit);
+ }
+
+ denom = BigInt.Div(denom, 10u);
+
+ while (!denom.IsZero())
+ {
+ uint digit = (uint)BigInt.Div(scratch, denom);
+ if (digit == 10) digit--;
+ scratch.Sub(BigInt.Mul(denom, digit));
+ output = String.Format("{0}{1}", output, digit);
+ denom = BigInt.Div(denom, 10u);
+ }
+
+ output = RoundString(output, RoundingDigits);
+ }
+
+ return output;
+ }
+
+ //**************** Special value handling for ops ***************
+
+ private void SetNaN()
+ {
+ exponent = Int32.MaxValue;
+ mantissa.SetBit(mantissa.Precision.NumBits - 1);
+ }
+
+ private void SetZero()
+ {
+ exponent = 0;
+ mantissa.Zero();
+ Sign = false;
+ }
+
+ private void SetInfPlus()
+ {
+ Sign = false;
+ exponent = Int32.MaxValue;
+ mantissa.Zero();
+ }
+
+ private void SetInfMinus()
+ {
+ Sign = true;
+ exponent = Int32.MaxValue;
+ mantissa.Zero();
+ }
+
+ private bool SpecialValueAddTest(BigFloat n2)
+ {
+ if (IsSpecialValue || n2.IsSpecialValue)
+ {
+ SpecialValueType s1 = SpecialValue;
+ SpecialValueType s2 = n2.SpecialValue;
+
+ if (s1 == SpecialValueType.NAN) return true;
+ if (s2 == SpecialValueType.NAN)
+ {
+ //Set NaN and return.
+ SetNaN();
+ return true;
+ }
+
+ if (s1 == SpecialValueType.INF_PLUS)
+ {
+ //INF+ + INF- = NAN
+ if (s2 == SpecialValueType.INF_MINUS)
+ {
+ SetNaN();
+ return true;
+ }
+
+ return true;
+ }
+
+ if (s1 == SpecialValueType.INF_MINUS)
+ {
+ //INF+ + INF- = NAN
+ if (s2 == SpecialValueType.INF_PLUS)
+ {
+ SetNaN();
+ return true;
+ }
+
+ return true;
+ }
+
+ if (s2 == SpecialValueType.ZERO)
+ {
+ return true;
+ }
+
+ if (s1 == SpecialValueType.ZERO)
+ {
+ Assign(n2);
+ return true;
+ }
+ }
+
+ return false;
+ }
+
+ private bool SpecialValueMulTest(BigFloat n2)
+ {
+ if (IsSpecialValue || n2.IsSpecialValue)
+ {
+ SpecialValueType s1 = SpecialValue;
+ SpecialValueType s2 = n2.SpecialValue;
+
+ if (s1 == SpecialValueType.NAN) return true;
+ if (s2 == SpecialValueType.NAN)
+ {
+ //Set NaN and return.
+ SetNaN();
+ return true;
+ }
+
+ if (s1 == SpecialValueType.INF_PLUS)
+ {
+ //Inf+ * Inf- = Inf-
+ if (s2 == SpecialValueType.INF_MINUS)
+ {
+ Assign(n2);
+ return true;
+ }
+
+ //Inf+ * 0 = NaN
+ if (s2 == SpecialValueType.ZERO)
+ {
+ //Set NaN and return.
+ SetNaN();
+ return true;
+ }
+
+ return true;
+ }
+
+ if (s1 == SpecialValueType.INF_MINUS)
+ {
+ //Inf- * Inf- = Inf+
+ if (s2 == SpecialValueType.INF_MINUS)
+ {
+ Sign = false;
+ return true;
+ }
+
+ //Inf- * 0 = NaN
+ if (s2 == SpecialValueType.ZERO)
+ {
+ //Set NaN and return.
+ SetNaN();
+ return true;
+ }
+
+ return true;
+ }
+
+ if (s2 == SpecialValueType.ZERO)
+ {
+ SetZero();
+ return true;
+ }
+
+ if (s1 == SpecialValueType.ZERO)
+ {
+ return true;
+ }
+ }
+
+ return false;
+ }
+
+ private bool SpecialValueDivTest(BigFloat n2)
+ {
+ if (IsSpecialValue || n2.IsSpecialValue)
+ {
+ SpecialValueType s1 = SpecialValue;
+ SpecialValueType s2 = n2.SpecialValue;
+
+ if (s1 == SpecialValueType.NAN) return true;
+ if (s2 == SpecialValueType.NAN)
+ {
+ //Set NaN and return.
+ SetNaN();
+ return true;
+ }
+
+ if ((s1 == SpecialValueType.INF_PLUS || s1 == SpecialValueType.INF_MINUS))
+ {
+ if (s2 == SpecialValueType.INF_PLUS || s2 == SpecialValueType.INF_MINUS)
+ {
+ //Set NaN and return.
+ SetNaN();
+ return true;
+ }
+
+ if (n2.Sign)
+ {
+ if (s1 == SpecialValueType.INF_PLUS)
+ {
+ SetInfMinus();
+ return true;
+ }
+
+ SetInfPlus();
+ return true;
+ }
+
+ //Keep inf
+ return true;
+ }
+
+ if (s2 == SpecialValueType.ZERO)
+ {
+ if (s1 == SpecialValueType.ZERO)
+ {
+ SetNaN();
+ return true;
+ }
+
+ if (Sign)
+ {
+ SetInfMinus();
+ return true;
+ }
+
+ SetInfPlus();
+ return true;
+ }
+ }
+
+ return false;
+ }
+
+ //****************** Internal helper functions *****************
+
+ ///
+ /// Used for fixed point speed-ups (where the extra precision is not required). Note that Denormalised
+ /// floats break the assumptions that underly Add() and Sub(), so they can only be used for multiplication
+ ///
+ ///
+ private void Denormalise(int targetExponent)
+ {
+ int diff = targetExponent - exponent;
+ if (diff <= 0) return;
+
+ //This only works to reduce the precision, so if the difference implies an increase, we can't do anything.
+ mantissa.RSH(diff);
+ exponent += diff;
+ }
+
+ ///
+ /// The binary logarithm, log2(x) - for precisions above 1000 bits, use Log() and convert the base.
+ ///
+ private void Log2()
+ {
+ if (scratch.Precision.NumBits != mantissa.Precision.NumBits)
+ {
+ scratch = new BigInt(mantissa.Precision);
+ }
+
+ int bits = mantissa.Precision.NumBits;
+ BigFloat temp = new BigFloat(this);
+ BigFloat result = new BigFloat(exponent, mantissa.Precision);
+ BigFloat pow2 = new BigFloat(1, mantissa.Precision);
+ temp.exponent = 0;
+ int bitsCalculated = 0;
+
+ while (bitsCalculated < bits)
+ {
+ int i;
+ for (i = 0; (temp.exponent == 0); i++)
+ {
+ temp.mantissa.SquareHiFast(scratch);
+ int shift = temp.mantissa.Normalise();
+ temp.exponent += 1 - shift;
+ if (i + bitsCalculated >= bits) break;
+ }
+
+ pow2.MulPow2(-i);
+ result.Add(pow2);
+ temp.exponent = 0;
+ bitsCalculated += i;
+ }
+
+ this.Assign(result);
+ }
+
+ ///
+ /// Tried the newton method for logs, but the exponential function is too slow to do it.
+ ///
+ private void LogNewton()
+ {
+ if (mantissa.IsZero() || mantissa.Sign)
+ {
+ return;
+ }
+
+ //Compute ln2.
+ if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
+ {
+ CalculateLog2(mantissa.Precision.NumBits);
+ }
+
+ int numBits = mantissa.Precision.NumBits;
+
+ //Use inverse exp function with Newton's method.
+ BigFloat xn = new BigFloat(this);
+ BigFloat oldExponent = new BigFloat(xn.exponent, mantissa.Precision);
+ xn.exponent = 0;
+ this.exponent = 0;
+ //Hack to subtract 1
+ xn.mantissa.ClearBit(numBits - 1);
+ //x0 = (x - 1) * log2 - this is a straight line fit between log(1) = 0 and log(2) = ln2
+ xn.Mul(ln2cache);
+ //x0 = (x - 1) * log2 + C - this corrects for minimum error over the range.
+ xn.Add(logNewtonConstant);
+ BigFloat term = new BigFloat(mantissa.Precision);
+ BigFloat one = new BigFloat(1, mantissa.Precision);
+
+ int precision = 32;
+ int normalPrecision = mantissa.Precision.NumBits;
+
+ int iterations = 0;
+
+ while (true)
+ {
+ term.Assign(xn);
+ term.mantissa.Sign = true;
+ term.Exp(precision);
+ term.Mul(this);
+ term.Sub(one);
+
+ iterations++;
+ if (term.exponent < -((precision >> 1) - 4))
+ {
+ if (precision == normalPrecision)
+ {
+ if (term.exponent < -(precision - 4)) break;
+ }
+ else
+ {
+ precision = precision << 1;
+ if (precision > normalPrecision) precision = normalPrecision;
+ }
+ }
+
+ xn.Add(term);
+ }
+
+ //log(2^n*s) = log(2^n) + log(s) = nlog(2) + log(s)
+ term.Assign(ln2cache);
+ term.Mul(oldExponent);
+
+ this.Assign(xn);
+ this.Add(term);
+ }
+
+ ///
+ /// Log(x) implemented as an Arithmetic-Geometric Mean. Fast for high precisions.
+ ///
+ private void LogAGM1()
+ {
+ if (mantissa.IsZero() || mantissa.Sign)
+ {
+ return;
+ }
+
+ //Compute ln2.
+ if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
+ {
+ CalculateLog2(mantissa.Precision.NumBits);
+ }
+
+ //Compute ln(x) using AGM formula
+
+ //1. Re-write the input as 2^n * (0.5 <= x < 1)
+ int power2 = exponent + 1;
+ exponent = -1;
+
+ //BigFloat res = new BigFloat(firstAGMcache);
+ BigFloat a0 = new BigFloat(1, mantissa.Precision);
+ BigFloat b0 = new BigFloat(pow10cache);
+ b0.Mul(this);
+
+ BigFloat r = R(a0, b0);
+
+ this.Assign(firstAGMcache);
+ this.Sub(r);
+
+ a0.Assign(ln2cache);
+ a0.Mul(new BigFloat(power2, mantissa.Precision));
+ this.Add(a0);
+ }
+
+ private void Exp(int numBits)
+ {
+ if (IsSpecialValue)
+ {
+ if (SpecialValue == SpecialValueType.ZERO)
+ {
+ //e^0 = 1
+ exponent = 0;
+ mantissa.SetHighDigit(0x80000000);
+ }
+ else if (SpecialValue == SpecialValueType.INF_MINUS)
+ {
+ //e^-inf = 0
+ SetZero();
+ }
+
+ return;
+ }
+
+ PrecisionSpec prec = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
+ numBits = prec.NumBits;
+
+ if (scratch.Precision.NumBits != prec.NumBits)
+ {
+ scratch = new BigInt(prec);
+ }
+
+ if (inverseFactorialCache == null || invFactorialCutoff < numBits)
+ {
+ CalculateFactorials(numBits);
+ }
+
+ //let x = 1 * 'this'.mantissa (i.e. 1 <= x < 2)
+ //exp(2^n * x) = e^(2^n * x) = (e^x)^2n = exp(x)^2n
+
+ int oldExponent = 0;
+
+ if (exponent > -4)
+ {
+ oldExponent = exponent + 4;
+ exponent = -4;
+ }
+
+ BigFloat thisSave = new BigFloat(this, prec);
+ BigFloat temp = new BigFloat(1, prec);
+ BigFloat temp2 = new BigFloat(this, prec);
+ BigFloat res = new BigFloat(1, prec);
+ int length = inverseFactorialCache.Length;
+
+ int iterations;
+ for (int i = 1; i < length; i++)
+ {
+ //temp = x^i
+ temp.Mul(thisSave);
+ temp2.Assign(inverseFactorialCache[i]);
+ temp2.Mul(temp);
+
+ if (temp2.exponent < -(numBits + 4)) { iterations = i; break; }
+
+ res.Add(temp2);
+ }
+
+ //res = exp(x)
+ //Now... x^(2^n) = (x^2)^(2^(n - 1))
+ for (int i = 0; i < oldExponent; i++)
+ {
+ res.mantissa.SquareHiFast(scratch);
+ int shift = res.mantissa.Normalise();
+ res.exponent = res.exponent << 1;
+ res.exponent += 1 - shift;
+ }
+
+ //Deal with +/- inf
+ if (res.exponent == Int32.MaxValue)
+ {
+ res.mantissa.Zero();
+ }
+
+ Assign(res);
+ }
+
+ ///
+ /// Calculates ln(2) and returns -10^(n/2 + a bit) for reuse, using the AGM method as described in
+ /// http://lacim.uqam.ca/~plouffe/articles/log2.pdf
+ ///
+ ///
+ ///
+ private static void CalculateLog2(int numBits)
+ {
+ //Use the AGM method formula to get log2 to N digits.
+ //R(a0, b0) = 1 / (1 - Sum(2^-n*(an^2 - bn^2)))
+ //log(1/2) = R(1, 10^-n) - R(1, 10^-n/2)
+ PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
+ PrecisionSpec extendedPres = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);
+ BigFloat a0 = new BigFloat(1, extendedPres);
+ BigFloat b0 = TenPow(-(int)((double)((numBits >> 1) + 2) * 0.302), extendedPres);
+ BigFloat pow10saved = new BigFloat(b0);
+ BigFloat firstAGMcacheSaved = new BigFloat(extendedPres);
+
+ //save power of 10 (in normal precision)
+ pow10cache = new BigFloat(b0, normalPres);
+
+ ln2cache = R(a0, b0);
+
+ //save the first half of the log calculation
+ firstAGMcache = new BigFloat(ln2cache, normalPres);
+ firstAGMcacheSaved.Assign(ln2cache);
+
+ b0.MulPow2(-1);
+ ln2cache.Sub(R(a0, b0));
+
+ //Convert to log(2)
+ ln2cache.mantissa.Sign = false;
+
+ //Save magic constant for newton log
+ //First guess in range 1 <= x < 2 is x0 = ln2 * (x - 1) + C
+ logNewtonConstant = new BigFloat(ln2cache);
+ logNewtonConstant.Mul(new BigFloat(3, extendedPres));
+ logNewtonConstant.exponent--;
+ logNewtonConstant.Sub(new BigFloat(1, extendedPres));
+ logNewtonConstant = new BigFloat(logNewtonConstant, normalPres);
+
+ //Save the inverse.
+ log2ecache = new BigFloat(ln2cache);
+ log2ecache = new BigFloat(log2ecache.Reciprocal(), normalPres);
+
+ //Now cache log10
+ //Because the log functions call this function to the precision to which they
+ //are called, we cannot call them without causing an infinite loop, so we need
+ //to inline the code.
+ log10recip = new BigFloat(10, extendedPres);
+
+ {
+ int power2 = log10recip.exponent + 1;
+ log10recip.exponent = -1;
+
+ //BigFloat res = new BigFloat(firstAGMcache);
+ BigFloat ax = new BigFloat(1, extendedPres);
+ BigFloat bx = new BigFloat(pow10saved);
+ bx.Mul(log10recip);
+
+ BigFloat r = R(ax, bx);
+
+ log10recip.Assign(firstAGMcacheSaved);
+ log10recip.Sub(r);
+
+ ax.Assign(ln2cache);
+ ax.Mul(new BigFloat(power2, log10recip.mantissa.Precision));
+ log10recip.Add(ax);
+ }
+
+ log10recip = log10recip.Reciprocal();
+ log10recip = new BigFloat(log10recip, normalPres);
+
+
+ //Trim to n bits
+ ln2cache = new BigFloat(ln2cache, normalPres);
+ }
+
+ private static BigFloat TenPow(int power, PrecisionSpec precision)
+ {
+ BigFloat acc = new BigFloat(1, precision);
+ BigFloat temp = new BigFloat(1, precision);
+
+ int powerTemp = power;
+
+ BigFloat multiplierToUse = new BigFloat(10, precision);
+
+ if (power < 0)
+ {
+ multiplierToUse = multiplierToUse.Reciprocal();
+ powerTemp = -power;
+ }
+
+ //Fast power function
+ while (powerTemp != 0)
+ {
+ temp.Mul(multiplierToUse);
+ multiplierToUse.Assign(temp);
+
+ if ((powerTemp & 1) != 0)
+ {
+ acc.Mul(temp);
+ }
+
+ powerTemp >>= 1;
+ }
+
+ return acc;
+ }
+
+ private static BigFloat R(BigFloat a0, BigFloat b0)
+ {
+ //Precision extend taken out.
+ int bits = a0.mantissa.Precision.NumBits;
+ PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);
+ BigFloat an = new BigFloat(a0, extendedPres);
+ BigFloat bn = new BigFloat(b0, extendedPres);
+ BigFloat sum = new BigFloat(extendedPres);
+ BigFloat term = new BigFloat(extendedPres);
+ BigFloat temp1 = new BigFloat(extendedPres);
+ BigFloat one = new BigFloat(1, extendedPres);
+
+ int iteration = 0;
+
+ for (iteration = 0; ; iteration++)
+ {
+ //Get the sum term for this iteration.
+ term.Assign(an);
+ term.Mul(an);
+ temp1.Assign(bn);
+ temp1.Mul(bn);
+ //term = an^2 - bn^2
+ term.Sub(temp1);
+ //term = 2^(n-1) * (an^2 - bn^2)
+ term.exponent += iteration - 1;
+ sum.Add(term);
+
+ if (term.exponent < -(bits - 8)) break;
+
+ //Calculate the new AGM estimates.
+ temp1.Assign(an);
+ an.Add(bn);
+ //a(n+1) = (an + bn) / 2
+ an.MulPow2(-1);
+
+ //b(n+1) = sqrt(an*bn)
+ bn.Mul(temp1);
+ bn.Sqrt();
+ }
+
+ one.Sub(sum);
+ one = one.Reciprocal();
+ return new BigFloat(one, a0.mantissa.Precision);
+ }
+
+ private static void CalculateFactorials(int numBits)
+ {
+ System.Collections.Generic.List list = new System.Collections.Generic.List(64);
+ System.Collections.Generic.List list2 = new System.Collections.Generic.List(64);
+
+ PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);
+ PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
+
+ BigFloat factorial = new BigFloat(1, extendedPrecision);
+ BigFloat reciprocal;
+
+ //Calculate e while we're at it
+ BigFloat e = new BigFloat(1, extendedPrecision);
+
+ list.Add(new BigFloat(factorial, normalPrecision));
+
+ for (int i = 1; i < Int32.MaxValue; i++)
+ {
+ BigFloat number = new BigFloat(i, extendedPrecision);
+ factorial.Mul(number);
+
+ if (factorial.exponent > numBits) break;
+
+ list2.Add(new BigFloat(factorial, normalPrecision));
+ reciprocal = factorial.Reciprocal();
+
+ e.Add(reciprocal);
+ list.Add(new BigFloat(reciprocal, normalPrecision));
+ }
+
+ //Set the cached static values.
+ inverseFactorialCache = list.ToArray();
+ factorialCache = list2.ToArray();
+ invFactorialCutoff = numBits;
+ eCache = new BigFloat(e, normalPrecision);
+ eRCPCache = new BigFloat(e.Reciprocal(), normalPrecision);
+ }
+
+ private static void CalculateEOnly(int numBits)
+ {
+ PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);
+ PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
+
+ int iExponent = (int)(Math.Sqrt(numBits));
+
+ BigFloat factorial = new BigFloat(1, extendedPrecision);
+ BigFloat constant = new BigFloat(1, extendedPrecision);
+ constant.exponent -= iExponent;
+ BigFloat numerator = new BigFloat(constant);
+ BigFloat reciprocal;
+
+ //Calculate the 2^iExponent th root of e
+ BigFloat e = new BigFloat(1, extendedPrecision);
+
+ int i;
+ for (i = 1; i < Int32.MaxValue; i++)
+ {
+ BigFloat number = new BigFloat(i, extendedPrecision);
+ factorial.Mul(number);
+ reciprocal = factorial.Reciprocal();
+ reciprocal.Mul(numerator);
+
+ if (-reciprocal.exponent > numBits) break;
+
+ e.Add(reciprocal);
+ numerator.Mul(constant);
+ System.GC.Collect();
+ }
+
+ for (i = 0; i < iExponent; i++)
+ {
+ numerator.Assign(e);
+ e.Mul(numerator);
+ }
+
+ //Set the cached static values.
+ eCache = new BigFloat(e, normalPrecision);
+ eRCPCache = new BigFloat(e.Reciprocal(), normalPrecision);
+ }
+
+ ///
+ /// Uses the Gauss-Legendre formula for pi
+ /// Taken from http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm
+ ///
+ ///
+ private static void CalculatePi(int numBits)
+ {
+ int bits = numBits + 32;
+ //Precision extend taken out.
+ PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
+ PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);
+
+ if (scratch.Precision.NumBits != bits)
+ {
+ scratch = new BigInt(extendedPres);
+ }
+
+ //a0 = 1
+ BigFloat an = new BigFloat(1, extendedPres);
+
+ //b0 = 1/sqrt(2)
+ BigFloat bn = new BigFloat(2, extendedPres);
+ bn.Sqrt();
+ bn.exponent--;
+
+ //to = 1/4
+ BigFloat tn = new BigFloat(1, extendedPres);
+ tn.exponent -= 2;
+
+ int pn = 0;
+
+ BigFloat anTemp = new BigFloat(extendedPres);
+
+ int iteration = 0;
+ int cutoffBits = numBits >> 5;
+
+ for (iteration = 0; ; iteration++)
+ {
+ //Save a(n)
+ anTemp.Assign(an);
+
+ //Calculate new an
+ an.Add(bn);
+ an.exponent--;
+
+ //Calculate new bn
+ bn.Mul(anTemp);
+ bn.Sqrt();
+
+ //Calculate new tn
+ anTemp.Sub(an);
+ anTemp.mantissa.SquareHiFast(scratch);
+ anTemp.exponent += anTemp.exponent + pn + 1 - anTemp.mantissa.Normalise();
+ tn.Sub(anTemp);
+
+ anTemp.Assign(an);
+ anTemp.Sub(bn);
+
+ if (anTemp.exponent < -(bits - cutoffBits)) break;
+
+ //New pn
+ pn++;
+ }
+
+ an.Add(bn);
+ an.mantissa.SquareHiFast(scratch);
+ an.exponent += an.exponent + 1 - an.mantissa.Normalise();
+ tn.exponent += 2;
+ an.Div(tn);
+
+ pi = new BigFloat(an, normalPres);
+ piBy2 = new BigFloat(pi);
+ piBy2.exponent--;
+ twoPi = new BigFloat(pi, normalPres);
+ twoPi.exponent++;
+ piRecip = new BigFloat(an.Reciprocal(), normalPres);
+ twoPiRecip = new BigFloat(piRecip);
+ twoPiRecip.exponent--;
+ //1/3 is going to be useful for sin.
+ threeRecip = new BigFloat((new BigFloat(3, extendedPres)).Reciprocal(), normalPres);
+ }
+
+ ///
+ /// Calculates the odd reciprocals of the natural numbers (for atan series)
+ ///
+ ///
+ ///
+ private static void CalculateReciprocals(int numBits, int terms)
+ {
+ int bits = numBits + 32;
+ PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);
+ PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
+
+ System.Collections.Generic.List list = new System.Collections.Generic.List(terms);
+
+ for (int i = 0; i < terms; i++)
+ {
+ BigFloat term = new BigFloat(i*2 + 1, extendedPres);
+ list.Add(new BigFloat(term.Reciprocal(), normalPres));
+ }
+
+ reciprocals = list.ToArray();
+ }
+
+ ///
+ /// Does decimal rounding, for numbers without E notation.
+ ///
+ ///
+ ///
+ ///
+ private static string RoundString(string input, int places)
+ {
+ if (places <= 0) return input;
+ string trim = input.Trim();
+ char[] digits = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'};
+
+ /*
+ for (int i = 1; i <= places; i++)
+ {
+ //Skip decimal points.
+ if (trim[trim.Length - i] == '.')
+ {
+ places++;
+ continue;
+ }
+
+ int index = Array.IndexOf(digits, trim[trim.Length - i]);
+
+ if (index < 0) return input;
+
+ value += ten * index;
+ ten *= 10;
+ }
+ * */
+
+ //Look for a decimal point
+ string decimalPoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;
+
+ int indexPoint = trim.LastIndexOf(decimalPoint);
+ if (indexPoint < 0)
+ {
+ //We can't modify a string which doesn't have a decimal point.
+ return trim;
+ }
+
+ int trimPoint = trim.Length - places;
+ if (trimPoint < indexPoint) trimPoint = indexPoint;
+
+ bool roundDown = false;
+
+ if (trim[trimPoint] == '.')
+ {
+ if (trimPoint + 1 >= trim.Length)
+ {
+ roundDown = true;
+ }
+ else
+ {
+ int digit = Array.IndexOf(digits, trim[trimPoint + 1]);
+ if (digit < 5) roundDown = true;
+ }
+ }
+ else
+ {
+ int digit = Array.IndexOf(digits, trim[trimPoint]);
+ if (digit < 5) roundDown = true;
+ }
+
+ string output;
+
+ //Round down - just return a new string without the extra digits.
+ if (roundDown)
+ {
+ if (RoundingMode == RoundingModeType.EXACT)
+ {
+ return trim.Substring(0, trimPoint);
+ }
+ else
+ {
+ char[] trimChars = { '0' };
+ output = trim.Substring(0, trimPoint).TrimEnd(trimChars);
+ trimChars[0] = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator[0];
+ return output.TrimEnd(trimChars);
+ }
+ }
+
+ //Round up - bit more complicated.
+ char [] arrayOutput = trim.ToCharArray();//0, trimPoint);
+
+ //Now, we round going from the back to the front.
+ int j;
+ for (j = trimPoint - 1; j >= 0; j--)
+ {
+ int index = Array.IndexOf(digits, arrayOutput[j]);
+
+ //Skip decimal points etc...
+ if (index < 0) continue;
+
+ if (index < 9)
+ {
+ arrayOutput[j] = digits[index + 1];
+ break;
+ }
+ else
+ {
+ arrayOutput[j] = digits[0];
+ }
+ }
+
+ output = new string(arrayOutput);
+
+ if (j < 0)
+ {
+ //Need to add a new digit.
+ output = String.Format("{0}{1}", "1", output);
+ }
+
+ if (RoundingMode == RoundingModeType.EXACT)
+ {
+ return output;
+ }
+ else
+ {
+ char[] trimChars = { '0' };
+ output = output.TrimEnd(trimChars);
+ trimChars[0] = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator[0];
+ return output.TrimEnd(trimChars);
+ }
+ }
+
+ //***************************** Data *****************************
+
+
+ //Side node - this way of doing things is far from optimal, both in terms of memory use and performance.
+ private ExponentAdaptor exponent;
+ private BigInt mantissa;
+
+ ///
+ /// Storage area for calculations.
+ ///
+ private static BigInt scratch;
+
+ private static BigFloat ln2cache; //Value of ln(2)
+ private static BigFloat log2ecache; //Value of log2(e) = 1/ln(2)
+ private static BigFloat pow10cache; //Cached power of 10 for AGM log calculation
+ private static BigFloat log10recip; //1/ln(10)
+ private static BigFloat firstAGMcache; //Cached half of AGM operation.
+ private static BigFloat[] factorialCache; //The values of n!
+ private static BigFloat[] inverseFactorialCache; //Values of 1/n! up to 2^-m where m = invFactorialCutoff (below)
+ private static int invFactorialCutoff; //The number of significant bits for the cutoff of the inverse factorials.
+ private static BigFloat eCache; //Value of e cached to invFactorialCutoff bits
+ private static BigFloat eRCPCache; //Reciprocal of e
+ private static BigFloat logNewtonConstant; //1.5*ln(2) - 1
+ private static BigFloat pi; //pi
+ private static BigFloat piBy2; //pi/2
+ private static BigFloat twoPi; //2*pi
+ private static BigFloat piRecip; //1/pi
+ private static BigFloat twoPiRecip; //1/2*pi
+ private static BigFloat threeRecip; //1/3
+ private static BigFloat[] reciprocals; //1/x
+
+ ///
+ /// The number of decimal digits to round the output of ToString() by
+ ///
+ public static int RoundingDigits { get; set; }
+
+ ///
+ /// The way in which ToString() should deal with insignificant trailing zeroes
+ ///
+ public static RoundingModeType RoundingMode { get; set; }
+ }
+}
\ No newline at end of file