psychlops silverlight
Rev. | 08bcb090f6a3e6ee38c712571f24b57a4dd32fb4 |
---|---|
Size | 127,679 bytes |
Time | 2012-03-29 01:46:02 |
Author | HOSOKAWA Kenchi |
Log Message | 5
|
// http://www.fractal-landscapes.co.uk/bigint.html
using System;
namespace BigNum
{
/// <summary>
/// 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.
/// </summary>
public class BigFloat
{
/// <summary>
/// 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
/// </summary>
public enum SpecialValueType
{
/// <summary>
/// Not a special value
/// </summary>
NONE = 0,
/// <summary>
/// Zero
/// </summary>
ZERO,
/// <summary>
/// Positive infinity
/// </summary>
INF_PLUS,
/// <summary>
/// Negative infinity
/// </summary>
INF_MINUS,
/// <summary>
/// Not a number
/// </summary>
NAN
}
/// <summary>
/// This affects the ToString() method.
///
/// With Trim rounding, all insignificant zero digits are drip
/// </summary>
public enum RoundingModeType
{
/// <summary>
/// Trim non-significant zeros from ToString output after rounding
/// </summary>
TRIM,
/// <summary>
/// Keep all non-significant zeroes in ToString output after rounding
/// </summary>
EXACT
}
/// <summary>
/// A wrapper for the signed exponent, avoiding overflow.
/// </summary>
protected struct ExponentAdaptor
{
/// <summary>
/// The 32-bit exponent
/// </summary>
public Int32 exponent
{
get { return expValue; }
set { expValue = value; }
}
/// <summary>
/// Implicit cast to Int32
/// </summary>
public static implicit operator Int32(ExponentAdaptor adaptor)
{
return adaptor.expValue;
}
/// <summary>
/// Implicit cast from Int32 to ExponentAdaptor
/// </summary>
/// <param name="value"></param>
/// <returns></returns>
public static implicit operator ExponentAdaptor(Int32 value)
{
ExponentAdaptor adaptor = new ExponentAdaptor();
adaptor.expValue = value;
return adaptor;
}
/// <summary>
/// Overloaded increment operator
/// </summary>
public static ExponentAdaptor operator ++(ExponentAdaptor adaptor)
{
adaptor = adaptor + 1;
return adaptor;
}
/// <summary>
/// Overloaded decrement operator
/// </summary>
public static ExponentAdaptor operator --(ExponentAdaptor adaptor)
{
adaptor = adaptor - 1;
return adaptor;
}
/// <summary>
/// Overloaded addition operator
/// </summary>
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;
}
/// <summary>
/// Overloaded subtraction operator
/// </summary>
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;
}
/// <summary>
/// Overloaded multiplication operator
/// </summary>
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;
}
/// <summary>
/// Overloaded division operator
/// </summary>
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;
}
/// <summary>
/// Overloaded right-shift operator
/// </summary>
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;
}
/// <summary>
/// Overloaded left-shift operator
/// </summary>
/// <param name="a1"></param>
/// <param name="shift"></param>
/// <returns></returns>
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 **************************
/// <summary>
/// Constructs a 128-bit BigFloat
///
/// Sets the value to zero
/// </summary>
static BigFloat()
{
RoundingDigits = 3;
RoundingMode = RoundingModeType.TRIM;
scratch = new BigInt(new PrecisionSpec(128, PrecisionSpec.BaseType.BIN));
}
/// <summary>
/// Constructs a BigFloat of the required precision
///
/// Sets the value to zero
/// </summary>
/// <param name="mantissaPrec"></param>
public BigFloat(PrecisionSpec mantissaPrec)
{
Init(mantissaPrec);
}
/// <summary>
/// Constructs a big float from a UInt32 to the required precision
/// </summary>
/// <param name="value"></param>
/// <param name="mantissaPrec"></param>
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;
}
/// <summary>
/// Constructs a BigFloat from an Int32 to the required precision
/// </summary>
/// <param name="value"></param>
/// <param name="mantissaPrec"></param>
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;
}
/// <summary>
/// Constructs a BigFloat from a 64-bit integer
/// </summary>
/// <param name="value"></param>
/// <param name="mantissaPrec"></param>
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;
}
/// <summary>
/// Constructs a BigFloat from a 64-bit unsigned integer
/// </summary>
/// <param name="value"></param>
/// <param name="mantissaPrec"></param>
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;
}
/// <summary>
/// Constructs a BigFloat from a BigInt, using the specified precision
/// </summary>
/// <param name="value"></param>
/// <param name="mantissaPrec"></param>
public BigFloat(BigInt value, PrecisionSpec mantissaPrec)
{
if (value.IsZero())
{
Init(mantissaPrec);
SetZero();
return;
}
mantissa = new BigInt(value, mantissaPrec);
exponent = BigInt.GetMSB(value);
mantissa.Normalise();
}
/// <summary>
/// Construct a BigFloat from a double-precision floating point number
/// </summary>
/// <param name="value"></param>
/// <param name="mantissaPrec"></param>
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;
}
/// <summary>
/// Copy constructor
/// </summary>
/// <param name="value"></param>
public BigFloat(BigFloat value)
{
Init(value.mantissa.Precision);
exponent = value.exponent;
mantissa.Assign(value.mantissa);
}
/// <summary>
/// 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.
/// </summary>
/// <param name="value"></param>
/// <param name="mantissaPrec"></param>
public BigFloat(BigFloat value, PrecisionSpec mantissaPrec)
{
Init(mantissaPrec);
exponent = value.exponent;
if (mantissa.AssignHigh(value.mantissa)) exponent++;
}
/// <summary>
/// Constructs a BigFloat from a string
/// </summary>
/// <param name="value"></param>
/// <param name="mantissaPrec"></param>
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 *************************
/// <summary>
/// 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.
/// </summary>
public PrecisionSpec Precision
{
get { return mantissa.Precision; }
}
/// <summary>
/// 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)
/// </summary>
public bool Sign
{
get { return mantissa.Sign; }
set { mantissa.Sign = value; }
}
/// <summary>
/// Read-only property.
/// True if the number is NAN, INF_PLUS, INF_MINUS or ZERO
/// False if the number has any other value.
/// </summary>
public bool IsSpecialValue
{
get
{
return (exponent == Int32.MaxValue || mantissa.IsZero());
}
}
/// <summary>
/// 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.
/// </summary>
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 *******************
/// <summary>
/// Gets pi to the indicated precision
/// </summary>
/// <param name="precision">The precision to perform the calculation to</param>
/// <returns>pi (the ratio of the area of a circle to its diameter)</returns>
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;
}
/// <summary>
/// Get e to the indicated precision
/// </summary>
/// <param name="precision">The preicision to perform the calculation to</param>
/// <returns>e (the number for which the d/dx(e^x) = e^x)</returns>
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 ********************
/// <summary>
/// Addition (this = this + n2)
/// </summary>
/// <param name="n2">The number to add</param>
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();
}
}
/// <summary>
/// Subtraction (this = this - n2)
/// </summary>
/// <param name="n2">The number to subtract from this</param>
public void Sub(BigFloat n2)
{
n2.mantissa.Sign = !n2.mantissa.Sign;
Add(n2);
n2.mantissa.Sign = !n2.mantissa.Sign;
}
/// <summary>
/// Multiplication (this = this * n2)
/// </summary>
/// <param name="n2">The number to multiply this by</param>
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;
}
/// <summary>
/// Division (this = this / n2)
/// </summary>
/// <param name="n2">The number to divide this by</param>
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);
}
}
/// <summary>
/// Multiply by a power of 2 (-ve implies division)
/// </summary>
/// <param name="pow2"></param>
public void MulPow2(int pow2)
{
exponent += pow2;
}
/// <summary>
/// Division-based reciprocal, fastest for small precisions up to 15,000 bits.
/// </summary>
/// <returns>The reciprocal 1/this</returns>
public BigFloat Reciprocal()
{
if (mantissa.Precision.NumBits >= 8192) return ReciprocalNewton();
BigFloat reciprocal = new BigFloat(1u, mantissa.Precision);
reciprocal.Div(this);
return reciprocal;
}
/// <summary>
/// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.
/// </summary>
/// <returns>The reciprocal 1/this</returns>
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;
}
/// <summary>
/// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.
/// </summary>
/// <returns>The reciprocal 1/this</returns>
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;
}
/// <summary>
/// Sets this equal to the input
/// </summary>
/// <param name="n2"></param>
public void Assign(BigFloat n2)
{
exponent = n2.exponent;
if (mantissa.AssignHigh(n2.mantissa)) exponent++;
}
//********************* Comparison Functions *******************
/// <summary>
/// Greater than comparison
/// </summary>
/// <param name="n2">the number to compare this to</param>
/// <returns>true iff this is greater than n2 (this > n2)</returns>
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);
}
/// <summary>
/// Less than comparison
/// </summary>
/// <param name="n2">the number to compare this to</param>
/// <returns>true iff this is less than n2 (this < n2)</returns>
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);
}
/// <summary>
/// Greater than comparison
/// </summary>
/// <param name="i">the number to compare this to</param>
/// <returns>true iff this is greater than n2 (this > n2)</returns>
public bool GreaterThan(int i)
{
BigFloat integer = new BigFloat(i, mantissa.Precision);
return GreaterThan(integer);
}
/// <summary>
/// Less than comparison
/// </summary>
/// <param name="i">the number to compare this to</param>
/// <returns>true iff this is less than n2 (this < n2)</returns>
public bool LessThan(int i)
{
BigFloat integer = new BigFloat(i, mantissa.Precision);
return LessThan(integer);
}
/// <summary>
/// Compare to zero
/// </summary>
/// <returns>true if this is zero (this == 0)</returns>
public bool IsZero()
{
return (mantissa.IsZero());
}
//******************** Mathematical Functions ******************
/// <summary>
/// Sets the number to the biggest integer numerically closer to zero, if possible.
/// </summary>
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));
}
/// <summary>
/// Sets the number to its fractional component (equivalent to 'this' - (int)'this')
/// </summary>
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();
}
/// <summary>
/// Calculates tan(x)
/// </summary>
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;
}
/// <summary>
/// Calculates Cos(x)
/// </summary>
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();
}
/// <summary>
/// Calculates Sin(x):
/// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].
/// </summary>
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;
}
/// <summary>
/// Hyperbolic Sin (sinh) function
/// </summary>
public void Sinh()
{
if (IsSpecialValue)
{
return;
}
Exp();
Sub(Reciprocal());
exponent--;
}
/// <summary>
/// Hyperbolic cosine (cosh) function
/// </summary>
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--;
}
/// <summary>
/// Hyperbolic tangent function (tanh)
/// </summary>
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);
}
/// <summary>
/// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)
/// </summary>
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;
}
}
/// <summary>
/// arccos(): the inverse function of cos(), range (0..pi)
/// </summary>
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;
}
}
}
/// <summary>
/// arctan(): the inverse function of sin(), range of (-pi/2..pi/2)
/// </summary>
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;
}
}
/// <summary>
/// Arcsinh(): the inverse sinh function
/// </summary>
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();
}
/// <summary>
/// Arccosh(): the inverse cosh() function
/// </summary>
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();
}
/// <summary>
/// Arctanh(): the inverse tanh function
/// </summary>
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();
}
}
}
/// <summary>
/// Two-variable iterative square root, taken from
/// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method
/// </summary>
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);
}
/// <summary>
/// The natural logarithm, ln(x)
/// </summary>
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);
}
/// <summary>
/// Log to the base 10
/// </summary>
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);
}
/// <summary>
/// The exponential function. Less accurate for high exponents, scales poorly with the number
/// of bits.
/// </summary>
public void Exp()
{
Exp(mantissa.Precision.NumBits);
}
/// <summary>
/// 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)
///
/// </summary>
/// <param name="power"></param>
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);
}
/// <summary>
/// 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.
/// </summary>
/// <param name="power"></param>
public void Pow(BigFloat power)
{
Log();
Mul(power);
Exp();
}
//******************** Static Math Functions *******************
/// <summary>
/// Returns the integer component of the input
/// </summary>
/// <param name="n1">The input number</param>
/// <remarks>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.</remarks>
public static BigFloat Floor(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Floor();
return n1;
}
/// <summary>
/// Returns the fractional (non-integer component of the input)
/// </summary>
/// <param name="n1">The input number</param>
public static BigFloat FPart(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.FPart();
return n1;
}
/// <summary>
/// Calculates tan(x)
/// </summary>
/// <param name="n1">The angle (in radians) to find the tangent of</param>
public static BigFloat Tan(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Tan();
return n1;
}
/// <summary>
/// Calculates Cos(x)
/// </summary>
/// <param name="n1">The angle (in radians) to find the cosine of</param>
/// <remarks>This is a reasonably fast function for smaller precisions, but
/// doesn't scale well for higher precision arguments</remarks>
public static BigFloat Cos(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Cos();
return n1;
}
/// <summary>
/// Calculates Sin(x):
/// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].
/// </summary>
/// <param name="n1">The angle to find the sine of (in radians)</param>
/// <remarks>This is a resonably fast function, for smaller precision arguments, but doesn't
/// scale very well with the number of bits in the input.</remarks>
public static BigFloat Sin(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Sin();
return n1;
}
/// <summary>
/// Hyperbolic Sin (sinh) function
/// </summary>
/// <param name="n1">The number to find the hyperbolic sine of</param>
public static BigFloat Sinh(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Sinh();
return n1;
}
/// <summary>
/// Hyperbolic cosine (cosh) function
/// </summary>
/// <param name="n1">The number to find the hyperbolic cosine of</param>
public static BigFloat Cosh(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Cosh();
return n1;
}
/// <summary>
/// Hyperbolic tangent function (tanh)
/// </summary>
/// <param name="n1">The number to find the hyperbolic tangent of</param>
public static BigFloat Tanh(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Tanh();
return n1;
}
/// <summary>
/// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)
/// </summary>
/// <param name="n1">The number to find the arcsine of (-pi/2..pi/2)</param>
/// <remarks>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.
/// </remarks>
public static BigFloat Arcsin(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Arcsin();
return n1;
}
/// <summary>
/// arccos(): the inverse function of cos(), input range (0..pi)
/// </summary>
/// <param name="n1">The number to find the arccosine of (0..pi)</param>
/// <remarks>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.
/// </remarks>
public static BigFloat Arccos(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Arccos();
return n1;
}
/// <summary>
/// arctan(): the inverse function of sin(), input range of (-pi/2..pi/2)
/// </summary>
/// <param name="n1">The number to find the arctangent of (-pi/2..pi/2)</param>
/// <remarks>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.
/// </remarks>
public static BigFloat Arctan(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Arctan();
return n1;
}
/// <summary>
/// Arcsinh(): the inverse sinh function
/// </summary>
/// <param name="n1">The number to find the inverse hyperbolic sine of</param>
public static BigFloat Arcsinh(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Arcsinh();
return n1;
}
/// <summary>
/// Arccosh(): the inverse cosh() function
/// </summary>
/// <param name="n1">The number to find the inverse hyperbolic cosine of</param>
public static BigFloat Arccosh(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Arccosh();
return n1;
}
/// <summary>
/// Arctanh(): the inverse tanh function
/// </summary>
/// <param name="n1">The number to fine the inverse hyperbolic tan of</param>
public static BigFloat Arctanh(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Arctanh();
return n1;
}
/// <summary>
/// Two-variable iterative square root, taken from
/// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method
/// </summary>
/// <remarks>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.
/// </remarks>
public static BigFloat Sqrt(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Sqrt();
return n1;
}
/// <summary>
/// The natural logarithm, ln(x) (log base e)
/// </summary>
/// <remarks>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.
/// </remarks>
/// <param name="n1">The number to find the natural logarithm of</param>
public static BigFloat Log(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Log();
return n1;
}
/// <summary>
/// Base 10 logarithm of a number
/// </summary>
/// <remarks>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.
/// </remarks>
/// <param name="n1">The number to find the base 10 logarithm of</param>
public static BigFloat Log10(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Log10();
return n1;
}
/// <summary>
/// The exponential function. Less accurate for high exponents, scales poorly with the number
/// of bits. This is quite fast for low-precision arguments.
/// </summary>
public static BigFloat Exp(BigFloat n1)
{
BigFloat res = new BigFloat(n1);
n1.Exp();
return n1;
}
/// <summary>
/// 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).
/// </summary>
/// <param name="n1">The number to raise to the power</param>
/// <param name="power">The power to raise it to</param>
public static BigFloat Pow(BigFloat n1, int power)
{
BigFloat res = new BigFloat(n1);
n1.Pow(power);
return n1;
}
/// <summary>
/// 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)
/// </summary>
/// <remarks>
/// 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.
/// </remarks>
/// <param name="n1">The number to raise to a power</param>
/// <param name="power">The power to raise it to</param>
public static BigFloat Pow(BigFloat n1, BigFloat power)
{
BigFloat res = new BigFloat(n1);
n1.Pow(power);
return n1;
}
//********************** Static functions **********************
/// <summary>
/// Adds two numbers and returns the result
/// </summary>
public static BigFloat Add(BigFloat n1, BigFloat n2)
{
BigFloat ret = new BigFloat(n1);
ret.Add(n2);
return ret;
}
/// <summary>
/// Subtracts two numbers and returns the result
/// </summary>
public static BigFloat Sub(BigFloat n1, BigFloat n2)
{
BigFloat ret = new BigFloat(n1);
ret.Sub(n2);
return ret;
}
/// <summary>
/// Multiplies two numbers and returns the result
/// </summary>
public static BigFloat Mul(BigFloat n1, BigFloat n2)
{
BigFloat ret = new BigFloat(n1);
ret.Mul(n2);
return ret;
}
/// <summary>
/// Divides two numbers and returns the result
/// </summary>
public static BigFloat Div(BigFloat n1, BigFloat n2)
{
BigFloat ret = new BigFloat(n1);
ret.Div(n2);
return ret;
}
/// <summary>
/// Tests whether n1 is greater than n2
/// </summary>
public static bool GreaterThan(BigFloat n1, BigFloat n2)
{
return n1.GreaterThan(n2);
}
/// <summary>
/// Tests whether n1 is less than n2
/// </summary>
public static bool LessThan(BigFloat n1, BigFloat n2)
{
return n1.LessThan(n2);
}
//******************* Fast static functions ********************
/// <summary>
/// Adds two numbers and assigns the result to res.
/// </summary>
/// <param name="res">a pre-existing BigFloat to take the result</param>
/// <param name="n1">the first number</param>
/// <param name="n2">the second number</param>
/// <returns>a handle to res</returns>
public static BigFloat Add(BigFloat res, BigFloat n1, BigFloat n2)
{
res.Assign(n1);
res.Add(n2);
return res;
}
/// <summary>
/// Subtracts two numbers and assigns the result to res.
/// </summary>
/// <param name="res">a pre-existing BigFloat to take the result</param>
/// <param name="n1">the first number</param>
/// <param name="n2">the second number</param>
/// <returns>a handle to res</returns>
public static BigFloat Sub(BigFloat res, BigFloat n1, BigFloat n2)
{
res.Assign(n1);
res.Sub(n2);
return res;
}
/// <summary>
/// Multiplies two numbers and assigns the result to res.
/// </summary>
/// <param name="res">a pre-existing BigFloat to take the result</param>
/// <param name="n1">the first number</param>
/// <param name="n2">the second number</param>
/// <returns>a handle to res</returns>
public static BigFloat Mul(BigFloat res, BigFloat n1, BigFloat n2)
{
res.Assign(n1);
res.Mul(n2);
return res;
}
/// <summary>
/// Divides two numbers and assigns the result to res.
/// </summary>
/// <param name="res">a pre-existing BigFloat to take the result</param>
/// <param name="n1">the first number</param>
/// <param name="n2">the second number</param>
/// <returns>a handle to res</returns>
public static BigFloat Div(BigFloat res, BigFloat n1, BigFloat n2)
{
res.Assign(n1);
res.Div(n2);
return res;
}
//************************* Operators **************************
/// <summary>
/// The addition operator
/// </summary>
public static BigFloat operator +(BigFloat n1, BigFloat n2)
{
return Add(n1, n2);
}
/// <summary>
/// The subtraction operator
/// </summary>
public static BigFloat operator -(BigFloat n1, BigFloat n2)
{
return Sub(n1, n2);
}
/// <summary>
/// The multiplication operator
/// </summary>
public static BigFloat operator *(BigFloat n1, BigFloat n2)
{
return Mul(n1, n2);
}
/// <summary>
/// The division operator
/// </summary>
public static BigFloat operator /(BigFloat n1, BigFloat n2)
{
return Div(n1, n2);
}
//************************** Conversions *************************
/// <summary>
/// Converts a BigFloat to an BigInt with the specified precision
/// </summary>
/// <param name="n1">The number to convert</param>
/// <param name="precision">The precision to convert it with</param>
/// <param name="round">Do we round the number if we are truncating the mantissa?</param>
/// <returns></returns>
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;
}
/// <summary>
/// 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
/// </summary>
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 *****************
/// <summary>
/// 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
/// </summary>
/// <param name="targetExponent"></param>
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;
}
/// <summary>
/// The binary logarithm, log2(x) - for precisions above 1000 bits, use Log() and convert the base.
/// </summary>
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);
}
/// <summary>
/// Tried the newton method for logs, but the exponential function is too slow to do it.
/// </summary>
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);
}
/// <summary>
/// Log(x) implemented as an Arithmetic-Geometric Mean. Fast for high precisions.
/// </summary>
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);
}
/// <summary>
/// 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
/// </summary>
/// <param name="numBits"></param>
/// <returns></returns>
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<BigFloat> list = new System.Collections.Generic.List<BigFloat>(64);
System.Collections.Generic.List<BigFloat> list2 = new System.Collections.Generic.List<BigFloat>(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);
}
/// <summary>
/// Uses the Gauss-Legendre formula for pi
/// Taken from http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm
/// </summary>
/// <param name="numBits"></param>
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);
}
/// <summary>
/// Calculates the odd reciprocals of the natural numbers (for atan series)
/// </summary>
/// <param name="numBits"></param>
/// <param name="terms"></param>
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<BigFloat> list = new System.Collections.Generic.List<BigFloat>(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();
}
/// <summary>
/// Does decimal rounding, for numbers without E notation.
/// </summary>
/// <param name="input"></param>
/// <param name="places"></param>
/// <returns></returns>
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;
/// <summary>
/// Storage area for calculations.
/// </summary>
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
/// <summary>
/// The number of decimal digits to round the output of ToString() by
/// </summary>
public static int RoundingDigits { get; set; }
/// <summary>
/// The way in which ToString() should deal with insignificant trailing zeroes
/// </summary>
public static RoundingModeType RoundingMode { get; set; }
}
}