• R/O
  • R/O (HTTP)
  • R/W (SSH)
  • R/W (HTTPS)

silverlight:

psychlops silverlight


File Info

Rev. 08bcb090f6a3e6ee38c712571f24b57a4dd32fb4
Size 127,679 bytes
Time 2012-03-29 01:46:02
Author HOSOKAWA Kenchi
Log Message

5

Content

  1. // http://www.fractal-landscapes.co.uk/bigint.html
  2. using System;
  3. namespace BigNum
  4. {
  5. /// <summary>
  6. /// An arbitrary-precision floating-point class
  7. ///
  8. /// Format:
  9. /// Each number is stored as an exponent (32-bit signed integer), and a mantissa
  10. /// (n-bit) BigInteger. The sign of the number is stored in the BigInteger
  11. ///
  12. /// Applicability and Performance:
  13. /// This class is designed to be used for small extended precisions. It may not be
  14. /// safe (and certainly won't be fast) to use it with mixed-precision arguments.
  15. /// It does support, but will not be efficient for, numbers over around 2048 bits.
  16. ///
  17. /// Notes:
  18. /// All conversions to and from strings are slow.
  19. ///
  20. /// Conversions from simple integer types Int32, Int64, UInt32, UInt64 are performed
  21. /// using the appropriate constructor, and are relatively fast.
  22. ///
  23. /// The class is written entirely in managed C# code, with not native or managed
  24. /// assembler. The use of native assembler would speed up the multiplication operations
  25. /// many times over, and therefore all higher-order operations too.
  26. /// </summary>
  27. public class BigFloat
  28. {
  29. /// <summary>
  30. /// Floats can have 4 special value types:
  31. ///
  32. /// NaN: Not a number (cannot be changed using any operations)
  33. /// Infinity: Positive infinity. Some operations e.g. Arctan() allow this input.
  34. /// -Infinity: Negative infinity. Some operations allow this input.
  35. /// Zero
  36. /// </summary>
  37. public enum SpecialValueType
  38. {
  39. /// <summary>
  40. /// Not a special value
  41. /// </summary>
  42. NONE = 0,
  43. /// <summary>
  44. /// Zero
  45. /// </summary>
  46. ZERO,
  47. /// <summary>
  48. /// Positive infinity
  49. /// </summary>
  50. INF_PLUS,
  51. /// <summary>
  52. /// Negative infinity
  53. /// </summary>
  54. INF_MINUS,
  55. /// <summary>
  56. /// Not a number
  57. /// </summary>
  58. NAN
  59. }
  60. /// <summary>
  61. /// This affects the ToString() method.
  62. ///
  63. /// With Trim rounding, all insignificant zero digits are drip
  64. /// </summary>
  65. public enum RoundingModeType
  66. {
  67. /// <summary>
  68. /// Trim non-significant zeros from ToString output after rounding
  69. /// </summary>
  70. TRIM,
  71. /// <summary>
  72. /// Keep all non-significant zeroes in ToString output after rounding
  73. /// </summary>
  74. EXACT
  75. }
  76. /// <summary>
  77. /// A wrapper for the signed exponent, avoiding overflow.
  78. /// </summary>
  79. protected struct ExponentAdaptor
  80. {
  81. /// <summary>
  82. /// The 32-bit exponent
  83. /// </summary>
  84. public Int32 exponent
  85. {
  86. get { return expValue; }
  87. set { expValue = value; }
  88. }
  89. /// <summary>
  90. /// Implicit cast to Int32
  91. /// </summary>
  92. public static implicit operator Int32(ExponentAdaptor adaptor)
  93. {
  94. return adaptor.expValue;
  95. }
  96. /// <summary>
  97. /// Implicit cast from Int32 to ExponentAdaptor
  98. /// </summary>
  99. /// <param name="value"></param>
  100. /// <returns></returns>
  101. public static implicit operator ExponentAdaptor(Int32 value)
  102. {
  103. ExponentAdaptor adaptor = new ExponentAdaptor();
  104. adaptor.expValue = value;
  105. return adaptor;
  106. }
  107. /// <summary>
  108. /// Overloaded increment operator
  109. /// </summary>
  110. public static ExponentAdaptor operator ++(ExponentAdaptor adaptor)
  111. {
  112. adaptor = adaptor + 1;
  113. return adaptor;
  114. }
  115. /// <summary>
  116. /// Overloaded decrement operator
  117. /// </summary>
  118. public static ExponentAdaptor operator --(ExponentAdaptor adaptor)
  119. {
  120. adaptor = adaptor - 1;
  121. return adaptor;
  122. }
  123. /// <summary>
  124. /// Overloaded addition operator
  125. /// </summary>
  126. public static ExponentAdaptor operator +(ExponentAdaptor a1, ExponentAdaptor a2)
  127. {
  128. if (a1.expValue == Int32.MaxValue) return a1;
  129. Int64 temp = (Int64)a1.expValue;
  130. temp += (Int64)(a2.expValue);
  131. if (temp > (Int64)Int32.MaxValue)
  132. {
  133. a1.expValue = Int32.MaxValue;
  134. }
  135. else if (temp < (Int64)Int32.MinValue)
  136. {
  137. a1.expValue = Int32.MinValue;
  138. }
  139. else
  140. {
  141. a1.expValue = (Int32)temp;
  142. }
  143. return a1;
  144. }
  145. /// <summary>
  146. /// Overloaded subtraction operator
  147. /// </summary>
  148. public static ExponentAdaptor operator -(ExponentAdaptor a1, ExponentAdaptor a2)
  149. {
  150. if (a1.expValue == Int32.MaxValue) return a1;
  151. Int64 temp = (Int64)a1.expValue;
  152. temp -= (Int64)(a2.expValue);
  153. if (temp > (Int64)Int32.MaxValue)
  154. {
  155. a1.expValue = Int32.MaxValue;
  156. }
  157. else if (temp < (Int64)Int32.MinValue)
  158. {
  159. a1.expValue = Int32.MinValue;
  160. }
  161. else
  162. {
  163. a1.expValue = (Int32)temp;
  164. }
  165. return a1;
  166. }
  167. /// <summary>
  168. /// Overloaded multiplication operator
  169. /// </summary>
  170. public static ExponentAdaptor operator *(ExponentAdaptor a1, ExponentAdaptor a2)
  171. {
  172. if (a1.expValue == Int32.MaxValue) return a1;
  173. Int64 temp = (Int64)a1.expValue;
  174. temp *= (Int64)a2.expValue;
  175. if (temp > (Int64)Int32.MaxValue)
  176. {
  177. a1.expValue = Int32.MaxValue;
  178. }
  179. else if (temp < (Int64)Int32.MinValue)
  180. {
  181. a1.expValue = Int32.MinValue;
  182. }
  183. else
  184. {
  185. a1.expValue = (Int32)temp;
  186. }
  187. return a1;
  188. }
  189. /// <summary>
  190. /// Overloaded division operator
  191. /// </summary>
  192. public static ExponentAdaptor operator /(ExponentAdaptor a1, ExponentAdaptor a2)
  193. {
  194. if (a1.expValue == Int32.MaxValue) return a1;
  195. ExponentAdaptor res = new ExponentAdaptor();
  196. res.expValue = a1.expValue / a2.expValue;
  197. return res;
  198. }
  199. /// <summary>
  200. /// Overloaded right-shift operator
  201. /// </summary>
  202. public static ExponentAdaptor operator >>(ExponentAdaptor a1, int shift)
  203. {
  204. if (a1.expValue == Int32.MaxValue) return a1;
  205. ExponentAdaptor res = new ExponentAdaptor();
  206. res.expValue = a1.expValue >> shift;
  207. return res;
  208. }
  209. /// <summary>
  210. /// Overloaded left-shift operator
  211. /// </summary>
  212. /// <param name="a1"></param>
  213. /// <param name="shift"></param>
  214. /// <returns></returns>
  215. public static ExponentAdaptor operator <<(ExponentAdaptor a1, int shift)
  216. {
  217. if (a1.expValue == 0) return a1;
  218. ExponentAdaptor res = new ExponentAdaptor();
  219. res.expValue = a1.expValue;
  220. if (shift > 31)
  221. {
  222. res.expValue = Int32.MaxValue;
  223. }
  224. else
  225. {
  226. Int64 temp = a1.expValue;
  227. temp = temp << shift;
  228. if (temp > (Int64)Int32.MaxValue)
  229. {
  230. res.expValue = Int32.MaxValue;
  231. }
  232. else if (temp < (Int64)Int32.MinValue)
  233. {
  234. res.expValue = Int32.MinValue;
  235. }
  236. else
  237. {
  238. res.expValue = (Int32)temp;
  239. }
  240. }
  241. return res;
  242. }
  243. private Int32 expValue;
  244. }
  245. //************************ Constructors **************************
  246. /// <summary>
  247. /// Constructs a 128-bit BigFloat
  248. ///
  249. /// Sets the value to zero
  250. /// </summary>
  251. static BigFloat()
  252. {
  253. RoundingDigits = 3;
  254. RoundingMode = RoundingModeType.TRIM;
  255. scratch = new BigInt(new PrecisionSpec(128, PrecisionSpec.BaseType.BIN));
  256. }
  257. /// <summary>
  258. /// Constructs a BigFloat of the required precision
  259. ///
  260. /// Sets the value to zero
  261. /// </summary>
  262. /// <param name="mantissaPrec"></param>
  263. public BigFloat(PrecisionSpec mantissaPrec)
  264. {
  265. Init(mantissaPrec);
  266. }
  267. /// <summary>
  268. /// Constructs a big float from a UInt32 to the required precision
  269. /// </summary>
  270. /// <param name="value"></param>
  271. /// <param name="mantissaPrec"></param>
  272. public BigFloat(UInt32 value, PrecisionSpec mantissaPrec)
  273. {
  274. int mbWords = ((mantissaPrec.NumBits) >> 5);
  275. if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
  276. int newManBits = mbWords << 5;
  277. //For efficiency, we just use a 32-bit exponent
  278. exponent = 0;
  279. mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  280. //scratch = new BigInt(mantissa.Precision);
  281. int bit = BigInt.GetMSB(value);
  282. if (bit == -1) return;
  283. int shift = mantissa.Precision.NumBits - (bit + 1);
  284. mantissa.LSH(shift);
  285. exponent = bit;
  286. }
  287. /// <summary>
  288. /// Constructs a BigFloat from an Int32 to the required precision
  289. /// </summary>
  290. /// <param name="value"></param>
  291. /// <param name="mantissaPrec"></param>
  292. public BigFloat(Int32 value, PrecisionSpec mantissaPrec)
  293. {
  294. int mbWords = ((mantissaPrec.NumBits) >> 5);
  295. if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
  296. int newManBits = mbWords << 5;
  297. //For efficiency, we just use a 32-bit exponent
  298. exponent = 0;
  299. UInt32 uValue;
  300. if (value < 0)
  301. {
  302. if (value == Int32.MinValue)
  303. {
  304. uValue = 0x80000000;
  305. }
  306. else
  307. {
  308. uValue = (UInt32)(-value);
  309. }
  310. }
  311. else
  312. {
  313. uValue = (UInt32)value;
  314. }
  315. mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  316. //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  317. int bit = BigInt.GetMSB(uValue);
  318. if (bit == -1) return;
  319. int shift = mantissa.Precision.NumBits - (bit + 1);
  320. mantissa.LSH(shift);
  321. exponent = bit;
  322. }
  323. /// <summary>
  324. /// Constructs a BigFloat from a 64-bit integer
  325. /// </summary>
  326. /// <param name="value"></param>
  327. /// <param name="mantissaPrec"></param>
  328. public BigFloat(Int64 value, PrecisionSpec mantissaPrec)
  329. {
  330. int mbWords = ((mantissaPrec.NumBits) >> 5);
  331. if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
  332. int newManBits = mbWords << 5;
  333. //For efficiency, we just use a 32-bit exponent
  334. exponent = 0;
  335. UInt64 uValue;
  336. if (value < 0)
  337. {
  338. if (value == Int64.MinValue)
  339. {
  340. uValue = 0x80000000;
  341. }
  342. else
  343. {
  344. uValue = (UInt64)(-value);
  345. }
  346. }
  347. else
  348. {
  349. uValue = (UInt64)value;
  350. }
  351. mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  352. //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  353. int bit = BigInt.GetMSB(uValue);
  354. if (bit == -1) return;
  355. int shift = mantissa.Precision.NumBits - (bit + 1);
  356. if (shift > 0)
  357. {
  358. mantissa.LSH(shift);
  359. }
  360. else
  361. {
  362. mantissa.SetHighDigit((uint)(uValue >> (-shift)));
  363. }
  364. exponent = bit;
  365. }
  366. /// <summary>
  367. /// Constructs a BigFloat from a 64-bit unsigned integer
  368. /// </summary>
  369. /// <param name="value"></param>
  370. /// <param name="mantissaPrec"></param>
  371. public BigFloat(UInt64 value, PrecisionSpec mantissaPrec)
  372. {
  373. int mbWords = ((mantissaPrec.NumBits) >> 5);
  374. if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
  375. int newManBits = mbWords << 5;
  376. //For efficiency, we just use a 32-bit exponent
  377. exponent = 0;
  378. int bit = BigInt.GetMSB(value);
  379. mantissa = new BigInt(value, new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  380. //scratch = new BigInt(mantissa.Precision);
  381. int shift = mantissa.Precision.NumBits - (bit + 1);
  382. if (shift > 0)
  383. {
  384. mantissa.LSH(shift);
  385. }
  386. else
  387. {
  388. mantissa.SetHighDigit((uint)(value >> (-shift)));
  389. }
  390. exponent = bit;
  391. }
  392. /// <summary>
  393. /// Constructs a BigFloat from a BigInt, using the specified precision
  394. /// </summary>
  395. /// <param name="value"></param>
  396. /// <param name="mantissaPrec"></param>
  397. public BigFloat(BigInt value, PrecisionSpec mantissaPrec)
  398. {
  399. if (value.IsZero())
  400. {
  401. Init(mantissaPrec);
  402. SetZero();
  403. return;
  404. }
  405. mantissa = new BigInt(value, mantissaPrec);
  406. exponent = BigInt.GetMSB(value);
  407. mantissa.Normalise();
  408. }
  409. /// <summary>
  410. /// Construct a BigFloat from a double-precision floating point number
  411. /// </summary>
  412. /// <param name="value"></param>
  413. /// <param name="mantissaPrec"></param>
  414. public BigFloat(double value, PrecisionSpec mantissaPrec)
  415. {
  416. if (value == 0.0)
  417. {
  418. Init(mantissaPrec);
  419. return;
  420. }
  421. bool sign = (value < 0) ? true : false;
  422. long bits = BitConverter.DoubleToInt64Bits(value);
  423. // Note that the shift is sign-extended, hence the test against -1 not 1
  424. int valueExponent = (int)((bits >> 52) & 0x7ffL);
  425. long valueMantissa = bits & 0xfffffffffffffL;
  426. //The mantissa is stored with the top bit implied.
  427. valueMantissa = valueMantissa | 0x10000000000000L;
  428. //The exponent is biased by 1023.
  429. exponent = valueExponent - 1023;
  430. //Round the number of bits to the nearest word.
  431. int mbWords = ((mantissaPrec.NumBits) >> 5);
  432. if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
  433. int newManBits = mbWords << 5;
  434. mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  435. //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  436. if (newManBits >= 64)
  437. {
  438. //The mantissa is 53 bits now, so add 11 to put it in the right place.
  439. mantissa.SetHighDigits(valueMantissa << 11);
  440. }
  441. else
  442. {
  443. //To get the top word of the mantissa, shift up by 11 and down by 32 = down by 21
  444. mantissa.SetHighDigit((uint)(valueMantissa >> 21));
  445. }
  446. mantissa.Sign = sign;
  447. }
  448. /// <summary>
  449. /// Copy constructor
  450. /// </summary>
  451. /// <param name="value"></param>
  452. public BigFloat(BigFloat value)
  453. {
  454. Init(value.mantissa.Precision);
  455. exponent = value.exponent;
  456. mantissa.Assign(value.mantissa);
  457. }
  458. /// <summary>
  459. /// Copy Constructor - constructs a new BigFloat with the specified precision, copying the old one.
  460. ///
  461. /// The value is rounded towards zero in the case where precision is decreased. The Round() function
  462. /// should be used beforehand if a correctly rounded result is required.
  463. /// </summary>
  464. /// <param name="value"></param>
  465. /// <param name="mantissaPrec"></param>
  466. public BigFloat(BigFloat value, PrecisionSpec mantissaPrec)
  467. {
  468. Init(mantissaPrec);
  469. exponent = value.exponent;
  470. if (mantissa.AssignHigh(value.mantissa)) exponent++;
  471. }
  472. /// <summary>
  473. /// Constructs a BigFloat from a string
  474. /// </summary>
  475. /// <param name="value"></param>
  476. /// <param name="mantissaPrec"></param>
  477. public BigFloat(string value, PrecisionSpec mantissaPrec)
  478. {
  479. Init(mantissaPrec);
  480. PrecisionSpec extendedPres = new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN);
  481. BigFloat ten = new BigFloat(10, extendedPres);
  482. BigFloat iPart = new BigFloat(extendedPres);
  483. BigFloat fPart = new BigFloat(extendedPres);
  484. BigFloat tenRCP = ten.Reciprocal();
  485. if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol))
  486. {
  487. SetNaN();
  488. return;
  489. }
  490. else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol))
  491. {
  492. SetInfPlus();
  493. return;
  494. }
  495. else if (value.Contains(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol))
  496. {
  497. SetInfMinus();
  498. return;
  499. }
  500. string decimalpoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;
  501. char[] digitChars = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ',', '.' };
  502. //Read in the integer part up the the decimal point.
  503. bool sign = false;
  504. value = value.Trim();
  505. int i = 0;
  506. if (value.Length > i && value[i] == '-')
  507. {
  508. sign = true;
  509. i++;
  510. }
  511. if (value.Length > i && value[i] == '+')
  512. {
  513. i++;
  514. }
  515. for ( ; i < value.Length; i++)
  516. {
  517. //break on decimal point
  518. if (value[i] == decimalpoint[0]) break;
  519. int digit = Array.IndexOf(digitChars, value[i]);
  520. if (digit < 0) break;
  521. //Ignore place separators (assumed either , or .)
  522. if (digit > 9) continue;
  523. if (i > 0) iPart.Mul(ten);
  524. iPart.Add(new BigFloat(digit, extendedPres));
  525. }
  526. //If we've run out of characters, assign everything and return
  527. if (i == value.Length)
  528. {
  529. iPart.mantissa.Sign = sign;
  530. exponent = iPart.exponent;
  531. if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
  532. return;
  533. }
  534. //Assign the characters after the decimal point to fPart
  535. if (value[i] == '.' && i < value.Length - 1)
  536. {
  537. BigFloat RecipToUse = new BigFloat(tenRCP);
  538. for (i++; i < value.Length; i++)
  539. {
  540. int digit = Array.IndexOf(digitChars, value[i]);
  541. if (digit < 0) break;
  542. BigFloat temp = new BigFloat(digit, extendedPres);
  543. temp.Mul(RecipToUse);
  544. RecipToUse.Mul(tenRCP);
  545. fPart.Add(temp);
  546. }
  547. }
  548. //If we're run out of characters, add fPart and iPart and return
  549. if (i == value.Length)
  550. {
  551. iPart.Add(fPart);
  552. iPart.mantissa.Sign = sign;
  553. exponent = iPart.exponent;
  554. if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
  555. return;
  556. }
  557. if (value[i] == '+' || value[i] == '-') i++;
  558. if (i == value.Length)
  559. {
  560. iPart.Add(fPart);
  561. iPart.mantissa.Sign = sign;
  562. exponent = iPart.exponent;
  563. if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
  564. return;
  565. }
  566. //Look for exponential notation.
  567. if ((value[i] == 'e' || value[i] == 'E') && i < value.Length - 1)
  568. {
  569. //Convert the exponent to an int.
  570. int exp;
  571. try
  572. {
  573. exp = System.Convert.ToInt32(new string(value.ToCharArray()));// i + 1, value.Length - (i + 1))));
  574. }
  575. catch (Exception)
  576. {
  577. iPart.Add(fPart);
  578. iPart.mantissa.Sign = sign;
  579. exponent = iPart.exponent;
  580. if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
  581. return;
  582. }
  583. //Raise or lower 10 to the power of the exponent
  584. BigFloat acc = new BigFloat(1, extendedPres);
  585. BigFloat temp = new BigFloat(1, extendedPres);
  586. int powerTemp = exp;
  587. BigFloat multiplierToUse;
  588. if (exp < 0)
  589. {
  590. multiplierToUse = new BigFloat(tenRCP);
  591. powerTemp = -exp;
  592. }
  593. else
  594. {
  595. multiplierToUse = new BigFloat(ten);
  596. }
  597. //Fast power function
  598. while (powerTemp != 0)
  599. {
  600. temp.Mul(multiplierToUse);
  601. multiplierToUse.Assign(temp);
  602. if ((powerTemp & 1) != 0)
  603. {
  604. acc.Mul(temp);
  605. }
  606. powerTemp >>= 1;
  607. }
  608. iPart.Add(fPart);
  609. iPart.Mul(acc);
  610. iPart.mantissa.Sign = sign;
  611. exponent = iPart.exponent;
  612. if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
  613. return;
  614. }
  615. iPart.Add(fPart);
  616. iPart.mantissa.Sign = sign;
  617. exponent = iPart.exponent;
  618. if (mantissa.AssignHigh(iPart.mantissa)) exponent++;
  619. }
  620. private void Init(PrecisionSpec mantissaPrec)
  621. {
  622. int mbWords = ((mantissaPrec.NumBits) >> 5);
  623. if ((mantissaPrec.NumBits & 31) != 0) mbWords++;
  624. int newManBits = mbWords << 5;
  625. //For efficiency, we just use a 32-bit exponent
  626. exponent = 0;
  627. mantissa = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  628. //scratch = new BigInt(new PrecisionSpec(newManBits, PrecisionSpec.BaseType.BIN));
  629. }
  630. //************************** Properties *************************
  631. /// <summary>
  632. /// Read-only property. Returns the precision specification of the mantissa.
  633. ///
  634. /// Floating point numbers are represented as 2^exponent * mantissa, where the
  635. /// mantissa and exponent are integers. Note that the exponent in this class is
  636. /// always a 32-bit integer. The precision therefore specifies how many bits
  637. /// the mantissa will have.
  638. /// </summary>
  639. public PrecisionSpec Precision
  640. {
  641. get { return mantissa.Precision; }
  642. }
  643. /// <summary>
  644. /// Writable property:
  645. /// true iff the number is negative or in some cases zero (&lt;0)
  646. /// false iff the number if positive or in some cases zero (&gt;0)
  647. /// </summary>
  648. public bool Sign
  649. {
  650. get { return mantissa.Sign; }
  651. set { mantissa.Sign = value; }
  652. }
  653. /// <summary>
  654. /// Read-only property.
  655. /// True if the number is NAN, INF_PLUS, INF_MINUS or ZERO
  656. /// False if the number has any other value.
  657. /// </summary>
  658. public bool IsSpecialValue
  659. {
  660. get
  661. {
  662. return (exponent == Int32.MaxValue || mantissa.IsZero());
  663. }
  664. }
  665. /// <summary>
  666. /// Read-only property, returns the type of number this is. Special values include:
  667. ///
  668. /// NONE - a regular number
  669. /// ZERO - zero
  670. /// NAN - Not a Number (some operations will return this if their inputs are out of range)
  671. /// INF_PLUS - Positive infinity, not really a number, but a valid input to and output of some functions.
  672. /// INF_MINUS - Negative infinity, not really a number, but a valid input to and output of some functions.
  673. /// </summary>
  674. public SpecialValueType SpecialValue
  675. {
  676. get
  677. {
  678. if (exponent == Int32.MaxValue)
  679. {
  680. if (mantissa.IsZero())
  681. {
  682. if (mantissa.Sign) return SpecialValueType.INF_MINUS;
  683. return SpecialValueType.INF_PLUS;
  684. }
  685. return SpecialValueType.NAN;
  686. }
  687. else
  688. {
  689. if (mantissa.IsZero()) return SpecialValueType.ZERO;
  690. return SpecialValueType.NONE;
  691. }
  692. }
  693. }
  694. //******************** Mathematical Constants *******************
  695. /// <summary>
  696. /// Gets pi to the indicated precision
  697. /// </summary>
  698. /// <param name="precision">The precision to perform the calculation to</param>
  699. /// <returns>pi (the ratio of the area of a circle to its diameter)</returns>
  700. public static BigFloat GetPi(PrecisionSpec precision)
  701. {
  702. if (pi == null || precision.NumBits <= pi.mantissa.Precision.NumBits)
  703. {
  704. CalculatePi(precision.NumBits);
  705. }
  706. BigFloat ret = new BigFloat (precision);
  707. ret.Assign(pi);
  708. return ret;
  709. }
  710. /// <summary>
  711. /// Get e to the indicated precision
  712. /// </summary>
  713. /// <param name="precision">The preicision to perform the calculation to</param>
  714. /// <returns>e (the number for which the d/dx(e^x) = e^x)</returns>
  715. public static BigFloat GetE(PrecisionSpec precision)
  716. {
  717. if (eCache == null || eCache.mantissa.Precision.NumBits < precision.NumBits)
  718. {
  719. CalculateEOnly(precision.NumBits);
  720. //CalculateFactorials(precision.NumBits);
  721. }
  722. BigFloat ret = new BigFloat(precision);
  723. ret.Assign(eCache);
  724. return ret;
  725. }
  726. //******************** Arithmetic Functions ********************
  727. /// <summary>
  728. /// Addition (this = this + n2)
  729. /// </summary>
  730. /// <param name="n2">The number to add</param>
  731. public void Add(BigFloat n2)
  732. {
  733. if (SpecialValueAddTest(n2)) return;
  734. if (scratch.Precision.NumBits != n2.mantissa.Precision.NumBits)
  735. {
  736. scratch = new BigInt(n2.mantissa.Precision);
  737. }
  738. if (exponent <= n2.exponent)
  739. {
  740. int diff = n2.exponent - exponent;
  741. exponent = n2.exponent;
  742. if (diff != 0)
  743. {
  744. mantissa.RSH(diff);
  745. }
  746. uint carry = mantissa.Add(n2.mantissa);
  747. if (carry != 0)
  748. {
  749. mantissa.RSH(1);
  750. mantissa.SetBit(mantissa.Precision.NumBits - 1);
  751. exponent++;
  752. }
  753. exponent -= mantissa.Normalise();
  754. }
  755. else
  756. {
  757. int diff = exponent - n2.exponent;
  758. scratch.Assign(n2.mantissa);
  759. scratch.RSH(diff);
  760. uint carry = scratch.Add(mantissa);
  761. if (carry != 0)
  762. {
  763. scratch.RSH(1);
  764. scratch.SetBit(mantissa.Precision.NumBits - 1);
  765. exponent++;
  766. }
  767. mantissa.Assign(scratch);
  768. exponent -= mantissa.Normalise();
  769. }
  770. }
  771. /// <summary>
  772. /// Subtraction (this = this - n2)
  773. /// </summary>
  774. /// <param name="n2">The number to subtract from this</param>
  775. public void Sub(BigFloat n2)
  776. {
  777. n2.mantissa.Sign = !n2.mantissa.Sign;
  778. Add(n2);
  779. n2.mantissa.Sign = !n2.mantissa.Sign;
  780. }
  781. /// <summary>
  782. /// Multiplication (this = this * n2)
  783. /// </summary>
  784. /// <param name="n2">The number to multiply this by</param>
  785. public void Mul(BigFloat n2)
  786. {
  787. if (SpecialValueMulTest(n2)) return;
  788. //Anything times 0 = 0
  789. if (n2.mantissa.IsZero())
  790. {
  791. mantissa.Assign(n2.mantissa);
  792. exponent = 0;
  793. return;
  794. }
  795. mantissa.MulHi(n2.mantissa);
  796. int shift = mantissa.Normalise();
  797. exponent = exponent + n2.exponent + 1 - shift;
  798. }
  799. /// <summary>
  800. /// Division (this = this / n2)
  801. /// </summary>
  802. /// <param name="n2">The number to divide this by</param>
  803. public void Div(BigFloat n2)
  804. {
  805. if (SpecialValueDivTest(n2)) return;
  806. if (mantissa.Precision.NumBits >= 8192)
  807. {
  808. BigFloat rcp = n2.Reciprocal();
  809. Mul(rcp);
  810. }
  811. else
  812. {
  813. int shift = mantissa.DivAndShift(n2.mantissa);
  814. exponent = exponent - (n2.exponent + shift);
  815. }
  816. }
  817. /// <summary>
  818. /// Multiply by a power of 2 (-ve implies division)
  819. /// </summary>
  820. /// <param name="pow2"></param>
  821. public void MulPow2(int pow2)
  822. {
  823. exponent += pow2;
  824. }
  825. /// <summary>
  826. /// Division-based reciprocal, fastest for small precisions up to 15,000 bits.
  827. /// </summary>
  828. /// <returns>The reciprocal 1/this</returns>
  829. public BigFloat Reciprocal()
  830. {
  831. if (mantissa.Precision.NumBits >= 8192) return ReciprocalNewton();
  832. BigFloat reciprocal = new BigFloat(1u, mantissa.Precision);
  833. reciprocal.Div(this);
  834. return reciprocal;
  835. }
  836. /// <summary>
  837. /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.
  838. /// </summary>
  839. /// <returns>The reciprocal 1/this</returns>
  840. public BigFloat ReciprocalNewton()
  841. {
  842. if (mantissa.IsZero())
  843. {
  844. exponent = Int32.MaxValue;
  845. return null;
  846. }
  847. bool oldSign = mantissa.Sign;
  848. int oldExponent = exponent;
  849. //Kill exponent for now (will re-institute later)
  850. exponent = 0;
  851. bool topBit = mantissa.IsTopBitOnlyBit();
  852. PrecisionSpec curPrec = new PrecisionSpec(32, PrecisionSpec.BaseType.BIN);
  853. BigFloat reciprocal = new BigFloat(curPrec);
  854. BigFloat constant2 = new BigFloat(curPrec);
  855. BigFloat temp = new BigFloat(curPrec);
  856. BigFloat thisPrec = new BigFloat(this, curPrec);
  857. reciprocal.exponent = 1;
  858. reciprocal.mantissa.SetHighDigit(3129112985u);
  859. constant2.exponent = 1;
  860. constant2.mantissa.SetHighDigit(0x80000000u);
  861. //D is deliberately left negative for all the following operations.
  862. thisPrec.mantissa.Sign = true;
  863. //Initial estimate.
  864. reciprocal.Add(thisPrec);
  865. //mantissa.Sign = false;
  866. //Shift down into 0.5 < this < 1 range
  867. thisPrec.mantissa.RSH(1);
  868. //Iteration.
  869. int accuracyBits = 2;
  870. int mantissaBits = mantissa.Precision.NumBits;
  871. //Each iteration is a pass of newton's method for RCP.
  872. //The is a substantial optimisation to be done here...
  873. //You can double the number of bits for the calculations
  874. //at each iteration, meaning that the whole process only
  875. //takes some constant multiplier of the time for the
  876. //full-scale multiplication.
  877. while (accuracyBits < mantissaBits)
  878. {
  879. //Increase the precision as needed
  880. if (accuracyBits >= curPrec.NumBits / 2)
  881. {
  882. int newBits = curPrec.NumBits * 2;
  883. if (newBits > mantissaBits) newBits = mantissaBits;
  884. curPrec = new PrecisionSpec(newBits, PrecisionSpec.BaseType.BIN);
  885. reciprocal = new BigFloat(reciprocal, curPrec);
  886. constant2 = new BigFloat(curPrec);
  887. constant2.exponent = 1;
  888. constant2.mantissa.SetHighDigit(0x80000000u);
  889. temp = new BigFloat(temp, curPrec);
  890. thisPrec = new BigFloat(this, curPrec);
  891. thisPrec.mantissa.Sign = true;
  892. thisPrec.mantissa.RSH(1);
  893. }
  894. //temp = Xn
  895. temp.exponent = reciprocal.exponent;
  896. temp.mantissa.Assign(reciprocal.mantissa);
  897. //temp = -Xn * D
  898. temp.Mul(thisPrec);
  899. //temp = -Xn * D + 2 (= 2 - Xn * D)
  900. temp.Add(constant2);
  901. //reciprocal = X(n+1) = Xn * (2 - Xn * D)
  902. reciprocal.Mul(temp);
  903. accuracyBits *= 2;
  904. }
  905. //'reciprocal' is now the reciprocal of the shifted down, zero-exponent mantissa of 'this'
  906. //Restore the mantissa.
  907. //mantissa.LSH(1);
  908. exponent = oldExponent;
  909. //mantissa.Sign = oldSign;
  910. if (topBit)
  911. {
  912. reciprocal.exponent = -(oldExponent);
  913. }
  914. else
  915. {
  916. reciprocal.exponent = -(oldExponent + 1);
  917. }
  918. reciprocal.mantissa.Sign = oldSign;
  919. return reciprocal;
  920. }
  921. /// <summary>
  922. /// Newton's method reciprocal, fastest for larger precisions over 15,000 bits.
  923. /// </summary>
  924. /// <returns>The reciprocal 1/this</returns>
  925. private BigFloat ReciprocalNewton2()
  926. {
  927. if (mantissa.IsZero())
  928. {
  929. exponent = Int32.MaxValue;
  930. return null;
  931. }
  932. bool oldSign = mantissa.Sign;
  933. int oldExponent = exponent;
  934. //Kill exponent for now (will re-institute later)
  935. exponent = 0;
  936. BigFloat reciprocal = new BigFloat(mantissa.Precision);
  937. BigFloat constant2 = new BigFloat(mantissa.Precision);
  938. BigFloat temp = new BigFloat(mantissa.Precision);
  939. reciprocal.exponent = 1;
  940. reciprocal.mantissa.SetHighDigit(3129112985u);
  941. constant2.exponent = 1;
  942. constant2.mantissa.SetHighDigit(0x80000000u);
  943. //D is deliberately left negative for all the following operations.
  944. mantissa.Sign = true;
  945. //Initial estimate.
  946. reciprocal.Add(this);
  947. //mantissa.Sign = false;
  948. //Shift down into 0.5 < this < 1 range
  949. mantissa.RSH(1);
  950. //Iteration.
  951. int accuracyBits = 2;
  952. int mantissaBits = mantissa.Precision.NumBits;
  953. //Each iteration is a pass of newton's method for RCP.
  954. //The is a substantial optimisation to be done here...
  955. //You can double the number of bits for the calculations
  956. //at each iteration, meaning that the whole process only
  957. //takes some constant multiplier of the time for the
  958. //full-scale multiplication.
  959. while (accuracyBits < mantissaBits)
  960. {
  961. //temp = Xn
  962. temp.exponent = reciprocal.exponent;
  963. temp.mantissa.Assign(reciprocal.mantissa);
  964. //temp = -Xn * D
  965. temp.Mul(this);
  966. //temp = -Xn * D + 2 (= 2 - Xn * D)
  967. temp.Add(constant2);
  968. //reciprocal = X(n+1) = Xn * (2 - Xn * D)
  969. reciprocal.Mul(temp);
  970. accuracyBits *= 2;
  971. }
  972. //'reciprocal' is now the reciprocal of the shifted down, zero-exponent mantissa of 'this'
  973. //Restore the mantissa.
  974. mantissa.LSH(1);
  975. exponent = oldExponent;
  976. mantissa.Sign = oldSign;
  977. reciprocal.exponent = -(oldExponent + 1);
  978. reciprocal.mantissa.Sign = oldSign;
  979. return reciprocal;
  980. }
  981. /// <summary>
  982. /// Sets this equal to the input
  983. /// </summary>
  984. /// <param name="n2"></param>
  985. public void Assign(BigFloat n2)
  986. {
  987. exponent = n2.exponent;
  988. if (mantissa.AssignHigh(n2.mantissa)) exponent++;
  989. }
  990. //********************* Comparison Functions *******************
  991. /// <summary>
  992. /// Greater than comparison
  993. /// </summary>
  994. /// <param name="n2">the number to compare this to</param>
  995. /// <returns>true iff this is greater than n2 (this &gt; n2)</returns>
  996. public bool GreaterThan(BigFloat n2)
  997. {
  998. if (IsSpecialValue || n2.IsSpecialValue)
  999. {
  1000. SpecialValueType s1 = SpecialValue;
  1001. SpecialValueType s2 = SpecialValue;
  1002. if (s1 == SpecialValueType.NAN || s2 == SpecialValueType.NAN) return false;
  1003. if (s1 == SpecialValueType.INF_MINUS) return false;
  1004. if (s2 == SpecialValueType.INF_PLUS) return false;
  1005. if (s1 == SpecialValueType.INF_PLUS) return true;
  1006. if (s2 == SpecialValueType.INF_MINUS) return true;
  1007. if (s1 == SpecialValueType.ZERO)
  1008. {
  1009. if (s2 != SpecialValueType.ZERO && n2.Sign)
  1010. {
  1011. return true;
  1012. }
  1013. else
  1014. {
  1015. return false;
  1016. }
  1017. }
  1018. if (s2 == SpecialValueType.ZERO)
  1019. {
  1020. return !Sign;
  1021. }
  1022. }
  1023. if (!mantissa.Sign && n2.mantissa.Sign) return true;
  1024. if (mantissa.Sign && !n2.mantissa.Sign) return false;
  1025. if (!mantissa.Sign)
  1026. {
  1027. if (exponent > n2.exponent) return true;
  1028. if (exponent < n2.exponent) return false;
  1029. }
  1030. if (mantissa.Sign)
  1031. {
  1032. if (exponent > n2.exponent) return false;
  1033. if (exponent < n2.exponent) return true;
  1034. }
  1035. return mantissa.GreaterThan(n2.mantissa);
  1036. }
  1037. /// <summary>
  1038. /// Less than comparison
  1039. /// </summary>
  1040. /// <param name="n2">the number to compare this to</param>
  1041. /// <returns>true iff this is less than n2 (this &lt; n2)</returns>
  1042. public bool LessThan(BigFloat n2)
  1043. {
  1044. if (IsSpecialValue || n2.IsSpecialValue)
  1045. {
  1046. SpecialValueType s1 = SpecialValue;
  1047. SpecialValueType s2 = SpecialValue;
  1048. if (s1 == SpecialValueType.NAN || s2 == SpecialValueType.NAN) return false;
  1049. if (s1 == SpecialValueType.INF_PLUS) return false;
  1050. if (s2 == SpecialValueType.INF_PLUS) return true;
  1051. if (s2 == SpecialValueType.INF_MINUS) return false;
  1052. if (s1 == SpecialValueType.INF_MINUS) return true;
  1053. if (s1 == SpecialValueType.ZERO)
  1054. {
  1055. if (s2 != SpecialValueType.ZERO && !n2.Sign)
  1056. {
  1057. return true;
  1058. }
  1059. else
  1060. {
  1061. return false;
  1062. }
  1063. }
  1064. if (s2 == SpecialValueType.ZERO)
  1065. {
  1066. return Sign;
  1067. }
  1068. }
  1069. if (!mantissa.Sign && n2.mantissa.Sign) return false;
  1070. if (mantissa.Sign && !n2.mantissa.Sign) return true;
  1071. if (!mantissa.Sign)
  1072. {
  1073. if (exponent > n2.exponent) return false;
  1074. if (exponent < n2.exponent) return true;
  1075. }
  1076. if (mantissa.Sign)
  1077. {
  1078. if (exponent > n2.exponent) return true;
  1079. if (exponent < n2.exponent) return false;
  1080. }
  1081. return mantissa.LessThan(n2.mantissa);
  1082. }
  1083. /// <summary>
  1084. /// Greater than comparison
  1085. /// </summary>
  1086. /// <param name="i">the number to compare this to</param>
  1087. /// <returns>true iff this is greater than n2 (this &gt; n2)</returns>
  1088. public bool GreaterThan(int i)
  1089. {
  1090. BigFloat integer = new BigFloat(i, mantissa.Precision);
  1091. return GreaterThan(integer);
  1092. }
  1093. /// <summary>
  1094. /// Less than comparison
  1095. /// </summary>
  1096. /// <param name="i">the number to compare this to</param>
  1097. /// <returns>true iff this is less than n2 (this &lt; n2)</returns>
  1098. public bool LessThan(int i)
  1099. {
  1100. BigFloat integer = new BigFloat(i, mantissa.Precision);
  1101. return LessThan(integer);
  1102. }
  1103. /// <summary>
  1104. /// Compare to zero
  1105. /// </summary>
  1106. /// <returns>true if this is zero (this == 0)</returns>
  1107. public bool IsZero()
  1108. {
  1109. return (mantissa.IsZero());
  1110. }
  1111. //******************** Mathematical Functions ******************
  1112. /// <summary>
  1113. /// Sets the number to the biggest integer numerically closer to zero, if possible.
  1114. /// </summary>
  1115. public void Floor()
  1116. {
  1117. //Already an integer.
  1118. if (exponent >= mantissa.Precision.NumBits) return;
  1119. if (exponent < 0)
  1120. {
  1121. mantissa.ZeroBits(mantissa.Precision.NumBits);
  1122. exponent = 0;
  1123. return;
  1124. }
  1125. mantissa.ZeroBits(mantissa.Precision.NumBits - (exponent + 1));
  1126. }
  1127. /// <summary>
  1128. /// Sets the number to its fractional component (equivalent to 'this' - (int)'this')
  1129. /// </summary>
  1130. public void FPart()
  1131. {
  1132. //Already fractional
  1133. if (exponent < 0)
  1134. {
  1135. return;
  1136. }
  1137. //Has no fractional part
  1138. if (exponent >= mantissa.Precision.NumBits)
  1139. {
  1140. mantissa.Zero();
  1141. exponent = 0;
  1142. return;
  1143. }
  1144. mantissa.ZeroBitsHigh(exponent + 1);
  1145. exponent -= mantissa.Normalise();
  1146. }
  1147. /// <summary>
  1148. /// Calculates tan(x)
  1149. /// </summary>
  1150. public void Tan()
  1151. {
  1152. if (IsSpecialValue)
  1153. {
  1154. //Tan(x) has no limit as x->inf
  1155. if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
  1156. {
  1157. SetNaN();
  1158. }
  1159. else if (SpecialValue == SpecialValueType.ZERO)
  1160. {
  1161. SetZero();
  1162. }
  1163. return;
  1164. }
  1165. if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
  1166. {
  1167. CalculatePi(mantissa.Precision.NumBits);
  1168. }
  1169. //Work out the sign change (involves replicating some rescaling).
  1170. bool sign = mantissa.Sign;
  1171. mantissa.Sign = false;
  1172. if (mantissa.IsZero())
  1173. {
  1174. return;
  1175. }
  1176. //Rescale into 0 <= x < pi
  1177. if (GreaterThan(pi))
  1178. {
  1179. //There will be an inherent loss of precision doing this.
  1180. BigFloat newAngle = new BigFloat(this);
  1181. newAngle.Mul(piRecip);
  1182. newAngle.FPart();
  1183. newAngle.Mul(pi);
  1184. Assign(newAngle);
  1185. }
  1186. //Rescale to -pi/2 <= x < pi/2
  1187. if (!LessThan(piBy2))
  1188. {
  1189. Sub(pi);
  1190. }
  1191. //Now the sign of the sin determines the sign of the tan.
  1192. //tan(x) = sin(x) / sqrt(1 - sin^2(x))
  1193. Sin();
  1194. BigFloat denom = new BigFloat(this);
  1195. denom.Mul(this);
  1196. denom.Sub(new BigFloat(1, mantissa.Precision));
  1197. denom.mantissa.Sign = !denom.mantissa.Sign;
  1198. if (denom.mantissa.Sign)
  1199. {
  1200. denom.SetZero();
  1201. }
  1202. denom.Sqrt();
  1203. Div(denom);
  1204. if (sign) mantissa.Sign = !mantissa.Sign;
  1205. }
  1206. /// <summary>
  1207. /// Calculates Cos(x)
  1208. /// </summary>
  1209. public void Cos()
  1210. {
  1211. if (IsSpecialValue)
  1212. {
  1213. //Cos(x) has no limit as x->inf
  1214. if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
  1215. {
  1216. SetNaN();
  1217. }
  1218. else if (SpecialValue == SpecialValueType.ZERO)
  1219. {
  1220. Assign(new BigFloat(1, mantissa.Precision));
  1221. }
  1222. return;
  1223. }
  1224. if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
  1225. {
  1226. CalculatePi(mantissa.Precision.NumBits);
  1227. }
  1228. Add(piBy2);
  1229. Sin();
  1230. }
  1231. /// <summary>
  1232. /// Calculates Sin(x):
  1233. /// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].
  1234. /// </summary>
  1235. public void Sin()
  1236. {
  1237. if (IsSpecialValue)
  1238. {
  1239. //Sin(x) has no limit as x->inf
  1240. if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
  1241. {
  1242. SetNaN();
  1243. }
  1244. return;
  1245. }
  1246. //Convert to positive range (0 <= x < inf)
  1247. bool sign = mantissa.Sign;
  1248. mantissa.Sign = false;
  1249. if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
  1250. {
  1251. CalculatePi(mantissa.Precision.NumBits);
  1252. }
  1253. if (inverseFactorialCache == null || invFactorialCutoff != mantissa.Precision.NumBits)
  1254. {
  1255. CalculateFactorials(mantissa.Precision.NumBits);
  1256. }
  1257. //Rescale into 0 <= x < 2*pi
  1258. if (GreaterThan(twoPi))
  1259. {
  1260. //There will be an inherent loss of precision doing this.
  1261. BigFloat newAngle = new BigFloat(this);
  1262. newAngle.Mul(twoPiRecip);
  1263. newAngle.FPart();
  1264. newAngle.Mul(twoPi);
  1265. Assign(newAngle);
  1266. }
  1267. //Rescale into range 0 <= x < pi
  1268. if (GreaterThan(pi))
  1269. {
  1270. //sin(pi + a) = sin(pi)cos(a) + sin(a)cos(pi) = 0 - sin(a) = -sin(a)
  1271. Sub(pi);
  1272. sign = !sign;
  1273. }
  1274. BigFloat temp = new BigFloat(mantissa.Precision);
  1275. //Rescale into range 0 <= x < pi/2
  1276. if (GreaterThan(piBy2))
  1277. {
  1278. temp.Assign(this);
  1279. Assign(pi);
  1280. Sub(temp);
  1281. }
  1282. //Rescale into range 0 <= x < pi/6 to accelerate convergence.
  1283. //This is done using sin(3x) = 3sin(x) - 4sin^3(x)
  1284. Mul(threeRecip);
  1285. if (mantissa.IsZero())
  1286. {
  1287. exponent = 0;
  1288. return;
  1289. }
  1290. BigFloat term = new BigFloat(this);
  1291. BigFloat square = new BigFloat(this);
  1292. square.Mul(term);
  1293. BigFloat sum = new BigFloat(this);
  1294. bool termSign = true;
  1295. int length = inverseFactorialCache.Length;
  1296. int numBits = mantissa.Precision.NumBits;
  1297. for (int i = 3; i < length; i += 2)
  1298. {
  1299. term.Mul(square);
  1300. temp.Assign(inverseFactorialCache[i]);
  1301. temp.Mul(term);
  1302. temp.mantissa.Sign = termSign;
  1303. termSign = !termSign;
  1304. if (temp.exponent < -numBits) break;
  1305. sum.Add(temp);
  1306. }
  1307. //Restore the triple-angle: sin(3x) = 3sin(x) - 4sin^3(x)
  1308. Assign(sum);
  1309. sum.Mul(this);
  1310. sum.Mul(this);
  1311. Mul(new BigFloat(3, mantissa.Precision));
  1312. sum.exponent += 2;
  1313. Sub(sum);
  1314. //Restore the sign
  1315. mantissa.Sign = sign;
  1316. }
  1317. /// <summary>
  1318. /// Hyperbolic Sin (sinh) function
  1319. /// </summary>
  1320. public void Sinh()
  1321. {
  1322. if (IsSpecialValue)
  1323. {
  1324. return;
  1325. }
  1326. Exp();
  1327. Sub(Reciprocal());
  1328. exponent--;
  1329. }
  1330. /// <summary>
  1331. /// Hyperbolic cosine (cosh) function
  1332. /// </summary>
  1333. public void Cosh()
  1334. {
  1335. if (IsSpecialValue)
  1336. {
  1337. if (SpecialValue == SpecialValueType.ZERO)
  1338. {
  1339. Assign(new BigFloat(1, mantissa.Precision));
  1340. }
  1341. else if (SpecialValue == SpecialValueType.INF_MINUS)
  1342. {
  1343. SetInfPlus();
  1344. }
  1345. return;
  1346. }
  1347. Exp();
  1348. Add(Reciprocal());
  1349. exponent--;
  1350. }
  1351. /// <summary>
  1352. /// Hyperbolic tangent function (tanh)
  1353. /// </summary>
  1354. public void Tanh()
  1355. {
  1356. if (IsSpecialValue)
  1357. {
  1358. if (SpecialValue == SpecialValueType.INF_MINUS)
  1359. {
  1360. Assign(new BigFloat(-1, mantissa.Precision));
  1361. }
  1362. else if (SpecialValue == SpecialValueType.INF_PLUS)
  1363. {
  1364. Assign(new BigFloat(1, mantissa.Precision));
  1365. }
  1366. return;
  1367. }
  1368. exponent++;
  1369. Exp();
  1370. BigFloat temp = new BigFloat(this);
  1371. BigFloat one = new BigFloat(1, mantissa.Precision);
  1372. temp.Add(one);
  1373. Sub(one);
  1374. Div(temp);
  1375. }
  1376. /// <summary>
  1377. /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)
  1378. /// </summary>
  1379. public void Arcsin()
  1380. {
  1381. if (IsSpecialValue)
  1382. {
  1383. if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)
  1384. {
  1385. SetNaN();
  1386. return;
  1387. }
  1388. return;
  1389. }
  1390. BigFloat one = new BigFloat(1, mantissa.Precision);
  1391. BigFloat plusABit = new BigFloat(1, mantissa.Precision);
  1392. plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));
  1393. BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);
  1394. onePlusABit.Add(plusABit);
  1395. bool sign = mantissa.Sign;
  1396. mantissa.Sign = false;
  1397. if (GreaterThan(onePlusABit))
  1398. {
  1399. SetNaN();
  1400. }
  1401. else if (LessThan(one))
  1402. {
  1403. BigFloat temp = new BigFloat(this);
  1404. temp.Mul(this);
  1405. temp.Sub(one);
  1406. temp.mantissa.Sign = !temp.mantissa.Sign;
  1407. temp.Sqrt();
  1408. temp.Add(one);
  1409. Div(temp);
  1410. Arctan();
  1411. exponent++;
  1412. mantissa.Sign = sign;
  1413. }
  1414. else
  1415. {
  1416. if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
  1417. {
  1418. CalculatePi(mantissa.Precision.NumBits);
  1419. }
  1420. Assign(piBy2);
  1421. if (sign) mantissa.Sign = true;
  1422. }
  1423. }
  1424. /// <summary>
  1425. /// arccos(): the inverse function of cos(), range (0..pi)
  1426. /// </summary>
  1427. public void Arccos()
  1428. {
  1429. if (IsSpecialValue)
  1430. {
  1431. if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS || SpecialValue == SpecialValueType.NAN)
  1432. {
  1433. SetNaN();
  1434. }
  1435. else if (SpecialValue == SpecialValueType.ZERO)
  1436. {
  1437. Assign(new BigFloat(1, mantissa.Precision));
  1438. exponent = 0;
  1439. Sign = false;
  1440. }
  1441. return;
  1442. }
  1443. BigFloat one = new BigFloat(1, mantissa.Precision);
  1444. BigFloat plusABit = new BigFloat(1, mantissa.Precision);
  1445. plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));
  1446. BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);
  1447. onePlusABit.Add(plusABit);
  1448. bool sign = mantissa.Sign;
  1449. mantissa.Sign = false;
  1450. if (GreaterThan(onePlusABit))
  1451. {
  1452. SetNaN();
  1453. }
  1454. else if (LessThan(one))
  1455. {
  1456. if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
  1457. {
  1458. CalculatePi(mantissa.Precision.NumBits);
  1459. }
  1460. mantissa.Sign = sign;
  1461. BigFloat temp = new BigFloat(this);
  1462. Mul(temp);
  1463. Sub(one);
  1464. mantissa.Sign = !mantissa.Sign;
  1465. Sqrt();
  1466. temp.Add(one);
  1467. Div(temp);
  1468. Arctan();
  1469. exponent++;
  1470. }
  1471. else
  1472. {
  1473. if (sign)
  1474. {
  1475. if (pi == null || pi.mantissa.Precision.NumBits != mantissa.Precision.NumBits)
  1476. {
  1477. CalculatePi(mantissa.Precision.NumBits);
  1478. }
  1479. Assign(pi);
  1480. }
  1481. else
  1482. {
  1483. mantissa.Zero();
  1484. exponent = 0;
  1485. }
  1486. }
  1487. }
  1488. /// <summary>
  1489. /// arctan(): the inverse function of sin(), range of (-pi/2..pi/2)
  1490. /// </summary>
  1491. public void Arctan()
  1492. {
  1493. //With 2 argument reductions, we increase precision by a minimum of 4 bits per term.
  1494. int numBits = mantissa.Precision.NumBits;
  1495. int maxTerms = numBits >> 2;
  1496. if (pi == null || pi.mantissa.Precision.NumBits != numBits)
  1497. {
  1498. CalculatePi(mantissa.Precision.NumBits);
  1499. }
  1500. //Make domain positive
  1501. bool sign = mantissa.Sign;
  1502. mantissa.Sign = false;
  1503. if (IsSpecialValue)
  1504. {
  1505. if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
  1506. {
  1507. Assign(piBy2);
  1508. mantissa.Sign = sign;
  1509. return;
  1510. }
  1511. return;
  1512. }
  1513. if (reciprocals == null || reciprocals[0].mantissa.Precision.NumBits != numBits || reciprocals.Length < maxTerms)
  1514. {
  1515. CalculateReciprocals(numBits, maxTerms);
  1516. }
  1517. bool invert = false;
  1518. BigFloat one = new BigFloat(1, mantissa.Precision);
  1519. //Invert if outside of convergence
  1520. if (GreaterThan(one))
  1521. {
  1522. invert = true;
  1523. Assign(Reciprocal());
  1524. }
  1525. //Reduce using half-angle formula:
  1526. //arctan(2x) = 2 arctan (x / (1 + sqrt(1 + x)))
  1527. //First reduction (guarantees 2 bits per iteration)
  1528. BigFloat temp = new BigFloat(this);
  1529. temp.Mul(this);
  1530. temp.Add(one);
  1531. temp.Sqrt();
  1532. temp.Add(one);
  1533. this.Div(temp);
  1534. //Second reduction (guarantees 4 bits per iteration)
  1535. temp.Assign(this);
  1536. temp.Mul(this);
  1537. temp.Add(one);
  1538. temp.Sqrt();
  1539. temp.Add(one);
  1540. this.Div(temp);
  1541. //Actual series calculation
  1542. int length = reciprocals.Length;
  1543. BigFloat term = new BigFloat(this);
  1544. //pow = x^2
  1545. BigFloat pow = new BigFloat(this);
  1546. pow.Mul(this);
  1547. BigFloat sum = new BigFloat(this);
  1548. for (int i = 1; i < length; i++)
  1549. {
  1550. //u(n) = u(n-1) * x^2
  1551. //t(n) = u(n) / (2n+1)
  1552. term.Mul(pow);
  1553. term.Sign = !term.Sign;
  1554. temp.Assign(term);
  1555. temp.Mul(reciprocals[i]);
  1556. if (temp.exponent < -numBits) break;
  1557. sum.Add(temp);
  1558. }
  1559. //Undo the reductions.
  1560. Assign(sum);
  1561. exponent += 2;
  1562. if (invert)
  1563. {
  1564. //Assign(Reciprocal());
  1565. mantissa.Sign = true;
  1566. Add(piBy2);
  1567. }
  1568. if (sign)
  1569. {
  1570. mantissa.Sign = sign;
  1571. }
  1572. }
  1573. /// <summary>
  1574. /// Arcsinh(): the inverse sinh function
  1575. /// </summary>
  1576. public void Arcsinh()
  1577. {
  1578. //Just let all special values fall through
  1579. if (IsSpecialValue)
  1580. {
  1581. return;
  1582. }
  1583. BigFloat temp = new BigFloat(this);
  1584. temp.Mul(this);
  1585. temp.Add(new BigFloat(1, mantissa.Precision));
  1586. temp.Sqrt();
  1587. Add(temp);
  1588. Log();
  1589. }
  1590. /// <summary>
  1591. /// Arccosh(): the inverse cosh() function
  1592. /// </summary>
  1593. public void Arccosh()
  1594. {
  1595. //acosh isn't defined for x < 1
  1596. if (IsSpecialValue)
  1597. {
  1598. if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.ZERO)
  1599. {
  1600. SetNaN();
  1601. return;
  1602. }
  1603. return;
  1604. }
  1605. BigFloat one = new BigFloat(1, mantissa.Precision);
  1606. if (LessThan(one))
  1607. {
  1608. SetNaN();
  1609. return;
  1610. }
  1611. BigFloat temp = new BigFloat(this);
  1612. temp.Mul(this);
  1613. temp.Sub(one);
  1614. temp.Sqrt();
  1615. Add(temp);
  1616. Log();
  1617. }
  1618. /// <summary>
  1619. /// Arctanh(): the inverse tanh function
  1620. /// </summary>
  1621. public void Arctanh()
  1622. {
  1623. //|x| <= 1 for a non-NaN output
  1624. if (IsSpecialValue)
  1625. {
  1626. if (SpecialValue == SpecialValueType.INF_MINUS || SpecialValue == SpecialValueType.INF_PLUS)
  1627. {
  1628. SetNaN();
  1629. return;
  1630. }
  1631. return;
  1632. }
  1633. BigFloat one = new BigFloat(1, mantissa.Precision);
  1634. BigFloat plusABit = new BigFloat(1, mantissa.Precision);
  1635. plusABit.exponent -= (mantissa.Precision.NumBits - (mantissa.Precision.NumBits >> 6));
  1636. BigFloat onePlusABit = new BigFloat(1, mantissa.Precision);
  1637. onePlusABit.Add(plusABit);
  1638. bool sign = mantissa.Sign;
  1639. mantissa.Sign = false;
  1640. if (GreaterThan(onePlusABit))
  1641. {
  1642. SetNaN();
  1643. }
  1644. else if (LessThan(one))
  1645. {
  1646. BigFloat temp = new BigFloat(this);
  1647. Add(one);
  1648. one.Sub(temp);
  1649. Div(one);
  1650. Log();
  1651. exponent--;
  1652. mantissa.Sign = sign;
  1653. }
  1654. else
  1655. {
  1656. if (sign)
  1657. {
  1658. SetInfMinus();
  1659. }
  1660. else
  1661. {
  1662. SetInfPlus();
  1663. }
  1664. }
  1665. }
  1666. /// <summary>
  1667. /// Two-variable iterative square root, taken from
  1668. /// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method
  1669. /// </summary>
  1670. public void Sqrt()
  1671. {
  1672. if (mantissa.Sign || IsSpecialValue)
  1673. {
  1674. if (SpecialValue == SpecialValueType.ZERO)
  1675. {
  1676. return;
  1677. }
  1678. if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)
  1679. {
  1680. SetNaN();
  1681. }
  1682. return;
  1683. }
  1684. BigFloat temp2;
  1685. BigFloat temp3 = new BigFloat(mantissa.Precision);
  1686. BigFloat three = new BigFloat(3, mantissa.Precision);
  1687. int exponentScale = 0;
  1688. //Rescale to 0.5 <= x < 2
  1689. if (exponent < -1)
  1690. {
  1691. int diff = -exponent;
  1692. if ((diff & 1) != 0)
  1693. {
  1694. diff--;
  1695. }
  1696. exponentScale = -diff;
  1697. exponent += diff;
  1698. }
  1699. else if (exponent > 0)
  1700. {
  1701. if ((exponent & 1) != 0)
  1702. {
  1703. exponentScale = exponent + 1;
  1704. exponent = -1;
  1705. }
  1706. else
  1707. {
  1708. exponentScale = exponent;
  1709. exponent = 0;
  1710. }
  1711. }
  1712. temp2 = new BigFloat(this);
  1713. temp2.Sub(new BigFloat(1, mantissa.Precision));
  1714. //if (temp2.mantissa.IsZero())
  1715. //{
  1716. // exponent += exponentScale;
  1717. // return;
  1718. //}
  1719. int numBits = mantissa.Precision.NumBits;
  1720. while ((exponent - temp2.exponent) < numBits && temp2.SpecialValue != SpecialValueType.ZERO)
  1721. {
  1722. //a(n+1) = an - an*cn / 2
  1723. temp3.Assign(this);
  1724. temp3.Mul(temp2);
  1725. temp3.MulPow2(-1);
  1726. this.Sub(temp3);
  1727. //c(n+1) = cn^2 * (cn - 3) / 4
  1728. temp3.Assign(temp2);
  1729. temp2.Sub(three);
  1730. temp2.Mul(temp3);
  1731. temp2.Mul(temp3);
  1732. temp2.MulPow2(-2);
  1733. }
  1734. exponent += (exponentScale >> 1);
  1735. }
  1736. /// <summary>
  1737. /// The natural logarithm, ln(x)
  1738. /// </summary>
  1739. public void Log()
  1740. {
  1741. if (IsSpecialValue || mantissa.Sign)
  1742. {
  1743. if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)
  1744. {
  1745. SetNaN();
  1746. }
  1747. else if (SpecialValue == SpecialValueType.ZERO)
  1748. {
  1749. SetInfMinus();
  1750. }
  1751. return;
  1752. }
  1753. if (mantissa.Precision.NumBits >= 512)
  1754. {
  1755. LogAGM1();
  1756. return;
  1757. }
  1758. //Compute ln2.
  1759. if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
  1760. {
  1761. CalculateLog2(mantissa.Precision.NumBits);
  1762. }
  1763. Log2();
  1764. Mul(ln2cache);
  1765. }
  1766. /// <summary>
  1767. /// Log to the base 10
  1768. /// </summary>
  1769. public void Log10()
  1770. {
  1771. if (IsSpecialValue || mantissa.Sign)
  1772. {
  1773. if (SpecialValue == SpecialValueType.INF_MINUS || mantissa.Sign)
  1774. {
  1775. SetNaN();
  1776. }
  1777. else if (SpecialValue == SpecialValueType.ZERO)
  1778. {
  1779. SetInfMinus();
  1780. }
  1781. return;
  1782. }
  1783. //Compute ln2.
  1784. if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
  1785. {
  1786. CalculateLog2(mantissa.Precision.NumBits);
  1787. }
  1788. Log();
  1789. Mul(log10recip);
  1790. }
  1791. /// <summary>
  1792. /// The exponential function. Less accurate for high exponents, scales poorly with the number
  1793. /// of bits.
  1794. /// </summary>
  1795. public void Exp()
  1796. {
  1797. Exp(mantissa.Precision.NumBits);
  1798. }
  1799. /// <summary>
  1800. /// Raises a number to an integer power (positive or negative). This is a very accurate and fast function,
  1801. /// comparable to or faster than division (although it is slightly slower for
  1802. /// negative powers, obviously)
  1803. ///
  1804. /// </summary>
  1805. /// <param name="power"></param>
  1806. public void Pow(int power)
  1807. {
  1808. BigFloat acc = new BigFloat(1, mantissa.Precision);
  1809. BigFloat temp = new BigFloat(1, mantissa.Precision);
  1810. int powerTemp = power;
  1811. if (power < 0)
  1812. {
  1813. Assign(Reciprocal());
  1814. powerTemp = -power;
  1815. }
  1816. //Fast power function
  1817. while (powerTemp != 0)
  1818. {
  1819. temp.Mul(this);
  1820. Assign(temp);
  1821. if ((powerTemp & 1) != 0)
  1822. {
  1823. acc.Mul(temp);
  1824. }
  1825. powerTemp >>= 1;
  1826. }
  1827. Assign(acc);
  1828. }
  1829. /// <summary>
  1830. /// Raises to an aribitrary power. This is both slow (uses Log) and inaccurate. If you need to
  1831. /// raise e^x use exp(). If you need an integer power, use the integer power function Pow(int)
  1832. /// Accuracy Note:
  1833. /// The function is only ever accurate to a maximum of 4 decimal digits
  1834. /// For every 10x larger (or smaller) the power gets, you lose an additional decimal digit
  1835. /// If you really need a precise result, do the calculation with an extra 32-bits and round
  1836. /// Domain Note:
  1837. /// This only works for powers of positive real numbers. Negative numbers will fail.
  1838. /// </summary>
  1839. /// <param name="power"></param>
  1840. public void Pow(BigFloat power)
  1841. {
  1842. Log();
  1843. Mul(power);
  1844. Exp();
  1845. }
  1846. //******************** Static Math Functions *******************
  1847. /// <summary>
  1848. /// Returns the integer component of the input
  1849. /// </summary>
  1850. /// <param name="n1">The input number</param>
  1851. /// <remarks>The integer component returned will always be numerically closer to zero
  1852. /// than the input: an input of -3.49 for instance would produce a value of 3.</remarks>
  1853. public static BigFloat Floor(BigFloat n1)
  1854. {
  1855. BigFloat res = new BigFloat(n1);
  1856. n1.Floor();
  1857. return n1;
  1858. }
  1859. /// <summary>
  1860. /// Returns the fractional (non-integer component of the input)
  1861. /// </summary>
  1862. /// <param name="n1">The input number</param>
  1863. public static BigFloat FPart(BigFloat n1)
  1864. {
  1865. BigFloat res = new BigFloat(n1);
  1866. n1.FPart();
  1867. return n1;
  1868. }
  1869. /// <summary>
  1870. /// Calculates tan(x)
  1871. /// </summary>
  1872. /// <param name="n1">The angle (in radians) to find the tangent of</param>
  1873. public static BigFloat Tan(BigFloat n1)
  1874. {
  1875. BigFloat res = new BigFloat(n1);
  1876. n1.Tan();
  1877. return n1;
  1878. }
  1879. /// <summary>
  1880. /// Calculates Cos(x)
  1881. /// </summary>
  1882. /// <param name="n1">The angle (in radians) to find the cosine of</param>
  1883. /// <remarks>This is a reasonably fast function for smaller precisions, but
  1884. /// doesn't scale well for higher precision arguments</remarks>
  1885. public static BigFloat Cos(BigFloat n1)
  1886. {
  1887. BigFloat res = new BigFloat(n1);
  1888. n1.Cos();
  1889. return n1;
  1890. }
  1891. /// <summary>
  1892. /// Calculates Sin(x):
  1893. /// This takes a little longer and is less accurate if the input is out of the range (-pi, pi].
  1894. /// </summary>
  1895. /// <param name="n1">The angle to find the sine of (in radians)</param>
  1896. /// <remarks>This is a resonably fast function, for smaller precision arguments, but doesn't
  1897. /// scale very well with the number of bits in the input.</remarks>
  1898. public static BigFloat Sin(BigFloat n1)
  1899. {
  1900. BigFloat res = new BigFloat(n1);
  1901. n1.Sin();
  1902. return n1;
  1903. }
  1904. /// <summary>
  1905. /// Hyperbolic Sin (sinh) function
  1906. /// </summary>
  1907. /// <param name="n1">The number to find the hyperbolic sine of</param>
  1908. public static BigFloat Sinh(BigFloat n1)
  1909. {
  1910. BigFloat res = new BigFloat(n1);
  1911. n1.Sinh();
  1912. return n1;
  1913. }
  1914. /// <summary>
  1915. /// Hyperbolic cosine (cosh) function
  1916. /// </summary>
  1917. /// <param name="n1">The number to find the hyperbolic cosine of</param>
  1918. public static BigFloat Cosh(BigFloat n1)
  1919. {
  1920. BigFloat res = new BigFloat(n1);
  1921. n1.Cosh();
  1922. return n1;
  1923. }
  1924. /// <summary>
  1925. /// Hyperbolic tangent function (tanh)
  1926. /// </summary>
  1927. /// <param name="n1">The number to find the hyperbolic tangent of</param>
  1928. public static BigFloat Tanh(BigFloat n1)
  1929. {
  1930. BigFloat res = new BigFloat(n1);
  1931. n1.Tanh();
  1932. return n1;
  1933. }
  1934. /// <summary>
  1935. /// arcsin(): the inverse function of sin(), range of (-pi/2..pi/2)
  1936. /// </summary>
  1937. /// <param name="n1">The number to find the arcsine of (-pi/2..pi/2)</param>
  1938. /// <remarks>Note that inverse trig functions are only defined within a specific range.
  1939. /// Values outside this range will return NaN, although some margin for error is assumed.
  1940. /// </remarks>
  1941. public static BigFloat Arcsin(BigFloat n1)
  1942. {
  1943. BigFloat res = new BigFloat(n1);
  1944. n1.Arcsin();
  1945. return n1;
  1946. }
  1947. /// <summary>
  1948. /// arccos(): the inverse function of cos(), input range (0..pi)
  1949. /// </summary>
  1950. /// <param name="n1">The number to find the arccosine of (0..pi)</param>
  1951. /// <remarks>Note that inverse trig functions are only defined within a specific range.
  1952. /// Values outside this range will return NaN, although some margin for error is assumed.
  1953. /// </remarks>
  1954. public static BigFloat Arccos(BigFloat n1)
  1955. {
  1956. BigFloat res = new BigFloat(n1);
  1957. n1.Arccos();
  1958. return n1;
  1959. }
  1960. /// <summary>
  1961. /// arctan(): the inverse function of sin(), input range of (-pi/2..pi/2)
  1962. /// </summary>
  1963. /// <param name="n1">The number to find the arctangent of (-pi/2..pi/2)</param>
  1964. /// <remarks>Note that inverse trig functions are only defined within a specific range.
  1965. /// Values outside this range will return NaN, although some margin for error is assumed.
  1966. /// </remarks>
  1967. public static BigFloat Arctan(BigFloat n1)
  1968. {
  1969. BigFloat res = new BigFloat(n1);
  1970. n1.Arctan();
  1971. return n1;
  1972. }
  1973. /// <summary>
  1974. /// Arcsinh(): the inverse sinh function
  1975. /// </summary>
  1976. /// <param name="n1">The number to find the inverse hyperbolic sine of</param>
  1977. public static BigFloat Arcsinh(BigFloat n1)
  1978. {
  1979. BigFloat res = new BigFloat(n1);
  1980. n1.Arcsinh();
  1981. return n1;
  1982. }
  1983. /// <summary>
  1984. /// Arccosh(): the inverse cosh() function
  1985. /// </summary>
  1986. /// <param name="n1">The number to find the inverse hyperbolic cosine of</param>
  1987. public static BigFloat Arccosh(BigFloat n1)
  1988. {
  1989. BigFloat res = new BigFloat(n1);
  1990. n1.Arccosh();
  1991. return n1;
  1992. }
  1993. /// <summary>
  1994. /// Arctanh(): the inverse tanh function
  1995. /// </summary>
  1996. /// <param name="n1">The number to fine the inverse hyperbolic tan of</param>
  1997. public static BigFloat Arctanh(BigFloat n1)
  1998. {
  1999. BigFloat res = new BigFloat(n1);
  2000. n1.Arctanh();
  2001. return n1;
  2002. }
  2003. /// <summary>
  2004. /// Two-variable iterative square root, taken from
  2005. /// http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#A_two-variable_iterative_method
  2006. /// </summary>
  2007. /// <remarks>This is quite a fast function, as elementary functions go. You can expect it to take
  2008. /// about twice as long as a floating-point division.
  2009. /// </remarks>
  2010. public static BigFloat Sqrt(BigFloat n1)
  2011. {
  2012. BigFloat res = new BigFloat(n1);
  2013. n1.Sqrt();
  2014. return n1;
  2015. }
  2016. /// <summary>
  2017. /// The natural logarithm, ln(x) (log base e)
  2018. /// </summary>
  2019. /// <remarks>This is a very slow function, despite repeated attempts at optimisation.
  2020. /// To make it any faster, different strategies would be needed for integer operations.
  2021. /// It does, however, scale well with the number of bits.
  2022. /// </remarks>
  2023. /// <param name="n1">The number to find the natural logarithm of</param>
  2024. public static BigFloat Log(BigFloat n1)
  2025. {
  2026. BigFloat res = new BigFloat(n1);
  2027. n1.Log();
  2028. return n1;
  2029. }
  2030. /// <summary>
  2031. /// Base 10 logarithm of a number
  2032. /// </summary>
  2033. /// <remarks>This is a very slow function, despite repeated attempts at optimisation.
  2034. /// To make it any faster, different strategies would be needed for integer operations.
  2035. /// It does, however, scale well with the number of bits.
  2036. /// </remarks>
  2037. /// <param name="n1">The number to find the base 10 logarithm of</param>
  2038. public static BigFloat Log10(BigFloat n1)
  2039. {
  2040. BigFloat res = new BigFloat(n1);
  2041. n1.Log10();
  2042. return n1;
  2043. }
  2044. /// <summary>
  2045. /// The exponential function. Less accurate for high exponents, scales poorly with the number
  2046. /// of bits. This is quite fast for low-precision arguments.
  2047. /// </summary>
  2048. public static BigFloat Exp(BigFloat n1)
  2049. {
  2050. BigFloat res = new BigFloat(n1);
  2051. n1.Exp();
  2052. return n1;
  2053. }
  2054. /// <summary>
  2055. /// Raises a number to an integer power (positive or negative). This is a very accurate and fast function,
  2056. /// comparable to or faster than division (although it is slightly slower for
  2057. /// negative powers, obviously).
  2058. /// </summary>
  2059. /// <param name="n1">The number to raise to the power</param>
  2060. /// <param name="power">The power to raise it to</param>
  2061. public static BigFloat Pow(BigFloat n1, int power)
  2062. {
  2063. BigFloat res = new BigFloat(n1);
  2064. n1.Pow(power);
  2065. return n1;
  2066. }
  2067. /// <summary>
  2068. /// Raises to an aribitrary power. This is both slow (uses Log) and inaccurate. If you need to
  2069. /// raise e^x use exp(). If you need an integer power, use the integer power function Pow(int)
  2070. /// </summary>
  2071. /// <remarks>
  2072. /// Accuracy Note:
  2073. /// The function is only ever accurate to a maximum of 4 decimal digits
  2074. /// For every 10x larger (or smaller) the power gets, you lose an additional decimal digit
  2075. /// If you really need a precise result, do the calculation with an extra 32-bits and round
  2076. ///
  2077. /// Domain Note:
  2078. /// This only works for powers of positive real numbers. Negative numbers will fail.
  2079. /// </remarks>
  2080. /// <param name="n1">The number to raise to a power</param>
  2081. /// <param name="power">The power to raise it to</param>
  2082. public static BigFloat Pow(BigFloat n1, BigFloat power)
  2083. {
  2084. BigFloat res = new BigFloat(n1);
  2085. n1.Pow(power);
  2086. return n1;
  2087. }
  2088. //********************** Static functions **********************
  2089. /// <summary>
  2090. /// Adds two numbers and returns the result
  2091. /// </summary>
  2092. public static BigFloat Add(BigFloat n1, BigFloat n2)
  2093. {
  2094. BigFloat ret = new BigFloat(n1);
  2095. ret.Add(n2);
  2096. return ret;
  2097. }
  2098. /// <summary>
  2099. /// Subtracts two numbers and returns the result
  2100. /// </summary>
  2101. public static BigFloat Sub(BigFloat n1, BigFloat n2)
  2102. {
  2103. BigFloat ret = new BigFloat(n1);
  2104. ret.Sub(n2);
  2105. return ret;
  2106. }
  2107. /// <summary>
  2108. /// Multiplies two numbers and returns the result
  2109. /// </summary>
  2110. public static BigFloat Mul(BigFloat n1, BigFloat n2)
  2111. {
  2112. BigFloat ret = new BigFloat(n1);
  2113. ret.Mul(n2);
  2114. return ret;
  2115. }
  2116. /// <summary>
  2117. /// Divides two numbers and returns the result
  2118. /// </summary>
  2119. public static BigFloat Div(BigFloat n1, BigFloat n2)
  2120. {
  2121. BigFloat ret = new BigFloat(n1);
  2122. ret.Div(n2);
  2123. return ret;
  2124. }
  2125. /// <summary>
  2126. /// Tests whether n1 is greater than n2
  2127. /// </summary>
  2128. public static bool GreaterThan(BigFloat n1, BigFloat n2)
  2129. {
  2130. return n1.GreaterThan(n2);
  2131. }
  2132. /// <summary>
  2133. /// Tests whether n1 is less than n2
  2134. /// </summary>
  2135. public static bool LessThan(BigFloat n1, BigFloat n2)
  2136. {
  2137. return n1.LessThan(n2);
  2138. }
  2139. //******************* Fast static functions ********************
  2140. /// <summary>
  2141. /// Adds two numbers and assigns the result to res.
  2142. /// </summary>
  2143. /// <param name="res">a pre-existing BigFloat to take the result</param>
  2144. /// <param name="n1">the first number</param>
  2145. /// <param name="n2">the second number</param>
  2146. /// <returns>a handle to res</returns>
  2147. public static BigFloat Add(BigFloat res, BigFloat n1, BigFloat n2)
  2148. {
  2149. res.Assign(n1);
  2150. res.Add(n2);
  2151. return res;
  2152. }
  2153. /// <summary>
  2154. /// Subtracts two numbers and assigns the result to res.
  2155. /// </summary>
  2156. /// <param name="res">a pre-existing BigFloat to take the result</param>
  2157. /// <param name="n1">the first number</param>
  2158. /// <param name="n2">the second number</param>
  2159. /// <returns>a handle to res</returns>
  2160. public static BigFloat Sub(BigFloat res, BigFloat n1, BigFloat n2)
  2161. {
  2162. res.Assign(n1);
  2163. res.Sub(n2);
  2164. return res;
  2165. }
  2166. /// <summary>
  2167. /// Multiplies two numbers and assigns the result to res.
  2168. /// </summary>
  2169. /// <param name="res">a pre-existing BigFloat to take the result</param>
  2170. /// <param name="n1">the first number</param>
  2171. /// <param name="n2">the second number</param>
  2172. /// <returns>a handle to res</returns>
  2173. public static BigFloat Mul(BigFloat res, BigFloat n1, BigFloat n2)
  2174. {
  2175. res.Assign(n1);
  2176. res.Mul(n2);
  2177. return res;
  2178. }
  2179. /// <summary>
  2180. /// Divides two numbers and assigns the result to res.
  2181. /// </summary>
  2182. /// <param name="res">a pre-existing BigFloat to take the result</param>
  2183. /// <param name="n1">the first number</param>
  2184. /// <param name="n2">the second number</param>
  2185. /// <returns>a handle to res</returns>
  2186. public static BigFloat Div(BigFloat res, BigFloat n1, BigFloat n2)
  2187. {
  2188. res.Assign(n1);
  2189. res.Div(n2);
  2190. return res;
  2191. }
  2192. //************************* Operators **************************
  2193. /// <summary>
  2194. /// The addition operator
  2195. /// </summary>
  2196. public static BigFloat operator +(BigFloat n1, BigFloat n2)
  2197. {
  2198. return Add(n1, n2);
  2199. }
  2200. /// <summary>
  2201. /// The subtraction operator
  2202. /// </summary>
  2203. public static BigFloat operator -(BigFloat n1, BigFloat n2)
  2204. {
  2205. return Sub(n1, n2);
  2206. }
  2207. /// <summary>
  2208. /// The multiplication operator
  2209. /// </summary>
  2210. public static BigFloat operator *(BigFloat n1, BigFloat n2)
  2211. {
  2212. return Mul(n1, n2);
  2213. }
  2214. /// <summary>
  2215. /// The division operator
  2216. /// </summary>
  2217. public static BigFloat operator /(BigFloat n1, BigFloat n2)
  2218. {
  2219. return Div(n1, n2);
  2220. }
  2221. //************************** Conversions *************************
  2222. /// <summary>
  2223. /// Converts a BigFloat to an BigInt with the specified precision
  2224. /// </summary>
  2225. /// <param name="n1">The number to convert</param>
  2226. /// <param name="precision">The precision to convert it with</param>
  2227. /// <param name="round">Do we round the number if we are truncating the mantissa?</param>
  2228. /// <returns></returns>
  2229. public static BigInt ConvertToInt(BigFloat n1, PrecisionSpec precision, bool round)
  2230. {
  2231. BigInt ret = new BigInt(precision);
  2232. int numBits = n1.mantissa.Precision.NumBits;
  2233. int shift = numBits - (n1.exponent + 1);
  2234. BigFloat copy = new BigFloat(n1);
  2235. bool inc = false;
  2236. //Rounding
  2237. if (copy.mantissa.Precision.NumBits > ret.Precision.NumBits)
  2238. {
  2239. inc = true;
  2240. for (int i = copy.exponent + 1; i <= ret.Precision.NumBits; i++)
  2241. {
  2242. if (copy.mantissa.GetBitFromTop(i) == 0)
  2243. {
  2244. inc = false;
  2245. break;
  2246. }
  2247. }
  2248. }
  2249. if (shift > 0)
  2250. {
  2251. copy.mantissa.RSH(shift);
  2252. }
  2253. else if (shift < 0)
  2254. {
  2255. copy.mantissa.LSH(-shift);
  2256. }
  2257. ret.Assign(copy.mantissa);
  2258. if (inc) ret.Increment();
  2259. return ret;
  2260. }
  2261. /// <summary>
  2262. /// Returns a base-10 string representing the number.
  2263. ///
  2264. /// Note: This is inefficient and possibly inaccurate. Please use with enough
  2265. /// rounding digits (set using the RoundingDigits property) to ensure accuracy
  2266. /// </summary>
  2267. public override string ToString()
  2268. {
  2269. if (IsSpecialValue)
  2270. {
  2271. SpecialValueType s = SpecialValue;
  2272. if (s == SpecialValueType.ZERO)
  2273. {
  2274. return String.Format("0{0}0", System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator);
  2275. }
  2276. else if (s == SpecialValueType.INF_PLUS)
  2277. {
  2278. return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.PositiveInfinitySymbol;
  2279. }
  2280. else if (s == SpecialValueType.INF_MINUS)
  2281. {
  2282. return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NegativeInfinitySymbol;
  2283. }
  2284. else if (s == SpecialValueType.NAN)
  2285. {
  2286. return System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NaNSymbol;
  2287. }
  2288. else
  2289. {
  2290. return "Unrecognised special type";
  2291. }
  2292. }
  2293. if (scratch.Precision.NumBits != mantissa.Precision.NumBits)
  2294. {
  2295. scratch = new BigInt(mantissa.Precision);
  2296. }
  2297. //The mantissa expresses 1.xxxxxxxxxxx
  2298. //The highest possible value for the mantissa without the implicit 1. is 0.9999999...
  2299. scratch.Assign(mantissa);
  2300. //scratch.Round(3);
  2301. scratch.Sign = false;
  2302. BigInt denom = new BigInt("0", mantissa.Precision);
  2303. denom.SetBit(mantissa.Precision.NumBits - 1);
  2304. bool useExponentialNotation = false;
  2305. int halfBits = mantissa.Precision.NumBits / 2;
  2306. if (halfBits > 60) halfBits = 60;
  2307. int precDec = 10;
  2308. if (exponent > 0)
  2309. {
  2310. if (exponent < halfBits)
  2311. {
  2312. denom.RSH(exponent);
  2313. }
  2314. else
  2315. {
  2316. useExponentialNotation = true;
  2317. }
  2318. }
  2319. else if (exponent < 0)
  2320. {
  2321. int shift = -(exponent);
  2322. if (shift < precDec)
  2323. {
  2324. scratch.RSH(shift);
  2325. }
  2326. else
  2327. {
  2328. useExponentialNotation = true;
  2329. }
  2330. }
  2331. string output;
  2332. if (useExponentialNotation)
  2333. {
  2334. int absExponent = exponent;
  2335. if (absExponent < 0) absExponent = -absExponent;
  2336. int powerOf10 = (int)((double)absExponent * Math.Log10(2.0));
  2337. //Use 1 extra digit of precision (this is actually 32 bits more, nb)
  2338. BigFloat thisFloat = new BigFloat(this, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
  2339. thisFloat.mantissa.Sign = false;
  2340. //Multiplicative correction factor to bring number into range.
  2341. BigFloat one = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
  2342. BigFloat ten = new BigFloat(10, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
  2343. BigFloat tenRCP = ten.Reciprocal();
  2344. //Accumulator for the power of 10 calculation.
  2345. BigFloat acc = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
  2346. BigFloat tenToUse;
  2347. if (exponent > 0)
  2348. {
  2349. tenToUse = new BigFloat(tenRCP, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
  2350. }
  2351. else
  2352. {
  2353. tenToUse = new BigFloat(ten, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
  2354. }
  2355. BigFloat tenToPower = new BigFloat(1, new PrecisionSpec(mantissa.Precision.NumBits + 1, PrecisionSpec.BaseType.BIN));
  2356. int powerTemp = powerOf10;
  2357. //Fast power function
  2358. while (powerTemp != 0)
  2359. {
  2360. tenToPower.Mul(tenToUse);
  2361. tenToUse.Assign(tenToPower);
  2362. if ((powerTemp & 1) != 0)
  2363. {
  2364. acc.Mul(tenToPower);
  2365. }
  2366. powerTemp >>= 1;
  2367. }
  2368. thisFloat.Mul(acc);
  2369. //If we are out of range, correct.
  2370. if (thisFloat.GreaterThan(ten))
  2371. {
  2372. thisFloat.Mul(tenRCP);
  2373. if (exponent > 0)
  2374. {
  2375. powerOf10++;
  2376. }
  2377. else
  2378. {
  2379. powerOf10--;
  2380. }
  2381. }
  2382. else if (thisFloat.LessThan(one))
  2383. {
  2384. thisFloat.Mul(ten);
  2385. if (exponent > 0)
  2386. {
  2387. powerOf10--;
  2388. }
  2389. else
  2390. {
  2391. powerOf10++;
  2392. }
  2393. }
  2394. //Restore the precision and the sign.
  2395. BigFloat printable = new BigFloat(thisFloat, mantissa.Precision);
  2396. printable.mantissa.Sign = mantissa.Sign;
  2397. output = printable.ToString();
  2398. if (exponent < 0) powerOf10 = -powerOf10;
  2399. output = String.Format("{0}E{1}", output, powerOf10);
  2400. }
  2401. else
  2402. {
  2403. BigInt bigDigit = BigInt.Div(scratch, denom);
  2404. bigDigit.Sign = false;
  2405. scratch.Sub(BigInt.Mul(denom, bigDigit));
  2406. if (mantissa.Sign)
  2407. {
  2408. output = String.Format("-{0}.", bigDigit);
  2409. }
  2410. else
  2411. {
  2412. output = String.Format("{0}.", bigDigit);
  2413. }
  2414. denom = BigInt.Div(denom, 10u);
  2415. while (!denom.IsZero())
  2416. {
  2417. uint digit = (uint)BigInt.Div(scratch, denom);
  2418. if (digit == 10) digit--;
  2419. scratch.Sub(BigInt.Mul(denom, digit));
  2420. output = String.Format("{0}{1}", output, digit);
  2421. denom = BigInt.Div(denom, 10u);
  2422. }
  2423. output = RoundString(output, RoundingDigits);
  2424. }
  2425. return output;
  2426. }
  2427. //**************** Special value handling for ops ***************
  2428. private void SetNaN()
  2429. {
  2430. exponent = Int32.MaxValue;
  2431. mantissa.SetBit(mantissa.Precision.NumBits - 1);
  2432. }
  2433. private void SetZero()
  2434. {
  2435. exponent = 0;
  2436. mantissa.Zero();
  2437. Sign = false;
  2438. }
  2439. private void SetInfPlus()
  2440. {
  2441. Sign = false;
  2442. exponent = Int32.MaxValue;
  2443. mantissa.Zero();
  2444. }
  2445. private void SetInfMinus()
  2446. {
  2447. Sign = true;
  2448. exponent = Int32.MaxValue;
  2449. mantissa.Zero();
  2450. }
  2451. private bool SpecialValueAddTest(BigFloat n2)
  2452. {
  2453. if (IsSpecialValue || n2.IsSpecialValue)
  2454. {
  2455. SpecialValueType s1 = SpecialValue;
  2456. SpecialValueType s2 = n2.SpecialValue;
  2457. if (s1 == SpecialValueType.NAN) return true;
  2458. if (s2 == SpecialValueType.NAN)
  2459. {
  2460. //Set NaN and return.
  2461. SetNaN();
  2462. return true;
  2463. }
  2464. if (s1 == SpecialValueType.INF_PLUS)
  2465. {
  2466. //INF+ + INF- = NAN
  2467. if (s2 == SpecialValueType.INF_MINUS)
  2468. {
  2469. SetNaN();
  2470. return true;
  2471. }
  2472. return true;
  2473. }
  2474. if (s1 == SpecialValueType.INF_MINUS)
  2475. {
  2476. //INF+ + INF- = NAN
  2477. if (s2 == SpecialValueType.INF_PLUS)
  2478. {
  2479. SetNaN();
  2480. return true;
  2481. }
  2482. return true;
  2483. }
  2484. if (s2 == SpecialValueType.ZERO)
  2485. {
  2486. return true;
  2487. }
  2488. if (s1 == SpecialValueType.ZERO)
  2489. {
  2490. Assign(n2);
  2491. return true;
  2492. }
  2493. }
  2494. return false;
  2495. }
  2496. private bool SpecialValueMulTest(BigFloat n2)
  2497. {
  2498. if (IsSpecialValue || n2.IsSpecialValue)
  2499. {
  2500. SpecialValueType s1 = SpecialValue;
  2501. SpecialValueType s2 = n2.SpecialValue;
  2502. if (s1 == SpecialValueType.NAN) return true;
  2503. if (s2 == SpecialValueType.NAN)
  2504. {
  2505. //Set NaN and return.
  2506. SetNaN();
  2507. return true;
  2508. }
  2509. if (s1 == SpecialValueType.INF_PLUS)
  2510. {
  2511. //Inf+ * Inf- = Inf-
  2512. if (s2 == SpecialValueType.INF_MINUS)
  2513. {
  2514. Assign(n2);
  2515. return true;
  2516. }
  2517. //Inf+ * 0 = NaN
  2518. if (s2 == SpecialValueType.ZERO)
  2519. {
  2520. //Set NaN and return.
  2521. SetNaN();
  2522. return true;
  2523. }
  2524. return true;
  2525. }
  2526. if (s1 == SpecialValueType.INF_MINUS)
  2527. {
  2528. //Inf- * Inf- = Inf+
  2529. if (s2 == SpecialValueType.INF_MINUS)
  2530. {
  2531. Sign = false;
  2532. return true;
  2533. }
  2534. //Inf- * 0 = NaN
  2535. if (s2 == SpecialValueType.ZERO)
  2536. {
  2537. //Set NaN and return.
  2538. SetNaN();
  2539. return true;
  2540. }
  2541. return true;
  2542. }
  2543. if (s2 == SpecialValueType.ZERO)
  2544. {
  2545. SetZero();
  2546. return true;
  2547. }
  2548. if (s1 == SpecialValueType.ZERO)
  2549. {
  2550. return true;
  2551. }
  2552. }
  2553. return false;
  2554. }
  2555. private bool SpecialValueDivTest(BigFloat n2)
  2556. {
  2557. if (IsSpecialValue || n2.IsSpecialValue)
  2558. {
  2559. SpecialValueType s1 = SpecialValue;
  2560. SpecialValueType s2 = n2.SpecialValue;
  2561. if (s1 == SpecialValueType.NAN) return true;
  2562. if (s2 == SpecialValueType.NAN)
  2563. {
  2564. //Set NaN and return.
  2565. SetNaN();
  2566. return true;
  2567. }
  2568. if ((s1 == SpecialValueType.INF_PLUS || s1 == SpecialValueType.INF_MINUS))
  2569. {
  2570. if (s2 == SpecialValueType.INF_PLUS || s2 == SpecialValueType.INF_MINUS)
  2571. {
  2572. //Set NaN and return.
  2573. SetNaN();
  2574. return true;
  2575. }
  2576. if (n2.Sign)
  2577. {
  2578. if (s1 == SpecialValueType.INF_PLUS)
  2579. {
  2580. SetInfMinus();
  2581. return true;
  2582. }
  2583. SetInfPlus();
  2584. return true;
  2585. }
  2586. //Keep inf
  2587. return true;
  2588. }
  2589. if (s2 == SpecialValueType.ZERO)
  2590. {
  2591. if (s1 == SpecialValueType.ZERO)
  2592. {
  2593. SetNaN();
  2594. return true;
  2595. }
  2596. if (Sign)
  2597. {
  2598. SetInfMinus();
  2599. return true;
  2600. }
  2601. SetInfPlus();
  2602. return true;
  2603. }
  2604. }
  2605. return false;
  2606. }
  2607. //****************** Internal helper functions *****************
  2608. /// <summary>
  2609. /// Used for fixed point speed-ups (where the extra precision is not required). Note that Denormalised
  2610. /// floats break the assumptions that underly Add() and Sub(), so they can only be used for multiplication
  2611. /// </summary>
  2612. /// <param name="targetExponent"></param>
  2613. private void Denormalise(int targetExponent)
  2614. {
  2615. int diff = targetExponent - exponent;
  2616. if (diff <= 0) return;
  2617. //This only works to reduce the precision, so if the difference implies an increase, we can't do anything.
  2618. mantissa.RSH(diff);
  2619. exponent += diff;
  2620. }
  2621. /// <summary>
  2622. /// The binary logarithm, log2(x) - for precisions above 1000 bits, use Log() and convert the base.
  2623. /// </summary>
  2624. private void Log2()
  2625. {
  2626. if (scratch.Precision.NumBits != mantissa.Precision.NumBits)
  2627. {
  2628. scratch = new BigInt(mantissa.Precision);
  2629. }
  2630. int bits = mantissa.Precision.NumBits;
  2631. BigFloat temp = new BigFloat(this);
  2632. BigFloat result = new BigFloat(exponent, mantissa.Precision);
  2633. BigFloat pow2 = new BigFloat(1, mantissa.Precision);
  2634. temp.exponent = 0;
  2635. int bitsCalculated = 0;
  2636. while (bitsCalculated < bits)
  2637. {
  2638. int i;
  2639. for (i = 0; (temp.exponent == 0); i++)
  2640. {
  2641. temp.mantissa.SquareHiFast(scratch);
  2642. int shift = temp.mantissa.Normalise();
  2643. temp.exponent += 1 - shift;
  2644. if (i + bitsCalculated >= bits) break;
  2645. }
  2646. pow2.MulPow2(-i);
  2647. result.Add(pow2);
  2648. temp.exponent = 0;
  2649. bitsCalculated += i;
  2650. }
  2651. this.Assign(result);
  2652. }
  2653. /// <summary>
  2654. /// Tried the newton method for logs, but the exponential function is too slow to do it.
  2655. /// </summary>
  2656. private void LogNewton()
  2657. {
  2658. if (mantissa.IsZero() || mantissa.Sign)
  2659. {
  2660. return;
  2661. }
  2662. //Compute ln2.
  2663. if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
  2664. {
  2665. CalculateLog2(mantissa.Precision.NumBits);
  2666. }
  2667. int numBits = mantissa.Precision.NumBits;
  2668. //Use inverse exp function with Newton's method.
  2669. BigFloat xn = new BigFloat(this);
  2670. BigFloat oldExponent = new BigFloat(xn.exponent, mantissa.Precision);
  2671. xn.exponent = 0;
  2672. this.exponent = 0;
  2673. //Hack to subtract 1
  2674. xn.mantissa.ClearBit(numBits - 1);
  2675. //x0 = (x - 1) * log2 - this is a straight line fit between log(1) = 0 and log(2) = ln2
  2676. xn.Mul(ln2cache);
  2677. //x0 = (x - 1) * log2 + C - this corrects for minimum error over the range.
  2678. xn.Add(logNewtonConstant);
  2679. BigFloat term = new BigFloat(mantissa.Precision);
  2680. BigFloat one = new BigFloat(1, mantissa.Precision);
  2681. int precision = 32;
  2682. int normalPrecision = mantissa.Precision.NumBits;
  2683. int iterations = 0;
  2684. while (true)
  2685. {
  2686. term.Assign(xn);
  2687. term.mantissa.Sign = true;
  2688. term.Exp(precision);
  2689. term.Mul(this);
  2690. term.Sub(one);
  2691. iterations++;
  2692. if (term.exponent < -((precision >> 1) - 4))
  2693. {
  2694. if (precision == normalPrecision)
  2695. {
  2696. if (term.exponent < -(precision - 4)) break;
  2697. }
  2698. else
  2699. {
  2700. precision = precision << 1;
  2701. if (precision > normalPrecision) precision = normalPrecision;
  2702. }
  2703. }
  2704. xn.Add(term);
  2705. }
  2706. //log(2^n*s) = log(2^n) + log(s) = nlog(2) + log(s)
  2707. term.Assign(ln2cache);
  2708. term.Mul(oldExponent);
  2709. this.Assign(xn);
  2710. this.Add(term);
  2711. }
  2712. /// <summary>
  2713. /// Log(x) implemented as an Arithmetic-Geometric Mean. Fast for high precisions.
  2714. /// </summary>
  2715. private void LogAGM1()
  2716. {
  2717. if (mantissa.IsZero() || mantissa.Sign)
  2718. {
  2719. return;
  2720. }
  2721. //Compute ln2.
  2722. if (ln2cache == null || mantissa.Precision.NumBits > ln2cache.mantissa.Precision.NumBits)
  2723. {
  2724. CalculateLog2(mantissa.Precision.NumBits);
  2725. }
  2726. //Compute ln(x) using AGM formula
  2727. //1. Re-write the input as 2^n * (0.5 <= x < 1)
  2728. int power2 = exponent + 1;
  2729. exponent = -1;
  2730. //BigFloat res = new BigFloat(firstAGMcache);
  2731. BigFloat a0 = new BigFloat(1, mantissa.Precision);
  2732. BigFloat b0 = new BigFloat(pow10cache);
  2733. b0.Mul(this);
  2734. BigFloat r = R(a0, b0);
  2735. this.Assign(firstAGMcache);
  2736. this.Sub(r);
  2737. a0.Assign(ln2cache);
  2738. a0.Mul(new BigFloat(power2, mantissa.Precision));
  2739. this.Add(a0);
  2740. }
  2741. private void Exp(int numBits)
  2742. {
  2743. if (IsSpecialValue)
  2744. {
  2745. if (SpecialValue == SpecialValueType.ZERO)
  2746. {
  2747. //e^0 = 1
  2748. exponent = 0;
  2749. mantissa.SetHighDigit(0x80000000);
  2750. }
  2751. else if (SpecialValue == SpecialValueType.INF_MINUS)
  2752. {
  2753. //e^-inf = 0
  2754. SetZero();
  2755. }
  2756. return;
  2757. }
  2758. PrecisionSpec prec = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
  2759. numBits = prec.NumBits;
  2760. if (scratch.Precision.NumBits != prec.NumBits)
  2761. {
  2762. scratch = new BigInt(prec);
  2763. }
  2764. if (inverseFactorialCache == null || invFactorialCutoff < numBits)
  2765. {
  2766. CalculateFactorials(numBits);
  2767. }
  2768. //let x = 1 * 'this'.mantissa (i.e. 1 <= x < 2)
  2769. //exp(2^n * x) = e^(2^n * x) = (e^x)^2n = exp(x)^2n
  2770. int oldExponent = 0;
  2771. if (exponent > -4)
  2772. {
  2773. oldExponent = exponent + 4;
  2774. exponent = -4;
  2775. }
  2776. BigFloat thisSave = new BigFloat(this, prec);
  2777. BigFloat temp = new BigFloat(1, prec);
  2778. BigFloat temp2 = new BigFloat(this, prec);
  2779. BigFloat res = new BigFloat(1, prec);
  2780. int length = inverseFactorialCache.Length;
  2781. int iterations;
  2782. for (int i = 1; i < length; i++)
  2783. {
  2784. //temp = x^i
  2785. temp.Mul(thisSave);
  2786. temp2.Assign(inverseFactorialCache[i]);
  2787. temp2.Mul(temp);
  2788. if (temp2.exponent < -(numBits + 4)) { iterations = i; break; }
  2789. res.Add(temp2);
  2790. }
  2791. //res = exp(x)
  2792. //Now... x^(2^n) = (x^2)^(2^(n - 1))
  2793. for (int i = 0; i < oldExponent; i++)
  2794. {
  2795. res.mantissa.SquareHiFast(scratch);
  2796. int shift = res.mantissa.Normalise();
  2797. res.exponent = res.exponent << 1;
  2798. res.exponent += 1 - shift;
  2799. }
  2800. //Deal with +/- inf
  2801. if (res.exponent == Int32.MaxValue)
  2802. {
  2803. res.mantissa.Zero();
  2804. }
  2805. Assign(res);
  2806. }
  2807. /// <summary>
  2808. /// Calculates ln(2) and returns -10^(n/2 + a bit) for reuse, using the AGM method as described in
  2809. /// http://lacim.uqam.ca/~plouffe/articles/log2.pdf
  2810. /// </summary>
  2811. /// <param name="numBits"></param>
  2812. /// <returns></returns>
  2813. private static void CalculateLog2(int numBits)
  2814. {
  2815. //Use the AGM method formula to get log2 to N digits.
  2816. //R(a0, b0) = 1 / (1 - Sum(2^-n*(an^2 - bn^2)))
  2817. //log(1/2) = R(1, 10^-n) - R(1, 10^-n/2)
  2818. PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
  2819. PrecisionSpec extendedPres = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);
  2820. BigFloat a0 = new BigFloat(1, extendedPres);
  2821. BigFloat b0 = TenPow(-(int)((double)((numBits >> 1) + 2) * 0.302), extendedPres);
  2822. BigFloat pow10saved = new BigFloat(b0);
  2823. BigFloat firstAGMcacheSaved = new BigFloat(extendedPres);
  2824. //save power of 10 (in normal precision)
  2825. pow10cache = new BigFloat(b0, normalPres);
  2826. ln2cache = R(a0, b0);
  2827. //save the first half of the log calculation
  2828. firstAGMcache = new BigFloat(ln2cache, normalPres);
  2829. firstAGMcacheSaved.Assign(ln2cache);
  2830. b0.MulPow2(-1);
  2831. ln2cache.Sub(R(a0, b0));
  2832. //Convert to log(2)
  2833. ln2cache.mantissa.Sign = false;
  2834. //Save magic constant for newton log
  2835. //First guess in range 1 <= x < 2 is x0 = ln2 * (x - 1) + C
  2836. logNewtonConstant = new BigFloat(ln2cache);
  2837. logNewtonConstant.Mul(new BigFloat(3, extendedPres));
  2838. logNewtonConstant.exponent--;
  2839. logNewtonConstant.Sub(new BigFloat(1, extendedPres));
  2840. logNewtonConstant = new BigFloat(logNewtonConstant, normalPres);
  2841. //Save the inverse.
  2842. log2ecache = new BigFloat(ln2cache);
  2843. log2ecache = new BigFloat(log2ecache.Reciprocal(), normalPres);
  2844. //Now cache log10
  2845. //Because the log functions call this function to the precision to which they
  2846. //are called, we cannot call them without causing an infinite loop, so we need
  2847. //to inline the code.
  2848. log10recip = new BigFloat(10, extendedPres);
  2849. {
  2850. int power2 = log10recip.exponent + 1;
  2851. log10recip.exponent = -1;
  2852. //BigFloat res = new BigFloat(firstAGMcache);
  2853. BigFloat ax = new BigFloat(1, extendedPres);
  2854. BigFloat bx = new BigFloat(pow10saved);
  2855. bx.Mul(log10recip);
  2856. BigFloat r = R(ax, bx);
  2857. log10recip.Assign(firstAGMcacheSaved);
  2858. log10recip.Sub(r);
  2859. ax.Assign(ln2cache);
  2860. ax.Mul(new BigFloat(power2, log10recip.mantissa.Precision));
  2861. log10recip.Add(ax);
  2862. }
  2863. log10recip = log10recip.Reciprocal();
  2864. log10recip = new BigFloat(log10recip, normalPres);
  2865. //Trim to n bits
  2866. ln2cache = new BigFloat(ln2cache, normalPres);
  2867. }
  2868. private static BigFloat TenPow(int power, PrecisionSpec precision)
  2869. {
  2870. BigFloat acc = new BigFloat(1, precision);
  2871. BigFloat temp = new BigFloat(1, precision);
  2872. int powerTemp = power;
  2873. BigFloat multiplierToUse = new BigFloat(10, precision);
  2874. if (power < 0)
  2875. {
  2876. multiplierToUse = multiplierToUse.Reciprocal();
  2877. powerTemp = -power;
  2878. }
  2879. //Fast power function
  2880. while (powerTemp != 0)
  2881. {
  2882. temp.Mul(multiplierToUse);
  2883. multiplierToUse.Assign(temp);
  2884. if ((powerTemp & 1) != 0)
  2885. {
  2886. acc.Mul(temp);
  2887. }
  2888. powerTemp >>= 1;
  2889. }
  2890. return acc;
  2891. }
  2892. private static BigFloat R(BigFloat a0, BigFloat b0)
  2893. {
  2894. //Precision extend taken out.
  2895. int bits = a0.mantissa.Precision.NumBits;
  2896. PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);
  2897. BigFloat an = new BigFloat(a0, extendedPres);
  2898. BigFloat bn = new BigFloat(b0, extendedPres);
  2899. BigFloat sum = new BigFloat(extendedPres);
  2900. BigFloat term = new BigFloat(extendedPres);
  2901. BigFloat temp1 = new BigFloat(extendedPres);
  2902. BigFloat one = new BigFloat(1, extendedPres);
  2903. int iteration = 0;
  2904. for (iteration = 0; ; iteration++)
  2905. {
  2906. //Get the sum term for this iteration.
  2907. term.Assign(an);
  2908. term.Mul(an);
  2909. temp1.Assign(bn);
  2910. temp1.Mul(bn);
  2911. //term = an^2 - bn^2
  2912. term.Sub(temp1);
  2913. //term = 2^(n-1) * (an^2 - bn^2)
  2914. term.exponent += iteration - 1;
  2915. sum.Add(term);
  2916. if (term.exponent < -(bits - 8)) break;
  2917. //Calculate the new AGM estimates.
  2918. temp1.Assign(an);
  2919. an.Add(bn);
  2920. //a(n+1) = (an + bn) / 2
  2921. an.MulPow2(-1);
  2922. //b(n+1) = sqrt(an*bn)
  2923. bn.Mul(temp1);
  2924. bn.Sqrt();
  2925. }
  2926. one.Sub(sum);
  2927. one = one.Reciprocal();
  2928. return new BigFloat(one, a0.mantissa.Precision);
  2929. }
  2930. private static void CalculateFactorials(int numBits)
  2931. {
  2932. System.Collections.Generic.List<BigFloat> list = new System.Collections.Generic.List<BigFloat>(64);
  2933. System.Collections.Generic.List<BigFloat> list2 = new System.Collections.Generic.List<BigFloat>(64);
  2934. PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);
  2935. PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
  2936. BigFloat factorial = new BigFloat(1, extendedPrecision);
  2937. BigFloat reciprocal;
  2938. //Calculate e while we're at it
  2939. BigFloat e = new BigFloat(1, extendedPrecision);
  2940. list.Add(new BigFloat(factorial, normalPrecision));
  2941. for (int i = 1; i < Int32.MaxValue; i++)
  2942. {
  2943. BigFloat number = new BigFloat(i, extendedPrecision);
  2944. factorial.Mul(number);
  2945. if (factorial.exponent > numBits) break;
  2946. list2.Add(new BigFloat(factorial, normalPrecision));
  2947. reciprocal = factorial.Reciprocal();
  2948. e.Add(reciprocal);
  2949. list.Add(new BigFloat(reciprocal, normalPrecision));
  2950. }
  2951. //Set the cached static values.
  2952. inverseFactorialCache = list.ToArray();
  2953. factorialCache = list2.ToArray();
  2954. invFactorialCutoff = numBits;
  2955. eCache = new BigFloat(e, normalPrecision);
  2956. eRCPCache = new BigFloat(e.Reciprocal(), normalPrecision);
  2957. }
  2958. private static void CalculateEOnly(int numBits)
  2959. {
  2960. PrecisionSpec extendedPrecision = new PrecisionSpec(numBits + 1, PrecisionSpec.BaseType.BIN);
  2961. PrecisionSpec normalPrecision = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
  2962. int iExponent = (int)(Math.Sqrt(numBits));
  2963. BigFloat factorial = new BigFloat(1, extendedPrecision);
  2964. BigFloat constant = new BigFloat(1, extendedPrecision);
  2965. constant.exponent -= iExponent;
  2966. BigFloat numerator = new BigFloat(constant);
  2967. BigFloat reciprocal;
  2968. //Calculate the 2^iExponent th root of e
  2969. BigFloat e = new BigFloat(1, extendedPrecision);
  2970. int i;
  2971. for (i = 1; i < Int32.MaxValue; i++)
  2972. {
  2973. BigFloat number = new BigFloat(i, extendedPrecision);
  2974. factorial.Mul(number);
  2975. reciprocal = factorial.Reciprocal();
  2976. reciprocal.Mul(numerator);
  2977. if (-reciprocal.exponent > numBits) break;
  2978. e.Add(reciprocal);
  2979. numerator.Mul(constant);
  2980. System.GC.Collect();
  2981. }
  2982. for (i = 0; i < iExponent; i++)
  2983. {
  2984. numerator.Assign(e);
  2985. e.Mul(numerator);
  2986. }
  2987. //Set the cached static values.
  2988. eCache = new BigFloat(e, normalPrecision);
  2989. eRCPCache = new BigFloat(e.Reciprocal(), normalPrecision);
  2990. }
  2991. /// <summary>
  2992. /// Uses the Gauss-Legendre formula for pi
  2993. /// Taken from http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm
  2994. /// </summary>
  2995. /// <param name="numBits"></param>
  2996. private static void CalculatePi(int numBits)
  2997. {
  2998. int bits = numBits + 32;
  2999. //Precision extend taken out.
  3000. PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
  3001. PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);
  3002. if (scratch.Precision.NumBits != bits)
  3003. {
  3004. scratch = new BigInt(extendedPres);
  3005. }
  3006. //a0 = 1
  3007. BigFloat an = new BigFloat(1, extendedPres);
  3008. //b0 = 1/sqrt(2)
  3009. BigFloat bn = new BigFloat(2, extendedPres);
  3010. bn.Sqrt();
  3011. bn.exponent--;
  3012. //to = 1/4
  3013. BigFloat tn = new BigFloat(1, extendedPres);
  3014. tn.exponent -= 2;
  3015. int pn = 0;
  3016. BigFloat anTemp = new BigFloat(extendedPres);
  3017. int iteration = 0;
  3018. int cutoffBits = numBits >> 5;
  3019. for (iteration = 0; ; iteration++)
  3020. {
  3021. //Save a(n)
  3022. anTemp.Assign(an);
  3023. //Calculate new an
  3024. an.Add(bn);
  3025. an.exponent--;
  3026. //Calculate new bn
  3027. bn.Mul(anTemp);
  3028. bn.Sqrt();
  3029. //Calculate new tn
  3030. anTemp.Sub(an);
  3031. anTemp.mantissa.SquareHiFast(scratch);
  3032. anTemp.exponent += anTemp.exponent + pn + 1 - anTemp.mantissa.Normalise();
  3033. tn.Sub(anTemp);
  3034. anTemp.Assign(an);
  3035. anTemp.Sub(bn);
  3036. if (anTemp.exponent < -(bits - cutoffBits)) break;
  3037. //New pn
  3038. pn++;
  3039. }
  3040. an.Add(bn);
  3041. an.mantissa.SquareHiFast(scratch);
  3042. an.exponent += an.exponent + 1 - an.mantissa.Normalise();
  3043. tn.exponent += 2;
  3044. an.Div(tn);
  3045. pi = new BigFloat(an, normalPres);
  3046. piBy2 = new BigFloat(pi);
  3047. piBy2.exponent--;
  3048. twoPi = new BigFloat(pi, normalPres);
  3049. twoPi.exponent++;
  3050. piRecip = new BigFloat(an.Reciprocal(), normalPres);
  3051. twoPiRecip = new BigFloat(piRecip);
  3052. twoPiRecip.exponent--;
  3053. //1/3 is going to be useful for sin.
  3054. threeRecip = new BigFloat((new BigFloat(3, extendedPres)).Reciprocal(), normalPres);
  3055. }
  3056. /// <summary>
  3057. /// Calculates the odd reciprocals of the natural numbers (for atan series)
  3058. /// </summary>
  3059. /// <param name="numBits"></param>
  3060. /// <param name="terms"></param>
  3061. private static void CalculateReciprocals(int numBits, int terms)
  3062. {
  3063. int bits = numBits + 32;
  3064. PrecisionSpec extendedPres = new PrecisionSpec(bits, PrecisionSpec.BaseType.BIN);
  3065. PrecisionSpec normalPres = new PrecisionSpec(numBits, PrecisionSpec.BaseType.BIN);
  3066. System.Collections.Generic.List<BigFloat> list = new System.Collections.Generic.List<BigFloat>(terms);
  3067. for (int i = 0; i < terms; i++)
  3068. {
  3069. BigFloat term = new BigFloat(i*2 + 1, extendedPres);
  3070. list.Add(new BigFloat(term.Reciprocal(), normalPres));
  3071. }
  3072. reciprocals = list.ToArray();
  3073. }
  3074. /// <summary>
  3075. /// Does decimal rounding, for numbers without E notation.
  3076. /// </summary>
  3077. /// <param name="input"></param>
  3078. /// <param name="places"></param>
  3079. /// <returns></returns>
  3080. private static string RoundString(string input, int places)
  3081. {
  3082. if (places <= 0) return input;
  3083. string trim = input.Trim();
  3084. char[] digits = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9'};
  3085. /*
  3086. for (int i = 1; i <= places; i++)
  3087. {
  3088. //Skip decimal points.
  3089. if (trim[trim.Length - i] == '.')
  3090. {
  3091. places++;
  3092. continue;
  3093. }
  3094. int index = Array.IndexOf(digits, trim[trim.Length - i]);
  3095. if (index < 0) return input;
  3096. value += ten * index;
  3097. ten *= 10;
  3098. }
  3099. * */
  3100. //Look for a decimal point
  3101. string decimalPoint = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator;
  3102. int indexPoint = trim.LastIndexOf(decimalPoint);
  3103. if (indexPoint < 0)
  3104. {
  3105. //We can't modify a string which doesn't have a decimal point.
  3106. return trim;
  3107. }
  3108. int trimPoint = trim.Length - places;
  3109. if (trimPoint < indexPoint) trimPoint = indexPoint;
  3110. bool roundDown = false;
  3111. if (trim[trimPoint] == '.')
  3112. {
  3113. if (trimPoint + 1 >= trim.Length)
  3114. {
  3115. roundDown = true;
  3116. }
  3117. else
  3118. {
  3119. int digit = Array.IndexOf(digits, trim[trimPoint + 1]);
  3120. if (digit < 5) roundDown = true;
  3121. }
  3122. }
  3123. else
  3124. {
  3125. int digit = Array.IndexOf(digits, trim[trimPoint]);
  3126. if (digit < 5) roundDown = true;
  3127. }
  3128. string output;
  3129. //Round down - just return a new string without the extra digits.
  3130. if (roundDown)
  3131. {
  3132. if (RoundingMode == RoundingModeType.EXACT)
  3133. {
  3134. return trim.Substring(0, trimPoint);
  3135. }
  3136. else
  3137. {
  3138. char[] trimChars = { '0' };
  3139. output = trim.Substring(0, trimPoint).TrimEnd(trimChars);
  3140. trimChars[0] = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator[0];
  3141. return output.TrimEnd(trimChars);
  3142. }
  3143. }
  3144. //Round up - bit more complicated.
  3145. char [] arrayOutput = trim.ToCharArray();//0, trimPoint);
  3146. //Now, we round going from the back to the front.
  3147. int j;
  3148. for (j = trimPoint - 1; j >= 0; j--)
  3149. {
  3150. int index = Array.IndexOf(digits, arrayOutput[j]);
  3151. //Skip decimal points etc...
  3152. if (index < 0) continue;
  3153. if (index < 9)
  3154. {
  3155. arrayOutput[j] = digits[index + 1];
  3156. break;
  3157. }
  3158. else
  3159. {
  3160. arrayOutput[j] = digits[0];
  3161. }
  3162. }
  3163. output = new string(arrayOutput);
  3164. if (j < 0)
  3165. {
  3166. //Need to add a new digit.
  3167. output = String.Format("{0}{1}", "1", output);
  3168. }
  3169. if (RoundingMode == RoundingModeType.EXACT)
  3170. {
  3171. return output;
  3172. }
  3173. else
  3174. {
  3175. char[] trimChars = { '0' };
  3176. output = output.TrimEnd(trimChars);
  3177. trimChars[0] = System.Globalization.CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator[0];
  3178. return output.TrimEnd(trimChars);
  3179. }
  3180. }
  3181. //***************************** Data *****************************
  3182. //Side node - this way of doing things is far from optimal, both in terms of memory use and performance.
  3183. private ExponentAdaptor exponent;
  3184. private BigInt mantissa;
  3185. /// <summary>
  3186. /// Storage area for calculations.
  3187. /// </summary>
  3188. private static BigInt scratch;
  3189. private static BigFloat ln2cache; //Value of ln(2)
  3190. private static BigFloat log2ecache; //Value of log2(e) = 1/ln(2)
  3191. private static BigFloat pow10cache; //Cached power of 10 for AGM log calculation
  3192. private static BigFloat log10recip; //1/ln(10)
  3193. private static BigFloat firstAGMcache; //Cached half of AGM operation.
  3194. private static BigFloat[] factorialCache; //The values of n!
  3195. private static BigFloat[] inverseFactorialCache; //Values of 1/n! up to 2^-m where m = invFactorialCutoff (below)
  3196. private static int invFactorialCutoff; //The number of significant bits for the cutoff of the inverse factorials.
  3197. private static BigFloat eCache; //Value of e cached to invFactorialCutoff bits
  3198. private static BigFloat eRCPCache; //Reciprocal of e
  3199. private static BigFloat logNewtonConstant; //1.5*ln(2) - 1
  3200. private static BigFloat pi; //pi
  3201. private static BigFloat piBy2; //pi/2
  3202. private static BigFloat twoPi; //2*pi
  3203. private static BigFloat piRecip; //1/pi
  3204. private static BigFloat twoPiRecip; //1/2*pi
  3205. private static BigFloat threeRecip; //1/3
  3206. private static BigFloat[] reciprocals; //1/x
  3207. /// <summary>
  3208. /// The number of decimal digits to round the output of ToString() by
  3209. /// </summary>
  3210. public static int RoundingDigits { get; set; }
  3211. /// <summary>
  3212. /// The way in which ToString() should deal with insignificant trailing zeroes
  3213. /// </summary>
  3214. public static RoundingModeType RoundingMode { get; set; }
  3215. }
  3216. }
旧リポジトリブラウザで表示