1 // http://www.fractal-landscapes.co.uk/bigint.html
\r
8 /// An arbitrary-precision floating-point class
\r
11 /// Each number is stored as an exponent (32-bit signed integer), and a mantissa
\r
12 /// (n-bit) BigInteger. The sign of the number is stored in the BigInteger
\r
14 /// Applicability and Performance:
\r
15 /// This class is designed to be used for small extended precisions. It may not be
\r
16 /// safe (and certainly won't be fast) to use it with mixed-precision arguments.
\r
17 /// It does support, but will not be efficient for, numbers over around 2048 bits.
\r
20 /// All conversions to and from strings are slow.
\r
22 /// Conversions from simple integer types Int32, Int64, UInt32, UInt64 are performed
\r
23 /// using the appropriate constructor, and are relatively fast.
\r
25 /// The class is written entirely in managed C# code, with not native or managed
\r
26 /// assembler. The use of native assembler would speed up the multiplication operations
\r
27 /// many times over, and therefore all higher-order operations too.
\r
29 public class BigFloat
\r
32 /// Floats can have 4 special value types:
\r
34 /// NaN: Not a number (cannot be changed using any operations)
\r
35 /// Infinity: Positive infinity. Some operations e.g. Arctan() allow this input.
\r
36 /// -Infinity: Negative infinity. Some operations allow this input.
\r
39 public enum SpecialValueType
\r
42 /// Not a special value
\r
50 /// Positive infinity
\r
54 /// Negative infinity
\r
64 /// This affects the ToString() method.
\r
66 /// With Trim rounding, all insignificant zero digits are drip
\r
68 public enum RoundingModeType
\r
71 /// Trim non-significant zeros from ToString output after rounding
\r
75 /// Keep all non-significant zeroes in ToString output after rounding
\r
81 /// A wrapper for the signed exponent, avoiding overflow.
\r
83 protected struct ExponentAdaptor
\r
86 /// The 32-bit exponent
\r
88 public Int32 exponent
\r
90 get { return expValue; }
\r
91 set { expValue = value; }
\r
95 /// Implicit cast to Int32
\r
97 public static implicit operator Int32(ExponentAdaptor adaptor)
\r
99 return adaptor.expValue;
\r
103 /// Implicit cast from Int32 to ExponentAdaptor
\r
105 /// <param name="value"></param>
\r
106 /// <returns></returns>
\r
107 public static implicit operator ExponentAdaptor(Int32 value)
\r
109 ExponentAdaptor adaptor = new ExponentAdaptor();
\r
110 adaptor.expValue = value;
\r
115 /// Overloaded increment operator
\r
117 public static ExponentAdaptor operator ++(ExponentAdaptor adaptor)
\r
119 adaptor = adaptor + 1;
\r
124 /// Overloaded decrement operator
\r
126 public static ExponentAdaptor operator --(ExponentAdaptor adaptor)
\r
128 adaptor = adaptor - 1;
\r
133 /// Overloaded addition operator
\r
135 public static ExponentAdaptor operator +(ExponentAdaptor a1, ExponentAdaptor a2)
\r
137 if (a1.expValue == Int32.MaxValue) return a1;
\r
139 Int64 temp = (Int64)a1.expValue;
\r
140 temp += (Int64)(a2.expValue);
\r
142 if (temp > (Int64)Int32.MaxValue)
\r
144 a1.expValue = Int32.MaxValue;
\r
146 else if (temp < (Int64)Int32.MinValue)
\r
148 a1.expValue = Int32.MinValue;
\r
152 a1.expValue = (Int32)temp;
\r
159 /// Overloaded subtraction operator
\r
161 public static ExponentAdaptor operator -(ExponentAdaptor a1, ExponentAdaptor a2)
\r
163 if (a1.expValue == Int32.MaxValue) return a1;
\r
165 Int64 temp = (Int64)a1.expValue;
\r
166 temp -= (Int64)(a2.expValue);
\r
168 if (temp > (Int64)Int32.MaxValue)
\r
170 a1.expValue = Int32.MaxValue;
\r
172 else if (temp < (Int64)Int32.MinValue)
\r
174 a1.expValue = Int32.MinValue;
\r
178 a1.expValue = (Int32)temp;
\r
185 /// Overloaded multiplication operator
\r
187 public static ExponentAdaptor operator *(ExponentAdaptor a1, ExponentAdaptor a2)
\r
189 if (a1.expValue == Int32.MaxValue) return a1;
\r
191 Int64 temp = (Int64)a1.expValue;
\r
192 temp *= (Int64)a2.expValue;
\r
194 if (temp > (Int64)Int32.MaxValue)
\r
196 a1.expValue = Int32.MaxValue;
\r
198 else if (temp < (Int64)Int32.MinValue)
\r
200 a1.expValue = Int32.MinValue;
\r
204 a1.expValue = (Int32)temp;
\r
211 /// Overloaded division operator
\r
213 public static ExponentAdaptor operator /(ExponentAdaptor a1, ExponentAdaptor a2)
\r
215 if (a1.expValue == Int32.MaxValue) return a1;
\r
217 ExponentAdaptor res = new ExponentAdaptor();
\r
218 res.expValue = a1.expValue / a2.expValue;
\r
223 /// Overloaded right-shift operator
\r
225 public static ExponentAdaptor operator >>(ExponentAdaptor a1, int shift)
\r
227 if (a1.expValue == Int32.MaxValue) return a1;
\r
229 ExponentAdaptor res = new ExponentAdaptor();
\r
230 res.expValue = a1.expValue >> shift;
\r
235 /// Overloaded left-shift operator
\r
237 /// <param name="a1"></param>
\r
238 /// <param name="shift"></param>
\r
239 /// <returns></returns>
\r
240 public static ExponentAdaptor operator <<(ExponentAdaptor a1, int shift)
\r
242 if (a1.expValue == 0) return a1;
\r
244 ExponentAdaptor res = new ExponentAdaptor();
\r
245 res.expValue = a1.expValue;
\r
249 res.expValue = Int32.MaxValue;
\r
253 Int64 temp = a1.expValue;
\r
254 temp = temp << shift;
\r
256 if (temp > (Int64)Int32.MaxValue)
\r
258 res.expValue = Int32.MaxValue;
\r
260 else if (temp < (Int64)Int32.MinValue)
\r
262 res.expValue = Int32.MinValue;
\r
266 res.expValue = (Int32)temp;
\r
273 private Int32 expValue;
\r
276 //************************ Constructors **************************
\r
279 /// Constructs a 128-bit BigFloat
\r
281 /// Sets the value to zero
\r
285 RoundingDigits = 3;
\r
286 RoundingMode = RoundingModeType.TRIM;
\r
287 scratch = new BigInt(new PrecisionSpec(128, PrecisionSpec.BaseType.BIN));
\r
291 /// Constructs a BigFloat of the required precision
\r
293 /// Sets the value to zero
\r
295 /// <param name="mantissaPrec"></param>
\r
296 public BigFloat(PrecisionSpec mantissaPrec)
\r
298 Init(mantissaPrec);
\r
302 /// Constructs a big float from a UInt32 to the required precision
\r
304 /// <param name="value"></param>
\r
305 /// <param name="mantissaPrec"></param>
\r
306 public BigFloat(UInt32 value, PrecisionSpec mantissaPrec)
\r
308 int mbWords = ((mantissaPrec.NumBits) >> 5);
\r
309 if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
\r
310 int newManBits = mbWords << 5;
\r
312 //For efficiency, we just use a 32-bit exponent
\r
315 mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
316 //scratch = new BigInt(mantissa.Precision);
\r
318 int bit = BigInt.GetMSB(value);
\r
319 if (bit == -1) return;
\r
321 int shift = mantissa.Precision.NumBits - (bit + 1);
\r
322 mantissa.LSH(shift);
\r
327 /// Constructs a BigFloat from an Int32 to the required precision
\r
329 /// <param name="value"></param>
\r
330 /// <param name="mantissaPrec"></param>
\r
331 public BigFloat(Int32 value, PrecisionSpec mantissaPrec)
\r
333 int mbWords = ((mantissaPrec.NumBits) >> 5);
\r
334 if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
\r
335 int newManBits = mbWords << 5;
\r
337 //For efficiency, we just use a 32-bit exponent
\r
343 if (value == Int32.MinValue)
\r
345 uValue = 0x80000000;
\r
349 uValue = (UInt32)(-value);
\r
354 uValue = (UInt32)value;
\r
357 mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
358 //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
360 int bit = BigInt.GetMSB(uValue);
\r
361 if (bit == -1) return;
\r
363 int shift = mantissa.Precision.NumBits - (bit + 1);
\r
364 mantissa.LSH(shift);
\r
369 /// Constructs a BigFloat from a 64-bit integer
\r
371 /// <param name="value"></param>
\r
372 /// <param name="mantissaPrec"></param>
\r
373 public BigFloat(Int64 value, PrecisionSpec mantissaPrec)
\r
375 int mbWords = ((mantissaPrec.NumBits) >> 5);
\r
376 if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
\r
377 int newManBits = mbWords << 5;
\r
379 //For efficiency, we just use a 32-bit exponent
\r
385 if (value == Int64.MinValue)
\r
387 uValue = 0x80000000;
\r
391 uValue = (UInt64)(-value);
\r
396 uValue = (UInt64)value;
\r
399 mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
400 //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
402 int bit = BigInt.GetMSB(uValue);
\r
403 if (bit == -1) return;
\r
405 int shift = mantissa.Precision.NumBits - (bit + 1);
\r
408 mantissa.LSH(shift);
\r
412 mantissa.SetHighDigit((uint)(uValue >> (-shift)));
\r
418 /// Constructs a BigFloat from a 64-bit unsigned integer
\r
420 /// <param name="value"></param>
\r
421 /// <param name="mantissaPrec"></param>
\r
422 public BigFloat(UInt64 value, PrecisionSpec mantissaPrec)
\r
424 int mbWords = ((mantissaPrec.NumBits) >> 5);
\r
425 if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
\r
426 int newManBits = mbWords << 5;
\r
428 //For efficiency, we just use a 32-bit exponent
\r
431 int bit = BigInt.GetMSB(value);
\r
433 mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
434 //scratch = new BigInt(mantissa.Precision);
\r
436 int shift = mantissa.Precision.NumBits - (bit + 1);
\r
439 mantissa.LSH(shift);
\r
443 mantissa.SetHighDigit((uint)(value >> (-shift)));
\r
449 /// Constructs a BigFloat from a BigInt, using the specified precision
\r
451 /// <param name="value"></param>
\r
452 /// <param name="mantissaPrec"></param>
\r
453 public BigFloat(BigInt value, PrecisionSpec mantissaPrec)
\r
455 if (value.IsZero())
\r
457 Init(mantissaPrec);
\r
462 mantissa = new BigInt(value, mantissaPrec);
\r
463 exponent = BigInt.GetMSB(value);
\r
464 mantissa.Normalise();
\r
468 /// Construct a BigFloat from a double-precision floating point number
\r
470 /// <param name="value"></param>
\r
471 /// <param name="mantissaPrec"></param>
\r
472 public BigFloat(double value, PrecisionSpec mantissaPrec)
\r
476 Init(mantissaPrec);
\r
480 bool sign = (value < 0) ? true : false;
\r
482 long bits = BitConverter.DoubleToInt64Bits(value);
\r
483 // Note that the shift is sign-extended, hence the test against -1 not 1
\r
484 int valueExponent = (int)((bits >> 52) & 0x7ffL);
\r
485 long valueMantissa = bits & 0xfffffffffffffL;
\r
487 //The mantissa is stored with the top bit implied.
\r
488 valueMantissa = valueMantissa | 0x10000000000000L;
\r
490 //The exponent is biased by 1023.
\r
491 exponent = valueExponent - 1023;
\r
493 //Round the number of bits to the nearest word.
\r
494 int mbWords = ((mantissaPrec.NumBits) >> 5);
\r
495 if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
\r
496 int newManBits = mbWords << 5;
\r
498 mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
499 //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
501 if (newManBits >= 64)
\r
503 //The mantissa is 53 bits now, so add 11 to put it in the right place.
\r
504 mantissa.SetHighDigits(valueMantissa << 11);
\r
508 //To get the top word of the mantissa, shift up by 11 and down by 32 = down by 21
\r
509 mantissa.SetHighDigit((uint)(valueMantissa >> 21));
\r
512 mantissa.Sign = sign;
\r
516 /// Copy constructor
\r
518 /// <param name="value"></param>
\r
519 public BigFloat(BigFloat value)
\r
521 Init(value.mantissa.Precision);
\r
522 exponent = value.exponent;
\r
523 mantissa.Assign(value.mantissa);
\r
527 /// Copy Constructor - constructs a new BigFloat with the specified precision, copying the old one.
\r
529 /// The value is rounded towards zero in the case where precision is decreased. The Round() function
\r
530 /// should be used beforehand if a correctly rounded result is required.
\r
532 /// <param name="value"></param>
\r
533 /// <param name="mantissaPrec"></param>
\r
534 public BigFloat(BigFloat value, PrecisionSpec mantissaPrec)
\r
536 Init(mantissaPrec);
\r
537 exponent = value.exponent;
\r
538 if (mantissa.AssignHigh(value.mantissa)) exponent++;
\r
542 /// Constructs a BigFloat from a string
\r
544 /// <param name="value"></param>
\r
545 /// <param name="mantissaPrec"></param>
\r
546 public BigFloat(string value, PrecisionSpec mantissaPrec)
\r
548 Init(mantissaPrec);
\r
550 PrecisionSpec extendedPres = new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN);
\r
551 BigFloat ten = new BigFloat(10, extendedPres);
\r
552 BigFloat iPart = new BigFloat(extendedPres);
\r
553 BigFloat fPart = new BigFloat(extendedPres);
\r
554 BigFloat tenRCP = ten.Reciprocal();
\r
556 if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol))
\r
561 else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol))
\r
566 else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol))
\r
572 string decimalpoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;
\r
574 char[] digitChars = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ',', '.' };
\r
576 //Read in the integer part up the the decimal point.
\r
578 value = value.Trim();
\r
582 if (value.Length > i && value[i] == '-')
\r
588 if (value.Length > i && value[i] == '+')
\r
593 for ( ; i < value.Length; i++)
\r
595 //break on decimal point
\r
596 if (value[i] == decimalpoint[0]) break;
\r
598 int digit = Array.IndexOf(digitChars, value[i]);
\r
599 if (digit < 0) break;
\r
601 //Ignore place separators (assumed either , or .)
\r
602 if (digit > 9) continue;
\r
604 if (i > 0) iPart.Mul(ten);
\r
605 iPart.Add(new BigFloat(digit, extendedPres));
\r
608 //If we've run out of characters, assign everything and return
\r
609 if (i == value.Length)
\r
611 iPart.mantissa.Sign = sign;
\r
612 exponent = iPart.exponent;
\r
613 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
\r
617 //Assign the characters after the decimal point to fPart
\r
618 if (value[i] == '.' && i < value.Length - 1)
\r
620 BigFloat RecipToUse = new BigFloat(tenRCP);
\r
622 for (i++; i < value.Length; i++)
\r
624 int digit = Array.IndexOf(digitChars, value[i]);
\r
625 if (digit < 0) break;
\r
626 BigFloat temp = new BigFloat(digit, extendedPres);
\r
627 temp.Mul(RecipToUse);
\r
628 RecipToUse.Mul(tenRCP);
\r
633 //If we're run out of characters, add fPart and iPart and return
\r
634 if (i == value.Length)
\r
637 iPart.mantissa.Sign = sign;
\r
638 exponent = iPart.exponent;
\r
639 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
\r
643 if (value[i] == '+' || value[i] == '-') i++;
\r
645 if (i == value.Length)
\r
648 iPart.mantissa.Sign = sign;
\r
649 exponent = iPart.exponent;
\r
650 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
\r
654 //Look for exponential notation.
\r
655 if ((value[i] == 'e' || value[i] == 'E') && i < value.Length - 1)
\r
657 //Convert the exponent to an int.
\r
662 exp = System.Convert.ToInt32(new string(value.ToCharArray()));// i + 1, value.Length - (i + 1))));
\r
667 iPart.mantissa.Sign = sign;
\r
668 exponent = iPart.exponent;
\r
669 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
\r
673 //Raise or lower 10 to the power of the exponent
\r
674 BigFloat acc = new BigFloat(1, extendedPres);
\r
675 BigFloat temp = new BigFloat(1, extendedPres);
\r
677 int powerTemp = exp;
\r
679 BigFloat multiplierToUse;
\r
683 multiplierToUse = new BigFloat(tenRCP);
\r
688 multiplierToUse = new BigFloat(ten);
\r
691 //Fast power function
\r
692 while (powerTemp != 0)
\r
694 temp.Mul(multiplierToUse);
\r
695 multiplierToUse.Assign(temp);
\r
697 if ((powerTemp & 1) != 0)
\r
707 iPart.mantissa.Sign = sign;
\r
708 exponent = iPart.exponent;
\r
709 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
\r
715 iPart.mantissa.Sign = sign;
\r
716 exponent = iPart.exponent;
\r
717 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
\r
721 private void Init(PrecisionSpec mantissaPrec)
\r
723 int mbWords = ((mantissaPrec.NumBits) >> 5);
\r
724 if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
\r
725 int newManBits = mbWords << 5;
\r
727 //For efficiency, we just use a 32-bit exponent
\r
729 mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
730 //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
\r
733 //************************** Properties *************************
\r
736 /// Read-only property. Returns the precision specification of the mantissa.
\r
738 /// Floating point numbers are represented as 2^exponent * mantissa, where the
\r
739 /// mantissa and exponent are integers. Note that the exponent in this class is
\r
740 /// always a 32-bit integer. The precision therefore specifies how many bits
\r
741 /// the mantissa will have.
\r
743 public PrecisionSpec Precision
\r
745 get { return mantissa.Precision; }
\r
749 /// Writable property:
\r
750 /// true iff the number is negative or in some cases zero (<0)
\r
751 /// false iff the number if positive or in some cases zero (>0)
\r
755 get { return mantissa.Sign; }
\r
756 set { mantissa.Sign = value; }
\r
760 /// Read-only property.
\r
761 /// True if the number is NAN, INF_PLUS, INF_MINUS or ZERO
\r
762 /// False if the number has any other value.
\r
764 public bool IsSpecialValue
\r
768 return (exponent == Int32.MaxValue || mantissa.IsZero());
\r
773 /// Read-only property, returns the type of number this is. Special values include:
\r
775 /// NONE - a regular number
\r
777 /// NAN - Not a Number (some operations will return this if their inputs are out of range)
\r
778 /// INF_PLUS - Positive infinity, not really a number, but a valid input to and output of some functions.
\r
779 /// INF_MINUS - Negative infinity, not really a number, but a valid input to and output of some functions.
\r
781 public SpecialValueType SpecialValue
\r
785 if (exponent == Int32.MaxValue)
\r
787 if (mantissa.IsZero())
\r
789 if (mantissa.Sign) return SpecialValueType.INF_MINUS;
\r
790 return SpecialValueType.INF_PLUS;
\r
793 return SpecialValueType.NAN;
\r
797 if (mantissa.IsZero()) return SpecialValueType.ZERO;
\r
798 return SpecialValueType.NONE;
\r
803 //******************** Mathematical Constants *******************
\r
806 /// Gets pi to the indicated precision
\r
808 /// <param name="precision">The precision to perform the calculation to</param>
\r
809 /// <returns>pi (the ratio of the area of a circle to its diameter)</returns>
\r
810 public static BigFloat GetPi(PrecisionSpec precision)
\r
812 if (pi == null || precision.NumBits <= pi.mantissa.Precision.NumBits)
\r
814 CalculatePi(precision.NumBits);
\r
817 BigFloat ret = new BigFloat (precision);
\r
824 /// Get e to the indicated precision
\r
826 /// <param name="precision">The preicision to perform the calculation to</param>
\r
827 /// <returns>e (the number for which the d/dx(e^x) = e^x)</returns>
\r
828 public static BigFloat GetE(PrecisionSpec precision)
\r
830 if (eCache == null || eCache.mantissa.Precision.NumBits < precision.NumBits)
\r
832 CalculateEOnly(precision.NumBits);
\r
833 //CalculateFactorials(precision.NumBits);
\r
836 BigFloat ret = new BigFloat(precision);
\r
837 ret.Assign(eCache);
\r
843 //******************** Arithmetic Functions ********************
\r
846 /// Addition (this = this + n2)
\r
848 /// <param name="n2">The number to add</param>
\r
849 public void Add(BigFloat n2)
\r
851 if (SpecialValueAddTest(n2)) return;
\r
853 if (scratch.Precision.NumBits != n2.mantissa.Precision.NumBits)
\r
855 scratch = new BigInt(n2.mantissa.Precision);
\r
858 if (exponent <= n2.exponent)
\r
860 int diff = n2.exponent - exponent;
\r
861 exponent = n2.exponent;
\r
865 mantissa.RSH(diff);
\r
868 uint carry = mantissa.Add(n2.mantissa);
\r
873 mantissa.SetBit(mantissa.Precision.NumBits - 1);
\r
877 exponent -= mantissa.Normalise();
\r
881 int diff = exponent - n2.exponent;
\r
883 scratch.Assign(n2.mantissa);
\r
886 uint carry = scratch.Add(mantissa);
\r
891 scratch.SetBit(mantissa.Precision.NumBits - 1);
\r
895 mantissa.Assign(scratch);
\r
897 exponent -= mantissa.Normalise();
\r
902 /// Subtraction (this = this - n2)
\r
904 /// <param name="n2">The number to subtract from this</param>
\r
905 public void Sub(BigFloat n2)
\r
907 n2.mantissa.Sign = !n2.mantissa.Sign;
\r
909 n2.mantissa.Sign = !n2.mantissa.Sign;
\r
913 /// Multiplication (this = this * n2)
\r
915 /// <param name="n2">The number to multiply this by</param>
\r
916 public void Mul(BigFloat n2)
\r
918 if (SpecialValueMulTest(n2)) return;
\r
920 //Anything times 0 = 0
\r
921 if (n2.mantissa.IsZero())
\r
923 mantissa.Assign(n2.mantissa);
\r
928 mantissa.MulHi(n2.mantissa);
\r
929 int shift = mantissa.Normalise();
\r
930 exponent = exponent + n2.exponent + 1 - shift;
\r
934 /// Division (this = this / n2)
\r
936 /// <param name="n2">The number to divide this by</param>
\r
937 public void Div(BigFloat n2)
\r
939 if (SpecialValueDivTest(n2)) return;
\r
941 if (mantissa.Precision.NumBits >= 8192)
\r
943 BigFloat rcp = n2.Reciprocal();
\r
948 int shift = mantissa.DivAndShift(n2.mantissa);
\r
949 exponent = exponent - (n2.exponent + shift);
\r
954 /// Multiply by a power of 2 (-ve implies division)
\r
956 /// <param name="pow2"></param>
\r
957 public void MulPow2(int pow2)
\r
963 /// Division-based reciprocal, fastest for small precisions up to 15,000 bits.
\r
965 /// <returns>The reciprocal 1/this</returns>
\r
966 public BigFloat Reciprocal()
\r
968 if (mantissa.Precision.NumBits >= 8192) return ReciprocalNewton();
\r
970 BigFloat reciprocal = new BigFloat(1u, mantissa.Precision);
\r
971 reciprocal.Div(this);
\r
976 /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.
\r
978 /// <returns>The reciprocal 1/this</returns>
\r
979 public BigFloat ReciprocalNewton()
\r
981 if (mantissa.IsZero())
\r
983 exponent = Int32.MaxValue;
\r
987 bool oldSign = mantissa.Sign;
\r
988 int oldExponent = exponent;
\r
990 //Kill exponent for now (will re-institute later)
\r
993 bool topBit = mantissa.IsTopBitOnlyBit();
\r
995 PrecisionSpec curPrec = new PrecisionSpec(32, PrecisionSpec.BaseType.BIN);
\r
997 BigFloat reciprocal = new BigFloat(curPrec);
\r
998 BigFloat constant2 = new BigFloat(curPrec);
\r
999 BigFloat temp = new BigFloat(curPrec);
\r
1000 BigFloat thisPrec = new BigFloat(this, curPrec);
\r
1002 reciprocal.exponent = 1;
\r
1003 reciprocal.mantissa.SetHighDigit(3129112985u);
\r
1005 constant2.exponent = 1;
\r
1006 constant2.mantissa.SetHighDigit(0x80000000u);
\r
1008 //D is deliberately left negative for all the following operations.
\r
1009 thisPrec.mantissa.Sign = true;
\r
1011 //Initial estimate.
\r
1012 reciprocal.Add(thisPrec);
\r
1014 //mantissa.Sign = false;
\r
1016 //Shift down into 0.5 < this < 1 range
\r
1017 thisPrec.mantissa.RSH(1);
\r
1020 int accuracyBits = 2;
\r
1021 int mantissaBits = mantissa.Precision.NumBits;
\r
1023 //Each iteration is a pass of newton's method for RCP.
\r
1024 //The is a substantial optimisation to be done here...
\r
1025 //You can double the number of bits for the calculations
\r
1026 //at each iteration, meaning that the whole process only
\r
1027 //takes some constant multiplier of the time for the
\r
1028 //full-scale multiplication.
\r
1029 while (accuracyBits < mantissaBits)
\r
1031 //Increase the precision as needed
\r
1032 if (accuracyBits >= curPrec.NumBits / 2)
\r
1034 int newBits = curPrec.NumBits * 2;
\r
1035 if (newBits > mantissaBits) newBits = mantissaBits;
\r
1036 curPrec = new PrecisionSpec(newBits, PrecisionSpec.BaseType.BIN);
\r
1038 reciprocal = new BigFloat(reciprocal, curPrec);
\r
1040 constant2 = new BigFloat(curPrec);
\r
1041 constant2.exponent = 1;
\r
1042 constant2.mantissa.SetHighDigit(0x80000000u);
\r
1044 temp = new BigFloat(temp, curPrec);
\r
1046 thisPrec = new BigFloat(this, curPrec);
\r
1047 thisPrec.mantissa.Sign = true;
\r
1048 thisPrec.mantissa.RSH(1);
\r
1052 temp.exponent = reciprocal.exponent;
\r
1053 temp.mantissa.Assign(reciprocal.mantissa);
\r
1055 temp.Mul(thisPrec);
\r
1056 //temp = -Xn * D + 2 (= 2 - Xn * D)
\r
1057 temp.Add(constant2);
\r
1058 //reciprocal = X(n+1) = Xn * (2 - Xn * D)
\r
1059 reciprocal.Mul(temp);
\r
1061 accuracyBits *= 2;
\r
1064 //'reciprocal' is now the reciprocal of the shifted down, zero-exponent mantissa of 'this'
\r
1065 //Restore the mantissa.
\r
1066 //mantissa.LSH(1);
\r
1067 exponent = oldExponent;
\r
1068 //mantissa.Sign = oldSign;
\r
1072 reciprocal.exponent = -(oldExponent);
\r
1076 reciprocal.exponent = -(oldExponent + 1);
\r
1078 reciprocal.mantissa.Sign = oldSign;
\r
1080 return reciprocal;
\r
1084 /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.
\r
1086 /// <returns>The reciprocal 1/this</returns>
\r
1087 private BigFloat ReciprocalNewton2()
\r
1089 if (mantissa.IsZero())
\r
1091 exponent = Int32.MaxValue;
\r
1095 bool oldSign = mantissa.Sign;
\r
1096 int oldExponent = exponent;
\r
1098 //Kill exponent for now (will re-institute later)
\r
1101 BigFloat reciprocal = new BigFloat(mantissa.Precision);
\r
1102 BigFloat constant2 = new BigFloat(mantissa.Precision);
\r
1103 BigFloat temp = new BigFloat(mantissa.Precision);
\r
1105 reciprocal.exponent = 1;
\r
1106 reciprocal.mantissa.SetHighDigit(3129112985u);
\r
1108 constant2.exponent = 1;
\r
1109 constant2.mantissa.SetHighDigit(0x80000000u);
\r
1111 //D is deliberately left negative for all the following operations.
\r
1112 mantissa.Sign = true;
\r
1114 //Initial estimate.
\r
1115 reciprocal.Add(this);
\r
1117 //mantissa.Sign = false;
\r
1119 //Shift down into 0.5 < this < 1 range
\r
1123 int accuracyBits = 2;
\r
1124 int mantissaBits = mantissa.Precision.NumBits;
\r
1126 //Each iteration is a pass of newton's method for RCP.
\r
1127 //The is a substantial optimisation to be done here...
\r
1128 //You can double the number of bits for the calculations
\r
1129 //at each iteration, meaning that the whole process only
\r
1130 //takes some constant multiplier of the time for the
\r
1131 //full-scale multiplication.
\r
1132 while (accuracyBits < mantissaBits)
\r
1135 temp.exponent = reciprocal.exponent;
\r
1136 temp.mantissa.Assign(reciprocal.mantissa);
\r
1139 //temp = -Xn * D + 2 (= 2 - Xn * D)
\r
1140 temp.Add(constant2);
\r
1141 //reciprocal = X(n+1) = Xn * (2 - Xn * D)
\r
1142 reciprocal.Mul(temp);
\r
1144 accuracyBits *= 2;
\r
1147 //'reciprocal' is now the reciprocal of the shifted down, zero-exponent mantissa of 'this'
\r
1148 //Restore the mantissa.
\r
1150 exponent = oldExponent;
\r
1151 mantissa.Sign = oldSign;
\r
1153 reciprocal.exponent = -(oldExponent + 1);
\r
1154 reciprocal.mantissa.Sign = oldSign;
\r
1156 return reciprocal;
\r
1160 /// Sets this equal to the input
\r
1162 /// <param name="n2"></param>
\r
1163 public void Assign(BigFloat n2)
\r
1165 exponent = n2.exponent;
\r
1166 if (mantissa.AssignHigh(n2.mantissa)) exponent++;
\r
1170 //********************* Comparison Functions *******************
\r
1173 /// Greater than comparison
\r
1175 /// <param name="n2">the number to compare this to</param>
\r
1176 /// <returns>true iff this is greater than n2 (this > n2)</returns>
\r
1177 public bool GreaterThan(BigFloat n2)
\r
1179 if (IsSpecialValue || n2.IsSpecialValue)
\r
1181 SpecialValueType s1 = SpecialValue;
\r
1182 SpecialValueType s2 = SpecialValue;
\r
1184 if (s1 == SpecialValueType.NAN || s2 == SpecialValueType.NAN) return false;
\r
1185 if (s1 == SpecialValueType.INF_MINUS) return false;
\r
1186 if (s2 == SpecialValueType.INF_PLUS) return false;
\r
1187 if (s1 == SpecialValueType.INF_PLUS) return true;
\r
1188 if (s2 == SpecialValueType.INF_MINUS) return true;
\r
1190 if (s1 == SpecialValueType.ZERO)
\r
1192 if (s2 != SpecialValueType.ZERO && n2.Sign)
\r
1202 if (s2 == SpecialValueType.ZERO)
\r
1208 if (!mantissa.Sign && n2.mantissa.Sign) return true;
\r
1209 if (mantissa.Sign && !n2.mantissa.Sign) return false;
\r
1210 if (!mantissa.Sign)
\r
1212 if (exponent > n2.exponent) return true;
\r
1213 if (exponent < n2.exponent) return false;
\r
1215 if (mantissa.Sign)
\r
1217 if (exponent > n2.exponent) return false;
\r
1218 if (exponent < n2.exponent) return true;
\r
1221 return mantissa.GreaterThan(n2.mantissa);
\r
1225 /// Less than comparison
\r
1227 /// <param name="n2">the number to compare this to</param>
\r
1228 /// <returns>true iff this is less than n2 (this < n2)</returns>
\r
1229 public bool LessThan(BigFloat n2)
\r
1231 if (IsSpecialValue || n2.IsSpecialValue)
\r
1233 SpecialValueType s1 = SpecialValue;
\r
1234 SpecialValueType s2 = SpecialValue;
\r
1236 if (s1 == SpecialValueType.NAN || s2 == SpecialValueType.NAN) return false;
\r
1237 if (s1 == SpecialValueType.INF_PLUS) return false;
\r
1238 if (s2 == SpecialValueType.INF_PLUS) return true;
\r
1239 if (s2 == SpecialValueType.INF_MINUS) return false;
\r
1240 if (s1 == SpecialValueType.INF_MINUS) return true;
\r
1242 if (s1 == SpecialValueType.ZERO)
\r
1244 if (s2 != SpecialValueType.ZERO && !n2.Sign)
\r
1254 if (s2 == SpecialValueType.ZERO)
\r
1260 if (!mantissa.Sign && n2.mantissa.Sign) return false;
\r
1261 if (mantissa.Sign && !n2.mantissa.Sign) return true;
\r
1262 if (!mantissa.Sign)
\r
1264 if (exponent > n2.exponent) return false;
\r
1265 if (exponent < n2.exponent) return true;
\r
1267 if (mantissa.Sign)
\r
1269 if (exponent > n2.exponent) return true;
\r
1270 if (exponent < n2.exponent) return false;
\r
1273 return mantissa.LessThan(n2.mantissa);
\r
1277 /// Greater than comparison
\r
1279 /// <param name="i">the number to compare this to</param>
\r
1280 /// <returns>true iff this is greater than n2 (this > n2)</returns>
\r
1281 public bool GreaterThan(int i)
\r
1283 BigFloat integer = new BigFloat(i, mantissa.Precision);
\r
1284 return GreaterThan(integer);
\r
1288 /// Less than comparison
\r
1290 /// <param name="i">the number to compare this to</param>
\r
1291 /// <returns>true iff this is less than n2 (this < n2)</returns>
\r
1292 public bool LessThan(int i)
\r
1294 BigFloat integer = new BigFloat(i, mantissa.Precision);
\r
1295 return LessThan(integer);
\r
1299 /// Compare to zero
\r
1301 /// <returns>true if this is zero (this == 0)</returns>
\r
1302 public bool IsZero()
\r
1304 return (mantissa.IsZero());
\r
1308 //******************** Mathematical Functions ******************
\r
1311 /// Sets the number to the biggest integer numerically closer to zero, if possible.
\r
1313 public void Floor()
\r
1315 //Already an integer.
\r
1316 if (exponent >= mantissa.Precision.NumBits) return;
\r
1320 mantissa.ZeroBits(mantissa.Precision.NumBits);
\r
1325 mantissa.ZeroBits(mantissa.Precision.NumBits - (exponent + 1));
\r
1329 /// Sets the number to its fractional component (equivalent to 'this' - (int)'this')
\r
1331 public void FPart()
\r
1333 //Already fractional
\r
1339 //Has no fractional part
\r
1340 if (exponent >= mantissa.Precision.NumBits)
\r
1347 mantissa.ZeroBitsHigh(exponent + 1);
\r
1348 exponent -= mantissa.Normalise();
\r
1352 /// Calculates tan(x)
\r
1356 if (IsSpecialValue)
\r
1358 //Tan(x) has no limit as x->inf
\r
1359 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
\r
1363 else if (SpecialValue == SpecialValueType.ZERO)
\r
1371 if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
\r
1373 CalculatePi(mantissa.Precision.NumBits);
\r
1376 //Work out the sign change (involves replicating some rescaling).
\r
1377 bool sign = mantissa.Sign;
\r
1378 mantissa.Sign = false;
\r
1380 if (mantissa.IsZero())
\r
1385 //Rescale into 0 <= x < pi
\r
1386 if (GreaterThan(pi))
\r
1388 //There will be an inherent loss of precision doing this.
\r
1389 BigFloat newAngle = new BigFloat(this);
\r
1390 newAngle.Mul(piRecip);
\r
1396 //Rescale to -pi/2 <= x < pi/2
\r
1397 if (!LessThan(piBy2))
\r
1402 //Now the sign of the sin determines the sign of the tan.
\r
1403 //tan(x) = sin(x) / sqrt(1 - sin^2(x))
\r
1405 BigFloat denom = new BigFloat(this);
\r
1407 denom.Sub(new BigFloat(1, mantissa.Precision));
\r
1408 denom.mantissa.Sign = !denom.mantissa.Sign;
\r
1410 if (denom.mantissa.Sign)
\r
1417 if (sign) mantissa.Sign = !mantissa.Sign;
\r
1421 /// Calculates Cos(x)
\r
1425 if (IsSpecialValue)
\r
1427 //Cos(x) has no limit as x->inf
\r
1428 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
\r
1432 else if (SpecialValue == SpecialValueType.ZERO)
\r
1434 Assign(new BigFloat(1, mantissa.Precision));
\r
1440 if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
\r
1442 CalculatePi(mantissa.Precision.NumBits);
\r
1450 /// Calculates Sin(x):
\r
1451 /// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].
\r
1455 if (IsSpecialValue)
\r
1457 //Sin(x) has no limit as x->inf
\r
1458 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
\r
1466 //Convert to positive range (0 <= x < inf)
\r
1467 bool sign = mantissa.Sign;
\r
1468 mantissa.Sign = false;
\r
1470 if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
\r
1472 CalculatePi(mantissa.Precision.NumBits);
\r
1475 if (inverseFactorialCache == null || invFactorialCutoff != mantissa.Precision.NumBits)
\r
1477 CalculateFactorials(mantissa.Precision.NumBits);
\r
1480 //Rescale into 0 <= x < 2*pi
\r
1481 if (GreaterThan(twoPi))
\r
1483 //There will be an inherent loss of precision doing this.
\r
1484 BigFloat newAngle = new BigFloat(this);
\r
1485 newAngle.Mul(twoPiRecip);
\r
1487 newAngle.Mul(twoPi);
\r
1491 //Rescale into range 0 <= x < pi
\r
1492 if (GreaterThan(pi))
\r
1494 //sin(pi + a) = sin(pi)cos(a) + sin(a)cos(pi) = 0 - sin(a) = -sin(a)
\r
1499 BigFloat temp = new BigFloat(mantissa.Precision);
\r
1501 //Rescale into range 0 <= x < pi/2
\r
1502 if (GreaterThan(piBy2))
\r
1504 temp.Assign(this);
\r
1509 //Rescale into range 0 <= x < pi/6 to accelerate convergence.
\r
1510 //This is done using sin(3x) = 3sin(x) - 4sin^3(x)
\r
1513 if (mantissa.IsZero())
\r
1519 BigFloat term = new BigFloat(this);
\r
1521 BigFloat square = new BigFloat(this);
\r
1524 BigFloat sum = new BigFloat(this);
\r
1526 bool termSign = true;
\r
1527 int length = inverseFactorialCache.Length;
\r
1528 int numBits = mantissa.Precision.NumBits;
\r
1530 for (int i = 3; i < length; i += 2)
\r
1533 temp.Assign(inverseFactorialCache[i]);
\r
1535 temp.mantissa.Sign = termSign;
\r
1536 termSign = !termSign;
\r
1538 if (temp.exponent < -numBits) break;
\r
1543 //Restore the triple-angle: sin(3x) = 3sin(x) - 4sin^3(x)
\r
1547 Mul(new BigFloat(3, mantissa.Precision));
\r
1548 sum.exponent += 2;
\r
1551 //Restore the sign
\r
1552 mantissa.Sign = sign;
\r
1556 /// Hyperbolic Sin (sinh) function
\r
1558 public void Sinh()
\r
1560 if (IsSpecialValue)
\r
1566 Sub(Reciprocal());
\r
1571 /// Hyperbolic cosine (cosh) function
\r
1573 public void Cosh()
\r
1575 if (IsSpecialValue)
\r
1577 if (SpecialValue == SpecialValueType.ZERO)
\r
1579 Assign(new BigFloat(1, mantissa.Precision));
\r
1581 else if (SpecialValue == SpecialValueType.INF_MINUS)
\r
1590 Add(Reciprocal());
\r
1595 /// Hyperbolic tangent function (tanh)
\r
1597 public void Tanh()
\r
1599 if (IsSpecialValue)
\r
1601 if (SpecialValue == SpecialValueType.INF_MINUS)
\r
1603 Assign(new BigFloat(-1, mantissa.Precision));
\r
1605 else if (SpecialValue == SpecialValueType.INF_PLUS)
\r
1607 Assign(new BigFloat(1, mantissa.Precision));
\r
1615 BigFloat temp = new BigFloat(this);
\r
1616 BigFloat one = new BigFloat(1, mantissa.Precision);
\r
1623 /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)
\r
1625 public void Arcsin()
\r
1627 if (IsSpecialValue)
\r
1629 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)
\r
1638 BigFloat one = new BigFloat(1, mantissa.Precision);
\r
1639 BigFloat plusABit = new BigFloat(1, mantissa.Precision);
\r
1640 plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));
\r
1641 BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);
\r
1642 onePlusABit.Add(plusABit);
\r
1644 bool sign = mantissa.Sign;
\r
1645 mantissa.Sign = false;
\r
1647 if (GreaterThan(onePlusABit))
\r
1651 else if (LessThan(one))
\r
1653 BigFloat temp = new BigFloat(this);
\r
1656 temp.mantissa.Sign = !temp.mantissa.Sign;
\r
1662 mantissa.Sign = sign;
\r
1666 if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
\r
1668 CalculatePi(mantissa.Precision.NumBits);
\r
1672 if (sign) mantissa.Sign = true;
\r
1677 /// arccos(): the inverse function of cos(), range (0..pi)
\r
1679 public void Arccos()
\r
1681 if (IsSpecialValue)
\r
1683 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)
\r
1687 else if (SpecialValue == SpecialValueType.ZERO)
\r
1689 Assign(new BigFloat(1, mantissa.Precision));
\r
1697 BigFloat one = new BigFloat(1, mantissa.Precision);
\r
1698 BigFloat plusABit = new BigFloat(1, mantissa.Precision);
\r
1699 plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));
\r
1700 BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);
\r
1701 onePlusABit.Add(plusABit);
\r
1703 bool sign = mantissa.Sign;
\r
1704 mantissa.Sign = false;
\r
1706 if (GreaterThan(onePlusABit))
\r
1710 else if (LessThan(one))
\r
1712 if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
\r
1714 CalculatePi(mantissa.Precision.NumBits);
\r
1717 mantissa.Sign = sign;
\r
1718 BigFloat temp = new BigFloat(this);
\r
1721 mantissa.Sign = !mantissa.Sign;
\r
1732 if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
\r
1734 CalculatePi(mantissa.Precision.NumBits);
\r
1748 /// arctan(): the inverse function of sin(), range of (-pi/2..pi/2)
\r
1750 public void Arctan()
\r
1752 //With 2 argument reductions, we increase precision by a minimum of 4 bits per term.
\r
1753 int numBits = mantissa.Precision.NumBits;
\r
1754 int maxTerms = numBits >> 2;
\r
1756 if (pi == null || pi.mantissa.Precision.NumBits != numBits)
\r
1758 CalculatePi(mantissa.Precision.NumBits);
\r
1761 //Make domain positive
\r
1762 bool sign = mantissa.Sign;
\r
1763 mantissa.Sign = false;
\r
1765 if (IsSpecialValue)
\r
1767 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
\r
1770 mantissa.Sign = sign;
\r
1777 if (reciprocals == null || reciprocals[0].mantissa.Precision.NumBits != numBits || reciprocals.Length < maxTerms)
\r
1779 CalculateReciprocals(numBits, maxTerms);
\r
1782 bool invert = false;
\r
1783 BigFloat one = new BigFloat(1, mantissa.Precision);
\r
1785 //Invert if outside of convergence
\r
1786 if (GreaterThan(one))
\r
1789 Assign(Reciprocal());
\r
1792 //Reduce using half-angle formula:
\r
1793 //arctan(2x) = 2 arctan (x / (1 + sqrt(1 + x)))
\r
1795 //First reduction (guarantees 2 bits per iteration)
\r
1796 BigFloat temp = new BigFloat(this);
\r
1803 //Second reduction (guarantees 4 bits per iteration)
\r
1804 temp.Assign(this);
\r
1811 //Actual series calculation
\r
1812 int length = reciprocals.Length;
\r
1813 BigFloat term = new BigFloat(this);
\r
1816 BigFloat pow = new BigFloat(this);
\r
1819 BigFloat sum = new BigFloat(this);
\r
1821 for (int i = 1; i < length; i++)
\r
1823 //u(n) = u(n-1) * x^2
\r
1824 //t(n) = u(n) / (2n+1)
\r
1826 term.Sign = !term.Sign;
\r
1827 temp.Assign(term);
\r
1828 temp.Mul(reciprocals[i]);
\r
1830 if (temp.exponent < -numBits) break;
\r
1835 //Undo the reductions.
\r
1841 //Assign(Reciprocal());
\r
1842 mantissa.Sign = true;
\r
1848 mantissa.Sign = sign;
\r
1853 /// Arcsinh(): the inverse sinh function
\r
1855 public void Arcsinh()
\r
1857 //Just let all special values fall through
\r
1858 if (IsSpecialValue)
\r
1863 BigFloat temp = new BigFloat(this);
\r
1865 temp.Add(new BigFloat(1, mantissa.Precision));
\r
1872 /// Arccosh(): the inverse cosh() function
\r
1874 public void Arccosh()
\r
1876 //acosh isn't defined for x < 1
\r
1877 if (IsSpecialValue)
\r
1879 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.ZERO)
\r
1888 BigFloat one = new BigFloat(1, mantissa.Precision);
\r
1889 if (LessThan(one))
\r
1895 BigFloat temp = new BigFloat(this);
\r
1904 /// Arctanh(): the inverse tanh function
\r
1906 public void Arctanh()
\r
1908 //|x| <= 1 for a non-NaN output
\r
1909 if (IsSpecialValue)
\r
1911 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
\r
1920 BigFloat one = new BigFloat(1, mantissa.Precision);
\r
1921 BigFloat plusABit = new BigFloat(1, mantissa.Precision);
\r
1922 plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));
\r
1923 BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);
\r
1924 onePlusABit.Add(plusABit);
\r
1926 bool sign = mantissa.Sign;
\r
1927 mantissa.Sign = false;
\r
1929 if (GreaterThan(onePlusABit))
\r
1933 else if (LessThan(one))
\r
1935 BigFloat temp = new BigFloat(this);
\r
1941 mantissa.Sign = sign;
\r
1957 /// Two-variable iterative square root, taken from
\r
1958 /// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method
\r
1960 public void Sqrt()
\r
1962 if (mantissa.Sign || IsSpecialValue)
\r
1964 if (SpecialValue == SpecialValueType.ZERO)
\r
1969 if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)
\r
1978 BigFloat temp3 = new BigFloat(mantissa.Precision);
\r
1979 BigFloat three = new BigFloat(3, mantissa.Precision);
\r
1981 int exponentScale = 0;
\r
1983 //Rescale to 0.5 <= x < 2
\r
1984 if (exponent < -1)
\r
1986 int diff = -exponent;
\r
1987 if ((diff & 1) != 0)
\r
1992 exponentScale = -diff;
\r
1995 else if (exponent > 0)
\r
1997 if ((exponent & 1) != 0)
\r
1999 exponentScale = exponent + 1;
\r
2004 exponentScale = exponent;
\r
2009 temp2 = new BigFloat(this);
\r
2010 temp2.Sub(new BigFloat(1, mantissa.Precision));
\r
2012 //if (temp2.mantissa.IsZero())
\r
2014 // exponent += exponentScale;
\r
2018 int numBits = mantissa.Precision.NumBits;
\r
2020 while ((exponent - temp2.exponent) < numBits && temp2.SpecialValue != SpecialValueType.ZERO)
\r
2022 //a(n+1) = an - an*cn / 2
\r
2023 temp3.Assign(this);
\r
2025 temp3.MulPow2(-1);
\r
2028 //c(n+1) = cn^2 * (cn - 3) / 4
\r
2029 temp3.Assign(temp2);
\r
2033 temp2.MulPow2(-2);
\r
2036 exponent += (exponentScale >> 1);
\r
2040 /// The natural logarithm, ln(x)
\r
2044 if (IsSpecialValue || mantissa.Sign)
\r
2046 if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)
\r
2050 else if (SpecialValue == SpecialValueType.ZERO)
\r
2058 if (mantissa.Precision.NumBits >= 512)
\r
2065 if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
\r
2067 CalculateLog2(mantissa.Precision.NumBits);
\r
2075 /// Log to the base 10
\r
2077 public void Log10()
\r
2079 if (IsSpecialValue || mantissa.Sign)
\r
2081 if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)
\r
2085 else if (SpecialValue == SpecialValueType.ZERO)
\r
2094 if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
\r
2096 CalculateLog2(mantissa.Precision.NumBits);
\r
2104 /// The exponential function. Less accurate for high exponents, scales poorly with the number
\r
2109 Exp(mantissa.Precision.NumBits);
\r
2113 /// Raises a number to an integer power (positive or negative). This is a very accurate and fast function,
\r
2114 /// comparable to or faster than division (although it is slightly slower for
\r
2115 /// negative powers, obviously)
\r
2118 /// <param name="power"></param>
\r
2119 public void Pow(int power)
\r
2121 BigFloat acc = new BigFloat(1, mantissa.Precision);
\r
2122 BigFloat temp = new BigFloat(1, mantissa.Precision);
\r
2124 int powerTemp = power;
\r
2128 Assign(Reciprocal());
\r
2129 powerTemp = -power;
\r
2132 //Fast power function
\r
2133 while (powerTemp != 0)
\r
2138 if ((powerTemp & 1) != 0)
\r
2150 /// Raises to an aribitrary power. This is both slow (uses Log) and inaccurate. If you need to
\r
2151 /// raise e^x use exp(). If you need an integer power, use the integer power function Pow(int)
\r
2152 /// Accuracy Note:
\r
2153 /// The function is only ever accurate to a maximum of 4 decimal digits
\r
2154 /// For every 10x larger (or smaller) the power gets, you lose an additional decimal digit
\r
2155 /// If you really need a precise result, do the calculation with an extra 32-bits and round
\r
2157 /// This only works for powers of positive real numbers. Negative numbers will fail.
\r
2159 /// <param name="power"></param>
\r
2160 public void Pow(BigFloat power)
\r
2168 //******************** Static Math Functions *******************
\r
2171 /// Returns the integer component of the input
\r
2173 /// <param name="n1">The input number</param>
\r
2174 /// <remarks>The integer component returned will always be numerically closer to zero
\r
2175 /// than the input: an input of -3.49 for instance would produce a value of 3.</remarks>
\r
2176 public static BigFloat Floor(BigFloat n1)
\r
2178 BigFloat res = new BigFloat(n1);
\r
2184 /// Returns the fractional (non-integer component of the input)
\r
2186 /// <param name="n1">The input number</param>
\r
2187 public static BigFloat FPart(BigFloat n1)
\r
2189 BigFloat res = new BigFloat(n1);
\r
2195 /// Calculates tan(x)
\r
2197 /// <param name="n1">The angle (in radians) to find the tangent of</param>
\r
2198 public static BigFloat Tan(BigFloat n1)
\r
2200 BigFloat res = new BigFloat(n1);
\r
2206 /// Calculates Cos(x)
\r
2208 /// <param name="n1">The angle (in radians) to find the cosine of</param>
\r
2209 /// <remarks>This is a reasonably fast function for smaller precisions, but
\r
2210 /// doesn't scale well for higher precision arguments</remarks>
\r
2211 public static BigFloat Cos(BigFloat n1)
\r
2213 BigFloat res = new BigFloat(n1);
\r
2219 /// Calculates Sin(x):
\r
2220 /// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].
\r
2222 /// <param name="n1">The angle to find the sine of (in radians)</param>
\r
2223 /// <remarks>This is a resonably fast function, for smaller precision arguments, but doesn't
\r
2224 /// scale very well with the number of bits in the input.</remarks>
\r
2225 public static BigFloat Sin(BigFloat n1)
\r
2227 BigFloat res = new BigFloat(n1);
\r
2233 /// Hyperbolic Sin (sinh) function
\r
2235 /// <param name="n1">The number to find the hyperbolic sine of</param>
\r
2236 public static BigFloat Sinh(BigFloat n1)
\r
2238 BigFloat res = new BigFloat(n1);
\r
2244 /// Hyperbolic cosine (cosh) function
\r
2246 /// <param name="n1">The number to find the hyperbolic cosine of</param>
\r
2247 public static BigFloat Cosh(BigFloat n1)
\r
2249 BigFloat res = new BigFloat(n1);
\r
2255 /// Hyperbolic tangent function (tanh)
\r
2257 /// <param name="n1">The number to find the hyperbolic tangent of</param>
\r
2258 public static BigFloat Tanh(BigFloat n1)
\r
2260 BigFloat res = new BigFloat(n1);
\r
2266 /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)
\r
2268 /// <param name="n1">The number to find the arcsine of (-pi/2..pi/2)</param>
\r
2269 /// <remarks>Note that inverse trig functions are only defined within a specific range.
\r
2270 /// Values outside this range will return NaN, although some margin for error is assumed.
\r
2272 public static BigFloat Arcsin(BigFloat n1)
\r
2274 BigFloat res = new BigFloat(n1);
\r
2280 /// arccos(): the inverse function of cos(), input range (0..pi)
\r
2282 /// <param name="n1">The number to find the arccosine of (0..pi)</param>
\r
2283 /// <remarks>Note that inverse trig functions are only defined within a specific range.
\r
2284 /// Values outside this range will return NaN, although some margin for error is assumed.
\r
2286 public static BigFloat Arccos(BigFloat n1)
\r
2288 BigFloat res = new BigFloat(n1);
\r
2294 /// arctan(): the inverse function of sin(), input range of (-pi/2..pi/2)
\r
2296 /// <param name="n1">The number to find the arctangent of (-pi/2..pi/2)</param>
\r
2297 /// <remarks>Note that inverse trig functions are only defined within a specific range.
\r
2298 /// Values outside this range will return NaN, although some margin for error is assumed.
\r
2300 public static BigFloat Arctan(BigFloat n1)
\r
2302 BigFloat res = new BigFloat(n1);
\r
2308 /// Arcsinh(): the inverse sinh function
\r
2310 /// <param name="n1">The number to find the inverse hyperbolic sine of</param>
\r
2311 public static BigFloat Arcsinh(BigFloat n1)
\r
2313 BigFloat res = new BigFloat(n1);
\r
2319 /// Arccosh(): the inverse cosh() function
\r
2321 /// <param name="n1">The number to find the inverse hyperbolic cosine of</param>
\r
2322 public static BigFloat Arccosh(BigFloat n1)
\r
2324 BigFloat res = new BigFloat(n1);
\r
2330 /// Arctanh(): the inverse tanh function
\r
2332 /// <param name="n1">The number to fine the inverse hyperbolic tan of</param>
\r
2333 public static BigFloat Arctanh(BigFloat n1)
\r
2335 BigFloat res = new BigFloat(n1);
\r
2341 /// Two-variable iterative square root, taken from
\r
2342 /// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method
\r
2344 /// <remarks>This is quite a fast function, as elementary functions go. You can expect it to take
\r
2345 /// about twice as long as a floating-point division.
\r
2347 public static BigFloat Sqrt(BigFloat n1)
\r
2349 BigFloat res = new BigFloat(n1);
\r
2355 /// The natural logarithm, ln(x) (log base e)
\r
2357 /// <remarks>This is a very slow function, despite repeated attempts at optimisation.
\r
2358 /// To make it any faster, different strategies would be needed for integer operations.
\r
2359 /// It does, however, scale well with the number of bits.
\r
2361 /// <param name="n1">The number to find the natural logarithm of</param>
\r
2362 public static BigFloat Log(BigFloat n1)
\r
2364 BigFloat res = new BigFloat(n1);
\r
2370 /// Base 10 logarithm of a number
\r
2372 /// <remarks>This is a very slow function, despite repeated attempts at optimisation.
\r
2373 /// To make it any faster, different strategies would be needed for integer operations.
\r
2374 /// It does, however, scale well with the number of bits.
\r
2376 /// <param name="n1">The number to find the base 10 logarithm of</param>
\r
2377 public static BigFloat Log10(BigFloat n1)
\r
2379 BigFloat res = new BigFloat(n1);
\r
2385 /// The exponential function. Less accurate for high exponents, scales poorly with the number
\r
2386 /// of bits. This is quite fast for low-precision arguments.
\r
2388 public static BigFloat Exp(BigFloat n1)
\r
2390 BigFloat res = new BigFloat(n1);
\r
2396 /// Raises a number to an integer power (positive or negative). This is a very accurate and fast function,
\r
2397 /// comparable to or faster than division (although it is slightly slower for
\r
2398 /// negative powers, obviously).
\r
2400 /// <param name="n1">The number to raise to the power</param>
\r
2401 /// <param name="power">The power to raise it to</param>
\r
2402 public static BigFloat Pow(BigFloat n1, int power)
\r
2404 BigFloat res = new BigFloat(n1);
\r
2410 /// Raises to an aribitrary power. This is both slow (uses Log) and inaccurate. If you need to
\r
2411 /// raise e^x use exp(). If you need an integer power, use the integer power function Pow(int)
\r
2414 /// Accuracy Note:
\r
2415 /// The function is only ever accurate to a maximum of 4 decimal digits
\r
2416 /// For every 10x larger (or smaller) the power gets, you lose an additional decimal digit
\r
2417 /// If you really need a precise result, do the calculation with an extra 32-bits and round
\r
2420 /// This only works for powers of positive real numbers. Negative numbers will fail.
\r
2422 /// <param name="n1">The number to raise to a power</param>
\r
2423 /// <param name="power">The power to raise it to</param>
\r
2424 public static BigFloat Pow(BigFloat n1, BigFloat power)
\r
2426 BigFloat res = new BigFloat(n1);
\r
2431 //********************** Static functions **********************
\r
2434 /// Adds two numbers and returns the result
\r
2436 public static BigFloat Add(BigFloat n1, BigFloat n2)
\r
2438 BigFloat ret = new BigFloat(n1);
\r
2444 /// Subtracts two numbers and returns the result
\r
2446 public static BigFloat Sub(BigFloat n1, BigFloat n2)
\r
2448 BigFloat ret = new BigFloat(n1);
\r
2454 /// Multiplies two numbers and returns the result
\r
2456 public static BigFloat Mul(BigFloat n1, BigFloat n2)
\r
2458 BigFloat ret = new BigFloat(n1);
\r
2464 /// Divides two numbers and returns the result
\r
2466 public static BigFloat Div(BigFloat n1, BigFloat n2)
\r
2468 BigFloat ret = new BigFloat(n1);
\r
2474 /// Tests whether n1 is greater than n2
\r
2476 public static bool GreaterThan(BigFloat n1, BigFloat n2)
\r
2478 return n1.GreaterThan(n2);
\r
2482 /// Tests whether n1 is less than n2
\r
2484 public static bool LessThan(BigFloat n1, BigFloat n2)
\r
2486 return n1.LessThan(n2);
\r
2490 //******************* Fast static functions ********************
\r
2493 /// Adds two numbers and assigns the result to res.
\r
2495 /// <param name="res">a pre-existing BigFloat to take the result</param>
\r
2496 /// <param name="n1">the first number</param>
\r
2497 /// <param name="n2">the second number</param>
\r
2498 /// <returns>a handle to res</returns>
\r
2499 public static BigFloat Add(BigFloat res, BigFloat n1, BigFloat n2)
\r
2507 /// Subtracts two numbers and assigns the result to res.
\r
2509 /// <param name="res">a pre-existing BigFloat to take the result</param>
\r
2510 /// <param name="n1">the first number</param>
\r
2511 /// <param name="n2">the second number</param>
\r
2512 /// <returns>a handle to res</returns>
\r
2513 public static BigFloat Sub(BigFloat res, BigFloat n1, BigFloat n2)
\r
2521 /// Multiplies two numbers and assigns the result to res.
\r
2523 /// <param name="res">a pre-existing BigFloat to take the result</param>
\r
2524 /// <param name="n1">the first number</param>
\r
2525 /// <param name="n2">the second number</param>
\r
2526 /// <returns>a handle to res</returns>
\r
2527 public static BigFloat Mul(BigFloat res, BigFloat n1, BigFloat n2)
\r
2535 /// Divides two numbers and assigns the result to res.
\r
2537 /// <param name="res">a pre-existing BigFloat to take the result</param>
\r
2538 /// <param name="n1">the first number</param>
\r
2539 /// <param name="n2">the second number</param>
\r
2540 /// <returns>a handle to res</returns>
\r
2541 public static BigFloat Div(BigFloat res, BigFloat n1, BigFloat n2)
\r
2549 //************************* Operators **************************
\r
2552 /// The addition operator
\r
2554 public static BigFloat operator +(BigFloat n1, BigFloat n2)
\r
2556 return Add(n1, n2);
\r
2560 /// The subtraction operator
\r
2562 public static BigFloat operator -(BigFloat n1, BigFloat n2)
\r
2564 return Sub(n1, n2);
\r
2568 /// The multiplication operator
\r
2570 public static BigFloat operator *(BigFloat n1, BigFloat n2)
\r
2572 return Mul(n1, n2);
\r
2576 /// The division operator
\r
2578 public static BigFloat operator /(BigFloat n1, BigFloat n2)
\r
2580 return Div(n1, n2);
\r
2583 //************************** Conversions *************************
\r
2586 /// Converts a BigFloat to an BigInt with the specified precision
\r
2588 /// <param name="n1">The number to convert</param>
\r
2589 /// <param name="precision">The precision to convert it with</param>
\r
2590 /// <param name="round">Do we round the number if we are truncating the mantissa?</param>
\r
2591 /// <returns></returns>
\r
2592 public static BigInt ConvertToInt(BigFloat n1, PrecisionSpec precision, bool round)
\r
2594 BigInt ret = new BigInt(precision);
\r
2596 int numBits = n1.mantissa.Precision.NumBits;
\r
2597 int shift = numBits - (n1.exponent + 1);
\r
2599 BigFloat copy = new BigFloat(n1);
\r
2603 if (copy.mantissa.Precision.NumBits > ret.Precision.NumBits)
\r
2607 for (int i = copy.exponent + 1; i <= ret.Precision.NumBits; i++)
\r
2609 if (copy.mantissa.GetBitFromTop(i) == 0)
\r
2619 copy.mantissa.RSH(shift);
\r
2621 else if (shift < 0)
\r
2623 copy.mantissa.LSH(-shift);
\r
2626 ret.Assign(copy.mantissa);
\r
2628 if (inc) ret.Increment();
\r
2634 /// Returns a base-10 string representing the number.
\r
2636 /// Note: This is inefficient and possibly inaccurate. Please use with enough
\r
2637 /// rounding digits (set using the RoundingDigits property) to ensure accuracy
\r
2639 public override string ToString()
\r
2641 if (IsSpecialValue)
\r
2643 SpecialValueType s = SpecialValue;
\r
2644 if (s == SpecialValueType.ZERO)
\r
2646 return String.Format("0{0}0", System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator);
\r
2648 else if (s == SpecialValueType.INF_PLUS)
\r
2650 return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol;
\r
2652 else if (s == SpecialValueType.INF_MINUS)
\r
2654 return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol;
\r
2656 else if (s == SpecialValueType.NAN)
\r
2658 return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol;
\r
2662 return "Unrecognised special type";
\r
2666 if (scratch.Precision.NumBits != mantissa.Precision.NumBits)
\r
2668 scratch = new BigInt(mantissa.Precision);
\r
2671 //The mantissa expresses 1.xxxxxxxxxxx
\r
2672 //The highest possible value for the mantissa without the implicit 1. is 0.9999999...
\r
2673 scratch.Assign(mantissa);
\r
2674 //scratch.Round(3);
\r
2675 scratch.Sign = false;
\r
2676 BigInt denom = new BigInt("0", mantissa.Precision);
\r
2677 denom.SetBit(mantissa.Precision.NumBits - 1);
\r
2679 bool useExponentialNotation = false;
\r
2680 int halfBits = mantissa.Precision.NumBits / 2;
\r
2681 if (halfBits > 60) halfBits = 60;
\r
2686 if (exponent < halfBits)
\r
2688 denom.RSH(exponent);
\r
2692 useExponentialNotation = true;
\r
2695 else if (exponent < 0)
\r
2697 int shift = -(exponent);
\r
2698 if (shift < precDec)
\r
2700 scratch.RSH(shift);
\r
2704 useExponentialNotation = true;
\r
2710 if (useExponentialNotation)
\r
2712 int absExponent = exponent;
\r
2713 if (absExponent < 0) absExponent = -absExponent;
\r
2714 int powerOf10 = (int)((double)absExponent * Math.Log10(2.0));
\r
2716 //Use 1 extra digit of precision (this is actually 32 bits more, nb)
\r
2717 BigFloat thisFloat = new BigFloat(this, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
\r
2718 thisFloat.mantissa.Sign = false;
\r
2720 //Multiplicative correction factor to bring number into range.
\r
2721 BigFloat one = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
\r
2722 BigFloat ten = new BigFloat(10, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
\r
2723 BigFloat tenRCP = ten.Reciprocal();
\r
2725 //Accumulator for the power of 10 calculation.
\r
2726 BigFloat acc = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
\r
2728 BigFloat tenToUse;
\r
2732 tenToUse = new BigFloat(tenRCP, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
\r
2736 tenToUse = new BigFloat(ten, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
\r
2739 BigFloat tenToPower = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
\r
2741 int powerTemp = powerOf10;
\r
2743 //Fast power function
\r
2744 while (powerTemp != 0)
\r
2746 tenToPower.Mul(tenToUse);
\r
2747 tenToUse.Assign(tenToPower);
\r
2749 if ((powerTemp & 1) != 0)
\r
2751 acc.Mul(tenToPower);
\r
2757 thisFloat.Mul(acc);
\r
2759 //If we are out of range, correct.
\r
2760 if (thisFloat.GreaterThan(ten))
\r
2762 thisFloat.Mul(tenRCP);
\r
2772 else if (thisFloat.LessThan(one))
\r
2774 thisFloat.Mul(ten);
\r
2785 //Restore the precision and the sign.
\r
2786 BigFloat printable = new BigFloat(thisFloat, mantissa.Precision);
\r
2787 printable.mantissa.Sign = mantissa.Sign;
\r
2788 output = printable.ToString();
\r
2790 if (exponent < 0) powerOf10 = -powerOf10;
\r
2792 output = String.Format("{0}E{1}", output, powerOf10);
\r
2796 BigInt bigDigit = BigInt.Div(scratch, denom);
\r
2797 bigDigit.Sign = false;
\r
2798 scratch.Sub(BigInt.Mul(denom, bigDigit));
\r
2800 if (mantissa.Sign)
\r
2802 output = String.Format("-{0}.", bigDigit);
\r
2806 output = String.Format("{0}.", bigDigit);
\r
2809 denom = BigInt.Div(denom, 10u);
\r
2811 while (!denom.IsZero())
\r
2813 uint digit = (uint)BigInt.Div(scratch, denom);
\r
2814 if (digit == 10) digit--;
\r
2815 scratch.Sub(BigInt.Mul(denom, digit));
\r
2816 output = String.Format("{0}{1}", output, digit);
\r
2817 denom = BigInt.Div(denom, 10u);
\r
2820 output = RoundString(output, RoundingDigits);
\r
2826 //**************** Special value handling for ops ***************
\r
2828 private void SetNaN()
\r
2830 exponent = Int32.MaxValue;
\r
2831 mantissa.SetBit(mantissa.Precision.NumBits - 1);
\r
2834 private void SetZero()
\r
2841 private void SetInfPlus()
\r
2844 exponent = Int32.MaxValue;
\r
2848 private void SetInfMinus()
\r
2851 exponent = Int32.MaxValue;
\r
2855 private bool SpecialValueAddTest(BigFloat n2)
\r
2857 if (IsSpecialValue || n2.IsSpecialValue)
\r
2859 SpecialValueType s1 = SpecialValue;
\r
2860 SpecialValueType s2 = n2.SpecialValue;
\r
2862 if (s1 == SpecialValueType.NAN) return true;
\r
2863 if (s2 == SpecialValueType.NAN)
\r
2865 //Set NaN and return.
\r
2870 if (s1 == SpecialValueType.INF_PLUS)
\r
2872 //INF+ + INF- = NAN
\r
2873 if (s2 == SpecialValueType.INF_MINUS)
\r
2882 if (s1 == SpecialValueType.INF_MINUS)
\r
2884 //INF+ + INF- = NAN
\r
2885 if (s2 == SpecialValueType.INF_PLUS)
\r
2894 if (s2 == SpecialValueType.ZERO)
\r
2899 if (s1 == SpecialValueType.ZERO)
\r
2909 private bool SpecialValueMulTest(BigFloat n2)
\r
2911 if (IsSpecialValue || n2.IsSpecialValue)
\r
2913 SpecialValueType s1 = SpecialValue;
\r
2914 SpecialValueType s2 = n2.SpecialValue;
\r
2916 if (s1 == SpecialValueType.NAN) return true;
\r
2917 if (s2 == SpecialValueType.NAN)
\r
2919 //Set NaN and return.
\r
2924 if (s1 == SpecialValueType.INF_PLUS)
\r
2926 //Inf+ * Inf- = Inf-
\r
2927 if (s2 == SpecialValueType.INF_MINUS)
\r
2934 if (s2 == SpecialValueType.ZERO)
\r
2936 //Set NaN and return.
\r
2944 if (s1 == SpecialValueType.INF_MINUS)
\r
2946 //Inf- * Inf- = Inf+
\r
2947 if (s2 == SpecialValueType.INF_MINUS)
\r
2954 if (s2 == SpecialValueType.ZERO)
\r
2956 //Set NaN and return.
\r
2964 if (s2 == SpecialValueType.ZERO)
\r
2970 if (s1 == SpecialValueType.ZERO)
\r
2979 private bool SpecialValueDivTest(BigFloat n2)
\r
2981 if (IsSpecialValue || n2.IsSpecialValue)
\r
2983 SpecialValueType s1 = SpecialValue;
\r
2984 SpecialValueType s2 = n2.SpecialValue;
\r
2986 if (s1 == SpecialValueType.NAN) return true;
\r
2987 if (s2 == SpecialValueType.NAN)
\r
2989 //Set NaN and return.
\r
2994 if ((s1 == SpecialValueType.INF_PLUS || s1 == SpecialValueType.INF_MINUS))
\r
2996 if (s2 == SpecialValueType.INF_PLUS || s2 == SpecialValueType.INF_MINUS)
\r
2998 //Set NaN and return.
\r
3005 if (s1 == SpecialValueType.INF_PLUS)
\r
3019 if (s2 == SpecialValueType.ZERO)
\r
3021 if (s1 == SpecialValueType.ZERO)
\r
3041 //****************** Internal helper functions *****************
\r
3044 /// Used for fixed point speed-ups (where the extra precision is not required). Note that Denormalised
\r
3045 /// floats break the assumptions that underly Add() and Sub(), so they can only be used for multiplication
\r
3047 /// <param name="targetExponent"></param>
\r
3048 private void Denormalise(int targetExponent)
\r
3050 int diff = targetExponent - exponent;
\r
3051 if (diff <= 0) return;
\r
3053 //This only works to reduce the precision, so if the difference implies an increase, we can't do anything.
\r
3054 mantissa.RSH(diff);
\r
3059 /// The binary logarithm, log2(x) - for precisions above 1000 bits, use Log() and convert the base.
\r
3061 private void Log2()
\r
3063 if (scratch.Precision.NumBits != mantissa.Precision.NumBits)
\r
3065 scratch = new BigInt(mantissa.Precision);
\r
3068 int bits = mantissa.Precision.NumBits;
\r
3069 BigFloat temp = new BigFloat(this);
\r
3070 BigFloat result = new BigFloat(exponent, mantissa.Precision);
\r
3071 BigFloat pow2 = new BigFloat(1, mantissa.Precision);
\r
3072 temp.exponent = 0;
\r
3073 int bitsCalculated = 0;
\r
3075 while (bitsCalculated < bits)
\r
3078 for (i = 0; (temp.exponent == 0); i++)
\r
3080 temp.mantissa.SquareHiFast(scratch);
\r
3081 int shift = temp.mantissa.Normalise();
\r
3082 temp.exponent += 1 - shift;
\r
3083 if (i + bitsCalculated >= bits) break;
\r
3088 temp.exponent = 0;
\r
3089 bitsCalculated += i;
\r
3092 this.Assign(result);
\r
3096 /// Tried the newton method for logs, but the exponential function is too slow to do it.
\r
3098 private void LogNewton()
\r
3100 if (mantissa.IsZero() || mantissa.Sign)
\r
3106 if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
\r
3108 CalculateLog2(mantissa.Precision.NumBits);
\r
3111 int numBits = mantissa.Precision.NumBits;
\r
3113 //Use inverse exp function with Newton's method.
\r
3114 BigFloat xn = new BigFloat(this);
\r
3115 BigFloat oldExponent = new BigFloat(xn.exponent, mantissa.Precision);
\r
3117 this.exponent = 0;
\r
3118 //Hack to subtract 1
\r
3119 xn.mantissa.ClearBit(numBits - 1);
\r
3120 //x0 = (x - 1) * log2 - this is a straight line fit between log(1) = 0 and log(2) = ln2
\r
3122 //x0 = (x - 1) * log2 + C - this corrects for minimum error over the range.
\r
3123 xn.Add(logNewtonConstant);
\r
3124 BigFloat term = new BigFloat(mantissa.Precision);
\r
3125 BigFloat one = new BigFloat(1, mantissa.Precision);
\r
3127 int precision = 32;
\r
3128 int normalPrecision = mantissa.Precision.NumBits;
\r
3130 int iterations = 0;
\r
3135 term.mantissa.Sign = true;
\r
3136 term.Exp(precision);
\r
3141 if (term.exponent < -((precision >> 1) - 4))
\r
3143 if (precision == normalPrecision)
\r
3145 if (term.exponent < -(precision - 4)) break;
\r
3149 precision = precision << 1;
\r
3150 if (precision > normalPrecision) precision = normalPrecision;
\r
3157 //log(2^n*s) = log(2^n) + log(s) = nlog(2) + log(s)
\r
3158 term.Assign(ln2cache);
\r
3159 term.Mul(oldExponent);
\r
3166 /// Log(x) implemented as an Arithmetic-Geometric Mean. Fast for high precisions.
\r
3168 private void LogAGM1()
\r
3170 if (mantissa.IsZero() || mantissa.Sign)
\r
3176 if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
\r
3178 CalculateLog2(mantissa.Precision.NumBits);
\r
3181 //Compute ln(x) using AGM formula
\r
3183 //1. Re-write the input as 2^n * (0.5 <= x < 1)
\r
3184 int power2 = exponent + 1;
\r
3187 //BigFloat res = new BigFloat(firstAGMcache);
\r
3188 BigFloat a0 = new BigFloat(1, mantissa.Precision);
\r
3189 BigFloat b0 = new BigFloat(pow10cache);
\r
3192 BigFloat r = R(a0, b0);
\r
3194 this.Assign(firstAGMcache);
\r
3197 a0.Assign(ln2cache);
\r
3198 a0.Mul(new BigFloat(power2, mantissa.Precision));
\r
3202 private void Exp(int numBits)
\r
3204 if (IsSpecialValue)
\r
3206 if (SpecialValue == SpecialValueType.ZERO)
\r
3210 mantissa.SetHighDigit(0x80000000);
\r
3212 else if (SpecialValue == SpecialValueType.INF_MINUS)
\r
3221 PrecisionSpec prec = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
\r
3222 numBits = prec.NumBits;
\r
3224 if (scratch.Precision.NumBits != prec.NumBits)
\r
3226 scratch = new BigInt(prec);
\r
3229 if (inverseFactorialCache == null || invFactorialCutoff < numBits)
\r
3231 CalculateFactorials(numBits);
\r
3234 //let x = 1 * 'this'.mantissa (i.e. 1 <= x < 2)
\r
3235 //exp(2^n * x) = e^(2^n * x) = (e^x)^2n = exp(x)^2n
\r
3237 int oldExponent = 0;
\r
3239 if (exponent > -4)
\r
3241 oldExponent = exponent + 4;
\r
3245 BigFloat thisSave = new BigFloat(this, prec);
\r
3246 BigFloat temp = new BigFloat(1, prec);
\r
3247 BigFloat temp2 = new BigFloat(this, prec);
\r
3248 BigFloat res = new BigFloat(1, prec);
\r
3249 int length = inverseFactorialCache.Length;
\r
3252 for (int i = 1; i < length; i++)
\r
3255 temp.Mul(thisSave);
\r
3256 temp2.Assign(inverseFactorialCache[i]);
\r
3259 if (temp2.exponent < -(numBits + 4)) { iterations = i; break; }
\r
3265 //Now... x^(2^n) = (x^2)^(2^(n - 1))
\r
3266 for (int i = 0; i < oldExponent; i++)
\r
3268 res.mantissa.SquareHiFast(scratch);
\r
3269 int shift = res.mantissa.Normalise();
\r
3270 res.exponent = res.exponent << 1;
\r
3271 res.exponent += 1 - shift;
\r
3274 //Deal with +/- inf
\r
3275 if (res.exponent == Int32.MaxValue)
\r
3277 res.mantissa.Zero();
\r
3284 /// Calculates ln(2) and returns -10^(n/2 + a bit) for reuse, using the AGM method as described in
\r
3285 /// http://lacim.uqam.ca/~plouffe/articles/log2.pdf
\r
3287 /// <param name="numBits"></param>
\r
3288 /// <returns></returns>
\r
3289 private static void CalculateLog2(int numBits)
\r
3291 //Use the AGM method formula to get log2 to N digits.
\r
3292 //R(a0, b0) = 1 / (1 - Sum(2^-n*(an^2 - bn^2)))
\r
3293 //log(1/2) = R(1, 10^-n) - R(1, 10^-n/2)
\r
3294 PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
\r
3295 PrecisionSpec extendedPres = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);
\r
3296 BigFloat a0 = new BigFloat(1, extendedPres);
\r
3297 BigFloat b0 = TenPow(-(int)((double)((numBits >> 1) + 2) * 0.302), extendedPres);
\r
3298 BigFloat pow10saved = new BigFloat(b0);
\r
3299 BigFloat firstAGMcacheSaved = new BigFloat(extendedPres);
\r
3301 //save power of 10 (in normal precision)
\r
3302 pow10cache = new BigFloat(b0, normalPres);
\r
3304 ln2cache = R(a0, b0);
\r
3306 //save the first half of the log calculation
\r
3307 firstAGMcache = new BigFloat(ln2cache, normalPres);
\r
3308 firstAGMcacheSaved.Assign(ln2cache);
\r
3311 ln2cache.Sub(R(a0, b0));
\r
3313 //Convert to log(2)
\r
3314 ln2cache.mantissa.Sign = false;
\r
3316 //Save magic constant for newton log
\r
3317 //First guess in range 1 <= x < 2 is x0 = ln2 * (x - 1) + C
\r
3318 logNewtonConstant = new BigFloat(ln2cache);
\r
3319 logNewtonConstant.Mul(new BigFloat(3, extendedPres));
\r
3320 logNewtonConstant.exponent--;
\r
3321 logNewtonConstant.Sub(new BigFloat(1, extendedPres));
\r
3322 logNewtonConstant = new BigFloat(logNewtonConstant, normalPres);
\r
3324 //Save the inverse.
\r
3325 log2ecache = new BigFloat(ln2cache);
\r
3326 log2ecache = new BigFloat(log2ecache.Reciprocal(), normalPres);
\r
3329 //Because the log functions call this function to the precision to which they
\r
3330 //are called, we cannot call them without causing an infinite loop, so we need
\r
3331 //to inline the code.
\r
3332 log10recip = new BigFloat(10, extendedPres);
\r
3335 int power2 = log10recip.exponent + 1;
\r
3336 log10recip.exponent = -1;
\r
3338 //BigFloat res = new BigFloat(firstAGMcache);
\r
3339 BigFloat ax = new BigFloat(1, extendedPres);
\r
3340 BigFloat bx = new BigFloat(pow10saved);
\r
3341 bx.Mul(log10recip);
\r
3343 BigFloat r = R(ax, bx);
\r
3345 log10recip.Assign(firstAGMcacheSaved);
\r
3346 log10recip.Sub(r);
\r
3348 ax.Assign(ln2cache);
\r
3349 ax.Mul(new BigFloat(power2, log10recip.mantissa.Precision));
\r
3350 log10recip.Add(ax);
\r
3353 log10recip = log10recip.Reciprocal();
\r
3354 log10recip = new BigFloat(log10recip, normalPres);
\r
3358 ln2cache = new BigFloat(ln2cache, normalPres);
\r
3361 private static BigFloat TenPow(int power, PrecisionSpec precision)
\r
3363 BigFloat acc = new BigFloat(1, precision);
\r
3364 BigFloat temp = new BigFloat(1, precision);
\r
3366 int powerTemp = power;
\r
3368 BigFloat multiplierToUse = new BigFloat(10, precision);
\r
3372 multiplierToUse = multiplierToUse.Reciprocal();
\r
3373 powerTemp = -power;
\r
3376 //Fast power function
\r
3377 while (powerTemp != 0)
\r
3379 temp.Mul(multiplierToUse);
\r
3380 multiplierToUse.Assign(temp);
\r
3382 if ((powerTemp & 1) != 0)
\r
3393 private static BigFloat R(BigFloat a0, BigFloat b0)
\r
3395 //Precision extend taken out.
\r
3396 int bits = a0.mantissa.Precision.NumBits;
\r
3397 PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);
\r
3398 BigFloat an = new BigFloat(a0, extendedPres);
\r
3399 BigFloat bn = new BigFloat(b0, extendedPres);
\r
3400 BigFloat sum = new BigFloat(extendedPres);
\r
3401 BigFloat term = new BigFloat(extendedPres);
\r
3402 BigFloat temp1 = new BigFloat(extendedPres);
\r
3403 BigFloat one = new BigFloat(1, extendedPres);
\r
3405 int iteration = 0;
\r
3407 for (iteration = 0; ; iteration++)
\r
3409 //Get the sum term for this iteration.
\r
3414 //term = an^2 - bn^2
\r
3416 //term = 2^(n-1) * (an^2 - bn^2)
\r
3417 term.exponent += iteration - 1;
\r
3420 if (term.exponent < -(bits - 8)) break;
\r
3422 //Calculate the new AGM estimates.
\r
3425 //a(n+1) = (an + bn) / 2
\r
3428 //b(n+1) = sqrt(an*bn)
\r
3434 one = one.Reciprocal();
\r
3435 return new BigFloat(one, a0.mantissa.Precision);
\r
3438 private static void CalculateFactorials(int numBits)
\r
3440 System.Collections.Generic.List<BigFloat> list = new System.Collections.Generic.List<BigFloat>(64);
\r
3441 System.Collections.Generic.List<BigFloat> list2 = new System.Collections.Generic.List<BigFloat>(64);
\r
3443 PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);
\r
3444 PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
\r
3446 BigFloat factorial = new BigFloat(1, extendedPrecision);
\r
3447 BigFloat reciprocal;
\r
3449 //Calculate e while we're at it
\r
3450 BigFloat e = new BigFloat(1, extendedPrecision);
\r
3452 list.Add(new BigFloat(factorial, normalPrecision));
\r
3454 for (int i = 1; i < Int32.MaxValue; i++)
\r
3456 BigFloat number = new BigFloat(i, extendedPrecision);
\r
3457 factorial.Mul(number);
\r
3459 if (factorial.exponent > numBits) break;
\r
3461 list2.Add(new BigFloat(factorial, normalPrecision));
\r
3462 reciprocal = factorial.Reciprocal();
\r
3464 e.Add(reciprocal);
\r
3465 list.Add(new BigFloat(reciprocal, normalPrecision));
\r
3468 //Set the cached static values.
\r
3469 inverseFactorialCache = list.ToArray();
\r
3470 factorialCache = list2.ToArray();
\r
3471 invFactorialCutoff = numBits;
\r
3472 eCache = new BigFloat(e, normalPrecision);
\r
3473 eRCPCache = new BigFloat(e.Reciprocal(), normalPrecision);
\r
3476 private static void CalculateEOnly(int numBits)
\r
3478 PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);
\r
3479 PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
\r
3481 int iExponent = (int)(Math.Sqrt(numBits));
\r
3483 BigFloat factorial = new BigFloat(1, extendedPrecision);
\r
3484 BigFloat constant = new BigFloat(1, extendedPrecision);
\r
3485 constant.exponent -= iExponent;
\r
3486 BigFloat numerator = new BigFloat(constant);
\r
3487 BigFloat reciprocal;
\r
3489 //Calculate the 2^iExponent th root of e
\r
3490 BigFloat e = new BigFloat(1, extendedPrecision);
\r
3493 for (i = 1; i < Int32.MaxValue; i++)
\r
3495 BigFloat number = new BigFloat(i, extendedPrecision);
\r
3496 factorial.Mul(number);
\r
3497 reciprocal = factorial.Reciprocal();
\r
3498 reciprocal.Mul(numerator);
\r
3500 if (-reciprocal.exponent > numBits) break;
\r
3502 e.Add(reciprocal);
\r
3503 numerator.Mul(constant);
\r
3504 System.GC.Collect();
\r
3507 for (i = 0; i < iExponent; i++)
\r
3509 numerator.Assign(e);
\r
3513 //Set the cached static values.
\r
3514 eCache = new BigFloat(e, normalPrecision);
\r
3515 eRCPCache = new BigFloat(e.Reciprocal(), normalPrecision);
\r
3519 /// Uses the Gauss-Legendre formula for pi
\r
3520 /// Taken from http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm
\r
3522 /// <param name="numBits"></param>
\r
3523 private static void CalculatePi(int numBits)
\r
3525 int bits = numBits + 32;
\r
3526 //Precision extend taken out.
\r
3527 PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
\r
3528 PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);
\r
3530 if (scratch.Precision.NumBits != bits)
\r
3532 scratch = new BigInt(extendedPres);
\r
3536 BigFloat an = new BigFloat(1, extendedPres);
\r
3539 BigFloat bn = new BigFloat(2, extendedPres);
\r
3544 BigFloat tn = new BigFloat(1, extendedPres);
\r
3549 BigFloat anTemp = new BigFloat(extendedPres);
\r
3551 int iteration = 0;
\r
3552 int cutoffBits = numBits >> 5;
\r
3554 for (iteration = 0; ; iteration++)
\r
3557 anTemp.Assign(an);
\r
3559 //Calculate new an
\r
3563 //Calculate new bn
\r
3567 //Calculate new tn
\r
3569 anTemp.mantissa.SquareHiFast(scratch);
\r
3570 anTemp.exponent += anTemp.exponent + pn + 1 - anTemp.mantissa.Normalise();
\r
3573 anTemp.Assign(an);
\r
3576 if (anTemp.exponent < -(bits - cutoffBits)) break;
\r
3583 an.mantissa.SquareHiFast(scratch);
\r
3584 an.exponent += an.exponent + 1 - an.mantissa.Normalise();
\r
3588 pi = new BigFloat(an, normalPres);
\r
3589 piBy2 = new BigFloat(pi);
\r
3591 twoPi = new BigFloat(pi, normalPres);
\r
3593 piRecip = new BigFloat(an.Reciprocal(), normalPres);
\r
3594 twoPiRecip = new BigFloat(piRecip);
\r
3595 twoPiRecip.exponent--;
\r
3596 //1/3 is going to be useful for sin.
\r
3597 threeRecip = new BigFloat((new BigFloat(3, extendedPres)).Reciprocal(), normalPres);
\r
3601 /// Calculates the odd reciprocals of the natural numbers (for atan series)
\r
3603 /// <param name="numBits"></param>
\r
3604 /// <param name="terms"></param>
\r
3605 private static void CalculateReciprocals(int numBits, int terms)
\r
3607 int bits = numBits + 32;
\r
3608 PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);
\r
3609 PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
\r
3611 System.Collections.Generic.List<BigFloat> list = new System.Collections.Generic.List<BigFloat>(terms);
\r
3613 for (int i = 0; i < terms; i++)
\r
3615 BigFloat term = new BigFloat(i*2 + 1, extendedPres);
\r
3616 list.Add(new BigFloat(term.Reciprocal(), normalPres));
\r
3619 reciprocals = list.ToArray();
\r
3623 /// Does decimal rounding, for numbers without E notation.
\r
3625 /// <param name="input"></param>
\r
3626 /// <param name="places"></param>
\r
3627 /// <returns></returns>
\r
3628 private static string RoundString(string input, int places)
\r
3630 if (places <= 0) return input;
\r
3631 string trim = input.Trim();
\r
3632 char[] digits = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'};
\r
3635 for (int i = 1; i <= places; i++)
\r
3637 //Skip decimal points.
\r
3638 if (trim[trim.Length - i] == '.')
\r
3644 int index = Array.IndexOf(digits, trim[trim.Length - i]);
\r
3646 if (index < 0) return input;
\r
3648 value += ten * index;
\r
3653 //Look for a decimal point
\r
3654 string decimalPoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;
\r
3656 int indexPoint = trim.LastIndexOf(decimalPoint);
\r
3657 if (indexPoint < 0)
\r
3659 //We can't modify a string which doesn't have a decimal point.
\r
3663 int trimPoint = trim.Length - places;
\r
3664 if (trimPoint < indexPoint) trimPoint = indexPoint;
\r
3666 bool roundDown = false;
\r
3668 if (trim[trimPoint] == '.')
\r
3670 if (trimPoint + 1 >= trim.Length)
\r
3676 int digit = Array.IndexOf(digits, trim[trimPoint + 1]);
\r
3677 if (digit < 5) roundDown = true;
\r
3682 int digit = Array.IndexOf(digits, trim[trimPoint]);
\r
3683 if (digit < 5) roundDown = true;
\r
3688 //Round down - just return a new string without the extra digits.
\r
3691 if (RoundingMode == RoundingModeType.EXACT)
\r
3693 return trim.Substring(0, trimPoint);
\r
3697 char[] trimChars = { '0' };
\r
3698 output = trim.Substring(0, trimPoint).TrimEnd(trimChars);
\r
3699 trimChars[0] = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator[0];
\r
3700 return output.TrimEnd(trimChars);
\r
3704 //Round up - bit more complicated.
\r
3705 char [] arrayOutput = trim.ToCharArray();//0, trimPoint);
\r
3707 //Now, we round going from the back to the front.
\r
3709 for (j = trimPoint - 1; j >= 0; j--)
\r
3711 int index = Array.IndexOf(digits, arrayOutput[j]);
\r
3713 //Skip decimal points etc...
\r
3714 if (index < 0) continue;
\r
3718 arrayOutput[j] = digits[index + 1];
\r
3723 arrayOutput[j] = digits[0];
\r
3727 output = new string(arrayOutput);
\r
3731 //Need to add a new digit.
\r
3732 output = String.Format("{0}{1}", "1", output);
\r
3735 if (RoundingMode == RoundingModeType.EXACT)
\r
3741 char[] trimChars = { '0' };
\r
3742 output = output.TrimEnd(trimChars);
\r
3743 trimChars[0] = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator[0];
\r
3744 return output.TrimEnd(trimChars);
\r
3748 //***************************** Data *****************************
\r
3751 //Side node - this way of doing things is far from optimal, both in terms of memory use and performance.
\r
3752 private ExponentAdaptor exponent;
\r
3753 private BigInt mantissa;
\r
3756 /// Storage area for calculations.
\r
3758 private static BigInt scratch;
\r
3760 private static BigFloat ln2cache; //Value of ln(2)
\r
3761 private static BigFloat log2ecache; //Value of log2(e) = 1/ln(2)
\r
3762 private static BigFloat pow10cache; //Cached power of 10 for AGM log calculation
\r
3763 private static BigFloat log10recip; //1/ln(10)
\r
3764 private static BigFloat firstAGMcache; //Cached half of AGM operation.
\r
3765 private static BigFloat[] factorialCache; //The values of n!
\r
3766 private static BigFloat[] inverseFactorialCache; //Values of 1/n! up to 2^-m where m = invFactorialCutoff (below)
\r
3767 private static int invFactorialCutoff; //The number of significant bits for the cutoff of the inverse factorials.
\r
3768 private static BigFloat eCache; //Value of e cached to invFactorialCutoff bits
\r
3769 private static BigFloat eRCPCache; //Reciprocal of e
\r
3770 private static BigFloat logNewtonConstant; //1.5*ln(2) - 1
\r
3771 private static BigFloat pi; //pi
\r
3772 private static BigFloat piBy2; //pi/2
\r
3773 private static BigFloat twoPi; //2*pi
\r
3774 private static BigFloat piRecip; //1/pi
\r
3775 private static BigFloat twoPiRecip; //1/2*pi
\r
3776 private static BigFloat threeRecip; //1/3
\r
3777 private static BigFloat[] reciprocals; //1/x
\r
3780 /// The number of decimal digits to round the output of ToString() by
\r
3782 public static int RoundingDigits { get; set; }
\r
3785 /// The way in which ToString() should deal with insignificant trailing zeroes
\r
3787 public static RoundingModeType RoundingMode { get; set; }
\r