dev4/psychlops/extention/math/BigFloat.cs

   1 // http://www.fractal-landscapes.co.uk/bigint.html\r
   2 \r
   3 using System;\r
   4 \r
   5 namespace BigNum\r
   6 {\r
   7     /// <summary>\r
   8     /// An arbitrary-precision floating-point class\r
   9     /// \r
  10     /// Format:\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
  13     /// \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
  18     /// \r
  19     /// Notes:\r
  20     /// All conversions to and from strings are slow.\r
  21     /// \r
  22     /// Conversions from simple integer types Int32, Int64, UInt32, UInt64 are performed\r
  23     /// using the appropriate constructor, and are relatively fast.\r
  24     /// \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
  28     /// </summary>\r
  29     public class BigFloat\r
  30     {\r
  31         /// <summary>\r
  32         /// Floats can have 4 special value types:\r
  33         /// \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
  37         /// Zero\r
  38         /// </summary>\r
  39         public enum SpecialValueType\r
  40         {\r
  41             /// <summary>\r
  42             /// Not a special value\r
  43             /// </summary>\r
  44             NONE = 0,\r
  45             /// <summary>\r
  46             /// Zero\r
  47             /// </summary>\r
  48             ZERO,\r
  49             /// <summary>\r
  50             /// Positive infinity\r
  51             /// </summary>\r
  52             INF_PLUS,\r
  53             /// <summary>\r
  54             /// Negative infinity\r
  55             /// </summary>\r
  56             INF_MINUS,\r
  57             /// <summary>\r
  58             /// Not a number\r
  59             /// </summary>\r
  60             NAN\r
  61         }\r
  62 \r
  63         /// <summary>\r
  64         /// This affects the ToString() method. \r
  65         /// \r
  66         /// With Trim rounding, all insignificant zero digits are drip\r
  67         /// </summary>\r
  68         public enum RoundingModeType\r
  69         {\r
  70             /// <summary>\r
  71             /// Trim non-significant zeros from ToString output after rounding\r
  72             /// </summary>\r
  73             TRIM,\r
  74             /// <summary>\r
  75             /// Keep all non-significant zeroes in ToString output after rounding\r
  76             /// </summary>\r
  77             EXACT\r
  78         }\r
  79 \r
  80         /// <summary>\r
  81         /// A wrapper for the signed exponent, avoiding overflow.\r
  82         /// </summary>\r
  83         protected struct ExponentAdaptor\r
  84         {\r
  85             /// <summary>\r
  86             /// The 32-bit exponent\r
  87             /// </summary>\r
  88             public Int32 exponent\r
  89             {\r
  90                 get { return expValue; }\r
  91                 set { expValue = value; }\r
  92             }\r
  93 \r
  94             /// <summary>\r
  95             /// Implicit cast to Int32\r
  96             /// </summary>\r
  97             public static implicit operator Int32(ExponentAdaptor adaptor)\r
  98             {\r
  99                 return adaptor.expValue;\r
 100             }\r
 101 \r
 102             /// <summary>\r
 103             /// Implicit cast from Int32 to ExponentAdaptor\r
 104             /// </summary>\r
 105             /// <param name="value"></param>\r
 106             /// <returns></returns>\r
 107             public static implicit operator ExponentAdaptor(Int32 value)\r
 108             {\r
 109                 ExponentAdaptor adaptor = new ExponentAdaptor();\r
 110                 adaptor.expValue = value;\r
 111                 return adaptor;\r
 112             }\r
 113 \r
 114             /// <summary>\r
 115             /// Overloaded increment operator\r
 116             /// </summary>\r
 117             public static ExponentAdaptor operator ++(ExponentAdaptor adaptor)\r
 118             {\r
 119                 adaptor = adaptor + 1;\r
 120                 return adaptor;\r
 121             }\r
 122 \r
 123             /// <summary>\r
 124             /// Overloaded decrement operator\r
 125             /// </summary>\r
 126             public static ExponentAdaptor operator --(ExponentAdaptor adaptor)\r
 127             {\r
 128                 adaptor = adaptor - 1;\r
 129                 return adaptor;\r
 130             }\r
 131 \r
 132             /// <summary>\r
 133             /// Overloaded addition operator\r
 134             /// </summary>\r
 135             public static ExponentAdaptor operator +(ExponentAdaptor a1, ExponentAdaptor a2)\r
 136             {\r
 137                 if (a1.expValue == Int32.MaxValue) return a1;\r
 138 \r
 139                 Int64 temp = (Int64)a1.expValue;\r
 140                 temp += (Int64)(a2.expValue);\r
 141 \r
 142                 if (temp > (Int64)Int32.MaxValue)\r
 143                 {\r
 144                     a1.expValue = Int32.MaxValue;\r
 145                 }\r
 146                 else if (temp < (Int64)Int32.MinValue)\r
 147                 {\r
 148                     a1.expValue = Int32.MinValue;\r
 149                 }\r
 150                 else\r
 151                 {\r
 152                     a1.expValue = (Int32)temp;\r
 153                 }\r
 154 \r
 155                 return a1;\r
 156             }\r
 157 \r
 158             /// <summary>\r
 159             /// Overloaded subtraction operator\r
 160             /// </summary>\r
 161             public static ExponentAdaptor operator -(ExponentAdaptor a1, ExponentAdaptor a2)\r
 162             {\r
 163                 if (a1.expValue == Int32.MaxValue) return a1;\r
 164 \r
 165                 Int64 temp = (Int64)a1.expValue;\r
 166                 temp -= (Int64)(a2.expValue);\r
 167 \r
 168                 if (temp > (Int64)Int32.MaxValue)\r
 169                 {\r
 170                     a1.expValue = Int32.MaxValue;\r
 171                 }\r
 172                 else if (temp < (Int64)Int32.MinValue)\r
 173                 {\r
 174                     a1.expValue = Int32.MinValue;\r
 175                 }\r
 176                 else\r
 177                 {\r
 178                     a1.expValue = (Int32)temp;\r
 179                 }\r
 180 \r
 181                 return a1;\r
 182             }\r
 183 \r
 184             /// <summary>\r
 185             /// Overloaded multiplication operator\r
 186             /// </summary>\r
 187             public static ExponentAdaptor operator *(ExponentAdaptor a1, ExponentAdaptor a2)\r
 188             {\r
 189                 if (a1.expValue == Int32.MaxValue) return a1;\r
 190 \r
 191                 Int64 temp = (Int64)a1.expValue;\r
 192                 temp *= (Int64)a2.expValue;\r
 193 \r
 194                 if (temp > (Int64)Int32.MaxValue)\r
 195                 {\r
 196                     a1.expValue = Int32.MaxValue;\r
 197                 }\r
 198                 else if (temp < (Int64)Int32.MinValue)\r
 199                 {\r
 200                     a1.expValue = Int32.MinValue;\r
 201                 }\r
 202                 else\r
 203                 {\r
 204                     a1.expValue = (Int32)temp;\r
 205                 }\r
 206 \r
 207                 return a1;\r
 208             }\r
 209 \r
 210             /// <summary>\r
 211             /// Overloaded division operator\r
 212             /// </summary>\r
 213             public static ExponentAdaptor operator /(ExponentAdaptor a1, ExponentAdaptor a2)\r
 214             {\r
 215                 if (a1.expValue == Int32.MaxValue) return a1;\r
 216 \r
 217                 ExponentAdaptor res = new ExponentAdaptor();\r
 218                 res.expValue = a1.expValue / a2.expValue;\r
 219                 return res;\r
 220             }\r
 221 \r
 222             /// <summary>\r
 223             /// Overloaded right-shift operator\r
 224             /// </summary>\r
 225             public static ExponentAdaptor operator >>(ExponentAdaptor a1, int shift)\r
 226             {\r
 227                 if (a1.expValue == Int32.MaxValue) return a1;\r
 228 \r
 229                 ExponentAdaptor res = new ExponentAdaptor();\r
 230                 res.expValue = a1.expValue >> shift;\r
 231                 return res;\r
 232             }\r
 233 \r
 234             /// <summary>\r
 235             /// Overloaded left-shift operator\r
 236             /// </summary>\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
 241             {\r
 242                 if (a1.expValue == 0) return a1;\r
 243 \r
 244                 ExponentAdaptor res = new ExponentAdaptor();\r
 245                 res.expValue = a1.expValue;\r
 246 \r
 247                 if (shift > 31)\r
 248                 {\r
 249                     res.expValue = Int32.MaxValue;\r
 250                 }\r
 251                 else\r
 252                 {\r
 253                     Int64 temp = a1.expValue;\r
 254                     temp = temp << shift;\r
 255 \r
 256                     if (temp > (Int64)Int32.MaxValue)\r
 257                     {\r
 258                         res.expValue = Int32.MaxValue;\r
 259                     }\r
 260                     else if (temp < (Int64)Int32.MinValue)\r
 261                     {\r
 262                         res.expValue = Int32.MinValue;\r
 263                     }\r
 264                     else\r
 265                     {\r
 266                         res.expValue = (Int32)temp;\r
 267                     }\r
 268                 }\r
 269 \r
 270                 return res;\r
 271             }\r
 272 \r
 273             private Int32 expValue;\r
 274         }\r
 275 \r
 276         //************************ Constructors **************************\r
 277 \r
 278         /// <summary>\r
 279         /// Constructs a 128-bit BigFloat\r
 280         /// \r
 281         /// Sets the value to zero\r
 282         /// </summary>\r
 283         static BigFloat()\r
 284         {\r
 285             RoundingDigits = 3;\r
 286             RoundingMode = RoundingModeType.TRIM;\r
 287             scratch = new BigInt(new PrecisionSpec(128, PrecisionSpec.BaseType.BIN));\r
 288         }\r
 289 \r
 290         /// <summary>\r
 291         /// Constructs a BigFloat of the required precision\r
 292         /// \r
 293         /// Sets the value to zero\r
 294         /// </summary>\r
 295         /// <param name="mantissaPrec"></param>\r
 296         public BigFloat(PrecisionSpec mantissaPrec)\r
 297         {\r
 298             Init(mantissaPrec);\r
 299         }\r
 300 \r
 301         /// <summary>\r
 302         /// Constructs a big float from a UInt32 to the required precision\r
 303         /// </summary>\r
 304         /// <param name="value"></param>\r
 305         /// <param name="mantissaPrec"></param>\r
 306         public BigFloat(UInt32 value, PrecisionSpec mantissaPrec)\r
 307         {\r
 308             int mbWords = ((mantissaPrec.NumBits) >> 5);\r
 309             if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
 310             int newManBits = mbWords << 5;\r
 311 \r
 312             //For efficiency, we just use a 32-bit exponent\r
 313             exponent = 0;\r
 314 \r
 315             mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 316             //scratch = new BigInt(mantissa.Precision);\r
 317 \r
 318             int bit = BigInt.GetMSB(value);\r
 319             if (bit == -1) return;\r
 320 \r
 321             int shift = mantissa.Precision.NumBits - (bit + 1);\r
 322             mantissa.LSH(shift);\r
 323             exponent = bit;\r
 324         }\r
 325 \r
 326         /// <summary>\r
 327         /// Constructs a BigFloat from an Int32 to the required precision\r
 328         /// </summary>\r
 329         /// <param name="value"></param>\r
 330         /// <param name="mantissaPrec"></param>\r
 331         public BigFloat(Int32 value, PrecisionSpec mantissaPrec)\r
 332         {\r
 333             int mbWords = ((mantissaPrec.NumBits) >> 5);\r
 334             if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
 335             int newManBits = mbWords << 5;\r
 336 \r
 337             //For efficiency, we just use a 32-bit exponent\r
 338             exponent = 0;\r
 339             UInt32 uValue;\r
 340             \r
 341             if (value < 0)\r
 342             {\r
 343                 if (value == Int32.MinValue)\r
 344                 {\r
 345                     uValue = 0x80000000;\r
 346                 }\r
 347                 else\r
 348                 {\r
 349                     uValue = (UInt32)(-value);\r
 350                 }\r
 351             }\r
 352             else\r
 353             {\r
 354                 uValue = (UInt32)value;\r
 355             }\r
 356 \r
 357             mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 358             //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 359 \r
 360             int bit = BigInt.GetMSB(uValue);\r
 361             if (bit == -1) return;\r
 362 \r
 363             int shift = mantissa.Precision.NumBits - (bit + 1);\r
 364             mantissa.LSH(shift);\r
 365             exponent = bit;\r
 366         }\r
 367 \r
 368         /// <summary>\r
 369         /// Constructs a BigFloat from a 64-bit integer\r
 370         /// </summary>\r
 371         /// <param name="value"></param>\r
 372         /// <param name="mantissaPrec"></param>\r
 373         public BigFloat(Int64 value, PrecisionSpec mantissaPrec)\r
 374         {\r
 375             int mbWords = ((mantissaPrec.NumBits) >> 5);\r
 376             if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
 377             int newManBits = mbWords << 5;\r
 378 \r
 379             //For efficiency, we just use a 32-bit exponent\r
 380             exponent = 0;\r
 381             UInt64 uValue;\r
 382 \r
 383             if (value < 0)\r
 384             {\r
 385                 if (value == Int64.MinValue)\r
 386                 {\r
 387                     uValue = 0x80000000;\r
 388                 }\r
 389                 else\r
 390                 {\r
 391                     uValue = (UInt64)(-value);\r
 392                 }\r
 393             }\r
 394             else\r
 395             {\r
 396                 uValue = (UInt64)value;\r
 397             }\r
 398 \r
 399             mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 400             //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 401 \r
 402             int bit = BigInt.GetMSB(uValue);\r
 403             if (bit == -1) return;\r
 404 \r
 405             int shift = mantissa.Precision.NumBits - (bit + 1);\r
 406             if (shift > 0)\r
 407             {\r
 408                 mantissa.LSH(shift);\r
 409             }\r
 410             else\r
 411             {\r
 412                 mantissa.SetHighDigit((uint)(uValue >> (-shift)));\r
 413             }\r
 414             exponent = bit;\r
 415         }\r
 416 \r
 417         /// <summary>\r
 418         /// Constructs a BigFloat from a 64-bit unsigned integer\r
 419         /// </summary>\r
 420         /// <param name="value"></param>\r
 421         /// <param name="mantissaPrec"></param>\r
 422         public BigFloat(UInt64 value, PrecisionSpec mantissaPrec)\r
 423         {\r
 424             int mbWords = ((mantissaPrec.NumBits) >> 5);\r
 425             if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
 426             int newManBits = mbWords << 5;\r
 427 \r
 428             //For efficiency, we just use a 32-bit exponent\r
 429             exponent = 0;\r
 430 \r
 431             int bit = BigInt.GetMSB(value);\r
 432 \r
 433             mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 434             //scratch = new BigInt(mantissa.Precision);\r
 435 \r
 436             int shift = mantissa.Precision.NumBits - (bit + 1);\r
 437             if (shift > 0)\r
 438             {\r
 439                 mantissa.LSH(shift);\r
 440             }\r
 441             else\r
 442             {\r
 443                 mantissa.SetHighDigit((uint)(value >> (-shift)));\r
 444             }\r
 445             exponent = bit;\r
 446         }\r
 447 \r
 448         /// <summary>\r
 449         /// Constructs a BigFloat from a BigInt, using the specified precision\r
 450         /// </summary>\r
 451         /// <param name="value"></param>\r
 452         /// <param name="mantissaPrec"></param>\r
 453         public BigFloat(BigInt value, PrecisionSpec mantissaPrec)\r
 454         {\r
 455             if (value.IsZero())\r
 456             {\r
 457                 Init(mantissaPrec);\r
 458                 SetZero();\r
 459                 return;\r
 460             }\r
 461 \r
 462             mantissa = new BigInt(value, mantissaPrec);\r
 463             exponent = BigInt.GetMSB(value);\r
 464             mantissa.Normalise();\r
 465         }\r
 466 \r
 467         /// <summary>\r
 468         /// Construct a BigFloat from a double-precision floating point number\r
 469         /// </summary>\r
 470         /// <param name="value"></param>\r
 471         /// <param name="mantissaPrec"></param>\r
 472         public BigFloat(double value, PrecisionSpec mantissaPrec)\r
 473         {\r
 474             if (value == 0.0)\r
 475             {\r
 476                 Init(mantissaPrec);\r
 477                 return;\r
 478             }\r
 479 \r
 480             bool sign = (value < 0) ? true : false;\r
 481 \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
 486 \r
 487             //The mantissa is stored with the top bit implied.\r
 488             valueMantissa = valueMantissa | 0x10000000000000L;\r
 489 \r
 490             //The exponent is biased by 1023.\r
 491             exponent = valueExponent - 1023;\r
 492 \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
 497 \r
 498             mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 499             //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 500 \r
 501             if (newManBits >= 64)\r
 502             {\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
 505             }\r
 506             else\r
 507             {\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
 510             }\r
 511 \r
 512             mantissa.Sign = sign;\r
 513         }\r
 514 \r
 515         /// <summary>\r
 516         /// Copy constructor\r
 517         /// </summary>\r
 518         /// <param name="value"></param>\r
 519         public BigFloat(BigFloat value)\r
 520         {\r
 521             Init(value.mantissa.Precision);\r
 522             exponent = value.exponent;\r
 523             mantissa.Assign(value.mantissa);\r
 524         }\r
 525 \r
 526         /// <summary>\r
 527         /// Copy Constructor - constructs a new BigFloat with the specified precision, copying the old one.\r
 528         /// \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
 531         /// </summary>\r
 532         /// <param name="value"></param>\r
 533         /// <param name="mantissaPrec"></param>\r
 534         public BigFloat(BigFloat value, PrecisionSpec mantissaPrec)\r
 535         {\r
 536             Init(mantissaPrec);\r
 537             exponent = value.exponent;\r
 538             if (mantissa.AssignHigh(value.mantissa)) exponent++;\r
 539         }\r
 540 \r
 541         /// <summary>\r
 542         /// Constructs a BigFloat from a string\r
 543         /// </summary>\r
 544         /// <param name="value"></param>\r
 545         /// <param name="mantissaPrec"></param>\r
 546         public BigFloat(string value, PrecisionSpec mantissaPrec)\r
 547         {\r
 548             Init(mantissaPrec);\r
 549 \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
 555 \r
 556             if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol))\r
 557             {\r
 558                 SetNaN();\r
 559                 return;\r
 560             }\r
 561             else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol))\r
 562             {\r
 563                 SetInfPlus();\r
 564                 return;\r
 565             }\r
 566             else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol))\r
 567             {\r
 568                 SetInfMinus();\r
 569                 return;\r
 570             }\r
 571 \r
 572             string decimalpoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;\r
 573 \r
 574             char[] digitChars = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ',', '.' };\r
 575 \r
 576             //Read in the integer part up the the decimal point.\r
 577             bool sign = false;\r
 578             value = value.Trim();\r
 579 \r
 580             int i = 0;\r
 581 \r
 582             if (value.Length > i && value[i] == '-')\r
 583             {\r
 584                 sign = true;\r
 585                 i++;\r
 586             }\r
 587 \r
 588             if (value.Length > i && value[i] == '+')\r
 589             {\r
 590                 i++;\r
 591             }\r
 592 \r
 593             for ( ; i < value.Length; i++)\r
 594             {\r
 595                 //break on decimal point\r
 596                 if (value[i] == decimalpoint[0]) break;\r
 597 \r
 598                 int digit = Array.IndexOf(digitChars, value[i]);\r
 599                 if (digit < 0) break;\r
 600 \r
 601                 //Ignore place separators (assumed either , or .)\r
 602                 if (digit > 9) continue;\r
 603 \r
 604                 if (i > 0) iPart.Mul(ten);\r
 605                 iPart.Add(new BigFloat(digit, extendedPres));\r
 606             }\r
 607 \r
 608             //If we've run out of characters, assign everything and return\r
 609             if (i == value.Length)\r
 610             {\r
 611                 iPart.mantissa.Sign = sign;\r
 612                 exponent = iPart.exponent;\r
 613                 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
 614                 return;\r
 615             }\r
 616 \r
 617             //Assign the characters after the decimal point to fPart\r
 618             if (value[i] == '.' && i < value.Length - 1)\r
 619             {\r
 620                 BigFloat RecipToUse = new BigFloat(tenRCP);\r
 621 \r
 622                 for (i++; i < value.Length; i++)\r
 623                 {\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
 629                     fPart.Add(temp);\r
 630                 }\r
 631             }\r
 632 \r
 633             //If we're run out of characters, add fPart and iPart and return\r
 634             if (i == value.Length)\r
 635             {\r
 636                 iPart.Add(fPart);\r
 637                 iPart.mantissa.Sign = sign;\r
 638                 exponent = iPart.exponent;\r
 639                 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
 640                 return;\r
 641             }\r
 642 \r
 643             if (value[i] == '+' || value[i] == '-') i++;\r
 644 \r
 645             if (i == value.Length)\r
 646             {\r
 647                 iPart.Add(fPart);\r
 648                 iPart.mantissa.Sign = sign;\r
 649                 exponent = iPart.exponent;\r
 650                 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
 651                 return;\r
 652             }\r
 653 \r
 654             //Look for exponential notation.\r
 655             if ((value[i] == 'e' || value[i] == 'E') && i < value.Length - 1)\r
 656             {\r
 657                 //Convert the exponent to an int.\r
 658                 int exp;\r
 659 \r
 660                 try\r
 661                 {\r
 662                                         exp = System.Convert.ToInt32(new string(value.ToCharArray()));// i + 1, value.Length - (i + 1))));\r
 663                 }\r
 664                 catch (Exception)\r
 665                 {\r
 666                     iPart.Add(fPart);\r
 667                     iPart.mantissa.Sign = sign;\r
 668                     exponent = iPart.exponent;\r
 669                     if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
 670                     return;\r
 671                 }\r
 672 \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
 676 \r
 677                 int powerTemp = exp;\r
 678 \r
 679                 BigFloat multiplierToUse;\r
 680 \r
 681                 if (exp < 0)\r
 682                 {\r
 683                     multiplierToUse = new BigFloat(tenRCP);\r
 684                     powerTemp = -exp;\r
 685                 }\r
 686                 else\r
 687                 {\r
 688                     multiplierToUse = new BigFloat(ten);\r
 689                 }\r
 690 \r
 691                 //Fast power function\r
 692                 while (powerTemp != 0)\r
 693                 {\r
 694                     temp.Mul(multiplierToUse);\r
 695                     multiplierToUse.Assign(temp);\r
 696 \r
 697                     if ((powerTemp & 1) != 0)\r
 698                     {\r
 699                         acc.Mul(temp);\r
 700                     }\r
 701 \r
 702                     powerTemp >>= 1;\r
 703                 }\r
 704 \r
 705                 iPart.Add(fPart);\r
 706                 iPart.Mul(acc);\r
 707                 iPart.mantissa.Sign = sign;\r
 708                 exponent = iPart.exponent;\r
 709                 if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
 710 \r
 711                 return;\r
 712             }\r
 713 \r
 714             iPart.Add(fPart);\r
 715             iPart.mantissa.Sign = sign;\r
 716             exponent = iPart.exponent;\r
 717             if (mantissa.AssignHigh(iPart.mantissa)) exponent++;\r
 718 \r
 719         }\r
 720 \r
 721         private void Init(PrecisionSpec mantissaPrec)\r
 722         {\r
 723             int mbWords = ((mantissaPrec.NumBits) >> 5);\r
 724             if ((mantissaPrec.NumBits & 31) != 0) mbWords++;\r
 725             int newManBits = mbWords << 5;\r
 726 \r
 727             //For efficiency, we just use a 32-bit exponent\r
 728             exponent = 0;\r
 729             mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 730             //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));\r
 731         }\r
 732 \r
 733         //************************** Properties *************************\r
 734 \r
 735         /// <summary>\r
 736         /// Read-only property. Returns the precision specification of the mantissa.\r
 737         /// \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
 742         /// </summary>\r
 743         public PrecisionSpec Precision\r
 744         {\r
 745             get { return mantissa.Precision; }\r
 746         }\r
 747 \r
 748         /// <summary>\r
 749         /// Writable property:\r
 750         ///     true iff the number is negative or in some cases zero (&lt;0)\r
 751         ///     false iff the number if positive or in some cases zero (&gt;0)\r
 752         /// </summary>\r
 753         public bool Sign \r
 754         { \r
 755             get { return mantissa.Sign; }\r
 756             set { mantissa.Sign = value; }\r
 757         }\r
 758 \r
 759         /// <summary>\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
 763         /// </summary>\r
 764         public bool IsSpecialValue\r
 765         {\r
 766             get\r
 767             {\r
 768                 return (exponent == Int32.MaxValue || mantissa.IsZero());\r
 769             }\r
 770         }\r
 771 \r
 772         /// <summary>\r
 773         /// Read-only property, returns the type of number this is. Special values include:\r
 774         /// \r
 775         /// NONE - a regular number\r
 776         /// ZERO - zero\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
 780         /// </summary>\r
 781         public SpecialValueType SpecialValue\r
 782         {\r
 783             get\r
 784             {\r
 785                 if (exponent == Int32.MaxValue)\r
 786                 {\r
 787                     if (mantissa.IsZero())\r
 788                     {\r
 789                         if (mantissa.Sign) return SpecialValueType.INF_MINUS;\r
 790                         return SpecialValueType.INF_PLUS;\r
 791                     }\r
 792 \r
 793                     return SpecialValueType.NAN;\r
 794                 }\r
 795                 else\r
 796                 {\r
 797                     if (mantissa.IsZero()) return SpecialValueType.ZERO;\r
 798                     return SpecialValueType.NONE;\r
 799                 }\r
 800             }\r
 801         }\r
 802 \r
 803         //******************** Mathematical Constants *******************\r
 804 \r
 805         /// <summary>\r
 806         /// Gets pi to the indicated precision\r
 807         /// </summary>\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
 811         {\r
 812             if (pi == null || precision.NumBits <= pi.mantissa.Precision.NumBits)\r
 813             {\r
 814                 CalculatePi(precision.NumBits);\r
 815             }\r
 816 \r
 817             BigFloat ret = new BigFloat (precision);\r
 818             ret.Assign(pi);\r
 819 \r
 820             return ret;\r
 821         }\r
 822 \r
 823         /// <summary>\r
 824         /// Get e to the indicated precision\r
 825         /// </summary>\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
 829         {\r
 830             if (eCache == null || eCache.mantissa.Precision.NumBits < precision.NumBits)\r
 831             {\r
 832                 CalculateEOnly(precision.NumBits);\r
 833                 //CalculateFactorials(precision.NumBits);\r
 834             }\r
 835 \r
 836             BigFloat ret = new BigFloat(precision);\r
 837             ret.Assign(eCache);\r
 838 \r
 839             return ret;\r
 840         }\r
 841 \r
 842 \r
 843         //******************** Arithmetic Functions ********************\r
 844 \r
 845         /// <summary>\r
 846         /// Addition (this = this + n2)\r
 847         /// </summary>\r
 848         /// <param name="n2">The number to add</param>\r
 849         public void Add(BigFloat n2)\r
 850         {\r
 851             if (SpecialValueAddTest(n2)) return;\r
 852 \r
 853             if (scratch.Precision.NumBits != n2.mantissa.Precision.NumBits)\r
 854             {\r
 855                 scratch = new BigInt(n2.mantissa.Precision);\r
 856             }\r
 857 \r
 858             if (exponent <= n2.exponent)\r
 859             {\r
 860                 int diff = n2.exponent - exponent;\r
 861                 exponent = n2.exponent;\r
 862 \r
 863                 if (diff != 0)\r
 864                 {\r
 865                     mantissa.RSH(diff);\r
 866                 }\r
 867 \r
 868                 uint carry = mantissa.Add(n2.mantissa);\r
 869 \r
 870                 if (carry != 0)\r
 871                 {\r
 872                     mantissa.RSH(1);\r
 873                     mantissa.SetBit(mantissa.Precision.NumBits - 1);\r
 874                     exponent++;\r
 875                 }\r
 876 \r
 877                 exponent -= mantissa.Normalise();\r
 878             }\r
 879             else\r
 880             {\r
 881                 int diff = exponent - n2.exponent;\r
 882 \r
 883                 scratch.Assign(n2.mantissa);\r
 884                 scratch.RSH(diff);\r
 885 \r
 886                 uint carry = scratch.Add(mantissa);\r
 887 \r
 888                 if (carry != 0)\r
 889                 {\r
 890                     scratch.RSH(1);\r
 891                     scratch.SetBit(mantissa.Precision.NumBits - 1);\r
 892                     exponent++;\r
 893                 }\r
 894 \r
 895                 mantissa.Assign(scratch);\r
 896 \r
 897                 exponent -= mantissa.Normalise();\r
 898             }\r
 899         }\r
 900 \r
 901         /// <summary>\r
 902         /// Subtraction (this = this - n2)\r
 903         /// </summary>\r
 904         /// <param name="n2">The number to subtract from this</param>\r
 905         public void Sub(BigFloat n2)\r
 906         {\r
 907             n2.mantissa.Sign = !n2.mantissa.Sign;\r
 908             Add(n2);\r
 909             n2.mantissa.Sign = !n2.mantissa.Sign;\r
 910         }\r
 911 \r
 912         /// <summary>\r
 913         /// Multiplication (this = this * n2)\r
 914         /// </summary>\r
 915         /// <param name="n2">The number to multiply this by</param>\r
 916         public void Mul(BigFloat n2)\r
 917         {\r
 918             if (SpecialValueMulTest(n2)) return;\r
 919 \r
 920             //Anything times 0 = 0\r
 921             if (n2.mantissa.IsZero())\r
 922             {\r
 923                 mantissa.Assign(n2.mantissa);\r
 924                 exponent = 0;\r
 925                 return;\r
 926             }\r
 927 \r
 928             mantissa.MulHi(n2.mantissa);\r
 929             int shift = mantissa.Normalise();\r
 930             exponent = exponent + n2.exponent + 1 - shift;\r
 931         }\r
 932 \r
 933         /// <summary>\r
 934         /// Division (this = this / n2)\r
 935         /// </summary>\r
 936         /// <param name="n2">The number to divide this by</param>\r
 937         public void Div(BigFloat n2)\r
 938         {\r
 939             if (SpecialValueDivTest(n2)) return;\r
 940 \r
 941             if (mantissa.Precision.NumBits >= 8192)\r
 942             {\r
 943                 BigFloat rcp = n2.Reciprocal();\r
 944                 Mul(rcp);\r
 945             }\r
 946             else\r
 947             {\r
 948                 int shift = mantissa.DivAndShift(n2.mantissa);\r
 949                 exponent = exponent - (n2.exponent + shift);\r
 950             }\r
 951         }\r
 952 \r
 953         /// <summary>\r
 954         /// Multiply by a power of 2 (-ve implies division)\r
 955         /// </summary>\r
 956         /// <param name="pow2"></param>\r
 957         public void MulPow2(int pow2)\r
 958         {\r
 959             exponent += pow2;\r
 960         }\r
 961 \r
 962         /// <summary>\r
 963         /// Division-based reciprocal, fastest for small precisions up to 15,000 bits.\r
 964         /// </summary>\r
 965         /// <returns>The reciprocal 1/this</returns>\r
 966         public BigFloat Reciprocal()\r
 967         {\r
 968             if (mantissa.Precision.NumBits >= 8192) return ReciprocalNewton();\r
 969 \r
 970             BigFloat reciprocal = new BigFloat(1u, mantissa.Precision);\r
 971             reciprocal.Div(this);\r
 972             return reciprocal;\r
 973         }\r
 974 \r
 975         /// <summary>\r
 976         /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.\r
 977         /// </summary>\r
 978         /// <returns>The reciprocal 1/this</returns>\r
 979         public BigFloat ReciprocalNewton()\r
 980         {\r
 981             if (mantissa.IsZero())\r
 982             {\r
 983                 exponent = Int32.MaxValue;\r
 984                 return null;\r
 985             }\r
 986 \r
 987             bool oldSign = mantissa.Sign;\r
 988             int oldExponent = exponent;\r
 989 \r
 990             //Kill exponent for now (will re-institute later)\r
 991             exponent = 0;\r
 992 \r
 993             bool topBit = mantissa.IsTopBitOnlyBit();\r
 994 \r
 995             PrecisionSpec curPrec = new PrecisionSpec(32, PrecisionSpec.BaseType.BIN);\r
 996 \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
1001 \r
1002             reciprocal.exponent = 1;\r
1003             reciprocal.mantissa.SetHighDigit(3129112985u);\r
1004 \r
1005             constant2.exponent = 1;\r
1006             constant2.mantissa.SetHighDigit(0x80000000u);\r
1007 \r
1008             //D is deliberately left negative for all the following operations.\r
1009             thisPrec.mantissa.Sign = true;\r
1010 \r
1011             //Initial estimate.\r
1012             reciprocal.Add(thisPrec);\r
1013 \r
1014             //mantissa.Sign = false;\r
1015 \r
1016             //Shift down into 0.5 < this < 1 range\r
1017             thisPrec.mantissa.RSH(1);\r
1018 \r
1019             //Iteration.\r
1020             int accuracyBits = 2;\r
1021             int mantissaBits = mantissa.Precision.NumBits;\r
1022 \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
1030             {\r
1031                 //Increase the precision as needed\r
1032                 if (accuracyBits >= curPrec.NumBits / 2)\r
1033                 {\r
1034                     int newBits = curPrec.NumBits * 2;\r
1035                     if (newBits > mantissaBits) newBits = mantissaBits;\r
1036                     curPrec = new PrecisionSpec(newBits, PrecisionSpec.BaseType.BIN);\r
1037 \r
1038                     reciprocal = new BigFloat(reciprocal, curPrec);\r
1039 \r
1040                     constant2 = new BigFloat(curPrec);\r
1041                     constant2.exponent = 1;\r
1042                     constant2.mantissa.SetHighDigit(0x80000000u);\r
1043 \r
1044                     temp = new BigFloat(temp, curPrec);\r
1045 \r
1046                     thisPrec = new BigFloat(this, curPrec);\r
1047                     thisPrec.mantissa.Sign = true;\r
1048                     thisPrec.mantissa.RSH(1);\r
1049                 }\r
1050 \r
1051                 //temp = Xn\r
1052                 temp.exponent = reciprocal.exponent;\r
1053                 temp.mantissa.Assign(reciprocal.mantissa);\r
1054                 //temp = -Xn * D\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
1060 \r
1061                 accuracyBits *= 2;\r
1062             }\r
1063 \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
1069 \r
1070             if (topBit)\r
1071             {\r
1072                 reciprocal.exponent = -(oldExponent);\r
1073             }\r
1074             else\r
1075             {\r
1076                 reciprocal.exponent = -(oldExponent + 1);\r
1077             }\r
1078             reciprocal.mantissa.Sign = oldSign;\r
1079 \r
1080             return reciprocal;\r
1081         }\r
1082 \r
1083         /// <summary>\r
1084         /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.\r
1085         /// </summary>\r
1086         /// <returns>The reciprocal 1/this</returns>\r
1087         private BigFloat ReciprocalNewton2()\r
1088         {\r
1089             if (mantissa.IsZero())\r
1090             {\r
1091                 exponent = Int32.MaxValue;\r
1092                 return null;\r
1093             }\r
1094 \r
1095             bool oldSign = mantissa.Sign;\r
1096             int oldExponent = exponent;\r
1097 \r
1098             //Kill exponent for now (will re-institute later)\r
1099             exponent = 0;\r
1100 \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
1104 \r
1105             reciprocal.exponent = 1;\r
1106             reciprocal.mantissa.SetHighDigit(3129112985u);\r
1107 \r
1108             constant2.exponent = 1;\r
1109             constant2.mantissa.SetHighDigit(0x80000000u);\r
1110 \r
1111             //D is deliberately left negative for all the following operations.\r
1112             mantissa.Sign = true;\r
1113 \r
1114             //Initial estimate.\r
1115             reciprocal.Add(this);\r
1116 \r
1117             //mantissa.Sign = false;\r
1118 \r
1119             //Shift down into 0.5 < this < 1 range\r
1120             mantissa.RSH(1);\r
1121             \r
1122             //Iteration.\r
1123             int accuracyBits = 2;\r
1124             int mantissaBits = mantissa.Precision.NumBits;\r
1125 \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
1133             {\r
1134                 //temp = Xn\r
1135                 temp.exponent = reciprocal.exponent;\r
1136                 temp.mantissa.Assign(reciprocal.mantissa);\r
1137                 //temp = -Xn * D\r
1138                 temp.Mul(this);\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
1143 \r
1144                 accuracyBits *= 2;\r
1145             }\r
1146 \r
1147             //'reciprocal' is now the reciprocal of the shifted down, zero-exponent mantissa of 'this'\r
1148             //Restore the mantissa.\r
1149             mantissa.LSH(1);\r
1150             exponent = oldExponent;\r
1151             mantissa.Sign = oldSign;\r
1152 \r
1153             reciprocal.exponent = -(oldExponent + 1);\r
1154             reciprocal.mantissa.Sign = oldSign;\r
1155 \r
1156             return reciprocal;\r
1157         }\r
1158 \r
1159         /// <summary>\r
1160         /// Sets this equal to the input\r
1161         /// </summary>\r
1162         /// <param name="n2"></param>\r
1163         public void Assign(BigFloat n2)\r
1164         {\r
1165             exponent = n2.exponent;\r
1166             if (mantissa.AssignHigh(n2.mantissa)) exponent++;\r
1167         }\r
1168 \r
1169 \r
1170         //********************* Comparison Functions *******************\r
1171 \r
1172         /// <summary>\r
1173         /// Greater than comparison\r
1174         /// </summary>\r
1175         /// <param name="n2">the number to compare this to</param>\r
1176         /// <returns>true iff this is greater than n2 (this &gt; n2)</returns>\r
1177         public bool GreaterThan(BigFloat n2)\r
1178         {\r
1179             if (IsSpecialValue || n2.IsSpecialValue)\r
1180             {\r
1181                 SpecialValueType s1 = SpecialValue;\r
1182                 SpecialValueType s2 = SpecialValue;\r
1183 \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
1189 \r
1190                 if (s1 == SpecialValueType.ZERO)\r
1191                 {\r
1192                     if (s2 != SpecialValueType.ZERO && n2.Sign)\r
1193                     {\r
1194                         return true;\r
1195                     }\r
1196                     else\r
1197                     {\r
1198                         return false;\r
1199                     }\r
1200                 }\r
1201 \r
1202                 if (s2 == SpecialValueType.ZERO)\r
1203                 {\r
1204                     return !Sign;\r
1205                 }\r
1206             }\r
1207 \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
1211             {\r
1212                 if (exponent > n2.exponent) return true;\r
1213                 if (exponent < n2.exponent) return false;\r
1214             }\r
1215             if (mantissa.Sign)\r
1216             {\r
1217                 if (exponent > n2.exponent) return false;\r
1218                 if (exponent < n2.exponent) return true;\r
1219             }\r
1220 \r
1221             return mantissa.GreaterThan(n2.mantissa);\r
1222         }\r
1223 \r
1224         /// <summary>\r
1225         /// Less than comparison\r
1226         /// </summary>\r
1227         /// <param name="n2">the number to compare this to</param>\r
1228         /// <returns>true iff this is less than n2 (this &lt; n2)</returns>\r
1229         public bool LessThan(BigFloat n2)\r
1230         {\r
1231             if (IsSpecialValue || n2.IsSpecialValue)\r
1232             {\r
1233                 SpecialValueType s1 = SpecialValue;\r
1234                 SpecialValueType s2 = SpecialValue;\r
1235 \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
1241 \r
1242                 if (s1 == SpecialValueType.ZERO)\r
1243                 {\r
1244                     if (s2 != SpecialValueType.ZERO && !n2.Sign)\r
1245                     {\r
1246                         return true;\r
1247                     }\r
1248                     else\r
1249                     {\r
1250                         return false;\r
1251                     }\r
1252                 }\r
1253 \r
1254                 if (s2 == SpecialValueType.ZERO)\r
1255                 {\r
1256                     return Sign;\r
1257                 }\r
1258             }\r
1259 \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
1263             {\r
1264                 if (exponent > n2.exponent) return false;\r
1265                 if (exponent < n2.exponent) return true;\r
1266             }\r
1267             if (mantissa.Sign)\r
1268             {\r
1269                 if (exponent > n2.exponent) return true;\r
1270                 if (exponent < n2.exponent) return false;\r
1271             }\r
1272 \r
1273             return mantissa.LessThan(n2.mantissa);\r
1274         }\r
1275 \r
1276         /// <summary>\r
1277         /// Greater than comparison\r
1278         /// </summary>\r
1279         /// <param name="i">the number to compare this to</param>\r
1280         /// <returns>true iff this is greater than n2 (this &gt; n2)</returns>\r
1281         public bool GreaterThan(int i)\r
1282         {\r
1283             BigFloat integer = new BigFloat(i, mantissa.Precision);\r
1284             return GreaterThan(integer);\r
1285         }\r
1286 \r
1287         /// <summary>\r
1288         /// Less than comparison\r
1289         /// </summary>\r
1290         /// <param name="i">the number to compare this to</param>\r
1291         /// <returns>true iff this is less than n2 (this &lt; n2)</returns>\r
1292         public bool LessThan(int i)\r
1293         {\r
1294             BigFloat integer = new BigFloat(i, mantissa.Precision);\r
1295             return LessThan(integer);\r
1296         }\r
1297 \r
1298         /// <summary>\r
1299         /// Compare to zero\r
1300         /// </summary>\r
1301         /// <returns>true if this is zero (this == 0)</returns>\r
1302         public bool IsZero()\r
1303         {\r
1304             return (mantissa.IsZero());\r
1305         }\r
1306 \r
1307 \r
1308         //******************** Mathematical Functions ******************\r
1309 \r
1310         /// <summary>\r
1311         /// Sets the number to the biggest integer numerically closer to zero, if possible.\r
1312         /// </summary>\r
1313         public void Floor()\r
1314         {\r
1315             //Already an integer.\r
1316             if (exponent >= mantissa.Precision.NumBits) return;\r
1317 \r
1318             if (exponent < 0)\r
1319             {\r
1320                 mantissa.ZeroBits(mantissa.Precision.NumBits);\r
1321                 exponent = 0;\r
1322                 return;\r
1323             }\r
1324 \r
1325             mantissa.ZeroBits(mantissa.Precision.NumBits - (exponent + 1));\r
1326         }\r
1327 \r
1328         /// <summary>\r
1329         /// Sets the number to its fractional component (equivalent to 'this' - (int)'this')\r
1330         /// </summary>\r
1331         public void FPart()\r
1332         {\r
1333             //Already fractional\r
1334             if (exponent < 0)\r
1335             {\r
1336                 return;\r
1337             }\r
1338 \r
1339             //Has no fractional part\r
1340             if (exponent >= mantissa.Precision.NumBits)\r
1341             {\r
1342                 mantissa.Zero();\r
1343                 exponent = 0;\r
1344                 return;\r
1345             }\r
1346 \r
1347             mantissa.ZeroBitsHigh(exponent + 1);\r
1348             exponent -= mantissa.Normalise();\r
1349         }\r
1350 \r
1351         /// <summary>\r
1352         /// Calculates tan(x)\r
1353         /// </summary>\r
1354         public void Tan()\r
1355         {\r
1356             if (IsSpecialValue)\r
1357             {\r
1358                 //Tan(x) has no limit as x->inf\r
1359                 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
1360                 {\r
1361                     SetNaN();\r
1362                 }\r
1363                 else if (SpecialValue == SpecialValueType.ZERO)\r
1364                 {\r
1365                     SetZero();\r
1366                 }\r
1367 \r
1368                 return;\r
1369             }\r
1370 \r
1371             if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
1372             {\r
1373                 CalculatePi(mantissa.Precision.NumBits);\r
1374             }\r
1375 \r
1376             //Work out the sign change (involves replicating some rescaling).\r
1377             bool sign = mantissa.Sign;\r
1378             mantissa.Sign = false;\r
1379 \r
1380             if (mantissa.IsZero())\r
1381             {\r
1382                 return;\r
1383             }\r
1384 \r
1385             //Rescale into 0 <= x < pi\r
1386             if (GreaterThan(pi))\r
1387             {\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
1391                 newAngle.FPart();\r
1392                 newAngle.Mul(pi);\r
1393                 Assign(newAngle);\r
1394             }\r
1395 \r
1396             //Rescale to -pi/2 <= x < pi/2\r
1397             if (!LessThan(piBy2))\r
1398             {\r
1399                 Sub(pi);\r
1400             }\r
1401 \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
1404             Sin();\r
1405             BigFloat denom = new BigFloat(this);\r
1406             denom.Mul(this);\r
1407             denom.Sub(new BigFloat(1, mantissa.Precision));\r
1408             denom.mantissa.Sign = !denom.mantissa.Sign;\r
1409 \r
1410             if (denom.mantissa.Sign)\r
1411             {\r
1412                 denom.SetZero();\r
1413             }\r
1414 \r
1415             denom.Sqrt();\r
1416             Div(denom);\r
1417             if (sign) mantissa.Sign = !mantissa.Sign;\r
1418         }\r
1419 \r
1420         /// <summary>\r
1421         /// Calculates Cos(x)\r
1422         /// </summary>\r
1423         public void Cos()\r
1424         {\r
1425             if (IsSpecialValue)\r
1426             {\r
1427                 //Cos(x) has no limit as x->inf\r
1428                 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
1429                 {\r
1430                     SetNaN();\r
1431                 }\r
1432                 else if (SpecialValue == SpecialValueType.ZERO)\r
1433                 {\r
1434                     Assign(new BigFloat(1, mantissa.Precision));\r
1435                 }\r
1436 \r
1437                 return;\r
1438             }\r
1439 \r
1440             if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
1441             {\r
1442                 CalculatePi(mantissa.Precision.NumBits);\r
1443             }\r
1444 \r
1445             Add(piBy2);\r
1446             Sin();\r
1447         }\r
1448 \r
1449         /// <summary>\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
1452         /// </summary>\r
1453         public void Sin()\r
1454         {\r
1455             if (IsSpecialValue)\r
1456             {\r
1457                 //Sin(x) has no limit as x->inf\r
1458                 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
1459                 {\r
1460                     SetNaN();\r
1461                 }\r
1462 \r
1463                 return;\r
1464             }\r
1465 \r
1466             //Convert to positive range (0 <= x < inf)\r
1467             bool sign = mantissa.Sign;\r
1468             mantissa.Sign = false;\r
1469 \r
1470             if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
1471             {\r
1472                 CalculatePi(mantissa.Precision.NumBits);\r
1473             }\r
1474 \r
1475             if (inverseFactorialCache == null || invFactorialCutoff != mantissa.Precision.NumBits)\r
1476             {\r
1477                 CalculateFactorials(mantissa.Precision.NumBits);\r
1478             }\r
1479 \r
1480             //Rescale into 0 <= x < 2*pi\r
1481             if (GreaterThan(twoPi))\r
1482             {\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
1486                 newAngle.FPart();\r
1487                 newAngle.Mul(twoPi);\r
1488                 Assign(newAngle);\r
1489             }\r
1490 \r
1491             //Rescale into range 0 <= x < pi\r
1492             if (GreaterThan(pi))\r
1493             {\r
1494                 //sin(pi + a) = sin(pi)cos(a) + sin(a)cos(pi) = 0 - sin(a) = -sin(a)\r
1495                 Sub(pi);\r
1496                 sign = !sign;\r
1497             }\r
1498 \r
1499             BigFloat temp = new BigFloat(mantissa.Precision);\r
1500 \r
1501             //Rescale into range 0 <= x < pi/2\r
1502             if (GreaterThan(piBy2))\r
1503             {\r
1504                 temp.Assign(this);\r
1505                 Assign(pi);\r
1506                 Sub(temp);\r
1507             }\r
1508 \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
1511             Mul(threeRecip);\r
1512 \r
1513             if (mantissa.IsZero())\r
1514             {\r
1515                 exponent = 0;\r
1516                 return;\r
1517             }\r
1518 \r
1519             BigFloat term = new BigFloat(this);\r
1520 \r
1521             BigFloat square = new BigFloat(this);\r
1522             square.Mul(term);\r
1523 \r
1524             BigFloat sum = new BigFloat(this);\r
1525 \r
1526             bool termSign = true;\r
1527             int length = inverseFactorialCache.Length;\r
1528             int numBits = mantissa.Precision.NumBits;\r
1529 \r
1530             for (int i = 3; i < length; i += 2)\r
1531             {\r
1532                 term.Mul(square);\r
1533                 temp.Assign(inverseFactorialCache[i]);\r
1534                 temp.Mul(term);\r
1535                 temp.mantissa.Sign = termSign;\r
1536                 termSign = !termSign;\r
1537 \r
1538                 if (temp.exponent < -numBits) break;\r
1539 \r
1540                 sum.Add(temp);\r
1541             }\r
1542 \r
1543             //Restore the triple-angle: sin(3x) = 3sin(x) - 4sin^3(x)\r
1544             Assign(sum);\r
1545             sum.Mul(this);\r
1546             sum.Mul(this);\r
1547             Mul(new BigFloat(3, mantissa.Precision));\r
1548             sum.exponent += 2;\r
1549             Sub(sum);\r
1550 \r
1551             //Restore the sign\r
1552             mantissa.Sign = sign;\r
1553         }\r
1554 \r
1555         /// <summary>\r
1556         /// Hyperbolic Sin (sinh) function\r
1557         /// </summary>\r
1558         public void Sinh()\r
1559         {\r
1560             if (IsSpecialValue)\r
1561             {\r
1562                 return;\r
1563             }\r
1564 \r
1565             Exp();\r
1566             Sub(Reciprocal());\r
1567             exponent--;\r
1568         }\r
1569 \r
1570         /// <summary>\r
1571         /// Hyperbolic cosine (cosh) function\r
1572         /// </summary>\r
1573         public void Cosh()\r
1574         {\r
1575             if (IsSpecialValue)\r
1576             {\r
1577                 if (SpecialValue == SpecialValueType.ZERO)\r
1578                 {\r
1579                     Assign(new BigFloat(1, mantissa.Precision));\r
1580                 }\r
1581                 else if (SpecialValue == SpecialValueType.INF_MINUS)\r
1582                 {\r
1583                     SetInfPlus();\r
1584                 }\r
1585 \r
1586                 return;\r
1587             }\r
1588 \r
1589             Exp();\r
1590             Add(Reciprocal());\r
1591             exponent--;\r
1592         }\r
1593 \r
1594         /// <summary>\r
1595         /// Hyperbolic tangent function (tanh)\r
1596         /// </summary>\r
1597         public void Tanh()\r
1598         {\r
1599             if (IsSpecialValue)\r
1600             {\r
1601                 if (SpecialValue == SpecialValueType.INF_MINUS)\r
1602                 {\r
1603                     Assign(new BigFloat(-1, mantissa.Precision));\r
1604                 }\r
1605                 else if (SpecialValue == SpecialValueType.INF_PLUS)\r
1606                 {\r
1607                     Assign(new BigFloat(1, mantissa.Precision));\r
1608                 }\r
1609 \r
1610                 return;\r
1611             }\r
1612 \r
1613             exponent++;\r
1614             Exp();\r
1615             BigFloat temp = new BigFloat(this);\r
1616             BigFloat one = new BigFloat(1, mantissa.Precision);\r
1617             temp.Add(one);\r
1618             Sub(one);\r
1619             Div(temp);\r
1620         }\r
1621 \r
1622         /// <summary>\r
1623         /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)\r
1624         /// </summary>\r
1625         public void Arcsin()\r
1626         {\r
1627             if (IsSpecialValue)\r
1628             {\r
1629                 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)\r
1630                 {\r
1631                     SetNaN();\r
1632                     return;\r
1633                 }\r
1634 \r
1635                 return;\r
1636             }\r
1637 \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
1643 \r
1644             bool sign = mantissa.Sign;\r
1645             mantissa.Sign = false;\r
1646 \r
1647             if (GreaterThan(onePlusABit))\r
1648             {\r
1649                 SetNaN();\r
1650             }\r
1651             else if (LessThan(one))\r
1652             {\r
1653                 BigFloat temp = new BigFloat(this);\r
1654                 temp.Mul(this);\r
1655                 temp.Sub(one);\r
1656                 temp.mantissa.Sign = !temp.mantissa.Sign;\r
1657                 temp.Sqrt();\r
1658                 temp.Add(one);\r
1659                 Div(temp);\r
1660                 Arctan();\r
1661                 exponent++;\r
1662                 mantissa.Sign = sign;\r
1663             }\r
1664             else\r
1665             {\r
1666                 if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
1667                 {\r
1668                     CalculatePi(mantissa.Precision.NumBits);\r
1669                 }\r
1670 \r
1671                 Assign(piBy2);\r
1672                 if (sign) mantissa.Sign = true;\r
1673             }\r
1674         }\r
1675 \r
1676         /// <summary>\r
1677         /// arccos(): the inverse function of cos(), range (0..pi)\r
1678         /// </summary>\r
1679         public void Arccos()\r
1680         {\r
1681             if (IsSpecialValue)\r
1682             {\r
1683                 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)\r
1684                 {\r
1685                     SetNaN();\r
1686                 }\r
1687                 else if (SpecialValue == SpecialValueType.ZERO)\r
1688                 {\r
1689                     Assign(new BigFloat(1, mantissa.Precision));\r
1690                     exponent = 0;\r
1691                     Sign = false;\r
1692                 }\r
1693 \r
1694                 return;\r
1695             }\r
1696 \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
1702 \r
1703             bool sign = mantissa.Sign;\r
1704             mantissa.Sign = false;\r
1705 \r
1706             if (GreaterThan(onePlusABit))\r
1707             {\r
1708                 SetNaN();\r
1709             }\r
1710             else if (LessThan(one))\r
1711             {\r
1712                 if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
1713                 {\r
1714                     CalculatePi(mantissa.Precision.NumBits);\r
1715                 }\r
1716 \r
1717                 mantissa.Sign = sign;\r
1718                 BigFloat temp = new BigFloat(this);\r
1719                 Mul(temp);\r
1720                 Sub(one);\r
1721                 mantissa.Sign = !mantissa.Sign;\r
1722                 Sqrt();\r
1723                 temp.Add(one);\r
1724                 Div(temp);\r
1725                 Arctan();\r
1726                 exponent++;\r
1727             }\r
1728             else\r
1729             {\r
1730                 if (sign)\r
1731                 {\r
1732                     if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)\r
1733                     {\r
1734                         CalculatePi(mantissa.Precision.NumBits);\r
1735                     }\r
1736 \r
1737                     Assign(pi);\r
1738                 }\r
1739                 else\r
1740                 {\r
1741                     mantissa.Zero();\r
1742                     exponent = 0;\r
1743                 }\r
1744             }\r
1745         }\r
1746 \r
1747         /// <summary>\r
1748         /// arctan(): the inverse function of sin(), range of (-pi/2..pi/2)\r
1749         /// </summary>\r
1750         public void Arctan()\r
1751         {\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
1755 \r
1756             if (pi == null || pi.mantissa.Precision.NumBits != numBits)\r
1757             {\r
1758                 CalculatePi(mantissa.Precision.NumBits);\r
1759             }\r
1760 \r
1761             //Make domain positive\r
1762             bool sign = mantissa.Sign;\r
1763             mantissa.Sign = false;\r
1764 \r
1765             if (IsSpecialValue)\r
1766             {\r
1767                 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
1768                 {\r
1769                     Assign(piBy2);\r
1770                     mantissa.Sign = sign;\r
1771                     return;\r
1772                 }\r
1773 \r
1774                 return;\r
1775             }\r
1776 \r
1777             if (reciprocals == null || reciprocals[0].mantissa.Precision.NumBits != numBits || reciprocals.Length < maxTerms)\r
1778             {\r
1779                 CalculateReciprocals(numBits, maxTerms);\r
1780             }\r
1781 \r
1782             bool invert = false;\r
1783             BigFloat one = new BigFloat(1, mantissa.Precision);\r
1784 \r
1785             //Invert if outside of convergence\r
1786             if (GreaterThan(one))\r
1787             {\r
1788                 invert = true;\r
1789                 Assign(Reciprocal());\r
1790             }\r
1791 \r
1792             //Reduce using half-angle formula:\r
1793             //arctan(2x) = 2 arctan (x / (1 + sqrt(1 + x)))\r
1794 \r
1795             //First reduction (guarantees 2 bits per iteration)\r
1796             BigFloat temp = new BigFloat(this);\r
1797             temp.Mul(this);\r
1798             temp.Add(one);\r
1799             temp.Sqrt();\r
1800             temp.Add(one);\r
1801             this.Div(temp);\r
1802 \r
1803             //Second reduction (guarantees 4 bits per iteration)\r
1804             temp.Assign(this);\r
1805             temp.Mul(this);\r
1806             temp.Add(one);\r
1807             temp.Sqrt();\r
1808             temp.Add(one);\r
1809             this.Div(temp);\r
1810 \r
1811             //Actual series calculation\r
1812             int length = reciprocals.Length;\r
1813             BigFloat term = new BigFloat(this);\r
1814 \r
1815             //pow = x^2\r
1816             BigFloat pow = new BigFloat(this);\r
1817             pow.Mul(this);\r
1818 \r
1819             BigFloat sum = new BigFloat(this);\r
1820 \r
1821             for (int i = 1; i < length; i++)\r
1822             {\r
1823                 //u(n) = u(n-1) * x^2\r
1824                 //t(n) = u(n) / (2n+1)\r
1825                 term.Mul(pow);\r
1826                 term.Sign = !term.Sign;\r
1827                 temp.Assign(term);\r
1828                 temp.Mul(reciprocals[i]);\r
1829 \r
1830                 if (temp.exponent < -numBits) break;\r
1831 \r
1832                 sum.Add(temp);\r
1833             }\r
1834 \r
1835             //Undo the reductions.\r
1836             Assign(sum);\r
1837             exponent += 2;\r
1838 \r
1839             if (invert)\r
1840             {\r
1841                 //Assign(Reciprocal());\r
1842                 mantissa.Sign = true;\r
1843                 Add(piBy2);\r
1844             }\r
1845 \r
1846             if (sign)\r
1847             {\r
1848                 mantissa.Sign = sign;\r
1849             }\r
1850         }\r
1851 \r
1852         /// <summary>\r
1853         /// Arcsinh(): the inverse sinh function\r
1854         /// </summary>\r
1855         public void Arcsinh()\r
1856         {\r
1857             //Just let all special values fall through\r
1858             if (IsSpecialValue)\r
1859             {\r
1860                 return;\r
1861             }\r
1862 \r
1863             BigFloat temp = new BigFloat(this);\r
1864             temp.Mul(this);\r
1865             temp.Add(new BigFloat(1, mantissa.Precision));\r
1866             temp.Sqrt();\r
1867             Add(temp);\r
1868             Log();\r
1869         }\r
1870 \r
1871         /// <summary>\r
1872         /// Arccosh(): the inverse cosh() function\r
1873         /// </summary>\r
1874         public void Arccosh()\r
1875         {\r
1876             //acosh isn't defined for x < 1\r
1877             if (IsSpecialValue)\r
1878             {\r
1879                 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.ZERO)\r
1880                 {\r
1881                     SetNaN();\r
1882                     return;\r
1883                 }\r
1884 \r
1885                 return;\r
1886             }\r
1887 \r
1888             BigFloat one = new BigFloat(1, mantissa.Precision);\r
1889             if (LessThan(one))\r
1890             {\r
1891                 SetNaN();\r
1892                 return;\r
1893             }\r
1894 \r
1895             BigFloat temp = new BigFloat(this);\r
1896             temp.Mul(this);\r
1897             temp.Sub(one);\r
1898             temp.Sqrt();\r
1899             Add(temp);\r
1900             Log();\r
1901         }\r
1902 \r
1903         /// <summary>\r
1904         /// Arctanh(): the inverse tanh function\r
1905         /// </summary>\r
1906         public void Arctanh()\r
1907         {\r
1908             //|x| <= 1 for a non-NaN output\r
1909             if (IsSpecialValue)\r
1910             {\r
1911                 if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)\r
1912                 {\r
1913                     SetNaN();\r
1914                     return;\r
1915                 }\r
1916 \r
1917                 return;\r
1918             }\r
1919 \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
1925 \r
1926             bool sign = mantissa.Sign;\r
1927             mantissa.Sign = false;\r
1928 \r
1929             if (GreaterThan(onePlusABit))\r
1930             {\r
1931                 SetNaN();\r
1932             }\r
1933             else if (LessThan(one))\r
1934             {\r
1935                 BigFloat temp = new BigFloat(this);\r
1936                 Add(one);\r
1937                 one.Sub(temp);\r
1938                 Div(one);\r
1939                 Log();\r
1940                 exponent--;\r
1941                 mantissa.Sign = sign;\r
1942             }\r
1943             else\r
1944             {\r
1945                 if (sign)\r
1946                 {\r
1947                     SetInfMinus();\r
1948                 }\r
1949                 else\r
1950                 {\r
1951                     SetInfPlus();\r
1952                 }\r
1953             }\r
1954         }\r
1955 \r
1956         /// <summary>\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
1959         /// </summary>\r
1960         public void Sqrt()\r
1961         {\r
1962             if (mantissa.Sign || IsSpecialValue)\r
1963             {\r
1964                 if (SpecialValue == SpecialValueType.ZERO)\r
1965                 {\r
1966                     return;\r
1967                 }\r
1968 \r
1969                 if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)\r
1970                 {\r
1971                     SetNaN();\r
1972                 }\r
1973 \r
1974                 return;\r
1975             }\r
1976 \r
1977             BigFloat temp2;\r
1978             BigFloat temp3 = new BigFloat(mantissa.Precision);\r
1979             BigFloat three = new BigFloat(3, mantissa.Precision);\r
1980 \r
1981             int exponentScale = 0;\r
1982 \r
1983             //Rescale to 0.5 <= x < 2\r
1984             if (exponent < -1)\r
1985             {\r
1986                 int diff = -exponent;\r
1987                 if ((diff & 1) != 0)\r
1988                 {\r
1989                     diff--;\r
1990                 }\r
1991 \r
1992                 exponentScale = -diff;\r
1993                 exponent += diff;\r
1994             }\r
1995             else if (exponent > 0)\r
1996             {\r
1997                 if ((exponent & 1) != 0)\r
1998                 {\r
1999                     exponentScale = exponent + 1;\r
2000                     exponent = -1;\r
2001                 }\r
2002                 else\r
2003                 {\r
2004                     exponentScale = exponent;\r
2005                     exponent = 0;\r
2006                 }\r
2007             }\r
2008 \r
2009             temp2 = new BigFloat(this);\r
2010             temp2.Sub(new BigFloat(1, mantissa.Precision));\r
2011 \r
2012             //if (temp2.mantissa.IsZero())\r
2013             //{\r
2014             //    exponent += exponentScale;\r
2015             //    return;\r
2016             //}\r
2017 \r
2018             int numBits = mantissa.Precision.NumBits;\r
2019 \r
2020             while ((exponent - temp2.exponent) < numBits && temp2.SpecialValue != SpecialValueType.ZERO)\r
2021             {\r
2022                 //a(n+1) = an - an*cn / 2\r
2023                 temp3.Assign(this);\r
2024                 temp3.Mul(temp2);\r
2025                 temp3.MulPow2(-1);\r
2026                 this.Sub(temp3);\r
2027 \r
2028                 //c(n+1) = cn^2 * (cn - 3) / 4\r
2029                 temp3.Assign(temp2);\r
2030                 temp2.Sub(three);\r
2031                 temp2.Mul(temp3);\r
2032                 temp2.Mul(temp3);\r
2033                 temp2.MulPow2(-2);\r
2034             }\r
2035 \r
2036             exponent += (exponentScale >> 1);\r
2037         }\r
2038 \r
2039         /// <summary>\r
2040         /// The natural logarithm, ln(x)\r
2041         /// </summary>\r
2042         public void Log()\r
2043         {\r
2044             if (IsSpecialValue || mantissa.Sign)\r
2045             {\r
2046                 if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)\r
2047                 {\r
2048                     SetNaN();\r
2049                 }\r
2050                 else if (SpecialValue == SpecialValueType.ZERO)\r
2051                 {\r
2052                     SetInfMinus();\r
2053                 }\r
2054 \r
2055                 return;\r
2056             }\r
2057 \r
2058             if (mantissa.Precision.NumBits >= 512)\r
2059             {\r
2060                 LogAGM1();\r
2061                 return;\r
2062             }\r
2063 \r
2064             //Compute ln2.\r
2065             if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)\r
2066             {\r
2067                 CalculateLog2(mantissa.Precision.NumBits);\r
2068             }\r
2069 \r
2070             Log2();\r
2071             Mul(ln2cache);\r
2072         }\r
2073 \r
2074         /// <summary>\r
2075         /// Log to the base 10\r
2076         /// </summary>\r
2077         public void Log10()\r
2078         {\r
2079             if (IsSpecialValue || mantissa.Sign)\r
2080             {\r
2081                 if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)\r
2082                 {\r
2083                     SetNaN();\r
2084                 }\r
2085                 else if (SpecialValue == SpecialValueType.ZERO)\r
2086                 {\r
2087                     SetInfMinus();\r
2088                 }\r
2089 \r
2090                 return;\r
2091             }\r
2092 \r
2093             //Compute ln2.\r
2094             if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)\r
2095             {\r
2096                 CalculateLog2(mantissa.Precision.NumBits);\r
2097             }\r
2098 \r
2099             Log();\r
2100             Mul(log10recip);\r
2101         }\r
2102 \r
2103         /// <summary>\r
2104         /// The exponential function. Less accurate for high exponents, scales poorly with the number\r
2105         /// of bits.\r
2106         /// </summary>\r
2107         public void Exp()\r
2108         {\r
2109             Exp(mantissa.Precision.NumBits);\r
2110         }\r
2111 \r
2112         /// <summary>\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
2116         /// \r
2117         /// </summary>\r
2118         /// <param name="power"></param>\r
2119         public void Pow(int power)\r
2120         {\r
2121             BigFloat acc = new BigFloat(1, mantissa.Precision);\r
2122             BigFloat temp = new BigFloat(1, mantissa.Precision);\r
2123 \r
2124             int powerTemp = power;\r
2125 \r
2126             if (power < 0)\r
2127             {\r
2128                 Assign(Reciprocal());\r
2129                 powerTemp = -power;\r
2130             }\r
2131 \r
2132             //Fast power function\r
2133             while (powerTemp != 0)\r
2134             {\r
2135                 temp.Mul(this);\r
2136                 Assign(temp);\r
2137 \r
2138                 if ((powerTemp & 1) != 0)\r
2139                 {\r
2140                     acc.Mul(temp);\r
2141                 }\r
2142 \r
2143                 powerTemp >>= 1;\r
2144             }\r
2145 \r
2146             Assign(acc);\r
2147         }\r
2148 \r
2149         /// <summary>\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
2156         /// Domain Note:\r
2157         ///    This only works for powers of positive real numbers. Negative numbers will fail.\r
2158         /// </summary>\r
2159         /// <param name="power"></param>\r
2160         public void Pow(BigFloat power)\r
2161         {\r
2162             Log();\r
2163             Mul(power);\r
2164             Exp();\r
2165         }\r
2166 \r
2167 \r
2168         //******************** Static Math Functions *******************\r
2169 \r
2170         /// <summary>\r
2171         /// Returns the integer component of the input\r
2172         /// </summary>\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
2177         {\r
2178             BigFloat res = new BigFloat(n1);\r
2179             n1.Floor();\r
2180             return n1;\r
2181         }\r
2182 \r
2183         /// <summary>\r
2184         /// Returns the fractional (non-integer component of the input)\r
2185         /// </summary>\r
2186         /// <param name="n1">The input number</param>\r
2187         public static BigFloat FPart(BigFloat n1)\r
2188         {\r
2189             BigFloat res = new BigFloat(n1);\r
2190             n1.FPart();\r
2191             return n1;\r
2192         }\r
2193 \r
2194         /// <summary>\r
2195         /// Calculates tan(x)\r
2196         /// </summary>\r
2197         /// <param name="n1">The angle (in radians) to find the tangent of</param>\r
2198         public static BigFloat Tan(BigFloat n1)\r
2199         {\r
2200             BigFloat res = new BigFloat(n1);\r
2201             n1.Tan();\r
2202             return n1;\r
2203         }\r
2204 \r
2205         /// <summary>\r
2206         /// Calculates Cos(x)\r
2207         /// </summary>\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
2212         {\r
2213             BigFloat res = new BigFloat(n1);\r
2214             n1.Cos();\r
2215             return n1;\r
2216         }\r
2217 \r
2218         /// <summary>\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
2221         /// </summary>\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
2226         {\r
2227             BigFloat res = new BigFloat(n1);\r
2228             n1.Sin();\r
2229             return n1;\r
2230         }\r
2231 \r
2232         /// <summary>\r
2233         /// Hyperbolic Sin (sinh) function\r
2234         /// </summary>\r
2235         /// <param name="n1">The number to find the hyperbolic sine of</param>\r
2236         public static BigFloat Sinh(BigFloat n1)\r
2237         {\r
2238             BigFloat res = new BigFloat(n1);\r
2239             n1.Sinh();\r
2240             return n1;\r
2241         }\r
2242 \r
2243         /// <summary>\r
2244         /// Hyperbolic cosine (cosh) function\r
2245         /// </summary>\r
2246         /// <param name="n1">The number to find the hyperbolic cosine of</param>\r
2247         public static BigFloat Cosh(BigFloat n1)\r
2248         {\r
2249             BigFloat res = new BigFloat(n1);\r
2250             n1.Cosh();\r
2251             return n1;\r
2252         }\r
2253 \r
2254         /// <summary>\r
2255         /// Hyperbolic tangent function (tanh)\r
2256         /// </summary>\r
2257         /// <param name="n1">The number to find the hyperbolic tangent of</param>\r
2258         public static BigFloat Tanh(BigFloat n1)\r
2259         {\r
2260             BigFloat res = new BigFloat(n1);\r
2261             n1.Tanh();\r
2262             return n1;\r
2263         }\r
2264 \r
2265         /// <summary>\r
2266         /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)\r
2267         /// </summary>\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
2271         /// </remarks>\r
2272         public static BigFloat Arcsin(BigFloat n1)\r
2273         {\r
2274             BigFloat res = new BigFloat(n1);\r
2275             n1.Arcsin();\r
2276             return n1;\r
2277         }\r
2278 \r
2279         /// <summary>\r
2280         /// arccos(): the inverse function of cos(), input range (0..pi)\r
2281         /// </summary>\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
2285         /// </remarks>\r
2286         public static BigFloat Arccos(BigFloat n1)\r
2287         {\r
2288             BigFloat res = new BigFloat(n1);\r
2289             n1.Arccos();\r
2290             return n1;\r
2291         }\r
2292 \r
2293         /// <summary>\r
2294         /// arctan(): the inverse function of sin(), input range of (-pi/2..pi/2)\r
2295         /// </summary>\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
2299         /// </remarks>\r
2300         public static BigFloat Arctan(BigFloat n1)\r
2301         {\r
2302             BigFloat res = new BigFloat(n1);\r
2303             n1.Arctan();\r
2304             return n1;\r
2305         }\r
2306 \r
2307         /// <summary>\r
2308         /// Arcsinh(): the inverse sinh function\r
2309         /// </summary>\r
2310         /// <param name="n1">The number to find the inverse hyperbolic sine of</param>\r
2311         public static BigFloat Arcsinh(BigFloat n1)\r
2312         {\r
2313             BigFloat res = new BigFloat(n1);\r
2314             n1.Arcsinh();\r
2315             return n1;\r
2316         }\r
2317 \r
2318         /// <summary>\r
2319         /// Arccosh(): the inverse cosh() function\r
2320         /// </summary>\r
2321         /// <param name="n1">The number to find the inverse hyperbolic cosine of</param>\r
2322         public static BigFloat Arccosh(BigFloat n1)\r
2323         {\r
2324             BigFloat res = new BigFloat(n1);\r
2325             n1.Arccosh();\r
2326             return n1;\r
2327         }\r
2328 \r
2329         /// <summary>\r
2330         /// Arctanh(): the inverse tanh function\r
2331         /// </summary>\r
2332         /// <param name="n1">The number to fine the inverse hyperbolic tan of</param>\r
2333         public static BigFloat Arctanh(BigFloat n1)\r
2334         {\r
2335             BigFloat res = new BigFloat(n1);\r
2336             n1.Arctanh();\r
2337             return n1;\r
2338         }\r
2339 \r
2340         /// <summary>\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
2343         /// </summary>\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
2346         /// </remarks>\r
2347         public static BigFloat Sqrt(BigFloat n1)\r
2348         {\r
2349             BigFloat res = new BigFloat(n1);\r
2350             n1.Sqrt();\r
2351             return n1;\r
2352         }\r
2353 \r
2354         /// <summary>\r
2355         /// The natural logarithm, ln(x) (log base e)\r
2356         /// </summary>\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
2360         /// </remarks>\r
2361         /// <param name="n1">The number to find the natural logarithm of</param>\r
2362         public static BigFloat Log(BigFloat n1)\r
2363         {\r
2364             BigFloat res = new BigFloat(n1);\r
2365             n1.Log();\r
2366             return n1;\r
2367         }\r
2368 \r
2369         /// <summary>\r
2370         /// Base 10 logarithm of a number\r
2371         /// </summary>\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
2375         /// </remarks>\r
2376         /// <param name="n1">The number to find the base 10 logarithm of</param>\r
2377         public static BigFloat Log10(BigFloat n1)\r
2378         {\r
2379             BigFloat res = new BigFloat(n1);\r
2380             n1.Log10();\r
2381             return n1;\r
2382         }\r
2383 \r
2384         /// <summary>\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
2387         /// </summary>\r
2388         public static BigFloat Exp(BigFloat n1)\r
2389         {\r
2390             BigFloat res = new BigFloat(n1);\r
2391             n1.Exp();\r
2392             return n1;\r
2393         }\r
2394 \r
2395         /// <summary>\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
2399         /// </summary>\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
2403         {\r
2404             BigFloat res = new BigFloat(n1);\r
2405             n1.Pow(power);\r
2406             return n1;\r
2407         }\r
2408 \r
2409         /// <summary>\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
2412         /// </summary>\r
2413         /// <remarks>\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
2418         ///    \r
2419         /// Domain Note:\r
2420         ///    This only works for powers of positive real numbers. Negative numbers will fail.\r
2421         /// </remarks>\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
2425         {\r
2426             BigFloat res = new BigFloat(n1);\r
2427             n1.Pow(power);\r
2428             return n1;\r
2429         }\r
2430 \r
2431         //********************** Static functions **********************\r
2432 \r
2433         /// <summary>\r
2434         /// Adds two numbers and returns the result\r
2435         /// </summary>\r
2436         public static BigFloat Add(BigFloat n1, BigFloat n2)\r
2437         {\r
2438             BigFloat ret = new BigFloat(n1);\r
2439             ret.Add(n2);\r
2440             return ret;\r
2441         }\r
2442 \r
2443         /// <summary>\r
2444         /// Subtracts two numbers and returns the result\r
2445         /// </summary>\r
2446         public static BigFloat Sub(BigFloat n1, BigFloat n2)\r
2447         {\r
2448             BigFloat ret = new BigFloat(n1);\r
2449             ret.Sub(n2);\r
2450             return ret;\r
2451         }\r
2452 \r
2453         /// <summary>\r
2454         /// Multiplies two numbers and returns the result\r
2455         /// </summary>\r
2456         public static BigFloat Mul(BigFloat n1, BigFloat n2)\r
2457         {\r
2458             BigFloat ret = new BigFloat(n1);\r
2459             ret.Mul(n2);\r
2460             return ret;\r
2461         }\r
2462 \r
2463         /// <summary>\r
2464         /// Divides two numbers and returns the result\r
2465         /// </summary>\r
2466         public static BigFloat Div(BigFloat n1, BigFloat n2)\r
2467         {\r
2468             BigFloat ret = new BigFloat(n1);\r
2469             ret.Div(n2);\r
2470             return ret;\r
2471         }\r
2472 \r
2473         /// <summary>\r
2474         /// Tests whether n1 is greater than n2\r
2475         /// </summary>\r
2476         public static bool GreaterThan(BigFloat n1, BigFloat n2)\r
2477         {\r
2478             return n1.GreaterThan(n2);\r
2479         }\r
2480 \r
2481         /// <summary>\r
2482         /// Tests whether n1 is less than n2\r
2483         /// </summary>\r
2484         public static bool LessThan(BigFloat n1, BigFloat n2)\r
2485         {\r
2486             return n1.LessThan(n2);\r
2487         }\r
2488 \r
2489 \r
2490         //******************* Fast static functions ********************\r
2491 \r
2492         /// <summary>\r
2493         /// Adds two numbers and assigns the result to res.\r
2494         /// </summary>\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
2500         {\r
2501             res.Assign(n1);\r
2502             res.Add(n2);\r
2503             return res;\r
2504         }\r
2505 \r
2506         /// <summary>\r
2507         /// Subtracts two numbers and assigns the result to res.\r
2508         /// </summary>\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
2514         {\r
2515             res.Assign(n1);\r
2516             res.Sub(n2);\r
2517             return res;\r
2518         }\r
2519 \r
2520         /// <summary>\r
2521         /// Multiplies two numbers and assigns the result to res.\r
2522         /// </summary>\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
2528         {\r
2529             res.Assign(n1);\r
2530             res.Mul(n2);\r
2531             return res;\r
2532         }\r
2533 \r
2534         /// <summary>\r
2535         /// Divides two numbers and assigns the result to res.\r
2536         /// </summary>\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
2542         {\r
2543             res.Assign(n1);\r
2544             res.Div(n2);\r
2545             return res;\r
2546         }\r
2547 \r
2548 \r
2549         //************************* Operators **************************\r
2550 \r
2551         /// <summary>\r
2552         /// The addition operator\r
2553         /// </summary>\r
2554         public static BigFloat operator +(BigFloat n1, BigFloat n2)\r
2555         {\r
2556             return Add(n1, n2);\r
2557         }\r
2558 \r
2559         /// <summary>\r
2560         /// The subtraction operator\r
2561         /// </summary>\r
2562         public static BigFloat operator -(BigFloat n1, BigFloat n2)\r
2563         {\r
2564             return Sub(n1, n2);\r
2565         }\r
2566 \r
2567         /// <summary>\r
2568         /// The multiplication operator\r
2569         /// </summary>\r
2570         public static BigFloat operator *(BigFloat n1, BigFloat n2)\r
2571         {\r
2572             return Mul(n1, n2);\r
2573         }\r
2574 \r
2575         /// <summary>\r
2576         /// The division operator\r
2577         /// </summary>\r
2578         public static BigFloat operator /(BigFloat n1, BigFloat n2)\r
2579         {\r
2580             return Div(n1, n2);\r
2581         }\r
2582 \r
2583         //************************** Conversions *************************\r
2584 \r
2585         /// <summary>\r
2586         /// Converts a BigFloat to an BigInt with the specified precision\r
2587         /// </summary>\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
2593         {\r
2594             BigInt ret = new BigInt(precision);\r
2595 \r
2596             int numBits = n1.mantissa.Precision.NumBits;\r
2597             int shift = numBits - (n1.exponent + 1);\r
2598 \r
2599             BigFloat copy = new BigFloat(n1);\r
2600             bool inc = false;\r
2601 \r
2602             //Rounding\r
2603             if (copy.mantissa.Precision.NumBits > ret.Precision.NumBits)\r
2604             {\r
2605                 inc = true;\r
2606 \r
2607                 for (int i = copy.exponent + 1; i <= ret.Precision.NumBits; i++)\r
2608                 {\r
2609                     if (copy.mantissa.GetBitFromTop(i) == 0)\r
2610                     {\r
2611                         inc = false;\r
2612                         break;\r
2613                     }\r
2614                 }\r
2615             }\r
2616 \r
2617             if (shift > 0)\r
2618             {\r
2619                 copy.mantissa.RSH(shift);\r
2620             }\r
2621             else if (shift < 0)\r
2622             {\r
2623                 copy.mantissa.LSH(-shift);\r
2624             }\r
2625 \r
2626             ret.Assign(copy.mantissa);\r
2627 \r
2628             if (inc) ret.Increment();\r
2629 \r
2630             return ret;\r
2631         }\r
2632 \r
2633         /// <summary>\r
2634         /// Returns a base-10 string representing the number.\r
2635         /// \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
2638         /// </summary>\r
2639         public override string ToString()\r
2640         {\r
2641             if (IsSpecialValue)\r
2642             {\r
2643                 SpecialValueType s = SpecialValue;\r
2644                 if (s == SpecialValueType.ZERO)\r
2645                 {\r
2646                     return String.Format("0{0}0", System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator);\r
2647                 }\r
2648                 else if (s == SpecialValueType.INF_PLUS)\r
2649                 {\r
2650                     return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol;\r
2651                 }\r
2652                 else if (s == SpecialValueType.INF_MINUS)\r
2653                 {\r
2654                     return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol;\r
2655                 }\r
2656                 else if (s == SpecialValueType.NAN)\r
2657                 {\r
2658                     return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol;\r
2659                 }\r
2660                 else\r
2661                 {\r
2662                     return "Unrecognised special type";\r
2663                 }\r
2664             }\r
2665 \r
2666             if (scratch.Precision.NumBits != mantissa.Precision.NumBits)\r
2667             {\r
2668                 scratch = new BigInt(mantissa.Precision);\r
2669             }\r
2670 \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
2678 \r
2679             bool useExponentialNotation = false;\r
2680             int halfBits = mantissa.Precision.NumBits / 2;\r
2681             if (halfBits > 60) halfBits = 60;\r
2682             int precDec = 10;\r
2683 \r
2684             if (exponent > 0)\r
2685             {\r
2686                 if (exponent < halfBits)\r
2687                 {\r
2688                     denom.RSH(exponent);\r
2689                 }\r
2690                 else\r
2691                 {\r
2692                     useExponentialNotation = true;\r
2693                 }\r
2694             }\r
2695             else if (exponent < 0)\r
2696             {\r
2697                 int shift = -(exponent);\r
2698                 if (shift < precDec)\r
2699                 {\r
2700                     scratch.RSH(shift);\r
2701                 }\r
2702                 else\r
2703                 {\r
2704                     useExponentialNotation = true;\r
2705                 }\r
2706             }\r
2707 \r
2708             string output;\r
2709 \r
2710             if (useExponentialNotation)\r
2711             {\r
2712                 int absExponent = exponent;\r
2713                 if (absExponent < 0) absExponent = -absExponent;\r
2714                 int powerOf10 = (int)((double)absExponent * Math.Log10(2.0));\r
2715 \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
2719 \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
2724 \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
2727 \r
2728                 BigFloat tenToUse;\r
2729 \r
2730                 if (exponent > 0)\r
2731                 {\r
2732                     tenToUse = new BigFloat(tenRCP, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
2733                 }\r
2734                 else\r
2735                 {\r
2736                     tenToUse = new BigFloat(ten, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
2737                 }\r
2738 \r
2739                 BigFloat tenToPower = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));\r
2740 \r
2741                 int powerTemp = powerOf10;\r
2742 \r
2743                 //Fast power function\r
2744                 while (powerTemp != 0)\r
2745                 {\r
2746                     tenToPower.Mul(tenToUse);\r
2747                     tenToUse.Assign(tenToPower);\r
2748 \r
2749                     if ((powerTemp & 1) != 0)\r
2750                     {\r
2751                         acc.Mul(tenToPower);\r
2752                     }\r
2753 \r
2754                     powerTemp >>= 1;\r
2755                 }\r
2756 \r
2757                 thisFloat.Mul(acc);\r
2758 \r
2759                 //If we are out of range, correct.           \r
2760                 if (thisFloat.GreaterThan(ten))\r
2761                 {\r
2762                     thisFloat.Mul(tenRCP);\r
2763                     if (exponent > 0)\r
2764                     {\r
2765                         powerOf10++;\r
2766                     }\r
2767                     else\r
2768                     {\r
2769                         powerOf10--;\r
2770                     }\r
2771                 }\r
2772                 else if (thisFloat.LessThan(one))\r
2773                 {\r
2774                     thisFloat.Mul(ten);\r
2775                     if (exponent > 0)\r
2776                     {\r
2777                         powerOf10--;\r
2778                     }\r
2779                     else\r
2780                     {\r
2781                         powerOf10++;\r
2782                     }\r
2783                 }\r
2784 \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
2789 \r
2790                 if (exponent < 0) powerOf10 = -powerOf10;\r
2791 \r
2792                 output = String.Format("{0}E{1}", output, powerOf10);\r
2793             }\r
2794             else\r
2795             {\r
2796                 BigInt bigDigit = BigInt.Div(scratch, denom);\r
2797                 bigDigit.Sign = false;\r
2798                 scratch.Sub(BigInt.Mul(denom, bigDigit));\r
2799 \r
2800                 if (mantissa.Sign)\r
2801                 {\r
2802                     output = String.Format("-{0}.", bigDigit);\r
2803                 }\r
2804                 else\r
2805                 {\r
2806                     output = String.Format("{0}.", bigDigit);\r
2807                 }\r
2808 \r
2809                 denom = BigInt.Div(denom, 10u);\r
2810 \r
2811                 while (!denom.IsZero())\r
2812                 {\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
2818                 }\r
2819 \r
2820                 output = RoundString(output, RoundingDigits);\r
2821             }\r
2822 \r
2823             return output;\r
2824         }\r
2825 \r
2826         //**************** Special value handling for ops ***************\r
2827 \r
2828         private void SetNaN()\r
2829         {\r
2830             exponent = Int32.MaxValue;\r
2831             mantissa.SetBit(mantissa.Precision.NumBits - 1);\r
2832         }\r
2833 \r
2834         private void SetZero()\r
2835         {\r
2836             exponent = 0;\r
2837             mantissa.Zero();\r
2838             Sign = false;\r
2839         }\r
2840 \r
2841         private void SetInfPlus()\r
2842         {\r
2843             Sign = false;\r
2844             exponent = Int32.MaxValue;\r
2845             mantissa.Zero();\r
2846         }\r
2847 \r
2848         private void SetInfMinus()\r
2849         {\r
2850             Sign = true;\r
2851             exponent = Int32.MaxValue;\r
2852             mantissa.Zero();\r
2853         }\r
2854 \r
2855         private bool SpecialValueAddTest(BigFloat n2)\r
2856         {\r
2857             if (IsSpecialValue || n2.IsSpecialValue)\r
2858             {\r
2859                 SpecialValueType s1 = SpecialValue;\r
2860                 SpecialValueType s2 = n2.SpecialValue;\r
2861 \r
2862                 if (s1 == SpecialValueType.NAN) return true;\r
2863                 if (s2 == SpecialValueType.NAN)\r
2864                 {\r
2865                     //Set NaN and return.\r
2866                     SetNaN();\r
2867                     return true;\r
2868                 }\r
2869 \r
2870                 if (s1 == SpecialValueType.INF_PLUS)\r
2871                 {\r
2872                     //INF+ + INF- = NAN\r
2873                     if (s2 == SpecialValueType.INF_MINUS)\r
2874                     {\r
2875                         SetNaN();\r
2876                         return true;\r
2877                     }\r
2878 \r
2879                     return true;\r
2880                 }\r
2881 \r
2882                 if (s1 == SpecialValueType.INF_MINUS)\r
2883                 {\r
2884                     //INF+ + INF- = NAN\r
2885                     if (s2 == SpecialValueType.INF_PLUS)\r
2886                     {\r
2887                         SetNaN();\r
2888                         return true;\r
2889                     }\r
2890 \r
2891                     return true;\r
2892                 }\r
2893 \r
2894                 if (s2 == SpecialValueType.ZERO)\r
2895                 {\r
2896                     return true;\r
2897                 }\r
2898 \r
2899                 if (s1 == SpecialValueType.ZERO)\r
2900                 {\r
2901                     Assign(n2);\r
2902                     return true;\r
2903                 }\r
2904             }\r
2905 \r
2906             return false;\r
2907         }\r
2908 \r
2909         private bool SpecialValueMulTest(BigFloat n2)\r
2910         {\r
2911             if (IsSpecialValue || n2.IsSpecialValue)\r
2912             {\r
2913                 SpecialValueType s1 = SpecialValue;\r
2914                 SpecialValueType s2 = n2.SpecialValue;\r
2915 \r
2916                 if (s1 == SpecialValueType.NAN) return true;\r
2917                 if (s2 == SpecialValueType.NAN)\r
2918                 {\r
2919                     //Set NaN and return.\r
2920                     SetNaN();\r
2921                     return true;\r
2922                 }\r
2923 \r
2924                 if (s1 == SpecialValueType.INF_PLUS)\r
2925                 {\r
2926                     //Inf+ * Inf- = Inf-\r
2927                     if (s2 == SpecialValueType.INF_MINUS)\r
2928                     {\r
2929                         Assign(n2);\r
2930                         return true;\r
2931                     }\r
2932 \r
2933                     //Inf+ * 0 = NaN\r
2934                     if (s2 == SpecialValueType.ZERO)\r
2935                     {\r
2936                         //Set NaN and return.\r
2937                         SetNaN();\r
2938                         return true;\r
2939                     }\r
2940 \r
2941                     return true;\r
2942                 }\r
2943 \r
2944                 if (s1 == SpecialValueType.INF_MINUS)\r
2945                 {\r
2946                     //Inf- * Inf- = Inf+\r
2947                     if (s2 == SpecialValueType.INF_MINUS)\r
2948                     {\r
2949                         Sign = false;\r
2950                         return true;\r
2951                     }\r
2952 \r
2953                     //Inf- * 0 = NaN\r
2954                     if (s2 == SpecialValueType.ZERO)\r
2955                     {\r
2956                         //Set NaN and return.\r
2957                         SetNaN();\r
2958                         return true;\r
2959                     }\r
2960 \r
2961                     return true;\r
2962                 }\r
2963 \r
2964                 if (s2 == SpecialValueType.ZERO)\r
2965                 {\r
2966                     SetZero();\r
2967                     return true;\r
2968                 }\r
2969 \r
2970                 if (s1 == SpecialValueType.ZERO)\r
2971                 {\r
2972                     return true;\r
2973                 }\r
2974             }\r
2975 \r
2976             return false;\r
2977         }\r
2978 \r
2979         private bool SpecialValueDivTest(BigFloat n2)\r
2980         {\r
2981             if (IsSpecialValue || n2.IsSpecialValue)\r
2982             {\r
2983                 SpecialValueType s1 = SpecialValue;\r
2984                 SpecialValueType s2 = n2.SpecialValue;\r
2985 \r
2986                 if (s1 == SpecialValueType.NAN) return true;\r
2987                 if (s2 == SpecialValueType.NAN)\r
2988                 {\r
2989                     //Set NaN and return.\r
2990                     SetNaN();\r
2991                     return true;\r
2992                 }\r
2993 \r
2994                 if ((s1 == SpecialValueType.INF_PLUS || s1 == SpecialValueType.INF_MINUS))\r
2995                 {\r
2996                     if (s2 == SpecialValueType.INF_PLUS || s2 == SpecialValueType.INF_MINUS)\r
2997                     {\r
2998                         //Set NaN and return.\r
2999                         SetNaN();\r
3000                         return true;\r
3001                     }\r
3002 \r
3003                     if (n2.Sign)\r
3004                     {\r
3005                         if (s1 == SpecialValueType.INF_PLUS)\r
3006                         {\r
3007                             SetInfMinus();\r
3008                             return true;\r
3009                         }\r
3010 \r
3011                         SetInfPlus();\r
3012                         return true;\r
3013                     }\r
3014 \r
3015                     //Keep inf\r
3016                     return true;\r
3017                 }\r
3018 \r
3019                 if (s2 == SpecialValueType.ZERO)\r
3020                 {\r
3021                     if (s1 == SpecialValueType.ZERO)\r
3022                     {\r
3023                         SetNaN();\r
3024                         return true;\r
3025                     }\r
3026 \r
3027                     if (Sign)\r
3028                     {\r
3029                         SetInfMinus();\r
3030                         return true;\r
3031                     }\r
3032 \r
3033                     SetInfPlus();\r
3034                     return true;\r
3035                 }\r
3036             }\r
3037 \r
3038             return false;\r
3039         }\r
3040 \r
3041         //****************** Internal helper functions *****************\r
3042 \r
3043         /// <summary>\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
3046         /// </summary>\r
3047         /// <param name="targetExponent"></param>\r
3048         private void Denormalise(int targetExponent)\r
3049         {\r
3050             int diff = targetExponent - exponent;\r
3051             if (diff <= 0) return;\r
3052 \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
3055             exponent += diff;\r
3056         }\r
3057 \r
3058         /// <summary>\r
3059         /// The binary logarithm, log2(x) - for precisions above 1000 bits, use Log() and convert the base.\r
3060         /// </summary>\r
3061         private void Log2()\r
3062         {\r
3063             if (scratch.Precision.NumBits != mantissa.Precision.NumBits)\r
3064             {\r
3065                 scratch = new BigInt(mantissa.Precision);\r
3066             }\r
3067 \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
3074 \r
3075             while (bitsCalculated < bits)\r
3076             {\r
3077                 int i;\r
3078                 for (i = 0; (temp.exponent == 0); i++)\r
3079                 {\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
3084                 }\r
3085 \r
3086                 pow2.MulPow2(-i);\r
3087                 result.Add(pow2);\r
3088                 temp.exponent = 0;\r
3089                 bitsCalculated += i;\r
3090             }\r
3091 \r
3092             this.Assign(result);\r
3093         }\r
3094 \r
3095         /// <summary>\r
3096         /// Tried the newton method for logs, but the exponential function is too slow to do it.\r
3097         /// </summary>\r
3098         private void LogNewton()\r
3099         {\r
3100             if (mantissa.IsZero() || mantissa.Sign)\r
3101             {\r
3102                 return;\r
3103             }\r
3104 \r
3105             //Compute ln2.\r
3106             if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)\r
3107             {\r
3108                 CalculateLog2(mantissa.Precision.NumBits);\r
3109             }\r
3110 \r
3111             int numBits = mantissa.Precision.NumBits;\r
3112 \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
3116             xn.exponent = 0;\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
3121             xn.Mul(ln2cache);\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
3126 \r
3127             int precision = 32;\r
3128             int normalPrecision = mantissa.Precision.NumBits;\r
3129 \r
3130             int iterations = 0;\r
3131 \r
3132             while (true)\r
3133             {\r
3134                 term.Assign(xn);\r
3135                 term.mantissa.Sign = true;\r
3136                 term.Exp(precision);\r
3137                 term.Mul(this);\r
3138                 term.Sub(one);\r
3139 \r
3140                 iterations++;\r
3141                 if (term.exponent < -((precision >> 1) - 4))\r
3142                 {\r
3143                     if (precision == normalPrecision)\r
3144                     {\r
3145                         if (term.exponent < -(precision - 4)) break;\r
3146                     }\r
3147                     else\r
3148                     {\r
3149                         precision = precision << 1;\r
3150                         if (precision > normalPrecision) precision = normalPrecision;\r
3151                     }\r
3152                 }\r
3153 \r
3154                 xn.Add(term);\r
3155             }\r
3156 \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
3160 \r
3161             this.Assign(xn);\r
3162             this.Add(term);\r
3163         }\r
3164 \r
3165         /// <summary>\r
3166         /// Log(x) implemented as an Arithmetic-Geometric Mean. Fast for high precisions.\r
3167         /// </summary>\r
3168         private void LogAGM1()\r
3169         {\r
3170             if (mantissa.IsZero() || mantissa.Sign)\r
3171             {\r
3172                 return;\r
3173             }\r
3174 \r
3175             //Compute ln2.\r
3176             if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)\r
3177             {\r
3178                 CalculateLog2(mantissa.Precision.NumBits);\r
3179             }\r
3180 \r
3181             //Compute ln(x) using AGM formula\r
3182 \r
3183             //1. Re-write the input as 2^n * (0.5 <= x < 1)\r
3184             int power2 = exponent + 1;\r
3185             exponent = -1;\r
3186 \r
3187             //BigFloat res = new BigFloat(firstAGMcache);\r
3188             BigFloat a0 = new BigFloat(1, mantissa.Precision);\r
3189             BigFloat b0 = new BigFloat(pow10cache);\r
3190             b0.Mul(this);\r
3191 \r
3192             BigFloat r = R(a0, b0);\r
3193 \r
3194             this.Assign(firstAGMcache);\r
3195             this.Sub(r);\r
3196 \r
3197             a0.Assign(ln2cache);\r
3198             a0.Mul(new BigFloat(power2, mantissa.Precision));\r
3199             this.Add(a0);\r
3200         }\r
3201 \r
3202         private void Exp(int numBits)\r
3203         {\r
3204             if (IsSpecialValue)\r
3205             {\r
3206                 if (SpecialValue == SpecialValueType.ZERO)\r
3207                 {\r
3208                     //e^0 = 1\r
3209                     exponent = 0;\r
3210                     mantissa.SetHighDigit(0x80000000);\r
3211                 }\r
3212                 else if (SpecialValue == SpecialValueType.INF_MINUS)\r
3213                 {\r
3214                     //e^-inf = 0\r
3215                     SetZero();\r
3216                 }\r
3217 \r
3218                 return;\r
3219             }\r
3220 \r
3221             PrecisionSpec prec = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);\r
3222             numBits = prec.NumBits;\r
3223 \r
3224             if (scratch.Precision.NumBits != prec.NumBits)\r
3225             {\r
3226                 scratch = new BigInt(prec);\r
3227             }\r
3228 \r
3229             if (inverseFactorialCache == null || invFactorialCutoff < numBits)\r
3230             {\r
3231                 CalculateFactorials(numBits);\r
3232             }\r
3233 \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
3236 \r
3237             int oldExponent = 0;\r
3238 \r
3239             if (exponent > -4)\r
3240             {\r
3241                 oldExponent = exponent + 4;\r
3242                 exponent = -4;\r
3243             }\r
3244 \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
3250 \r
3251             int iterations;\r
3252             for (int i = 1; i < length; i++)\r
3253             {\r
3254                 //temp = x^i\r
3255                 temp.Mul(thisSave);\r
3256                 temp2.Assign(inverseFactorialCache[i]);\r
3257                 temp2.Mul(temp);\r
3258 \r
3259                 if (temp2.exponent < -(numBits + 4)) { iterations = i; break; }\r
3260 \r
3261                 res.Add(temp2);\r
3262             }\r
3263 \r
3264             //res = exp(x)\r
3265             //Now... x^(2^n) = (x^2)^(2^(n - 1))\r
3266             for (int i = 0; i < oldExponent; i++)\r
3267             {\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
3272             }\r
3273 \r
3274             //Deal with +/- inf\r
3275             if (res.exponent == Int32.MaxValue)\r
3276             {\r
3277                 res.mantissa.Zero();\r
3278             }\r
3279 \r
3280             Assign(res);\r
3281         }\r
3282 \r
3283         /// <summary>\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
3286         /// </summary>\r
3287         /// <param name="numBits"></param>\r
3288         /// <returns></returns>\r
3289         private static void CalculateLog2(int numBits)\r
3290         {\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
3300 \r
3301             //save power of 10 (in normal precision)\r
3302             pow10cache = new BigFloat(b0, normalPres);\r
3303 \r
3304             ln2cache = R(a0, b0);\r
3305 \r
3306             //save the first half of the log calculation\r
3307             firstAGMcache = new BigFloat(ln2cache, normalPres);\r
3308             firstAGMcacheSaved.Assign(ln2cache);\r
3309 \r
3310             b0.MulPow2(-1);\r
3311             ln2cache.Sub(R(a0, b0));\r
3312 \r
3313             //Convert to log(2)\r
3314             ln2cache.mantissa.Sign = false;\r
3315 \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
3323 \r
3324             //Save the inverse.\r
3325             log2ecache = new BigFloat(ln2cache);\r
3326             log2ecache = new BigFloat(log2ecache.Reciprocal(), normalPres);\r
3327 \r
3328             //Now cache log10\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
3333 \r
3334             {\r
3335                 int power2 = log10recip.exponent + 1;\r
3336                 log10recip.exponent = -1;\r
3337 \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
3342 \r
3343                 BigFloat r = R(ax, bx);\r
3344 \r
3345                 log10recip.Assign(firstAGMcacheSaved);\r
3346                 log10recip.Sub(r);\r
3347 \r
3348                 ax.Assign(ln2cache);\r
3349                 ax.Mul(new BigFloat(power2, log10recip.mantissa.Precision));\r
3350                 log10recip.Add(ax);\r
3351             }\r
3352 \r
3353             log10recip = log10recip.Reciprocal();\r
3354             log10recip = new BigFloat(log10recip, normalPres);\r
3355 \r
3356 \r
3357             //Trim to n bits\r
3358             ln2cache = new BigFloat(ln2cache, normalPres);\r
3359         }\r
3360 \r
3361         private static BigFloat TenPow(int power, PrecisionSpec precision)\r
3362         {\r
3363             BigFloat acc = new BigFloat(1, precision);\r
3364             BigFloat temp = new BigFloat(1, precision);\r
3365 \r
3366             int powerTemp = power;\r
3367 \r
3368             BigFloat multiplierToUse = new BigFloat(10, precision);\r
3369 \r
3370             if (power < 0)\r
3371             {\r
3372                 multiplierToUse = multiplierToUse.Reciprocal();\r
3373                 powerTemp = -power;\r
3374             }\r
3375 \r
3376             //Fast power function\r
3377             while (powerTemp != 0)\r
3378             {\r
3379                 temp.Mul(multiplierToUse);\r
3380                 multiplierToUse.Assign(temp);\r
3381 \r
3382                 if ((powerTemp & 1) != 0)\r
3383                 {\r
3384                     acc.Mul(temp);\r
3385                 }\r
3386 \r
3387                 powerTemp >>= 1;\r
3388             }\r
3389 \r
3390             return acc;\r
3391         }\r
3392 \r
3393         private static BigFloat R(BigFloat a0, BigFloat b0)\r
3394         {\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
3404 \r
3405             int iteration = 0;\r
3406 \r
3407             for (iteration = 0; ; iteration++)\r
3408             {\r
3409                 //Get the sum term for this iteration.\r
3410                 term.Assign(an);\r
3411                 term.Mul(an);\r
3412                 temp1.Assign(bn);\r
3413                 temp1.Mul(bn);\r
3414                 //term = an^2 - bn^2\r
3415                 term.Sub(temp1);\r
3416                 //term = 2^(n-1) * (an^2 - bn^2)\r
3417                 term.exponent += iteration - 1;\r
3418                 sum.Add(term);\r
3419 \r
3420                 if (term.exponent < -(bits - 8)) break;\r
3421 \r
3422                 //Calculate the new AGM estimates.\r
3423                 temp1.Assign(an);\r
3424                 an.Add(bn);\r
3425                 //a(n+1) = (an + bn) / 2\r
3426                 an.MulPow2(-1);\r
3427 \r
3428                 //b(n+1) = sqrt(an*bn)\r
3429                 bn.Mul(temp1);\r
3430                 bn.Sqrt();\r
3431             }\r
3432 \r
3433             one.Sub(sum);\r
3434             one = one.Reciprocal();\r
3435             return new BigFloat(one, a0.mantissa.Precision);\r
3436         }\r
3437 \r
3438         private static void CalculateFactorials(int numBits)\r
3439         {\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
3442 \r
3443             PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);\r
3444             PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);\r
3445 \r
3446             BigFloat factorial = new BigFloat(1, extendedPrecision);\r
3447             BigFloat reciprocal;\r
3448 \r
3449             //Calculate e while we're at it\r
3450             BigFloat e = new BigFloat(1, extendedPrecision);\r
3451 \r
3452             list.Add(new BigFloat(factorial, normalPrecision));\r
3453 \r
3454             for (int i = 1; i < Int32.MaxValue; i++)\r
3455             {\r
3456                 BigFloat number = new BigFloat(i, extendedPrecision);\r
3457                 factorial.Mul(number);\r
3458 \r
3459                 if (factorial.exponent > numBits) break;\r
3460 \r
3461                 list2.Add(new BigFloat(factorial, normalPrecision));\r
3462                 reciprocal = factorial.Reciprocal();\r
3463 \r
3464                 e.Add(reciprocal);\r
3465                 list.Add(new BigFloat(reciprocal, normalPrecision));\r
3466             }\r
3467 \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
3474         }\r
3475 \r
3476         private static void CalculateEOnly(int numBits)\r
3477         {\r
3478             PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);\r
3479             PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);\r
3480 \r
3481             int iExponent = (int)(Math.Sqrt(numBits));\r
3482 \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
3488 \r
3489             //Calculate the 2^iExponent th root of e\r
3490             BigFloat e = new BigFloat(1, extendedPrecision);\r
3491 \r
3492             int i;\r
3493             for (i = 1; i < Int32.MaxValue; i++)\r
3494             {\r
3495                 BigFloat number = new BigFloat(i, extendedPrecision);\r
3496                 factorial.Mul(number);\r
3497                 reciprocal = factorial.Reciprocal();\r
3498                 reciprocal.Mul(numerator);\r
3499 \r
3500                 if (-reciprocal.exponent > numBits) break;\r
3501 \r
3502                 e.Add(reciprocal);\r
3503                 numerator.Mul(constant);\r
3504                 System.GC.Collect();\r
3505             }\r
3506 \r
3507             for (i = 0; i < iExponent; i++)\r
3508             {\r
3509                 numerator.Assign(e);\r
3510                 e.Mul(numerator);\r
3511             }\r
3512 \r
3513             //Set the cached static values.\r
3514             eCache = new BigFloat(e, normalPrecision);\r
3515             eRCPCache = new BigFloat(e.Reciprocal(), normalPrecision);\r
3516         }\r
3517 \r
3518         /// <summary>\r
3519         /// Uses the Gauss-Legendre formula for pi\r
3520         /// Taken from http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm\r
3521         /// </summary>\r
3522         /// <param name="numBits"></param>\r
3523         private static void CalculatePi(int numBits)\r
3524         {\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
3529 \r
3530             if (scratch.Precision.NumBits != bits)\r
3531             {\r
3532                 scratch = new BigInt(extendedPres);\r
3533             }\r
3534 \r
3535             //a0 = 1\r
3536             BigFloat an = new BigFloat(1, extendedPres);\r
3537 \r
3538             //b0 = 1/sqrt(2)\r
3539             BigFloat bn = new BigFloat(2, extendedPres);\r
3540             bn.Sqrt();\r
3541             bn.exponent--;\r
3542 \r
3543             //to = 1/4\r
3544             BigFloat tn = new BigFloat(1, extendedPres);\r
3545             tn.exponent -= 2;\r
3546 \r
3547             int pn = 0;\r
3548 \r
3549             BigFloat anTemp = new BigFloat(extendedPres);\r
3550 \r
3551             int iteration = 0;\r
3552             int cutoffBits = numBits >> 5;\r
3553 \r
3554             for (iteration = 0; ; iteration++)\r
3555             {\r
3556                 //Save a(n)\r
3557                 anTemp.Assign(an);\r
3558 \r
3559                 //Calculate new an\r
3560                 an.Add(bn);\r
3561                 an.exponent--;\r
3562 \r
3563                 //Calculate new bn\r
3564                 bn.Mul(anTemp);\r
3565                 bn.Sqrt();\r
3566 \r
3567                 //Calculate new tn\r
3568                 anTemp.Sub(an);\r
3569                 anTemp.mantissa.SquareHiFast(scratch);\r
3570                 anTemp.exponent += anTemp.exponent + pn + 1 - anTemp.mantissa.Normalise();\r
3571                 tn.Sub(anTemp);\r
3572 \r
3573                 anTemp.Assign(an);\r
3574                 anTemp.Sub(bn);\r
3575 \r
3576                 if (anTemp.exponent < -(bits - cutoffBits)) break;\r
3577 \r
3578                 //New pn\r
3579                 pn++;\r
3580             }\r
3581 \r
3582             an.Add(bn);\r
3583             an.mantissa.SquareHiFast(scratch);\r
3584             an.exponent += an.exponent + 1 - an.mantissa.Normalise();\r
3585             tn.exponent += 2;\r
3586             an.Div(tn);\r
3587 \r
3588             pi = new BigFloat(an, normalPres);\r
3589             piBy2 = new BigFloat(pi);\r
3590             piBy2.exponent--;\r
3591             twoPi = new BigFloat(pi, normalPres);\r
3592             twoPi.exponent++;\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
3598         }\r
3599 \r
3600         /// <summary>\r
3601         /// Calculates the odd reciprocals of the natural numbers (for atan series)\r
3602         /// </summary>\r
3603         /// <param name="numBits"></param>\r
3604         /// <param name="terms"></param>\r
3605         private static void CalculateReciprocals(int numBits, int terms)\r
3606         {\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
3610 \r
3611             System.Collections.Generic.List<BigFloat> list = new System.Collections.Generic.List<BigFloat>(terms);\r
3612 \r
3613             for (int i = 0; i < terms; i++)\r
3614             {\r
3615                 BigFloat term = new BigFloat(i*2 + 1, extendedPres);\r
3616                 list.Add(new BigFloat(term.Reciprocal(), normalPres));\r
3617             }\r
3618 \r
3619             reciprocals = list.ToArray();\r
3620         }\r
3621 \r
3622         /// <summary>\r
3623         /// Does decimal rounding, for numbers without E notation.\r
3624         /// </summary>\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
3629         {\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
3633 \r
3634             /*\r
3635             for (int i = 1; i <= places; i++)\r
3636             {\r
3637                 //Skip decimal points.\r
3638                 if (trim[trim.Length - i] == '.')\r
3639                 {\r
3640                     places++;\r
3641                     continue;\r
3642                 }\r
3643 \r
3644                 int index = Array.IndexOf(digits, trim[trim.Length - i]);\r
3645 \r
3646                 if (index < 0) return input;\r
3647 \r
3648                 value += ten * index;\r
3649                 ten *= 10;\r
3650             }\r
3651              * */\r
3652 \r
3653             //Look for a decimal point\r
3654             string decimalPoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;\r
3655 \r
3656             int indexPoint = trim.LastIndexOf(decimalPoint);\r
3657             if (indexPoint < 0)\r
3658             {\r
3659                 //We can't modify a string which doesn't have a decimal point.\r
3660                 return trim;\r
3661             }\r
3662 \r
3663             int trimPoint = trim.Length - places;\r
3664             if (trimPoint < indexPoint) trimPoint = indexPoint;\r
3665 \r
3666             bool roundDown = false;\r
3667 \r
3668             if (trim[trimPoint] == '.')\r
3669             {\r
3670                 if (trimPoint + 1 >= trim.Length)\r
3671                 {\r
3672                     roundDown = true;\r
3673                 }\r
3674                 else\r
3675                 {\r
3676                     int digit = Array.IndexOf(digits, trim[trimPoint + 1]);\r
3677                     if (digit < 5) roundDown = true;\r
3678                 }\r
3679             }\r
3680             else\r
3681             {\r
3682                 int digit = Array.IndexOf(digits, trim[trimPoint]);\r
3683                 if (digit < 5) roundDown = true;\r
3684             }\r
3685 \r
3686             string output;\r
3687 \r
3688             //Round down - just return a new string without the extra digits.\r
3689             if (roundDown)\r
3690             {\r
3691                 if (RoundingMode == RoundingModeType.EXACT)\r
3692                 {\r
3693                     return trim.Substring(0, trimPoint);\r
3694                 }\r
3695                 else\r
3696                 {\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
3701                 }\r
3702             }\r
3703             \r
3704             //Round up - bit more complicated.\r
3705             char [] arrayOutput = trim.ToCharArray();//0, trimPoint);\r
3706 \r
3707             //Now, we round going from the back to the front.\r
3708             int j;\r
3709             for (j = trimPoint - 1; j >= 0; j--)\r
3710             {\r
3711                 int index = Array.IndexOf(digits, arrayOutput[j]);\r
3712 \r
3713                 //Skip decimal points etc...\r
3714                 if (index < 0) continue;\r
3715 \r
3716                 if (index < 9)\r
3717                 {\r
3718                     arrayOutput[j] = digits[index + 1];\r
3719                     break;\r
3720                 }\r
3721                 else\r
3722                 {\r
3723                     arrayOutput[j] = digits[0];\r
3724                 }\r
3725             }\r
3726 \r
3727             output = new string(arrayOutput);\r
3728 \r
3729             if (j < 0)\r
3730             {\r
3731                 //Need to add a new digit.\r
3732                 output = String.Format("{0}{1}", "1", output);\r
3733             }\r
3734 \r
3735             if (RoundingMode == RoundingModeType.EXACT)\r
3736             {\r
3737                 return output;\r
3738             }\r
3739             else\r
3740             {\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
3745             }\r
3746         }\r
3747 \r
3748         //***************************** Data *****************************\r
3749 \r
3750 \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
3754 \r
3755         /// <summary>\r
3756         /// Storage area for calculations.\r
3757         /// </summary>\r
3758         private static BigInt scratch;\r
3759 \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
3778         \r
3779         /// <summary>\r
3780         /// The number of decimal digits to round the output of ToString() by\r
3781         /// </summary>\r
3782         public static int RoundingDigits { get; set; }\r
3783 \r
3784         /// <summary>\r
3785         /// The way in which ToString() should deal with insignificant trailing zeroes\r
3786         /// </summary>\r
3787         public static RoundingModeType RoundingMode { get; set; }\r
3788     }\r
3789 }