The Labs \ Source Viewer \ SSCLI \ System.Xml.Xsl \ BigNumber

  1. //------------------------------------------------------------------------------
  2. // <copyright file="XPathConvert.cs" company="Microsoft">
  3. //
  4. // Copyright (c) 2006 Microsoft Corporation. All rights reserved.
  5. //
  6. // The use and distribution terms for this software are contained in the file
  7. // named license.txt, which can be found in the root of this distribution.
  8. // By using this software in any fashion, you are agreeing to be bound by the
  9. // terms of this license.
  10. //
  11. // You must not remove this notice, or any other, from this software.
  12. //
  13. // </copyright>
  14. //------------------------------------------------------------------------------
  15. /**
  16. * Routines used to manipulate IEEE 754 double-precision numbers, taken from JScript codebase.
  17. *
  18. * Define NOPARSE if you do not need FloatingDecimal -> double conversions
  19. *
  20. */
  21. using System.Diagnostics;
  22. using System.Globalization;
  23. namespace System.Xml.Xsl
  24. {
  25.    
  26. /**
  27.     * Converts according to XPath/XSLT rules.
  28.     */   
  29.     static internal class XPathConvert
  30.     {
  31.        
  32.         public static uint DblHi(double dbl)
  33.         {
  34.             return (uint)(BitConverter.DoubleToInt64Bits(dbl) >> 32);
  35.         }
  36.        
  37.         public static uint DblLo(double dbl)
  38.         {
  39.             return (uint)BitConverter.DoubleToInt64Bits(dbl);
  40.         }
  41.        
  42.         // Returns true if value is infinite or NaN (exponent bits are all ones)
  43.         public static bool IsSpecial(double dbl)
  44.         {
  45.             return 0 == (~DblHi(dbl) & 2146435072);
  46.         }
  47.        
  48.         #if DEBUG
  49.         // Returns the next representable neighbor of x in the direction toward y
  50.         public static double NextAfter(double x, double y)
  51.         {
  52.             long bits;
  53.            
  54.             if (Double.IsNaN(x)) {
  55.                 return x;
  56.             }
  57.             if (Double.IsNaN(y)) {
  58.                 return y;
  59.             }
  60.             if (x == y) {
  61.                 return y;
  62.             }
  63.             if (x == 0) {
  64.                 bits = BitConverter.DoubleToInt64Bits(y) & 1l << 63;
  65.                 return BitConverter.Int64BitsToDouble(bits | 1);
  66.             }
  67.            
  68.             // At this point x!=y, and x!=0. x can be treated as a 64bit
  69.             // integer in sign/magnitude representation. To get the next
  70.             // representable neighbor we add or subtract one from this
  71.             // integer.
  72.            
  73.             bits = BitConverter.DoubleToInt64Bits(x);
  74.             if (0 < x && x < y || 0 > x && x > y) {
  75.                 bits++;
  76.             }
  77.             else {
  78.                 bits--;
  79.             }
  80.             return BitConverter.Int64BitsToDouble(bits);
  81.         }
  82.        
  83.         public static double Succ(double x)
  84.         {
  85.             return NextAfter(x, Double.PositiveInfinity);
  86.         }
  87.        
  88.         public static double Pred(double x)
  89.         {
  90.             return NextAfter(x, Double.NegativeInfinity);
  91.         }
  92.         #endif
  93.        
  94.         // Small powers of ten. These are all the powers of ten that have an exact
  95.         // representation in IEEE double precision format.
  96.         public static readonly double[] C10toN = {1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0, 1000000.0, 10000000.0, 100000000.0, 1000000000.0,
  97.         10000000000.0, 100000000000.0, 1000000000000.0, 10000000000000.0, 100000000000000.0, 1E+15, 1E+16, 1E+17, 1E+18, 1E+19,
  98.         1E+20, 1E+21, 1E+22};
  99.        
  100.         // Returns 1 if argument is non-zero, and 0 otherwise
  101.         public static uint NotZero(uint u)
  102.         {
  103.             return 0 != u ? 1u : 0u;
  104.         }
  105.        
  106. /*  ----------------------------------------------------------------------------
  107.             AddU()
  108.             Add two unsigned ints and return the carry bit.
  109.         */       
  110.         public static uint AddU(ref uint u1, uint u2)
  111.         {
  112.             u1 += u2;
  113.             return u1 < u2 ? 1u : 0u;
  114.         }
  115.        
  116. /*  ----------------------------------------------------------------------------
  117.             MulU()
  118.             Multiply two unsigned ints. Return the low uint and fill uHi with
  119.             the high uint.
  120.         */       
  121.         public static uint MulU(uint u1, uint u2, out uint uHi)
  122.         {
  123.             ulong result = (ulong)u1 * u2;
  124.             uHi = (uint)(result >> 32);
  125.             return (uint)result;
  126.         }
  127.        
  128. /*  ----------------------------------------------------------------------------
  129.             CbitZeroLeft()
  130.             Return a count of the number of leading 0 bits in u.
  131.         */       
  132.         public static int CbitZeroLeft(uint u)
  133.         {
  134.             int cbit = 0;
  135.            
  136.             if (0 == (u & 4294901760u)) {
  137.                 cbit += 16;
  138.                 u <<= 16;
  139.             }
  140.             if (0 == (u & 4278190080u)) {
  141.                 cbit += 8;
  142.                 u <<= 8;
  143.             }
  144.             if (0 == (u & 4026531840u)) {
  145.                 cbit += 4;
  146.                 u <<= 4;
  147.             }
  148.             if (0 == (u & 3221225472u)) {
  149.                 cbit += 2;
  150.                 u <<= 2;
  151.             }
  152.             if (0 == (u & 2147483648u)) {
  153.                 cbit += 1;
  154.                 u <<= 1;
  155.             }
  156.             Debug.Assert(0 != (u & 2147483648u));
  157.            
  158.             return cbit;
  159.         }
  160.        
  161. /*  ----------------------------------------------------------------------------
  162.             IsInteger()
  163.             If dbl is a whole number in the range of INT_MIN to INT_MAX, return true
  164.             and the integer in value.  Otherwise, return false.
  165.         */       
  166.         public static bool IsInteger(double dbl, out int value)
  167.         {
  168.             if (!IsSpecial(dbl)) {
  169.                 int i = (int)dbl;
  170.                 double dblRound = (double)i;
  171.                
  172.                 if (dbl == dblRound) {
  173.                     value = i;
  174.                     return true;
  175.                 }
  176.             }
  177.             value = 0;
  178.             return false;
  179.         }
  180.        
  181. /**
  182.         * Implementation of a big floating point number used to ensure adequate
  183.         * precision when performing calculations.
  184.         *
  185.         * Hungarian: num
  186.         *
  187.         */       
  188.         private struct BigNumber
  189.         {
  190.             private uint u0;
  191.             private uint u1;
  192.             private uint u2;
  193.             private int exp;
  194.            
  195.             // This is a bound on the absolute value of the error. It is based at
  196.             // one bit before the least significant bit of u0.
  197.             private uint error;
  198.            
  199.             public uint Error {
  200.                 get { return error; }
  201.             }
  202.            
  203.             public BigNumber(uint u0, uint u1, uint u2, int exp, uint error)
  204.             {
  205.                 this.u0 = u0;
  206.                 this.u1 = u1;
  207.                 this.u2 = u2;
  208.                 this.exp = exp;
  209.                 this.error = error;
  210.             }
  211.            
  212.             #if !NOPARSE || DEBUG
  213.             // Set the value from floating decimal
  214.             public BigNumber(FloatingDecimal dec)
  215.             {
  216.                 Debug.Assert(dec.MantissaSize > 0);
  217.                
  218.                 int ib;
  219.                 int exponent;
  220.                 int mantissaSize;
  221.                 int wT;
  222.                 uint uExtra;
  223.                 BigNumber[] TenPowers;
  224.                
  225.                 ib = 0;
  226.                 exponent = dec.Exponent;
  227.                 mantissaSize = dec.MantissaSize;
  228.                
  229.                 // Record the first digit
  230.                 Debug.Assert(dec[ib] > 0 && dec[ib] <= 9);
  231.                 this.u2 = (uint)(dec[ib]) << 28;
  232.                 this.u1 = 0;
  233.                 this.u0 = 0;
  234.                 this.exp = 4;
  235.                 this.error = 0;
  236.                 exponent--;
  237.                 Normalize();
  238.                
  239.                 while (++ib < mantissaSize) {
  240.                     Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
  241.                     uExtra = MulTenAdd(dec[ib]);
  242.                     exponent--;
  243.                     if (0 != uExtra) {
  244.                         // We've filled up our precision.
  245.                         Round(uExtra);
  246.                         if (ib < mantissaSize - 1) {
  247.                             // There are more digits, so add another error bit just for
  248.                             // safety's sake.
  249.                             this.error++;
  250.                         }
  251.                         break;
  252.                     }
  253.                 }
  254.                
  255.                 // Now multiply by 10^exp
  256.                 if (0 == exponent) {
  257.                     return;
  258.                 }
  259.                
  260.                 if (exponent < 0) {
  261.                     TenPowers = TenPowersNeg;
  262.                     exponent = -exponent;
  263.                 }
  264.                 else {
  265.                     TenPowers = TenPowersPos;
  266.                 }
  267.                
  268.                 Debug.Assert(exponent > 0 && exponent < 512);
  269.                 wT = exponent & 31;
  270.                 if (wT > 0) {
  271.                     Mul(ref TenPowers[wT - 1]);
  272.                 }
  273.                
  274.                 wT = (exponent >> 5) & 15;
  275.                 if (wT > 0) {
  276.                     Mul(ref TenPowers[wT + 30]);
  277.                 }
  278.             }
  279.            
  280.             // Multiply by ten and add a base 10 digit.
  281.             unsafe private uint MulTenAdd(uint digit)
  282.             {
  283.                 Debug.Assert(digit <= 9);
  284.                 Debug.Assert(0 != (this.u2 & 2147483648u));
  285.                
  286.                 // First "multipy" by eight
  287.                 this.exp += 3;
  288.                 Debug.Assert(this.exp >= 4);
  289.                
  290.                 // Initialize the carry values based on digit and exp.
  291.                 uint* rgu = stackalloc uint[5];
  292.                 for (int i = 0; i < 5; i++) {
  293.                     rgu[i] = 0;
  294.                 }
  295.                 if (0 != digit) {
  296.                     int idx = 3 - (this.exp >> 5);
  297.                     if (idx < 0) {
  298.                         rgu[0] = 1;
  299.                     }
  300.                     else {
  301.                         int ibit = this.exp & 31;
  302.                         if (ibit < 4) {
  303.                             rgu[idx + 1] = digit >> ibit;
  304.                             if (ibit > 0) {
  305.                                 rgu[idx] = digit << (32 - ibit);
  306.                             }
  307.                         }
  308.                         else {
  309.                             rgu[idx] = digit << (32 - ibit);
  310.                         }
  311.                     }
  312.                 }
  313.                
  314.                 // Shift and add to multiply by ten.
  315.                 rgu[1] += AddU(ref rgu[0], this.u0 << 30);
  316.                 rgu[2] += AddU(ref this.u0, (this.u0 >> 2) + (this.u1 << 30));
  317.                 if (0 != rgu[1]) {
  318.                     rgu[2] += AddU(ref this.u0, rgu[1]);
  319.                 }
  320.                 rgu[3] += AddU(ref this.u1, (this.u1 >> 2) + (this.u2 << 30));
  321.                 if (0 != rgu[2]) {
  322.                     rgu[3] += AddU(ref this.u1, rgu[2]);
  323.                 }
  324.                 rgu[4] = AddU(ref this.u2, (this.u2 >> 2) + rgu[3]);
  325.                
  326.                 // Handle the final carry.
  327.                 if (0 != rgu[4]) {
  328.                     Debug.Assert(1 == rgu[4]);
  329.                     rgu[0] = (rgu[0] >> 1) | (rgu[0] & 1) | (this.u0 << 31);
  330.                     this.u0 = (this.u0 >> 1) | (this.u1 << 31);
  331.                     this.u1 = (this.u1 >> 1) | (this.u2 << 31);
  332.                     this.u2 = (this.u2 >> 1) | 2147483648u;
  333.                     this.exp++;
  334.                 }
  335.                
  336.                 return rgu[0];
  337.             }
  338.            
  339.             // Round based on the given extra data using IEEE round to nearest rule.
  340.             private void Round(uint uExtra)
  341.             {
  342.                 if (0 == (uExtra & 2147483648u) || 0 == (uExtra & 2147483647) && 0 == (this.u0 & 1)) {
  343.                     if (0 != uExtra) {
  344.                         this.error++;
  345.                     }
  346.                     return;
  347.                 }
  348.                 this.error++;
  349.                
  350.                 // Round up.
  351.                 if (0 != AddU(ref this.u0, 1) && 0 != AddU(ref this.u1, 1) && 0 != AddU(ref this.u2, 1)) {
  352.                     Debug.Assert(this.IsZero);
  353.                     this.u2 = 2147483648u;
  354.                     this.exp++;
  355.                 }
  356.             }
  357.             #endif
  358.            
  359.             // Test to see if the num is zero. This works even if we're not normalized.
  360.             private bool IsZero {
  361.                 get { return (0 == this.u2) && (0 == this.u1) && (0 == this.u0); }
  362.             }
  363.            
  364. /*  ----------------------------------------------------------------------------
  365.                 Normalize()
  366.                 Normalize the big number - make sure the high bit is 1 or everything is zero
  367.                 (including the exponent).
  368.             */           
  369.             private void Normalize()
  370.             {
  371.                 int w1;
  372.                 int w2;
  373.                
  374.                 // Normalize mantissa
  375.                 if (0 == this.u2) {
  376.                     if (0 == this.u1) {
  377.                         if (0 == this.u0) {
  378.                             this.exp = 0;
  379.                             return;
  380.                         }
  381.                         this.u2 = this.u0;
  382.                         this.u0 = 0;
  383.                         this.exp -= 64;
  384.                     }
  385.                     else {
  386.                         this.u2 = this.u1;
  387.                         this.u1 = this.u0;
  388.                         this.u0 = 0;
  389.                         this.exp -= 32;
  390.                     }
  391.                 }
  392.                
  393.                 if (0 != (w1 = CbitZeroLeft(this.u2))) {
  394.                     w2 = 32 - w1;
  395.                     this.u2 = (this.u2 << w1) | (this.u1 >> w2);
  396.                     this.u1 = (this.u1 << w1) | (this.u0 >> w2);
  397.                     this.u0 = (this.u0 << w1);
  398.                     this.exp -= w1;
  399.                 }
  400.             }
  401.            
  402. /*  ----------------------------------------------------------------------------
  403.                 Mul()
  404.                 Multiply this big number by another big number.
  405.             */           
  406.             private void Mul(ref BigNumber numOp)
  407.             {
  408.                 Debug.Assert(0 != (this.u2 & 2147483648u));
  409.                 Debug.Assert(0 != (numOp.u2 & 2147483648u));
  410.                
  411.                 //uint *rgu = stackalloc uint[6];
  412.                 uint rgu0 = 0;
  413.                 uint rgu1 = 0;
  414.                 uint rgu2 = 0;
  415.                 uint rgu3 = 0;
  416.                 uint rgu4 = 0;
  417.                 uint rgu5 = 0;
  418.                 uint uLo;
  419.                 uint uHi;
  420.                 uint uT;
  421.                 uint wCarry;
  422.                
  423.                 if (0 != (uT = this.u0)) {
  424.                     uLo = MulU(uT, numOp.u0, out uHi);
  425.                     rgu0 = uLo;
  426.                     rgu1 = uHi;
  427.                    
  428.                     uLo = MulU(uT, numOp.u1, out uHi);
  429.                     Debug.Assert(uHi < 4294967295u);
  430.                     wCarry = AddU(ref rgu1, uLo);
  431.                     AddU(ref rgu2, uHi + wCarry);
  432.                    
  433.                     uLo = MulU(uT, numOp.u2, out uHi);
  434.                     Debug.Assert(uHi < 4294967295u);
  435.                     wCarry = AddU(ref rgu2, uLo);
  436.                     AddU(ref rgu3, uHi + wCarry);
  437.                 }
  438.                
  439.                 if (0 != (uT = this.u1)) {
  440.                     uLo = MulU(uT, numOp.u0, out uHi);
  441.                     Debug.Assert(uHi < 4294967295u);
  442.                     wCarry = AddU(ref rgu1, uLo);
  443.                     wCarry = AddU(ref rgu2, uHi + wCarry);
  444.                     if (0 != wCarry && 0 != AddU(ref rgu3, 1)) {
  445.                         AddU(ref rgu4, 1);
  446.                     }
  447.                     uLo = MulU(uT, numOp.u1, out uHi);
  448.                     Debug.Assert(uHi < 4294967295u);
  449.                     wCarry = AddU(ref rgu2, uLo);
  450.                     wCarry = AddU(ref rgu3, uHi + wCarry);
  451.                     if (0 != wCarry) {
  452.                         AddU(ref rgu4, 1);
  453.                     }
  454.                     uLo = MulU(uT, numOp.u2, out uHi);
  455.                     Debug.Assert(uHi < 4294967295u);
  456.                     wCarry = AddU(ref rgu3, uLo);
  457.                     AddU(ref rgu4, uHi + wCarry);
  458.                 }
  459.                
  460.                 uT = this.u2;
  461.                 Debug.Assert(0 != uT);
  462.                 uLo = MulU(uT, numOp.u0, out uHi);
  463.                 Debug.Assert(uHi < 4294967295u);
  464.                 wCarry = AddU(ref rgu2, uLo);
  465.                 wCarry = AddU(ref rgu3, uHi + wCarry);
  466.                 if (0 != wCarry && 0 != AddU(ref rgu4, 1)) {
  467.                     AddU(ref rgu5, 1);
  468.                 }
  469.                 uLo = MulU(uT, numOp.u1, out uHi);
  470.                 Debug.Assert(uHi < 4294967295u);
  471.                 wCarry = AddU(ref rgu3, uLo);
  472.                 wCarry = AddU(ref rgu4, uHi + wCarry);
  473.                 if (0 != wCarry) {
  474.                     AddU(ref rgu5, 1);
  475.                 }
  476.                 uLo = MulU(uT, numOp.u2, out uHi);
  477.                 Debug.Assert(uHi < 4294967295u);
  478.                 wCarry = AddU(ref rgu4, uLo);
  479.                 AddU(ref rgu5, uHi + wCarry);
  480.                
  481.                 // Compute the new exponent
  482.                 this.exp += numOp.exp;
  483.                
  484.                 // Accumulate the error. Adding doesn't necessarily give an accurate
  485.                 // bound if both of the errors are bigger than 2.
  486.                 Debug.Assert(this.error <= 2 || numOp.error <= 2);
  487.                 this.error += numOp.error;
  488.                
  489.                 // Handle rounding and normalize.
  490.                 if (0 == (rgu5 & 2147483648u)) {
  491.                     if (0 != (rgu2 & 1073741824) && (0 != (rgu2 & 3221225471u) || 0 != rgu1 || 0 != rgu0)) {
  492.                         // Round up by 1
  493.                         if (0 != AddU(ref rgu2, 1073741824) && 0 != AddU(ref rgu3, 1) && 0 != AddU(ref rgu4, 1)) {
  494.                             AddU(ref rgu5, 1);
  495.                             if (0 != (rgu5 & 2147483648u)) {
  496.                                 goto LNormalized;
  497.                             }
  498.                         }
  499.                     }
  500.                    
  501.                     // have to shift by one
  502.                     Debug.Assert(0 != (rgu5 & 1073741824));
  503.                     this.u2 = (rgu5 << 1) | (rgu4 >> 31);
  504.                     this.u1 = (rgu4 << 1) | (rgu3 >> 31);
  505.                     this.u0 = (rgu3 << 1) | (rgu2 >> 31);
  506.                     this.exp--;
  507.                     this.error <<= 1;
  508.                    
  509.                     // Add one for the error.
  510.                     if (0 != (rgu2 & 2147483647) || 0 != rgu1 || 0 != rgu0) {
  511.                         this.error++;
  512.                     }
  513.                     return;
  514.                 }
  515.                 else {
  516.                     if (0 != (rgu2 & 2147483648u) && (0 != (rgu3 & 1) || 0 != (rgu2 & 2147483647) || 0 != rgu1 || 0 != rgu0)) {
  517.                         // Round up by 1
  518.                         if (0 != AddU(ref rgu3, 1) && 0 != AddU(ref rgu4, 1) && 0 != AddU(ref rgu5, 1)) {
  519.                             Debug.Assert(0 == rgu3);
  520.                             Debug.Assert(0 == rgu4);
  521.                             Debug.Assert(0 == rgu5);
  522.                             rgu5 = 2147483648u;
  523.                             this.exp++;
  524.                         }
  525.                     }
  526.                 }
  527.                 LNormalized:
  528.                
  529.                 this.u2 = rgu5;
  530.                 this.u1 = rgu4;
  531.                 this.u0 = rgu3;
  532.                
  533.                 // Add one for the error.
  534.                 if (0 != rgu2 || 0 != rgu1 || 0 != rgu0) {
  535.                     this.error++;
  536.                 }
  537.             }
  538.            
  539.             // Get the double value.
  540.             public static explicit operator double(BigNumber bn)
  541.             {
  542.                 uint uEx;
  543.                 int exp;
  544.                 uint dblHi;
  545.                 uint dblLo;
  546.                
  547.                 Debug.Assert(0 != (bn.u2 & 2147483648u));
  548.                
  549.                 exp = bn.exp + 1022;
  550.                 if (exp >= 2047) {
  551.                     return Double.PositiveInfinity;
  552.                 }
  553.                
  554.                 // Round after filling in the bits. In the extra uint, we set the low bit
  555.                 // if there are any extra non-zero bits. This is for breaking the tie when
  556.                 // deciding whether to round up or down.
  557.                 if (exp > 0) {
  558.                     // Normalized.
  559.                     dblHi = ((uint)exp << 20) | ((bn.u2 & 2147483647) >> 11);
  560.                     dblLo = bn.u2 << 21 | bn.u1 >> 11;
  561.                     uEx = bn.u1 << 21 | NotZero(bn.u0);
  562.                 }
  563.                 else if (exp > -20) {
  564.                     // Denormal with some high bits.
  565.                     int wT = 12 - exp;
  566.                     Debug.Assert(wT >= 12 && wT < 32);
  567.                    
  568.                     dblHi = bn.u2 >> wT;
  569.                     dblLo = (bn.u2 << (32 - wT)) | (bn.u1 >> wT);
  570.                     uEx = (bn.u1 << (32 - wT)) | NotZero(bn.u0);
  571.                 }
  572.                 else if (exp == -20) {
  573.                     // Denormal with no high bits.
  574.                     dblHi = 0;
  575.                     dblLo = bn.u2;
  576.                     uEx = bn.u1 | (uint)(0 != bn.u0 ? 1 : 0);
  577.                 }
  578.                 else if (exp > -52) {
  579.                     // Denormal with no high bits.
  580.                     int wT = -exp - 20;
  581.                     Debug.Assert(wT > 0 && wT < 32);
  582.                    
  583.                     dblHi = 0;
  584.                     dblLo = bn.u2 >> wT;
  585.                     uEx = bn.u2 << (32 - wT) | NotZero(bn.u1) | NotZero(bn.u0);
  586.                 }
  587.                 else if (exp == -52) {
  588.                     // Zero unless we round up below.
  589.                     dblHi = 0;
  590.                     dblLo = 0;
  591.                     uEx = bn.u2 | NotZero(bn.u1) | NotZero(bn.u0);
  592.                 }
  593.                 else {
  594.                     return 0.0;
  595.                 }
  596.                
  597.                 // Handle rounding
  598.                 if (0 != (uEx & 2147483648u) && (0 != (uEx & 2147483647) || 0 != (dblLo & 1))) {
  599.                     // Round up. Note that this works even when we overflow into the
  600.                     // exponent.
  601.                     if (0 != AddU(ref dblLo, 1)) {
  602.                         AddU(ref dblHi, 1);
  603.                     }
  604.                 }
  605.                 return BitConverter.Int64BitsToDouble((long)dblHi << 32 | dblLo);
  606.             }
  607.            
  608.             // Lop off the integer part and return it.
  609.             private uint UMod1()
  610.             {
  611.                 if (this.exp <= 0) {
  612.                     return 0;
  613.                 }
  614.                 Debug.Assert(this.exp <= 32);
  615.                 uint uT = this.u2 >> (32 - this.exp);
  616.                 this.u2 &= (uint)2147483647 >> (this.exp - 1);
  617.                 Normalize();
  618.                 return uT;
  619.             }
  620.            
  621.             // If error is not zero, add it and set error to zero.
  622.             public void MakeUpperBound()
  623.             {
  624.                 Debug.Assert(this.error < 4294967295u);
  625.                 uint uT = (this.error + 1) >> 1;
  626.                
  627.                 if (0 != uT && 0 != AddU(ref this.u0, uT) && 0 != AddU(ref this.u1, 1) && 0 != AddU(ref this.u2, 1)) {
  628.                     Debug.Assert(0 == this.u2 && 0 == this.u1);
  629.                     this.u2 = 2147483648u;
  630.                     this.u0 = (this.u0 >> 1) + (this.u0 & 1);
  631.                     this.exp++;
  632.                 }
  633.                 this.error = 0;
  634.             }
  635.            
  636.             // If error is not zero, subtract it and set error to zero.
  637.             public void MakeLowerBound()
  638.             {
  639.                 Debug.Assert(this.error < 4294967295u);
  640.                 uint uT = (this.error + 1) >> 1;
  641.                
  642.                 if (0 != uT && 0 == AddU(ref this.u0, (uint)-(int)uT) && 0 == AddU(ref this.u1, 4294967295u)) {
  643.                     AddU(ref this.u2, 4294967295u);
  644.                     if (0 == (2147483648u & this.u2)) {
  645.                         Normalize();
  646.                     }
  647.                 }
  648.                 this.error = 0;
  649.             }
  650.            
  651. /*  ----------------------------------------------------------------------------
  652.                 DblToRgbFast()
  653.                 Get mantissa bytes (BCD).
  654.             */           
  655.             public static bool DblToRgbFast(double dbl, byte[] mantissa, out int exponent, out int mantissaSize)
  656.             {
  657.                 BigNumber numHH;
  658.                 BigNumber numHL;
  659.                 BigNumber numLH;
  660.                 BigNumber numLL;
  661.                 BigNumber numBase;
  662.                 BigNumber tenPower;
  663.                 int ib;
  664.                 int iT;
  665.                 uint uT;
  666.                 uint uScale;
  667.                 byte bHH;
  668.                 byte bHL;
  669.                 byte bLH;
  670.                 byte bLL;
  671.                 uint uHH;
  672.                 uint uHL;
  673.                 uint uLH;
  674.                 uint uLL;
  675.                 int wExp2;
  676.                 int wExp10 = 0;
  677.                 double dblInt;
  678.                 uint dblHi;
  679.                 uint dblLo;
  680.                
  681.                 dblHi = DblHi(dbl);
  682.                 dblLo = DblLo(dbl);
  683.                
  684.                 // Caller should take care of 0, negative and non-finite values.
  685.                 Debug.Assert(!IsSpecial(dbl));
  686.                 Debug.Assert(0 < dbl);
  687.                
  688.                 // Get numHH and numLL such that numLL < dbl < numHH and the
  689.                 // difference between adjacent values is half the distance to the next
  690.                 // representable value (in a double).
  691.                 wExp2 = (int)((dblHi >> 20) & 2047);
  692.                 if (wExp2 > 0) {
  693.                     // See if dbl is a small integer.
  694.                     if (wExp2 >= 1023 && wExp2 <= 1075 && dbl == Math.Floor(dbl)) {
  695.                         goto LSmallInt;
  696.                     }
  697.                    
  698.                     // Normalized
  699.                     numBase.u2 = 2147483648u | ((dblHi & 16777215) << 11) | (dblLo >> 21);
  700.                     numBase.u1 = dblLo << 11;
  701.                     numBase.u0 = 0;
  702.                     numBase.exp = wExp2 - 1022;
  703.                     numBase.error = 0;
  704.                    
  705.                     // Get the upper bound
  706.                     numHH = numBase;
  707.                     numHH.u1 |= (1 << 10);
  708.                    
  709.                     // Get the lower bound. A power of 2 must be special cased.
  710.                     numLL = numBase;
  711.                     if (2147483648u == numLL.u2 && 0 == numLL.u1) {
  712.                         // Subtract (0x00000000, 0x00000200, 0x00000000). Same as adding
  713.                         // (0xFFFFFFFF, 0xFFFFFE00, 0x00000000)
  714.                         uT = 4294966784u;
  715.                     }
  716.                     else {
  717.                         // Subtract (0x00000000, 0x00000400, 0x00000000). Same as adding
  718.                         // (0xFFFFFFFF, 0xFFFFFC00, 0x00000000)
  719.                         uT = 4294966272u;
  720.                     }
  721.                     if (0 == AddU(ref numLL.u1, uT)) {
  722.                         AddU(ref numLL.u2, 4294967295u);
  723.                         if (0 == (2147483648u & numLL.u2)) {
  724.                             numLL.Normalize();
  725.                         }
  726.                     }
  727.                 }
  728.                 else {
  729.                     // Denormal
  730.                     numBase.u2 = dblHi & 1048575;
  731.                     numBase.u1 = dblLo;
  732.                     numBase.u0 = 0;
  733.                     numBase.exp = -1010;
  734.                     numBase.error = 0;
  735.                    
  736.                     // Get the upper bound
  737.                     numHH = numBase;
  738.                     numHH.u0 = 2147483648u;
  739.                    
  740.                     // Get the lower bound
  741.                     numLL = numHH;
  742.                     if (0 == AddU(ref numLL.u1, 4294967295u)) {
  743.                         AddU(ref numLL.u2, 4294967295u);
  744.                     }
  745.                    
  746.                     numBase.Normalize();
  747.                     numHH.Normalize();
  748.                     numLL.Normalize();
  749.                 }
  750.                
  751.                 // Multiply by powers of ten until 0 < numHH.exp < 32.
  752.                 if (numHH.exp >= 32) {
  753.                     iT = (numHH.exp - 25) * 15 / -TenPowersNeg[45].exp;
  754.                     Debug.Assert(iT >= 0 && iT < 16);
  755.                     if (iT > 0) {
  756.                         tenPower = TenPowersNeg[30 + iT];
  757.                         Debug.Assert(numHH.exp + tenPower.exp > 1);
  758.                         numHH.Mul(ref tenPower);
  759.                         numLL.Mul(ref tenPower);
  760.                         wExp10 += iT * 32;
  761.                     }
  762.                    
  763.                     if (numHH.exp >= 32) {
  764.                         iT = (numHH.exp - 25) * 32 / -TenPowersNeg[31].exp;
  765.                         Debug.Assert(iT > 0 && iT <= 32);
  766.                         tenPower = TenPowersNeg[iT - 1];
  767.                         Debug.Assert(numHH.exp + tenPower.exp > 1);
  768.                         numHH.Mul(ref tenPower);
  769.                         numLL.Mul(ref tenPower);
  770.                         wExp10 += iT;
  771.                     }
  772.                 }
  773.                 else if (numHH.exp < 1) {
  774.                     iT = (25 - numHH.exp) * 15 / TenPowersPos[45].exp;
  775.                     Debug.Assert(iT >= 0 && iT < 16);
  776.                     if (iT > 0) {
  777.                         tenPower = TenPowersPos[30 + iT];
  778.                         Debug.Assert(numHH.exp + tenPower.exp <= 32);
  779.                         numHH.Mul(ref tenPower);
  780.                         numLL.Mul(ref tenPower);
  781.                         wExp10 -= iT * 32;
  782.                     }
  783.                    
  784.                     if (numHH.exp < 1) {
  785.                         iT = (25 - numHH.exp) * 32 / TenPowersPos[31].exp;
  786.                         Debug.Assert(iT > 0 && iT <= 32);
  787.                         tenPower = TenPowersPos[iT - 1];
  788.                         Debug.Assert(numHH.exp + tenPower.exp <= 32);
  789.                         numHH.Mul(ref tenPower);
  790.                         numLL.Mul(ref tenPower);
  791.                         wExp10 -= iT;
  792.                     }
  793.                 }
  794.                
  795.                 Debug.Assert(numHH.exp > 0 && numHH.exp < 32);
  796.                
  797.                 // Get the upper and lower bounds for these.
  798.                 numHL = numHH;
  799.                 numHH.MakeUpperBound();
  800.                 numHL.MakeLowerBound();
  801.                 uHH = numHH.UMod1();
  802.                 uHL = numHL.UMod1();
  803.                 numLH = numLL;
  804.                 numLH.MakeUpperBound();
  805.                 numLL.MakeLowerBound();
  806.                 uLH = numLH.UMod1();
  807.                 uLL = numLL.UMod1();
  808.                 Debug.Assert(uLL <= uLH && uLH <= uHL && uHL <= uHH);
  809.                
  810.                 // Find the starting scale
  811.                 uScale = 1;
  812.                 if (uHH >= 100000000) {
  813.                     uScale = 100000000;
  814.                     wExp10 += 8;
  815.                 }
  816.                 else {
  817.                     if (uHH >= 10000) {
  818.                         uScale = 10000;
  819.                         wExp10 += 4;
  820.                     }
  821.                     if (uHH >= 100 * uScale) {
  822.                         uScale *= 100;
  823.                         wExp10 += 2;
  824.                     }
  825.                 }
  826.                 if (uHH >= 10 * uScale) {
  827.                     uScale *= 10;
  828.                     wExp10++;
  829.                 }
  830.                 wExp10++;
  831.                 Debug.Assert(uHH >= uScale && uHH / uScale < 10);
  832.                
  833.                 for (ib = 0;;) {
  834.                     Debug.Assert(uLL <= uHH);
  835.                     bHH = (byte)(uHH / uScale);
  836.                     uHH %= uScale;
  837.                     bLL = (byte)(uLL / uScale);
  838.                     uLL %= uScale;
  839.                    
  840.                     if (bHH != bLL) {
  841.                         break;
  842.                     }
  843.                    
  844.                     Debug.Assert(0 != uHH || !numHH.IsZero);
  845.                     mantissa[ib++] = bHH;
  846.                    
  847.                     if (1 == uScale) {
  848.                         // Multiply by 10^8.
  849.                         uScale = 10000000;
  850.                        
  851.                         numHH.Mul(ref TenPowersPos[7]);
  852.                         numHH.MakeUpperBound();
  853.                         uHH = numHH.UMod1();
  854.                         if (uHH >= 100000000) {
  855.                             goto LFail;
  856.                         }
  857.                         numHL.Mul(ref TenPowersPos[7]);
  858.                         numHL.MakeLowerBound();
  859.                         uHL = numHL.UMod1();
  860.                        
  861.                         numLH.Mul(ref TenPowersPos[7]);
  862.                         numLH.MakeUpperBound();
  863.                         uLH = numLH.UMod1();
  864.                         numLL.Mul(ref TenPowersPos[7]);
  865.                         numLL.MakeLowerBound();
  866.                         uLL = numLL.UMod1();
  867.                     }
  868.                     else {
  869.                         uScale /= 10;
  870.                     }
  871.                 }
  872.                
  873.                 // LL and HH diverged. Get the digit values for LH and HL.
  874.                 Debug.Assert(0 <= bLL && bLL < bHH && bHH <= 9);
  875.                 bLH = (byte)((uLH / uScale) % 10);
  876.                 uLH %= uScale;
  877.                 bHL = (byte)((uHL / uScale) % 10);
  878.                 uHL %= uScale;
  879.                
  880.                 if (bLH >= bHL) {
  881.                     goto LFail;
  882.                 }
  883.                
  884.                 // LH and HL also diverged.
  885.                
  886.                 // We can get by with one fewer digit if: LL == LH and bLH is zero
  887.                 // and the current value of LH is zero and the least significant bit of
  888.                 // the double is zero. In this case, we have exactly the digit sequence
  889.                 // for the original numLL and IEEE and will rounds numLL up to the double.
  890.                 if (0 == bLH && 0 == uLH && numLH.IsZero && 0 == (dblLo & 1)) {
  891.                 }
  892.                 else if (bHL - bLH > 1) {
  893.                     // HL and LH differ by at least two in this digit, so split
  894.                     // the difference.
  895.                     mantissa[ib++] = (byte)((bHL + bLH + 1) / 2);
  896.                 }
  897.                 else if (0 != uHL || !numHL.IsZero || 0 == (dblLo & 1)) {
  898.                     // We can just use bHL because this guarantees that we're bigger than
  899.                     // LH and less than HL, so must convert to the double.
  900.                     mantissa[ib++] = bHL;
  901.                 }
  902.                 else {
  903.                     goto LFail;
  904.                 }
  905.                
  906.                 exponent = wExp10;
  907.                 mantissaSize = ib;
  908.                 goto LSucceed;
  909.                 LSmallInt:
  910.                
  911.                 // dbl should be an integer from 1 to (2^53 - 1).
  912.                 dblInt = dbl;
  913.                 Debug.Assert(dblInt == Math.Floor(dblInt) && 1 <= dblInt && dblInt <= 9.00719925474099E+15);
  914.                
  915.                 iT = 0;
  916.                 if (dblInt >= C10toN[iT + 8]) {
  917.                     iT += 8;
  918.                 }
  919.                 if (dblInt >= C10toN[iT + 4]) {
  920.                     iT += 4;
  921.                 }
  922.                 if (dblInt >= C10toN[iT + 2]) {
  923.                     iT += 2;
  924.                 }
  925.                 if (dblInt >= C10toN[iT + 1]) {
  926.                     iT += 1;
  927.                 }
  928.                 Debug.Assert(iT >= 0 && iT <= 15);
  929.                 Debug.Assert(dblInt >= C10toN[iT] && dblInt < C10toN[iT + 1]);
  930.                 exponent = iT + 1;
  931.                
  932.                 for (ib = 0; 0 != dblInt; iT--) {
  933.                     Debug.Assert(iT >= 0);
  934.                     bHH = (byte)(dblInt / C10toN[iT]);
  935.                     dblInt -= bHH * C10toN[iT];
  936.                     Debug.Assert(dblInt == Math.Floor(dblInt) && 0 <= dblInt && dblInt < C10toN[iT]);
  937.                     mantissa[ib++] = bHH;
  938.                 }
  939.                 mantissaSize = ib;
  940.                 LSucceed:
  941.                
  942.                
  943.                 #if DEBUG
  944.                 // Verify that precise is working and gives the same answer
  945.                 if (mantissaSize > 0) {
  946.                     byte[] mantissaPrec = new byte[20];
  947.                     int exponentPrec;
  948.                     int mantissaSizePrec;
  949.                     int idx;
  950.                    
  951.                     DblToRgbPrecise(dbl, mantissaPrec, out exponentPrec, out mantissaSizePrec);
  952.                     Debug.Assert(exponent == exponentPrec && mantissaSize == mantissaSizePrec);
  953.                     // Assert(!memcmp(mantissa, mantissaPrec, mantissaSizePrec - 1));
  954.                     bool equal = true;
  955.                     for (idx = 0; idx < mantissaSize; idx++) {
  956.                         equal &= ((mantissa[idx] == mantissaPrec[idx]) || (idx == mantissaSize - 1) && Math.Abs(mantissa[idx] - mantissaPrec[idx]) <= 1);
  957.                     }
  958.                     Debug.Assert(equal, "DblToRgbFast and DblToRgbPrecise should give the same result");
  959.                 }
  960.                 #endif
  961.                
  962.                 return true;
  963.                 LFail:
  964.                
  965.                 exponent = mantissaSize = 0;
  966.                 return false;
  967.             }
  968.            
  969. /*  ----------------------------------------------------------------------------
  970.                 DblToRgbPrecise()
  971.                 Uses big integer arithmetic to get the sequence of digits.
  972.             */           
  973.             public static void DblToRgbPrecise(double dbl, byte[] mantissa, out int exponent, out int mantissaSize)
  974.             {
  975.                 BigInteger biNum = new BigInteger();
  976.                 BigInteger biDen = new BigInteger();
  977.                 BigInteger biHi = new BigInteger();
  978.                 BigInteger biLo = new BigInteger();
  979.                 BigInteger biT = new BigInteger();
  980.                 BigInteger biHiLo;
  981.                 byte bT;
  982.                 bool fPow2;
  983.                 int ib;
  984.                 int cu;
  985.                 int wExp10;
  986.                 int wExp2;
  987.                 int w1;
  988.                 int w2;
  989.                 int c2Num;
  990.                 int c2Den;
  991.                 int c5Num;
  992.                 int c5Den;
  993.                 double dblT;
  994.                 //uint *rgu = stackalloc uint[2];
  995.                 uint rgu0;
  996.                 uint rgu1;
  997.                 uint dblHi;
  998.                 uint dblLo;
  999.                
  1000.                 dblHi = DblHi(dbl);
  1001.                 dblLo = DblLo(dbl);
  1002.                
  1003.                 // Caller should take care of 0, negative and non-finite values.
  1004.                 Debug.Assert(!IsSpecial(dbl));
  1005.                 Debug.Assert(0 < dbl);
  1006.                
  1007.                 // Init the Denominator, Hi error and Lo error bigints.
  1008.                 biDen.InitFromDigits(1, 0, 1);
  1009.                 biHi.InitFromDigits(1, 0, 1);
  1010.                
  1011.                 wExp2 = (int)(((dblHi & 2146435072) >> 20) - 1075);
  1012.                 rgu1 = dblHi & 1048575;
  1013.                 rgu0 = dblLo;
  1014.                 cu = 2;
  1015.                 fPow2 = false;
  1016.                 if (wExp2 == -1075) {
  1017.                     // dbl is denormalized.
  1018.                     Debug.Assert(0 == (dblHi & 2146435072));
  1019.                     if (0 == rgu1) {
  1020.                         cu = 1;
  1021.                     }
  1022.                    
  1023.                     // Get dblT such that dbl / dblT is a power of 2 and 1 <= dblT < 2.
  1024.                     // First multiply by a power of 2 to get a normalized value.
  1025.                     dblT = BitConverter.Int64BitsToDouble(1341128704l << 32);
  1026.                     dblT *= dbl;
  1027.                     Debug.Assert(0 != (DblHi(dblT) & 2146435072));
  1028.                    
  1029.                     // This is the power of 2.
  1030.                     w1 = (int)((DblHi(dblT) & 2146435072) >> 20) - (256 + 1023);
  1031.                    
  1032.                     dblHi = DblHi(dblT);
  1033.                     dblHi &= 1048575;
  1034.                     dblHi |= 1072693248;
  1035.                     dblT = BitConverter.Int64BitsToDouble((long)dblHi << 32 | DblLo(dblT));
  1036.                    
  1037.                     // Adjust wExp2 because we don't have the implicit bit.
  1038.                     wExp2++;
  1039.                 }
  1040.                 else {
  1041.                     // Get dblT such that dbl / dblT is a power of 2 and 1 <= dblT < 2.
  1042.                     // First multiply by a power of 2 to get a normalized value.
  1043.                     dblHi &= 1048575;
  1044.                     dblHi |= 1072693248;
  1045.                     dblT = BitConverter.Int64BitsToDouble((long)dblHi << 32 | dblLo);
  1046.                    
  1047.                     // This is the power of 2.
  1048.                     w1 = wExp2 + 52;
  1049.                    
  1050.                     if (0 == rgu0 && 0 == rgu1 && wExp2 > -1074) {
  1051.                         // Power of 2 bigger than smallest normal. The next smaller
  1052.                         // representable value is closer than the next larger value.
  1053.                         rgu1 = 2097152;
  1054.                         wExp2--;
  1055.                         fPow2 = true;
  1056.                     }
  1057.                     else {
  1058.                         // Normalized and not a power of 2 or the smallest normal. The
  1059.                         // representable values on either side are the same distance away.
  1060.                         rgu1 |= 1048576;
  1061.                     }
  1062.                 }
  1063.                
  1064.                 // Compute an approximation to the base 10 log. This is borrowed from
  1065.                 // David Gay's paper.
  1066.                 Debug.Assert(1 <= dblT && dblT < 2);
  1067.                 dblT = (dblT - 1.5) * 0.289529654602168 + 0.1760912590558 + w1 * 0.301029995663981;
  1068.                 wExp10 = (int)dblT;
  1069.                 if (dblT < 0 && dblT != wExp10) {
  1070.                     wExp10--;
  1071.                 }
  1072.                
  1073.                 if (wExp2 >= 0) {
  1074.                     c2Num = wExp2;
  1075.                     c2Den = 0;
  1076.                 }
  1077.                 else {
  1078.                     c2Num = 0;
  1079.                     c2Den = -wExp2;
  1080.                 }
  1081.                
  1082.                 if (wExp10 >= 0) {
  1083.                     c5Num = 0;
  1084.                     c5Den = wExp10;
  1085.                     c2Den += wExp10;
  1086.                 }
  1087.                 else {
  1088.                     c2Num -= wExp10;
  1089.                     c5Num = -wExp10;
  1090.                     c5Den = 0;
  1091.                 }
  1092.                
  1093.                 if (c2Num > 0 && c2Den > 0) {
  1094.                     w1 = c2Num < c2Den ? c2Num : c2Den;
  1095.                     c2Num -= w1;
  1096.                     c2Den -= w1;
  1097.                 }
  1098.                 // We need a bit for the Hi and Lo values.
  1099.                 c2Num++;
  1100.                 c2Den++;
  1101.                
  1102.                 // Initialize biNum and multiply by powers of 5.
  1103.                 if (c5Num > 0) {
  1104.                     Debug.Assert(0 == c5Den);
  1105.                     biHi.MulPow5(c5Num);
  1106.                     biNum.InitFromBigint(biHi);
  1107.                     if (1 == cu) {
  1108.                         biNum.MulAdd(rgu0, 0);
  1109.                     }
  1110.                     else {
  1111.                         biNum.MulAdd(rgu1, 0);
  1112.                         biNum.ShiftLeft(32);
  1113.                         if (0 != rgu0) {
  1114.                             biT.InitFromBigint(biHi);
  1115.                             biT.MulAdd(rgu0, 0);
  1116.                             biNum.Add(biT);
  1117.                         }
  1118.                     }
  1119.                 }
  1120.                 else {
  1121.                     Debug.Assert(cu <= 2);
  1122.                     biNum.InitFromDigits(rgu0, rgu1, cu);
  1123.                     if (c5Den > 0) {
  1124.                         biDen.MulPow5(c5Den);
  1125.                     }
  1126.                 }
  1127.                
  1128.                 // BigInteger.DivRem only works if the 4 high bits of the divisor are 0.
  1129.                 // It works most efficiently if there are exactly 4 zero high bits.
  1130.                 // Adjust c2Den and c2Num to guarantee this.
  1131.                 w1 = CbitZeroLeft(biDen[biDen.Length - 1]);
  1132.                 w1 = (w1 + 28 - c2Den) & 31;
  1133.                 c2Num += w1;
  1134.                 c2Den += w1;
  1135.                
  1136.                 // Multiply by powers of 2.
  1137.                 Debug.Assert(c2Num > 0 && c2Den > 0);
  1138.                 biNum.ShiftLeft(c2Num);
  1139.                 if (c2Num > 1) {
  1140.                     biHi.ShiftLeft(c2Num - 1);
  1141.                 }
  1142.                 biDen.ShiftLeft(c2Den);
  1143.                 Debug.Assert(0 == (biDen[biDen.Length - 1] & 4026531840u));
  1144.                 Debug.Assert(0 != (biDen[biDen.Length - 1] & 134217728));
  1145.                
  1146.                 // Get biHiLo and handle the power of 2 case where biHi needs to be doubled.
  1147.                 if (fPow2) {
  1148.                     biHiLo = biLo;
  1149.                     biHiLo.InitFromBigint(biHi);
  1150.                     biHi.ShiftLeft(1);
  1151.                 }
  1152.                 else {
  1153.                     biHiLo = biHi;
  1154.                 }
  1155.                
  1156.                 for (ib = 0;;) {
  1157.                     bT = (byte)biNum.DivRem(biDen);
  1158.                     if (0 == ib && 0 == bT) {
  1159.                         // Our estimate of wExp10 was too big. Oh well.
  1160.                         wExp10--;
  1161.                         goto LSkip;
  1162.                     }
  1163.                    
  1164.                     // w1 = sign(biNum - biHiLo).
  1165.                     w1 = biNum.CompareTo(biHiLo);
  1166.                    
  1167.                     // w2 = sign(biNum + biHi - biDen).
  1168.                     if (biDen.CompareTo(biHi) < 0) {
  1169.                         w2 = 1;
  1170.                     }
  1171.                     else {
  1172.                         biT.InitFromBigint(biDen);
  1173.                         biT.Subtract(biHi);
  1174.                         w2 = biNum.CompareTo(biT);
  1175.                     }
  1176.                    
  1177.                     // if (biNum + biHi == biDen && even)
  1178.                     if (0 == w2 && 0 == (dblLo & 1)) {
  1179.                         // Rounding up this digit produces exactly (biNum + biHi) which
  1180.                         // StrToDbl will round down to dbl.
  1181.                         if (bT == 9) {
  1182.                             goto LRoundUp9;
  1183.                         }
  1184.                         if (w1 > 0) {
  1185.                             bT++;
  1186.                         }
  1187.                         mantissa[ib++] = bT;
  1188.                         break;
  1189.                     }
  1190.                    
  1191.                     // if (biNum < biHiLo || biNum == biHiLo && even)
  1192.                     if (w1 < 0 || 0 == w1 && 0 == (dblLo & 1)) {
  1193.                         // if (biNum + biHi > biDen)
  1194.                         if (w2 > 0) {
  1195.                             // Decide whether to round up.
  1196.                             biNum.ShiftLeft(1);
  1197.                             w2 = biNum.CompareTo(biDen);
  1198.                             if ((w2 > 0 || 0 == w2 && (0 != (bT & 1))) && bT++ == 9) {
  1199.                                 goto LRoundUp9;
  1200.                             }
  1201.                         }
  1202.                         mantissa[ib++] = bT;
  1203.                         break;
  1204.                     }
  1205.                    
  1206.                     // if (biNum + biHi > biDen)
  1207.                     if (w2 > 0) {
  1208.                         // Round up and be done with it.
  1209.                         if (bT != 9) {
  1210.                             mantissa[ib++] = (byte)(bT + 1);
  1211.                             break;
  1212.                         }
  1213.                         goto LRoundUp9;
  1214.                     }
  1215.                    
  1216.                     // Save the digit.
  1217.                     mantissa[ib++] = bT;
  1218.                     LSkip:
  1219.                    
  1220.                     biNum.MulAdd(10, 0);
  1221.                     biHi.MulAdd(10, 0);
  1222.                     if ((object)biHiLo != (object)biHi) {
  1223.                         biHiLo.MulAdd(10, 0);
  1224.                     }
  1225.                     continue;
  1226.                     LRoundUp9:
  1227.                    
  1228.                     while (ib > 0) {
  1229.                         if (mantissa[--ib] != 9) {
  1230.                             mantissa[ib++]++;
  1231.                             goto LReturn;
  1232.                         }
  1233.                     }
  1234.                     wExp10++;
  1235.                     mantissa[ib++] = 1;
  1236.                     break;
  1237.                 }
  1238.                 LReturn:
  1239.                
  1240.                 exponent = wExp10 + 1;
  1241.                 mantissaSize = ib;
  1242.             }
  1243.            
  1244.             #region Powers of ten
  1245.             private static readonly BigNumber[] TenPowersPos = new BigNumber[46] {new BigNumber(0, 0, 2684354560u, 4, 0), new BigNumber(0, 0, 3355443200u, 7, 0), new BigNumber(0, 0, 4194304000u, 10, 0), new BigNumber(0, 0, 2621440000u, 14, 0), new BigNumber(0, 0, 3276800000u, 17, 0), new BigNumber(0, 0, 4096000000u, 20, 0), new BigNumber(0, 0, 2560000000u, 24, 0), new BigNumber(0, 0, 3200000000u, 27, 0), new BigNumber(0, 0, 4000000000u, 30, 0), new BigNumber(0, 0, 2500000000u, 34, 0),
  1246.             new BigNumber(0, 0, 3125000000u, 37, 0), new BigNumber(0, 0, 3906250000u, 40, 0), new BigNumber(0, 0, 2441406250u, 44, 0), new BigNumber(0, 2147483648u, 3051757812u, 47, 0), new BigNumber(0, 2684354560u, 3814697265u, 50, 0), new BigNumber(0, 67108864, 2384185791u, 54, 0), new BigNumber(0, 3305111552u, 2980232238u, 57, 0), new BigNumber(0, 1983905792, 3725290298u, 60, 0), new BigNumber(0, 2313682944u, 2328306436u, 64, 0), new BigNumber(0, 2892103680u, 2910383045u, 67, 0),
  1247.             new BigNumber(0, 393904128, 3637978807u, 70, 0), new BigNumber(0, 1856802816, 2273736754u, 74, 0), new BigNumber(0, 173519872, 2842170943u, 77, 0), new BigNumber(0, 3438125312u, 3552713678u, 80, 0), new BigNumber(0, 1075086496, 2220446049u, 84, 0), new BigNumber(0, 2417599944u, 2775557561u, 87, 0), new BigNumber(0, 4095741754u, 3469446951u, 90, 0), new BigNumber(1073741824, 4170451332u, 2168404344u, 94, 0), new BigNumber(1342177280, 918096869, 2710505431u, 97, 0), new BigNumber(2751463424u, 73879262, 3388131789u, 100, 0),
  1248.             new BigNumber(1291845632, 1166090902, 4235164736u, 103, 0), new BigNumber(4028628992u, 728806813, 2646977960u, 107, 0), new BigNumber(1019177842, 4291798741u, 3262652233u, 213, 1), new BigNumber(3318737231u, 3315274914u, 4021529366u, 319, 1), new BigNumber(3329176428u, 2162789599u, 2478458825u, 426, 1), new BigNumber(1467717739, 2145785770, 3054936363u, 532, 1), new BigNumber(2243682900u, 958879082, 3765499789u, 638, 1), new BigNumber(2193451889u, 3812411695u, 2320668415u, 745, 1), new BigNumber(3720056860u, 2650398349u, 2860444667u, 851, 1), new BigNumber(1937977068, 1550462860, 3525770265u, 957, 1),
  1249.                 // Positive powers of 10 to 96 bits precision.
  1250.                 // 10^1
  1251.                 // 10^2
  1252.                 // 10^3
  1253.                 // 10^4
  1254.                 // 10^5
  1255.                 // 10^6
  1256.                 // 10^7
  1257.                 // 10^8
  1258.                 // 10^9
  1259.                 // 10^10
  1260.                 // 10^11
  1261.                 // 10^12
  1262.                 // 10^13
  1263.                 // 10^14
  1264.                 // 10^15
  1265.                 // 10^16
  1266.                 // 10^17
  1267.                 // 10^18
  1268.                 // 10^19
  1269.                 // 10^20
  1270.                 // 10^21
  1271.                 // 10^22
  1272.                 // 10^23
  1273.                 // 10^24
  1274.                 // 10^25
  1275.                 // 10^26
  1276.                 // 10^27
  1277.                 // 10^28
  1278.                 // 10^29
  1279.                 // 10^30
  1280.                 // 10^31
  1281.                 // 10^32
  1282.                 // 10^64 (rounded up)
  1283.                 // 10^96 (rounded up)
  1284.                 // 10^128
  1285.                 // 10^160
  1286.                 // 10^192 (rounded up)
  1287.                 // 10^224 (rounded up)
  1288.                 // 10^256 (rounded up)
  1289.                 // 10^288
  1290.                 // 10^320
  1291.                 // 10^352 (rounded up)
  1292.                 // 10^384
  1293.                 // 10^416
  1294.                 // 10^448
  1295.                 // 10^480
  1296.             new BigNumber(3869316483u, 4073513845u, 2172923689u, 1064, 1), new BigNumber(1589582007, 3683650258u, 2678335232u, 1170, 1), new BigNumber(271056885, 2935532055u, 3301303056u, 1276, 1), new BigNumber(3051704177u, 3920665688u, 4069170183u, 1382, 1), new BigNumber(2817170568u, 3958895571u, 2507819745u, 1489, 1), new BigNumber(2113145460, 127246946, 3091126492u, 1595, 1)};
  1297.            
  1298.             private static readonly BigNumber[] TenPowersNeg = new BigNumber[46] {new BigNumber(3435973837u, 3435973836u, 3435973836u, -3, 1), new BigNumber(1030792151, 1889785610, 2748779069u, -6, 1), new BigNumber(1683627180, 2370821947u, 2199023255u, -9, 1), new BigNumber(3552796947u, 3793315115u, 3518437208u, -13, 1), new BigNumber(265257180, 457671715, 2814749767u, -16, 1), new BigNumber(2789186122u, 2943117749u, 2251799813u, -19, 1), new BigNumber(1026723958, 3849994940u, 3602879701u, -23, 1), new BigNumber(4257353003u, 2221002492u, 2882303761u, -26, 1), new BigNumber(828902025, 917808535, 2305843009u, -29, 1), new BigNumber(3044230158u, 3186480574u, 3689348814u, -33, 1),
  1299.             new BigNumber(4153371045u, 3408177918u, 2951479051u, -36, 1), new BigNumber(4181690295u, 1867548875, 2361183241u, -39, 1), new BigNumber(677750258, 1270091283, 3777893186u, -43, 1), new BigNumber(1401193666, 157079567, 3022314549u, -46, 1), new BigNumber(261961473, 984657113, 2417851639u, -49, 1), new BigNumber(1278131816, 3293438299u, 3868562622u, -53, 1), new BigNumber(163511994, 916763721, 3094850098u, -56, 1), new BigNumber(989803054, 2451397895u, 2475880078u, -59, 1), new BigNumber(724691428, 3063243173u, 3961408125u, -63, 1), new BigNumber(2297740061u, 2450594538u, 3169126500u, -66, 1),
  1300.             new BigNumber(3556178967u, 1960475630, 2535301200u, -69, 1), new BigNumber(1394919051, 3136761009u, 4056481920u, -73, 1), new BigNumber(1974928700, 2509408807u, 3245185536u, -76, 1), new BigNumber(3297929878u, 1148533586, 2596148429u, -79, 1), new BigNumber(981720510, 3555640657u, 4153837486u, -83, 1), new BigNumber(2503363326u, 1985519066, 3323069989u, -86, 1), new BigNumber(2002690661, 2447408712u, 2658455991u, -89, 1), new BigNumber(2345311598u, 2197867021u, 4253529586u, -93, 1), new BigNumber(158262360, 899300158, 3402823669u, -96, 1), new BigNumber(2703590266u, 1578433585, 2722258935u, -99, 1),
  1301.             new BigNumber(2162872213u, 1262746868, 2177807148u, -102, 1), new BigNumber(1742608622, 1161401530, 3484491437u, -106, 1), new BigNumber(1059297495, 2772036005u, 2826955303u, -212, 1), new BigNumber(299617026, 4252324763u, 2293498615u, -318, 1), new BigNumber(2893853687u, 1690100896, 3721414268u, -425, 1), new BigNumber(1508712807, 3681788051u, 3019169939u, -531, 1), new BigNumber(2070087331, 1411632134, 2449441655u, -637, 1), new BigNumber(2767765334u, 1244745405, 3974446316u, -744, 1), new BigNumber(4203811158u, 1668946233, 3224453925u, -850, 1), new BigNumber(1323526137, 2204812663u, 2615987810u, -956, 1),
  1302.                 // Negative powers of 10 to 96 bits precision.
  1303.                 // 10^-1 (rounded up)
  1304.                 // 10^-2
  1305.                 // 10^-3
  1306.                 // 10^-4
  1307.                 // 10^-5
  1308.                 // 10^-6 (rounded up)
  1309.                 // 10^-7
  1310.                 // 10^-8
  1311.                 // 10^-9 (rounded up)
  1312.                 // 10^-10
  1313.                 // 10^-11
  1314.                 // 10^-12
  1315.                 // 10^-13
  1316.                 // 10^-14 (rounded up)
  1317.                 // 10^-15
  1318.                 // 10^-16
  1319.                 // 10^-17 (rounded up)
  1320.                 // 10^-18
  1321.                 // 10^-19 (rounded up)
  1322.                 // 10^-20 (rounded up)
  1323.                 // 10^-21 (rounded up)
  1324.                 // 10^-22
  1325.                 // 10^-23
  1326.                 // 10^-24
  1327.                 // 10^-25 (rounded up)
  1328.                 // 10^-26
  1329.                 // 10^-27 (rounded up)
  1330.                 // 10^-28
  1331.                 // 10^-29
  1332.                 // 10^-30 (rounded up)
  1333.                 // 10^-31 (rounded up)
  1334.                 // 10^-32 (rounded up)
  1335.                 // 10^-64
  1336.                 // 10^-96
  1337.                 // 10^-128 (rounded up)
  1338.                 // 10^-160
  1339.                 // 10^-192
  1340.                 // 10^-224
  1341.                 // 10^-256 (rounded up)
  1342.                 // 10^-288
  1343.                 // 10^-320 (rounded up)
  1344.                 // 10^-352 (rounded up)
  1345.                 // 10^-384 (rounded up)
  1346.                 // 10^-416 (rounded up)
  1347.                 // 10^-448 (rounded up)
  1348.                 // 10^-480
  1349.             new BigNumber(2300620953u, 1199716560, 4244682903u, -1063, 1), new BigNumber(9598332, 1190350717, 3443695891u, -1169, 1), new BigNumber(2296094720u, 2971338839u, 2793858024u, -1275, 1), new BigNumber(441364487, 1073506470, 2266646913u, -1381, 1), new BigNumber(2227594191u, 3053929028u, 3677844889u, -1488, 1), new BigNumber(1642812130, 2030073654, 2983822260u, -1594, 1)};
  1350.             #endregion
  1351.            
  1352.         }
  1353.        
  1354. /**
  1355.         * Implementation of very large variable-precision non-negative integers.
  1356.         *
  1357.         * Hungarian: bi
  1358.         *
  1359.         */       
  1360.         private class BigInteger : IComparable
  1361.         {
  1362.             // Make this big enough that we rarely have to reallocate.
  1363.             private const int InitCapacity = 30;
  1364.            
  1365.             private int capacity;
  1366.             private int length;
  1367.             private uint[] digits;
  1368.            
  1369.             public BigInteger()
  1370.             {
  1371.                 capacity = InitCapacity;
  1372.                 length = 0;
  1373.                 digits = new uint[InitCapacity];
  1374.                 AssertValid();
  1375.             }
  1376.            
  1377.             public int Length {
  1378.                 get { return length; }
  1379.             }
  1380.            
  1381.             public uint this[int idx]
  1382.             {
  1383.                 get {
  1384.                     AssertValid();
  1385.                     Debug.Assert(0 <= idx && idx < length);
  1386.                     return digits[idx];
  1387.                 }
  1388.             }
  1389.            
  1390.             [Conditional("DEBUG")]
  1391.             private void AssertValidNoVal()
  1392.             {
  1393.                 Debug.Assert(capacity >= InitCapacity);
  1394.                 Debug.Assert(length >= 0 && length <= capacity);
  1395.             }
  1396.            
  1397.             [Conditional("DEBUG")]
  1398.             private void AssertValid()
  1399.             {
  1400.                 AssertValidNoVal();
  1401.                 Debug.Assert(0 == length || 0 != digits[length - 1]);
  1402.             }
  1403.            
  1404.             private void Ensure(int cu)
  1405.             {
  1406.                 AssertValidNoVal();
  1407.                
  1408.                 if (cu <= capacity) {
  1409.                     return;
  1410.                 }
  1411.                
  1412.                 cu += cu;
  1413.                 uint[] newDigits = new uint[cu];
  1414.                 digits.CopyTo(newDigits, 0);
  1415.                 digits = newDigits;
  1416.                 capacity = cu;
  1417.                
  1418.                 AssertValidNoVal();
  1419.             }
  1420.            
  1421. /*  ----------------------------------------------------------------------------
  1422.                 InitFromRgu()
  1423.                 Initialize this big integer from an array of uint values.
  1424.             */           
  1425.             public void InitFromRgu(uint[] rgu, int cu)
  1426.             {
  1427.                 AssertValid();
  1428.                 Debug.Assert(cu >= 0);
  1429.                
  1430.                 Ensure(cu);
  1431.                 length = cu;
  1432.                 for (int i = 0; i < cu; i++) {
  1433.                     digits[i] = rgu[i];
  1434.                 }
  1435.                 AssertValid();
  1436.             }
  1437.            
  1438. /*  ----------------------------------------------------------------------------
  1439.                 InitFromRgu()
  1440.                 Initialize this big integer from 0, 1, or 2 uint values.
  1441.             */           
  1442.             public void InitFromDigits(uint u0, uint u1, int cu)
  1443.             {
  1444.                 AssertValid();
  1445.                 Debug.Assert(2 <= capacity);
  1446.                
  1447.                 length = cu;
  1448.                 digits[0] = u0;
  1449.                 digits[1] = u1;
  1450.                 AssertValid();
  1451.             }
  1452.            
  1453. /*  ----------------------------------------------------------------------------
  1454.                 InitFromBigint()
  1455.                 Initialize this big integer from another BigInteger object.
  1456.             */           
  1457.             public void InitFromBigint(BigInteger biSrc)
  1458.             {
  1459.                 AssertValid();
  1460.                 biSrc.AssertValid();
  1461.                 Debug.Assert((object)this != (object)biSrc);
  1462.                
  1463.                 InitFromRgu(biSrc.digits, biSrc.length);
  1464.             }
  1465.            
  1466.             #if !NOPARSE || DEBUG
  1467. /*  ----------------------------------------------------------------------------
  1468.                 InitFromFloatingDecimal()
  1469.                 Initialize this big integer from a FloatingDecimal object.
  1470.             */           
  1471.             public void InitFromFloatingDecimal(FloatingDecimal dec)
  1472.             {
  1473.                 AssertValid();
  1474.                 Debug.Assert(dec.MantissaSize >= 0);
  1475.                
  1476.                 uint uAdd;
  1477.                 uint uMul;
  1478.                 int cu = (dec.MantissaSize + 8) / 9;
  1479.                 int mantissaSize = dec.MantissaSize;
  1480.                
  1481.                 Ensure(cu);
  1482.                 length = 0;
  1483.                
  1484.                 uAdd = 0;
  1485.                 uMul = 1;
  1486.                 for (int ib = 0; ib < mantissaSize; ib++) {
  1487.                     Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
  1488.                     if (1000000000 == uMul) {
  1489.                         MulAdd(uMul, uAdd);
  1490.                         uMul = 1;
  1491.                         uAdd = 0;
  1492.                     }
  1493.                     uMul *= 10;
  1494.                     uAdd = uAdd * 10 + dec[ib];
  1495.                 }
  1496.                 Debug.Assert(1 < uMul);
  1497.                 MulAdd(uMul, uAdd);
  1498.                
  1499.                 AssertValid();
  1500.             }
  1501.             #endif
  1502.            
  1503.             public void MulAdd(uint uMul, uint uAdd)
  1504.             {
  1505.                 AssertValid();
  1506.                 Debug.Assert(0 != uMul);
  1507.                
  1508.                 for (int i = 0; i < length; i++) {
  1509.                     uint d;
  1510.                     uint uT;
  1511.                     d = MulU(digits[i], uMul, out uT);
  1512.                     if (0 != uAdd) {
  1513.                         uT += AddU(ref d, uAdd);
  1514.                     }
  1515.                     digits[i] = d;
  1516.                     uAdd = uT;
  1517.                 }
  1518.                 if (0 != uAdd) {
  1519.                     Ensure(length + 1);
  1520.                     digits[length++] = uAdd;
  1521.                 }
  1522.                 AssertValid();
  1523.             }
  1524.            
  1525.             public void MulPow5(int c5)
  1526.             {
  1527.                 AssertValid();
  1528.                 Debug.Assert(c5 >= 0);
  1529.                
  1530.                 const uint C5to13 = 1220703125;
  1531.                 int cu = (c5 + 12) / 13;
  1532.                
  1533.                 if (0 == length || 0 == c5) {
  1534.                     return;
  1535.                 }
  1536.                
  1537.                 Ensure(length + cu);
  1538.                
  1539.                 for (; c5 >= 13; c5 -= 13) {
  1540.                     MulAdd(C5to13, 0);
  1541.                 }
  1542.                
  1543.                 if (c5 > 0) {
  1544.                     uint uT;
  1545.                     for (uT = 5; --c5 > 0;) {
  1546.                         uT *= 5;
  1547.                     }
  1548.                     MulAdd(uT, 0);
  1549.                 }
  1550.                 AssertValid();
  1551.             }
  1552.            
  1553.             public void ShiftLeft(int cbit)
  1554.             {
  1555.                 AssertValid();
  1556.                 Debug.Assert(cbit >= 0);
  1557.                
  1558.                 int idx;
  1559.                 int cu;
  1560.                 uint uExtra;
  1561.                
  1562.                 if (0 == cbit || 0 == length) {
  1563.                     return;
  1564.                 }
  1565.                
  1566.                 cu = cbit >> 5;
  1567.                 cbit &= 31;
  1568.                
  1569.                 if (cbit > 0) {
  1570.                     idx = length - 1;
  1571.                     uExtra = digits[idx] >> (32 - cbit);
  1572.                    
  1573.                     for (;; idx--) {
  1574.                         digits[idx] <<= cbit;
  1575.                         if (0 == idx) {
  1576.                             break;
  1577.                         }
  1578.                         digits[idx] |= digits[idx - 1] >> (32 - cbit);
  1579.                     }
  1580.                 }
  1581.                 else {
  1582.                     uExtra = 0;
  1583.                 }
  1584.                
  1585.                 if (cu > 0 || 0 != uExtra) {
  1586.                     // Make sure there's enough room.
  1587.                     idx = length + (0 != uExtra ? 1 : 0) + cu;
  1588.                     Ensure(idx);
  1589.                    
  1590.                     if (cu > 0) {
  1591.                         for (int i = length; 0 != i--;) {
  1592.                             digits[cu + i] = digits[i];
  1593.                         }
  1594.                         for (int i = 0; i < cu; i++) {
  1595.                             digits[i] = 0;
  1596.                         }
  1597.                         length += cu;
  1598.                     }
  1599.                    
  1600.                     // Throw on the extra one.
  1601.                     if (0 != uExtra) {
  1602.                         digits[length++] = uExtra;
  1603.                     }
  1604.                 }
  1605.                 AssertValid();
  1606.             }
  1607.            
  1608.             public void ShiftUsRight(int cu)
  1609.             {
  1610.                 AssertValid();
  1611.                 Debug.Assert(cu >= 0);
  1612.                
  1613.                 if (cu >= length) {
  1614.                     length = 0;
  1615.                 }
  1616.                 else if (cu > 0) {
  1617.                     for (int i = 0; i < length - cu; i++) {
  1618.                         digits[i] = digits[cu + i];
  1619.                     }
  1620.                     length -= cu;
  1621.                 }
  1622.                 AssertValid();
  1623.             }
  1624.            
  1625.             public void ShiftRight(int cbit)
  1626.             {
  1627.                 AssertValid();
  1628.                 Debug.Assert(cbit >= 0);
  1629.                
  1630.                 int idx;
  1631.                 int cu = cbit >> 5;
  1632.                 cbit &= 31;
  1633.                
  1634.                 if (cu > 0) {
  1635.                     ShiftUsRight(cu);
  1636.                 }
  1637.                
  1638.                 if (0 == cbit || 0 == length) {
  1639.                     AssertValid();
  1640.                     return;
  1641.                 }
  1642.                
  1643.                 for (idx = 0;;) {
  1644.                     digits[idx] >>= cbit;
  1645.                     if (++idx >= length) {
  1646.                         // Last one.
  1647.                         if (0 == digits[idx - 1]) {
  1648.                             length--;
  1649.                         }
  1650.                         break;
  1651.                     }
  1652.                     digits[idx - 1] |= digits[idx] << (32 - cbit);
  1653.                 }
  1654.                 AssertValid();
  1655.             }
  1656.            
  1657.             public int CompareTo(object obj)
  1658.             {
  1659.                 BigInteger bi = (BigInteger)obj;
  1660.                 AssertValid();
  1661.                 bi.AssertValid();
  1662.                
  1663.                 if (length > bi.length) {
  1664.                     return 1;
  1665.                 }
  1666.                 else if (length < bi.length) {
  1667.                     return -1;
  1668.                 }
  1669.                 else if (0 == length) {
  1670.                     return 0;
  1671.                 }
  1672.                
  1673.                 int idx;
  1674.                
  1675.                 for (idx = length - 1; digits[idx] == bi.digits[idx]; idx--) {
  1676.                     if (0 == idx) {
  1677.                         return 0;
  1678.                     }
  1679.                 }
  1680.                 Debug.Assert(idx >= 0 && idx < length);
  1681.                 Debug.Assert(digits[idx] != bi.digits[idx]);
  1682.                
  1683.                 return (digits[idx] > bi.digits[idx]) ? 1 : -1;
  1684.             }
  1685.            
  1686.             public void Add(BigInteger bi)
  1687.             {
  1688.                 AssertValid();
  1689.                 bi.AssertValid();
  1690.                 Debug.Assert((object)this != (object)bi);
  1691.                
  1692.                 int idx;
  1693.                 int cuMax;
  1694.                 int cuMin;
  1695.                 uint wCarry;
  1696.                
  1697.                 if ((cuMax = length) < (cuMin = bi.length)) {
  1698.                     cuMax = bi.length;
  1699.                     cuMin = length;
  1700.                     Ensure(cuMax + 1);
  1701.                 }
  1702.                
  1703.                 wCarry = 0;
  1704.                 for (idx = 0; idx < cuMin; idx++) {
  1705.                     if (0 != wCarry) {
  1706.                         wCarry = AddU(ref digits[idx], wCarry);
  1707.                     }
  1708.                     wCarry += AddU(ref digits[idx], bi.digits[idx]);
  1709.                 }
  1710.                
  1711.                 if (length < bi.length) {
  1712.                     for (; idx < cuMax; idx++) {
  1713.                         digits[idx] = bi.digits[idx];
  1714.                         if (0 != wCarry) {
  1715.                             wCarry = AddU(ref digits[idx], wCarry);
  1716.                         }
  1717.                     }
  1718.                     length = cuMax;
  1719.                 }
  1720.                 else {
  1721.                     for (; 0 != wCarry && idx < cuMax; idx++) {
  1722.                         wCarry = AddU(ref digits[idx], wCarry);
  1723.                     }
  1724.                 }
  1725.                
  1726.                 if (0 != wCarry) {
  1727.                     Ensure(length + 1);
  1728.                     digits[length++] = wCarry;
  1729.                 }
  1730.                 AssertValid();
  1731.             }
  1732.            
  1733.             public void Subtract(BigInteger bi)
  1734.             {
  1735.                 AssertValid();
  1736.                 bi.AssertValid();
  1737.                 Debug.Assert((object)this != (object)bi);
  1738.                
  1739.                 int idx;
  1740.                 uint wCarry;
  1741.                 uint uT;
  1742.                
  1743.                 if (length < bi.length) {
  1744.                     goto LNegative;
  1745.                 }
  1746.                
  1747.                 wCarry = 1;
  1748.                 for (idx = 0; idx < bi.length; idx++) {
  1749.                     Debug.Assert(0 == wCarry || 1 == wCarry);
  1750.                     uT = bi.digits[idx];
  1751.                    
  1752.                     // NOTE: We should really do:
  1753.                     // wCarry = AddU(ref digits[idx], wCarry);
  1754.                     // wCarry += AddU(ref digits[idx], ~uT);
  1755.                     // The only case where this is different than
  1756.                     // wCarry = AddU(ref digits[idx], ~uT + wCarry);
  1757.                     // is when 0 == uT and 1 == wCarry, in which case we don't
  1758.                     // need to add anything and wCarry should still be 1, so we can
  1759.                     // just skip the operations.
  1760.                    
  1761.                     if (0 != uT || 0 == wCarry) {
  1762.                         wCarry = AddU(ref digits[idx], ~uT + wCarry);
  1763.                     }
  1764.                 }
  1765.                 while (0 == wCarry && idx < length) {
  1766.                     wCarry = AddU(ref digits[idx], 4294967295u);
  1767.                 }
  1768.                
  1769.                 if (0 != wCarry) {
  1770.                     if (idx == length) {
  1771.                         // Trim off zeros.
  1772.                         while (--idx >= 0 && 0 == digits[idx]) {
  1773.                         }
  1774.                         length = idx + 1;
  1775.                     }
  1776.                     AssertValid();
  1777.                     return;
  1778.                 }
  1779.                 LNegative:
  1780.                
  1781.                 // bi was bigger than this.
  1782.                 Debug.Assert(false, "Who's subtracting to negative?");
  1783.                 length = 0;
  1784.                 AssertValid();
  1785.             }
  1786.            
  1787.             public uint DivRem(BigInteger bi)
  1788.             {
  1789.                 AssertValid();
  1790.                 bi.AssertValid();
  1791.                 Debug.Assert((object)this != (object)bi);
  1792.                
  1793.                 int idx;
  1794.                 int cu;
  1795.                 uint uQuo;
  1796.                 uint wCarry;
  1797.                 int wT;
  1798.                 uint uT;
  1799.                 uint uHi;
  1800.                 uint uLo;
  1801.                
  1802.                 cu = bi.length;
  1803.                 Debug.Assert(length <= cu);
  1804.                 if (length < cu) {
  1805.                     return 0;
  1806.                 }
  1807.                
  1808.                 // Get a lower bound on the quotient.
  1809.                 uQuo = (uint)(digits[cu - 1] / (bi.digits[cu - 1] + 1));
  1810.                 Debug.Assert(uQuo >= 0 && uQuo <= 9);
  1811.                
  1812.                 // Handle 0 and 1 as special cases.
  1813.                 switch (uQuo) {
  1814.                     case 0:
  1815.                         break;
  1816.                     case 1:
  1817.                         Subtract(bi);
  1818.                         break;
  1819.                     default:
  1820.                         uHi = 0;
  1821.                         wCarry = 1;
  1822.                         for (idx = 0; idx < cu; idx++) {
  1823.                             Debug.Assert(0 == wCarry || 1 == wCarry);
  1824.                            
  1825.                             // Compute the product.
  1826.                             uLo = MulU(uQuo, bi.digits[idx], out uT);
  1827.                             uHi = uT + AddU(ref uLo, uHi);
  1828.                            
  1829.                             // Subtract the product. See note in BigInteger.Subtract.
  1830.                             if (0 != uLo || 0 == wCarry) {
  1831.                                 wCarry = AddU(ref digits[idx], ~uLo + wCarry);
  1832.                             }
  1833.                         }
  1834.                         Debug.Assert(1 == wCarry);
  1835.                         Debug.Assert(idx == cu);
  1836.                        
  1837.                         // Trim off zeros.
  1838.                         while (--idx >= 0 && 0 == digits[idx]) {
  1839.                         }
  1840.                         length = idx + 1;
  1841.                         break;
  1842.                 }
  1843.                
  1844.                 if (uQuo < 9 && (wT = CompareTo(bi)) >= 0) {
  1845.                     // Quotient was off too small (by one).
  1846.                     uQuo++;
  1847.                     if (0 == wT) {
  1848.                         length = 0;
  1849.                     }
  1850.                     else {
  1851.                         Subtract(bi);
  1852.                     }
  1853.                 }
  1854.                 Debug.Assert(CompareTo(bi) < 0);
  1855.                
  1856.                 return uQuo;
  1857.             }
  1858.            
  1859.             #if NEVER
  1860.             public static explicit operator double(BigInteger bi)
  1861.             {
  1862.                 uint uHi;
  1863.                 uint uLo;
  1864.                 uint u1;
  1865.                 uint u2;
  1866.                 uint u3;
  1867.                 int idx;
  1868.                 int cbit;
  1869.                 uint dblHi;
  1870.                 uint dblLo;
  1871.                
  1872.                 switch (bi.length) {
  1873.                     case 0:
  1874.                         return 0;
  1875.                     case 1:
  1876.                         return bi.digits[0];
  1877.                     case 2:
  1878.                         return (ulong)bi.digits[1] << 32 | bi.digits[0];
  1879.                 }
  1880.                
  1881.                 Debug.Assert(3 <= bi.length);
  1882.                 if (bi.length > 32) {
  1883.                     // Result is infinite.
  1884.                     return BitConverter.Int64BitsToDouble(2146435072l << 32);
  1885.                 }
  1886.                
  1887.                 u1 = bi.digits[bi.length - 1];
  1888.                 u2 = bi.digits[bi.length - 2];
  1889.                 u3 = bi.digits[bi.length - 3];
  1890.                 Debug.Assert(0 != u1);
  1891.                 cbit = 31 - CbitZeroLeft(u1);
  1892.                
  1893.                 if (0 == cbit) {
  1894.                     uHi = u2;
  1895.                     uLo = u3;
  1896.                 }
  1897.                 else {
  1898.                     uHi = (u1 << (32 - cbit)) | (u2 >> cbit);
  1899.                     // Or 1 if there are any remaining nonzero bits in u3, so we take
  1900.                     // them into account when rounding.
  1901.                     uLo = (u2 << (32 - cbit)) | (u3 >> cbit) | NotZero(u3 << (32 - cbit));
  1902.                 }
  1903.                
  1904.                 // Set the mantissa bits.
  1905.                 dblHi = uHi >> 12;
  1906.                 dblLo = (uHi << 20) | (uLo >> 12);
  1907.                
  1908.                 // Set the exponent field.
  1909.                 dblHi |= (uint)(1023 + cbit + (bi.length - 1) * 32) << 20;
  1910.                
  1911.                 // Do IEEE rounding.
  1912.                 if (0 != (uLo & 2048)) {
  1913.                     if (0 != (uLo & 2047) || 0 != (dblLo & 1)) {
  1914.                         if (0 == ++dblLo) {
  1915.                             ++dblHi;
  1916.                         }
  1917.                     }
  1918.                     else {
  1919.                         // If there are any non-zero bits in digits from 0 to length - 4,
  1920.                         // round up.
  1921.                         for (idx = bi.length - 4; idx >= 0; idx--) {
  1922.                             if (0 != bi.digits[idx]) {
  1923.                                 if (0 == ++dblLo) {
  1924.                                     ++dblHi;
  1925.                                 }
  1926.                                 break;
  1927.                             }
  1928.                         }
  1929.                     }
  1930.                 }
  1931.                 return BitConverter.Int64BitsToDouble((long)dblHi << 32 | dblLo);
  1932.             }
  1933.             #endif
  1934.         }
  1935.        
  1936. /**
  1937.         * Floating point number represented in base-10.
  1938.         */       
  1939.         private class FloatingDecimal
  1940.         {
  1941.             public const int MaxDigits = 50;
  1942.             private const int MaxExp10 = 310;
  1943.             // Upper bound on base 10 exponent
  1944.             private const int MinExp10 = -325;
  1945.             // Lower bound on base 10 exponent
  1946.             private int exponent;
  1947.             // Base-10 scaling factor (0 means decimal point immediately precedes first digit)
  1948.             private int sign;
  1949.             // Sign is -1 or 1, depending on sign of number
  1950.             private int mantissaSize;
  1951.             // Size of mantissa
  1952.             private byte[] mantissa = new byte[MaxDigits];
  1953.             // Array of base-10 digits
  1954.             public int Exponent {
  1955.                 get { return exponent; }
  1956.                 set { exponent = value; }
  1957.             }
  1958.             public int Sign {
  1959.                 get { return sign; }
  1960.                 set { sign = value; }
  1961.             }
  1962.             public byte[] Mantissa {
  1963.                 get { return mantissa; }
  1964.             }
  1965.            
  1966.             public int MantissaSize {
  1967.                 get { return mantissaSize; }
  1968.                 set {
  1969.                     Debug.Assert(value <= MaxDigits);
  1970.                     mantissaSize = value;
  1971.                 }
  1972.             }
  1973.            
  1974.             public byte this[int ib]
  1975.             {
  1976.                 get {
  1977.                     Debug.Assert(0 <= ib && ib < mantissaSize);
  1978.                     return mantissa[ib];
  1979.                 }
  1980.             }
  1981.            
  1982.             public FloatingDecimal()
  1983.             {
  1984.                 exponent = 0;
  1985.                 sign = 1;
  1986.                 mantissaSize = 0;
  1987.             }
  1988.            
  1989.             public FloatingDecimal(double dbl)
  1990.             {
  1991.                 InitFromDouble(dbl);
  1992.                
  1993.                 #if DEBUG
  1994.                 if (0 != mantissaSize) {
  1995.                     Debug.Assert(dbl == (double)this);
  1996.                    
  1997.                     FloatingDecimal decAfter = new FloatingDecimal();
  1998.                     decAfter.InitFromDouble(Succ(dbl));
  1999.                     // Assert(memcmp(this, &decAfter, sizeof(*this) - MaxDigits + mantissaSize));
  2000.                     Debug.Assert(!this.Equals(decAfter));
  2001.                    
  2002.                     FloatingDecimal decBefore = new FloatingDecimal();
  2003.                     decBefore.InitFromDouble(Pred(dbl));
  2004.                     // Assert(memcmp(this, &decBefore, sizeof(*this) - MaxDigits + mantissaSize));
  2005.                     Debug.Assert(!this.Equals(decBefore));
  2006.                 }
  2007.                 #endif
  2008.             }
  2009.            
  2010.             #if DEBUG
  2011.             private bool Equals(FloatingDecimal other)
  2012.             {
  2013.                 if (exponent != other.exponent || sign != other.sign || mantissaSize != other.mantissaSize) {
  2014.                     return false;
  2015.                 }
  2016.                 for (int idx = 0; idx < mantissaSize; idx++) {
  2017.                     if (mantissa[idx] != other.mantissa[idx]) {
  2018.                         return false;
  2019.                     }
  2020.                 }
  2021.                 return true;
  2022.             }
  2023.             #endif
  2024.            
  2025.             #if NEVER
  2026. /*  ----------------------------------------------------------------------------
  2027.                 RoundTo()
  2028.                 Rounds off the BCD representation of a number to a specified number of digits.
  2029.                 This may result in the exponent being incremented (e.g. if digits were 999).
  2030.             */           
  2031.             public void RoundTo(int sizeMantissa)
  2032.             {
  2033.                 if (sizeMantissa >= mantissaSize) {
  2034.                     // No change required
  2035.                     return;
  2036.                 }
  2037.                
  2038.                 if (sizeMantissa >= 0) {
  2039.                     bool fRoundUp = mantissa[sizeMantissa] >= 5;
  2040.                     mantissaSize = sizeMantissa;
  2041.                    
  2042.                     // Round up if necessary and trim trailing zeros
  2043.                     for (int idx = mantissaSize - 1; idx >= 0; idx--) {
  2044.                         if (fRoundUp) {
  2045.                             if (++(mantissa[idx]) <= 9) {
  2046.                                 // Trailing digit is non-zero, so break
  2047.                                 fRoundUp = false;
  2048.                                 break;
  2049.                             }
  2050.                         }
  2051.                         else if (mantissa[idx] > 0) {
  2052.                             // Trailing digit is non-zero, so break
  2053.                             break;
  2054.                         }
  2055.                        
  2056.                         // Trim trailing zeros
  2057.                         mantissaSize--;
  2058.                     }
  2059.                    
  2060.                     if (fRoundUp) {
  2061.                         // Number consisted only of 9's
  2062.                         Debug.Assert(0 == mantissaSize);
  2063.                         mantissa[0] = 1;
  2064.                         mantissaSize = 1;
  2065.                         exponent++;
  2066.                     }
  2067.                 }
  2068.                 else {
  2069.                     // Number was rounded past any significant digits (e.g. 0.001 rounded to 1 fractional place), so round to 0.0
  2070.                     mantissaSize = 0;
  2071.                 }
  2072.                
  2073.                 if (0 == mantissaSize) {
  2074.                     // 0.0
  2075.                     sign = 1;
  2076.                     exponent = 0;
  2077.                 }
  2078.             }
  2079.             #endif
  2080.            
  2081.             #if !NOPARSE || DEBUG
  2082. /*  ----------------------------------------------------------------------------
  2083.                 explicit operator double()
  2084.                 Returns the double value of this floating decimal.
  2085.             */           
  2086.             public static explicit operator double(FloatingDecimal dec)
  2087.             {
  2088.                 BigNumber num;
  2089.                 BigNumber numHi;
  2090.                 BigNumber numLo;
  2091.                 uint ul;
  2092.                 int scale;
  2093.                 double dbl;
  2094.                 double dblLowPrec;
  2095.                 double dblLo;
  2096.                 int mantissaSize = dec.mantissaSize;
  2097.                
  2098.                 // Verify that there are no leading or trailing zeros in the mantissa
  2099.                 Debug.Assert(0 != mantissaSize && 0 != dec[0] && 0 != dec[mantissaSize - 1]);
  2100.                
  2101.                 // See if we can just use IEEE double arithmetic.
  2102.                 scale = dec.exponent - mantissaSize;
  2103.                 if (mantissaSize <= 15 && scale >= -22 && dec.exponent <= 37) {
  2104.                     // These calculations are all exact since mantissaSize <= 15.
  2105.                     if (mantissaSize <= 9) {
  2106.                         // Can use the ALU to perform fast integer arithmetic
  2107.                         ul = 0;
  2108.                         for (int ib = 0; ib < mantissaSize; ib++) {
  2109.                             Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
  2110.                             ul = ul * 10 + dec[ib];
  2111.                         }
  2112.                         dbl = ul;
  2113.                     }
  2114.                     else {
  2115.                         // Use floating point arithmetic
  2116.                         dbl = 0.0;
  2117.                         for (int ib = 0; ib < mantissaSize; ib++) {
  2118.                             Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
  2119.                             dbl = dbl * 10.0 + dec[ib];
  2120.                         }
  2121.                     }
  2122.                    
  2123.                     // This is the only (potential) rounding operation and we assume
  2124.                     // the compiler does the correct IEEE rounding.
  2125.                     if (scale > 0) {
  2126.                         // Need to scale upwards by powers of 10
  2127.                         if (scale > 22) {
  2128.                             // This one is exact. We're using the fact that mantissaSize < 15
  2129.                             // to handle exponents bigger than 22.
  2130.                             dbl *= C10toN[scale - 22];
  2131.                             dbl *= C10toN[22];
  2132.                         }
  2133.                         else {
  2134.                             dbl *= C10toN[scale];
  2135.                         }
  2136.                     }
  2137.                     else if (scale < 0) {
  2138.                         // Scale number by negative power of 10
  2139.                         dbl /= C10toN[-scale];
  2140.                     }
  2141.                    
  2142.                     #if DEBUG
  2143.                     // In the debug version, execute the high precision code also and
  2144.                     // verify that the results are the same.
  2145.                     dblLowPrec = dbl;
  2146.                 }
  2147.                 else {
  2148.                     dblLowPrec = Double.NaN;
  2149.                 }
  2150.                 #else
  2151.                 goto LDone;
  2152.             }
  2153.         }
  2154.     }
  2155. }
  2156. #endif

Developer Fusion