mcs/class/referencesource/System.Data.SqlXml/System/Xml/Xsl/XPathConvert.cs

   1 //------------------------------------------------------------------------------
   2 // <copyright file="XPathConvert.cs" company="Microsoft">
   3 //     Copyright (c) Microsoft Corporation.  All rights reserved.
   4 // </copyright>
   5 // <owner current="true" primary="true">Microsoft</owner>
   6 //------------------------------------------------------------------------------
   7
   8 /**
   9 * Routines used to manipulate IEEE 754 double-precision numbers, taken from JScript codebase.
  10 *
  11 * Define NOPARSE if you do not need FloatingDecimal -> double conversions
  12 *
  13 */
  14
  15 using System.Diagnostics;
  16 using System.Globalization;
  17
  18 namespace System.Xml.Xsl {
  19
  20     /**
  21     * Converts according to XPath/XSLT rules.
  22     */
  23     internal static class XPathConvert {
  24
  25         public static uint DblHi(double dbl) {
  26             return (uint)(BitConverter.DoubleToInt64Bits(dbl) >> 32);
  27         }
  28
  29         public static uint DblLo(double dbl) {
  30             return (uint)BitConverter.DoubleToInt64Bits(dbl);
  31         }
  32
  33         // Returns true if value is infinite or NaN (exponent bits are all ones)
  34         public static bool IsSpecial(double dbl) {
  35             return 0 == (~DblHi(dbl) & 0x7FF00000);
  36         }
  37
  38     #if DEBUG
  39         // Returns the next representable neighbor of x in the direction toward y
  40         public static double NextAfter(double x, double y) {
  41             long bits;
  42
  43             if (Double.IsNaN(x)) {
  44                 return x;
  45             }
  46             if (Double.IsNaN(y)) {
  47                 return y;
  48             }
  49             if (x == y) {
  50                 return y;
  51             }
  52             if (x == 0) {
  53                 bits = BitConverter.DoubleToInt64Bits(y) & 1L<<63;
  54                 return BitConverter.Int64BitsToDouble(bits | 1);
  55             }
  56
  57             // At this point x!=y, and x!=0. x can be treated as a 64bit
  58             // integer in sign/magnitude representation. To get the next
  59             // representable neighbor we add or subtract one from this
  60             // integer.
  61
  62             bits = BitConverter.DoubleToInt64Bits(x);
  63             if (0 < x && x < y || 0 > x && x > y) {
  64                 bits++;
  65             } else {
  66                 bits--;
  67             }
  68             return BitConverter.Int64BitsToDouble(bits);
  69         }
  70
  71         public static double Succ(double x) {
  72             return NextAfter(x, Double.PositiveInfinity);
  73         }
  74
  75         public static double Pred(double x) {
  76             return NextAfter(x, Double.NegativeInfinity);
  77         }
  78     #endif
  79
  80         // Small powers of ten. These are all the powers of ten that have an exact
  81         // representation in IEEE double precision format.
  82         public static readonly double[] C10toN = {
  83             1e00, 1e01, 1e02, 1e03, 1e04, 1e05, 1e06, 1e07, 1e08, 1e09,
  84             1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
  85             1e20, 1e21, 1e22,
  86         };
  87
  88         // Returns 1 if argument is non-zero, and 0 otherwise
  89         public static uint NotZero(uint u) {
  90             return 0 != u ? 1u : 0u;
  91         }
  92
  93         /*  ----------------------------------------------------------------------------
  94             AddU()
  95
  96             Add two unsigned ints and return the carry bit.
  97         */
  98         public static uint AddU(ref uint u1, uint u2) {
  99             u1 += u2;
 100             return u1 < u2 ? 1u : 0u;
 101         }
 102
 103         /*  ----------------------------------------------------------------------------
 104             MulU()
 105
 106             Multiply two unsigned ints. Return the low uint and fill uHi with
 107             the high uint.
 108         */
 109         public static uint MulU(uint u1, uint u2, out uint uHi) {
 110             ulong result = (ulong)u1 * u2;
 111             uHi = (uint)(result >> 32);
 112             return (uint)result;
 113         }
 114
 115         /*  ----------------------------------------------------------------------------
 116             CbitZeroLeft()
 117
 118             Return a count of the number of leading 0 bits in u.
 119         */
 120         public static int CbitZeroLeft(uint u) {
 121             int cbit = 0;
 122
 123             if (0 == (u & 0xFFFF0000)) {
 124                 cbit += 16;
 125                 u <<= 16;
 126             }
 127             if (0 == (u & 0xFF000000)) {
 128                 cbit += 8;
 129                 u <<= 8;
 130             }
 131             if (0 == (u & 0xF0000000)) {
 132                 cbit += 4;
 133                 u <<= 4;
 134             }
 135             if (0 == (u & 0xC0000000)) {
 136                 cbit += 2;
 137                 u <<= 2;
 138             }
 139             if (0 == (u & 0x80000000)) {
 140                 cbit += 1;
 141                 u <<= 1;
 142             }
 143             Debug.Assert(0 != (u & 0x80000000));
 144
 145             return cbit;
 146         }
 147
 148         /*  ----------------------------------------------------------------------------
 149             IsInteger()
 150
 151             If dbl is a whole number in the range of INT_MIN to INT_MAX, return true
 152             and the integer in value.  Otherwise, return false.
 153         */
 154         public static bool IsInteger(double dbl, out int value) {
 155             if (!IsSpecial(dbl)) {
 156                 int i = (int) dbl;
 157                 double dblRound = (double) i;
 158
 159                 if (dbl == dblRound) {
 160                     value = i;
 161                     return true;
 162                 }
 163             }
 164             value = 0;
 165             return false;
 166         }
 167
 168         /**
 169         * Implementation of a big floating point number used to ensure adequate
 170         * precision when performing calculations.
 171         *
 172         * Hungarian: num
 173         *
 174         */
 175         private struct BigNumber {
 176             private uint u0;
 177             private uint u1;
 178             private uint u2;
 179             private int  exp;
 180
 181             // This is a bound on the absolute value of the error. It is based at
 182             // one bit before the least significant bit of u0.
 183             private uint error;
 184
 185             public uint Error { get { return error; } }
 186
 187             public BigNumber(uint u0, uint u1, uint u2, int exp, uint error) {
 188                 this.u0 = u0;
 189                 this.u1 = u1;
 190                 this.u2 = u2;
 191                 this.exp = exp;
 192                 this.error = error;
 193             }
 194
 195         #if !NOPARSE || DEBUG
 196             // Set the value from floating decimal
 197             public BigNumber(FloatingDecimal dec) {
 198                 Debug.Assert(dec.MantissaSize > 0);
 199
 200                 int  ib, exponent, mantissaSize, wT;
 201                 uint uExtra;
 202                 BigNumber[] TenPowers;
 203
 204                 ib = 0;
 205                 exponent = dec.Exponent;
 206                 mantissaSize = dec.MantissaSize;
 207
 208                 // Record the first digit
 209                 Debug.Assert(dec[ib] > 0 && dec[ib] <= 9);
 210                 this.u2 = (uint)(dec[ib]) << 28;
 211                 this.u1 = 0;
 212                 this.u0 = 0;
 213                 this.exp = 4;
 214                 this.error = 0;
 215                 exponent--;
 216                 Normalize();
 217
 218                 while (++ib < mantissaSize) {
 219                     Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
 220                     uExtra = MulTenAdd(dec[ib]);
 221                     exponent--;
 222                     if (0 != uExtra) {
 223                         // We've filled up our precision.
 224                         Round(uExtra);
 225                         if (ib < mantissaSize - 1) {
 226                             // There are more digits, so add another error bit just for
 227                             // safety's sake.
 228                             this.error++;
 229                         }
 230                         break;
 231                     }
 232                 }
 233
 234                 // Now multiply by 10^exp
 235                 if (0 == exponent) {
 236                     return;
 237                 }
 238
 239                 if (exponent < 0) {
 240                     TenPowers = TenPowersNeg;
 241                     exponent = -exponent;
 242                 } else {
 243                     TenPowers = TenPowersPos;
 244                 }
 245
 246                 Debug.Assert(exponent > 0 && exponent < 512);
 247                 wT = exponent & 0x1F;
 248                 if (wT > 0) {
 249                     Mul(ref TenPowers[wT - 1]);
 250                 }
 251
 252                 wT = (exponent >> 5) & 0x0F;
 253                 if (wT > 0) {
 254                     Mul(ref TenPowers[wT + 30]);
 255                 }
 256             }
 257
 258             // Multiply by ten and add a base 10 digit.
 259             private unsafe uint MulTenAdd(uint digit) {
 260                 Debug.Assert(digit <= 9);
 261                 Debug.Assert(0 != (this.u2 & 0x80000000));
 262
 263                 // First "multipy" by eight
 264                 this.exp += 3;
 265                 Debug.Assert(this.exp >= 4);
 266
 267                 // Initialize the carry values based on digit and exp.
 268                 uint *rgu = stackalloc uint[5];
 269                 for (int i = 0; i < 5; i++) {
 270                     rgu[i] = 0;
 271                 }
 272                 if (0 != digit) {
 273                     int idx = 3 - (this.exp >> 5);
 274                     if (idx < 0) {
 275                         rgu[0] = 1;
 276                     } else {
 277                         int ibit = this.exp & 0x1F;
 278                         if (ibit < 4) {
 279                             rgu[idx + 1] = digit >> ibit;
 280                             if (ibit > 0) {
 281                                 rgu[idx] = digit << (32 - ibit);
 282                             }
 283                         } else {
 284                             rgu[idx] = digit << (32 - ibit);
 285                         }
 286                     }
 287                 }
 288
 289                 // Shift and add to multiply by ten.
 290                 rgu[1] += AddU(ref rgu[0], this.u0 << 30);
 291                 rgu[2] += AddU(ref this.u0, (this.u0 >> 2) + (this.u1 << 30));
 292                 if (0 != rgu[1]) {
 293                     rgu[2] += AddU(ref this.u0, rgu[1]);
 294                 }
 295                 rgu[3] += AddU(ref this.u1, (this.u1 >> 2) + (this.u2 << 30));
 296                 if (0 != rgu[2]) {
 297                     rgu[3] += AddU(ref this.u1, rgu[2]);
 298                 }
 299                 rgu[4] = AddU(ref this.u2, (this.u2 >> 2) + rgu[3]);
 300
 301                 // Handle the final carry.
 302                 if (0 != rgu[4]) {
 303                     Debug.Assert(1 == rgu[4]);
 304                     rgu[0] = (rgu[0] >> 1) | (rgu[0] & 1) | (this.u0 << 31);
 305                     this.u0 = (this.u0 >> 1) | (this.u1 << 31);
 306                     this.u1 = (this.u1 >> 1) | (this.u2 << 31);
 307                     this.u2 = (this.u2 >> 1) | 0x80000000;
 308                     this.exp++;
 309                 }
 310
 311                 return rgu[0];
 312             }
 313
 314             // Round based on the given extra data using IEEE round to nearest rule.
 315             private void Round(uint uExtra) {
 316                 if (0 == (uExtra & 0x80000000) || 0 == (uExtra & 0x7FFFFFFF) && 0 == (this.u0 & 1)) {
 317                     if (0 != uExtra) {
 318                         this.error++;
 319                     }
 320                     return;
 321                 }
 322                 this.error++;
 323
 324                 // Round up.
 325                 if (0 != AddU(ref this.u0, 1) && 0 != AddU(ref this.u1, 1) && 0 != AddU(ref this.u2, 1)) {
 326                     Debug.Assert(this.IsZero);
 327                     this.u2 = 0x80000000;
 328                     this.exp++;
 329                 }
 330             }
 331         #endif
 332
 333             // Test to see if the num is zero. This works even if we're not normalized.
 334             private bool IsZero {
 335                 get {
 336                     return (0 == this.u2) && (0 == this.u1) && (0 == this.u0);
 337                 }
 338             }
 339
 340             /*  ----------------------------------------------------------------------------
 341                 Normalize()
 342
 343                 Normalize the big number - make sure the high bit is 1 or everything is zero
 344                 (including the exponent).
 345             */
 346             private void Normalize() {
 347                 int w1, w2;
 348
 349                 // Normalize mantissa
 350                 if (0 == this.u2) {
 351                     if (0 == this.u1) {
 352                         if (0 == this.u0) {
 353                             this.exp = 0;
 354                             return;
 355                         }
 356                         this.u2 = this.u0;
 357                         this.u0 = 0;
 358                         this.exp -= 64;
 359                     } else {
 360                         this.u2 = this.u1;
 361                         this.u1 = this.u0;
 362                         this.u0 = 0;
 363                         this.exp -= 32;
 364                     }
 365                 }
 366
 367                 if (0 != (w1 = CbitZeroLeft(this.u2))) {
 368                     w2 = 32 - w1;
 369                     this.u2 = (this.u2 << w1) | (this.u1 >> w2);
 370                     this.u1 = (this.u1 << w1) | (this.u0 >> w2);
 371                     this.u0 = (this.u0 << w1);
 372                     this.exp -= w1;
 373                 }
 374             }
 375
 376             /*  ----------------------------------------------------------------------------
 377                 Mul()
 378
 379                 Multiply this big number by another big number.
 380             */
 381             private void Mul(ref BigNumber numOp) {
 382                 Debug.Assert(0 != (this.u2 & 0x80000000));
 383                 Debug.Assert(0 != (numOp.u2 & 0x80000000));
 384
 385                 //uint *rgu = stackalloc uint[6];
 386                 uint rgu0 = 0, rgu1 = 0, rgu2 = 0, rgu3 = 0, rgu4 = 0, rgu5 = 0;
 387                 uint uLo, uHi, uT;
 388                 uint wCarry;
 389
 390                 if (0 != (uT = this.u0)) {
 391                     uLo = MulU(uT, numOp.u0, out uHi);
 392                     rgu0 = uLo;
 393                     rgu1 = uHi;
 394
 395                     uLo = MulU(uT, numOp.u1, out uHi);
 396                     Debug.Assert(uHi < 0xFFFFFFFF);
 397                     wCarry = AddU(ref rgu1, uLo);
 398                     AddU(ref rgu2, uHi + wCarry);
 399
 400                     uLo = MulU(uT, numOp.u2, out uHi);
 401                     Debug.Assert(uHi < 0xFFFFFFFF);
 402                     wCarry = AddU(ref rgu2, uLo);
 403                     AddU(ref rgu3, uHi + wCarry);
 404                 }
 405
 406                 if (0 != (uT = this.u1)) {
 407                     uLo = MulU(uT, numOp.u0, out uHi);
 408                     Debug.Assert(uHi < 0xFFFFFFFF);
 409                     wCarry = AddU(ref rgu1, uLo);
 410                     wCarry = AddU(ref rgu2, uHi + wCarry);
 411                     if (0 != wCarry && 0 != AddU(ref rgu3, 1)) {
 412                         AddU(ref rgu4, 1);
 413                     }
 414                     uLo = MulU(uT, numOp.u1, out uHi);
 415                     Debug.Assert(uHi < 0xFFFFFFFF);
 416                     wCarry = AddU(ref rgu2, uLo);
 417                     wCarry = AddU(ref rgu3, uHi + wCarry);
 418                     if (0 != wCarry) {
 419                         AddU(ref rgu4, 1);
 420                     }
 421                     uLo = MulU(uT, numOp.u2, out uHi);
 422                     Debug.Assert(uHi < 0xFFFFFFFF);
 423                     wCarry = AddU(ref rgu3, uLo);
 424                     AddU(ref rgu4, uHi + wCarry);
 425                 }
 426
 427                 uT = this.u2;
 428                 Debug.Assert(0 != uT);
 429                 uLo = MulU(uT, numOp.u0, out uHi);
 430                 Debug.Assert(uHi < 0xFFFFFFFF);
 431                 wCarry = AddU(ref rgu2, uLo);
 432                 wCarry = AddU(ref rgu3, uHi + wCarry);
 433                 if (0 != wCarry && 0 != AddU(ref rgu4, 1)) {
 434                     AddU(ref rgu5, 1);
 435                 }
 436                 uLo = MulU(uT, numOp.u1, out uHi);
 437                 Debug.Assert(uHi < 0xFFFFFFFF);
 438                 wCarry = AddU(ref rgu3, uLo);
 439                 wCarry = AddU(ref rgu4, uHi + wCarry);
 440                 if (0 != wCarry) {
 441                     AddU(ref rgu5, 1);
 442                 }
 443                 uLo = MulU(uT, numOp.u2, out uHi);
 444                 Debug.Assert(uHi < 0xFFFFFFFF);
 445                 wCarry = AddU(ref rgu4, uLo);
 446                 AddU(ref rgu5, uHi + wCarry);
 447
 448                 // Compute the new exponent
 449                 this.exp += numOp.exp;
 450
 451                 // Accumulate the error. Adding doesn't necessarily give an accurate
 452                 // bound if both of the errors are bigger than 2.
 453                 Debug.Assert(this.error <= 2 || numOp.error <= 2);
 454                 this.error += numOp.error;
 455
 456                 // Handle rounding and normalize.
 457                 if (0 == (rgu5 & 0x80000000)) {
 458                     if (0 != (rgu2 & 0x40000000) &&
 459                             (0 != (rgu2 & 0xBFFFFFFF) || 0 != rgu1 || 0 != rgu0)) {
 460                         // Round up by 1
 461                         if (0 != AddU(ref rgu2, 0x40000000)
 462                             && 0 != AddU(ref rgu3, 1)
 463                             && 0 != AddU(ref rgu4, 1)
 464                         ) {
 465                             AddU(ref rgu5, 1);
 466                             if (0 != (rgu5 & 0x80000000)) {
 467                                 goto LNormalized;
 468                             }
 469                         }
 470                     }
 471
 472                     // have to shift by one
 473                     Debug.Assert(0 != (rgu5 & 0x40000000));
 474                     this.u2 = (rgu5 << 1) | (rgu4 >> 31);
 475                     this.u1 = (rgu4 << 1) | (rgu3 >> 31);
 476                     this.u0 = (rgu3 << 1) | (rgu2 >> 31);
 477                     this.exp--;
 478                     this.error <<= 1;
 479
 480                     // Add one for the error.
 481                     if (0 != (rgu2 & 0x7FFFFFFF) || 0 != rgu1 || 0 != rgu0) {
 482                         this.error++;
 483                     }
 484                     return;
 485                 } else {
 486                     if (0 != (rgu2 & 0x80000000) &&
 487                         (0 != (rgu3 & 1) || 0 != (rgu2 & 0x7FFFFFFF) ||
 488                             0 != rgu1 || 0 != rgu0)) {
 489                         // Round up by 1
 490                         if (0 != AddU(ref rgu3, 1) && 0 != AddU(ref rgu4, 1) && 0 != AddU(ref rgu5, 1)) {
 491                             Debug.Assert(0 == rgu3);
 492                             Debug.Assert(0 == rgu4);
 493                             Debug.Assert(0 == rgu5);
 494                             rgu5 = 0x80000000;
 495                             this.exp++;
 496                         }
 497                     }
 498                 }
 499
 500             LNormalized:
 501                 this.u2 = rgu5;
 502                 this.u1 = rgu4;
 503                 this.u0 = rgu3;
 504
 505                 // Add one for the error.
 506                 if (0 != rgu2 || 0 != rgu1 || 0 != rgu0) {
 507                     this.error++;
 508                 }
 509             }
 510
 511             // Get the double value.
 512             public static explicit operator double(BigNumber bn) {
 513                 uint uEx;
 514                 int exp;
 515                 uint dblHi, dblLo;
 516
 517                 Debug.Assert(0 != (bn.u2 & 0x80000000));
 518
 519                 exp = bn.exp + 1022;
 520                 if (exp >= 2047) {
 521                     return Double.PositiveInfinity;
 522                 }
 523
 524                 // Round after filling in the bits. In the extra uint, we set the low bit
 525                 // if there are any extra non-zero bits. This is for breaking the tie when
 526                 // deciding whether to round up or down.
 527                 if (exp > 0) {
 528                     // Normalized.
 529                     dblHi = ((uint)exp << 20) | ((bn.u2 & 0x7FFFFFFF) >> 11);
 530                     dblLo = bn.u2 << 21 | bn.u1 >> 11;
 531                     uEx = bn.u1 << 21 | NotZero(bn.u0);
 532                 } else if (exp > -20) {
 533                     // Denormal with some high bits.
 534                     int wT = 12 - exp;
 535                     Debug.Assert(wT >= 12 && wT < 32);
 536
 537                     dblHi = bn.u2 >> wT;
 538                     dblLo = (bn.u2 << (32 - wT)) | (bn.u1 >> wT);
 539                     uEx = (bn.u1 << (32 - wT)) | NotZero(bn.u0);
 540                 } else if (exp == -20) {
 541                     // Denormal with no high bits.
 542                     dblHi = 0;
 543                     dblLo = bn.u2;
 544                     uEx = bn.u1 | (uint)(0 != bn.u0 ? 1 : 0);
 545                 } else if (exp > -52) {
 546                     // Denormal with no high bits.
 547                     int wT = -exp - 20;
 548                     Debug.Assert(wT > 0 && wT < 32);
 549
 550                     dblHi = 0;
 551                     dblLo = bn.u2 >> wT;
 552                     uEx = bn.u2 << (32 - wT) | NotZero(bn.u1) | NotZero(bn.u0);
 553                 } else if (exp == -52) {
 554                     // Zero unless we round up below.
 555                     dblHi = 0;
 556                     dblLo = 0;
 557                     uEx = bn.u2 | NotZero(bn.u1) | NotZero(bn.u0);
 558                 } else {
 559                     return 0.0;
 560                 }
 561
 562                 // Handle rounding
 563                 if (0 != (uEx & 0x80000000) && (0 != (uEx & 0x7FFFFFFF) || 0 != (dblLo & 1))) {
 564                     // Round up. Note that this works even when we overflow into the
 565                     // exponent.
 566                     if (0 != AddU(ref dblLo, 1)) {
 567                         AddU(ref dblHi, 1);
 568                     }
 569                 }
 570                 return BitConverter.Int64BitsToDouble((long)dblHi << 32 | dblLo);
 571             }
 572
 573             // Lop off the integer part and return it.
 574             private uint UMod1() {
 575                 if (this.exp <= 0) {
 576                     return 0;
 577                 }
 578                 Debug.Assert(this.exp <= 32);
 579                 uint uT = this.u2 >> (32 - this.exp);
 580                 this.u2 &= (uint)0x7FFFFFFF >> (this.exp - 1);
 581                 Normalize();
 582                 return uT;
 583             }
 584
 585             // If error is not zero, add it and set error to zero.
 586             public void MakeUpperBound() {
 587                 Debug.Assert(this.error < 0xFFFFFFFF);
 588                 uint uT = (this.error + 1) >> 1;
 589
 590                 if (0 != uT && 0 != AddU(ref this.u0, uT) && 0 != AddU(ref this.u1, 1) && 0 != AddU(ref this.u2, 1)) {
 591                     Debug.Assert(0 == this.u2 && 0 == this.u1);
 592                     this.u2 = 0x80000000;
 593                     this.u0 = (this.u0 >> 1) + (this.u0 & 1);
 594                     this.exp++;
 595                 }
 596                 this.error = 0;
 597             }
 598
 599             // If error is not zero, subtract it and set error to zero.
 600             public void MakeLowerBound() {
 601                 Debug.Assert(this.error < 0xFFFFFFFF);
 602                 uint uT = (this.error + 1) >> 1;
 603
 604                 if (0 != uT && 0 == AddU(ref this.u0, (uint)-(int)uT) && 0 == AddU(ref this.u1, 0xFFFFFFFF)) {
 605                     AddU(ref this.u2, 0xFFFFFFFF);
 606                     if (0 == (0x80000000 & this.u2)) {
 607                         Normalize();
 608                     }
 609                 }
 610                 this.error = 0;
 611             }
 612
 613             /*  ----------------------------------------------------------------------------
 614                 DblToRgbFast()
 615
 616                 Get mantissa bytes (BCD).
 617             */
 618             public static bool DblToRgbFast(double dbl, byte[] mantissa, out int exponent, out int mantissaSize) {
 619                 BigNumber numHH, numHL, numLH, numLL;
 620                 BigNumber numBase;
 621                 BigNumber tenPower;
 622                 int ib, iT;
 623                 uint uT, uScale;
 624                 byte bHH, bHL, bLH, bLL;
 625                 uint uHH, uHL, uLH, uLL;
 626                 int wExp2, wExp10 = 0;
 627                 double dblInt;
 628                 uint dblHi, dblLo;
 629
 630                 dblHi = DblHi(dbl);
 631                 dblLo = DblLo(dbl);
 632
 633                 // Caller should take care of 0, negative and non-finite values.
 634                 Debug.Assert(!IsSpecial(dbl));
 635                 Debug.Assert(0 < dbl);
 636
 637                 // Get numHH and numLL such that numLL < dbl < numHH and the
 638                 // difference between adjacent values is half the distance to the next
 639                 // representable value (in a double).
 640                 wExp2 = (int)((dblHi >> 20) & 0x07FF);
 641                 if (wExp2 > 0) {
 642                     // See if dbl is a small integer.
 643                     if (wExp2 >= 1023 && wExp2 <= 1075 && dbl == Math.Floor(dbl)) {
 644                         goto LSmallInt;
 645                     }
 646
 647                     // Normalized
 648                     numBase.u2 = 0x80000000 | ((dblHi & 0x000FFFFFF) << 11) | (dblLo >> 21);
 649                     numBase.u1 = dblLo << 11;
 650                     numBase.u0 = 0;
 651                     numBase.exp = wExp2 - 1022;
 652                     numBase.error = 0;
 653
 654                     // Get the upper bound
 655                     numHH = numBase;
 656                     numHH.u1 |= (1 << 10);
 657
 658                     // Get the lower bound. A power of 2 must be special cased.
 659                     numLL = numBase;
 660                     if (0x80000000 == numLL.u2 && 0 == numLL.u1) {
 661                         // Subtract (0x00000000, 0x00000200, 0x00000000). Same as adding
 662                         // (0xFFFFFFFF, 0xFFFFFE00, 0x00000000)
 663                         uT = 0xFFFFFE00;
 664                     } else {
 665                         // Subtract (0x00000000, 0x00000400, 0x00000000). Same as adding
 666                         // (0xFFFFFFFF, 0xFFFFFC00, 0x00000000)
 667                         uT = 0xFFFFFC00;
 668                     }
 669                     if (0 == AddU(ref numLL.u1, uT)) {
 670                         AddU(ref numLL.u2, 0xFFFFFFFF);
 671                         if (0 == (0x80000000 & numLL.u2)) {
 672                             numLL.Normalize();
 673                         }
 674                     }
 675                 } else {
 676                     // Denormal
 677                     numBase.u2 = dblHi & 0x000FFFFF;
 678                     numBase.u1 = dblLo;
 679                     numBase.u0 = 0;
 680                     numBase.exp = -1010;
 681                     numBase.error = 0;
 682
 683                     // Get the upper bound
 684                     numHH = numBase;
 685                     numHH.u0 = 0x80000000;
 686
 687                     // Get the lower bound
 688                     numLL = numHH;
 689                     if (0 == AddU(ref numLL.u1, 0xFFFFFFFF)) {
 690                         AddU(ref numLL.u2, 0xFFFFFFFF);
 691                     }
 692
 693                     numBase.Normalize();
 694                     numHH.Normalize();
 695                     numLL.Normalize();
 696                 }
 697
 698                 // Multiply by powers of ten until 0 < numHH.exp < 32.
 699                 if (numHH.exp >= 32) {
 700                     iT = (numHH.exp - 25) * 15 / -TenPowersNeg[45].exp;
 701                     Debug.Assert(iT >= 0 && iT < 16);
 702                     if (iT > 0) {
 703                         tenPower = TenPowersNeg[30 + iT];
 704                         Debug.Assert(numHH.exp + tenPower.exp > 1);
 705                         numHH.Mul(ref tenPower);
 706                         numLL.Mul(ref tenPower);
 707                         wExp10 += iT * 32;
 708                     }
 709
 710                     if (numHH.exp >= 32) {
 711                         iT = (numHH.exp - 25) * 32 / -TenPowersNeg[31].exp;
 712                         Debug.Assert(iT > 0 && iT <= 32);
 713                         tenPower = TenPowersNeg[iT - 1];
 714                         Debug.Assert(numHH.exp + tenPower.exp > 1);
 715                         numHH.Mul(ref tenPower);
 716                         numLL.Mul(ref tenPower);
 717                         wExp10 += iT;
 718                     }
 719                 } else if (numHH.exp < 1) {
 720                     iT = (25 - numHH.exp) * 15 / TenPowersPos[45].exp;
 721                     Debug.Assert(iT >= 0 && iT < 16);
 722                     if (iT > 0) {
 723                         tenPower = TenPowersPos[30 + iT];
 724                         Debug.Assert(numHH.exp + tenPower.exp <= 32);
 725                         numHH.Mul(ref tenPower);
 726                         numLL.Mul(ref tenPower);
 727                         wExp10 -= iT * 32;
 728                     }
 729
 730                     if (numHH.exp < 1) {
 731                         iT = (25 - numHH.exp) * 32 / TenPowersPos[31].exp;
 732                         Debug.Assert(iT > 0 && iT <= 32);
 733                         tenPower = TenPowersPos[iT - 1];
 734                         Debug.Assert(numHH.exp + tenPower.exp <= 32);
 735                         numHH.Mul(ref tenPower);
 736                         numLL.Mul(ref tenPower);
 737                         wExp10 -= iT;
 738                     }
 739                 }
 740
 741                 Debug.Assert(numHH.exp > 0 && numHH.exp < 32);
 742
 743                 // Get the upper and lower bounds for these.
 744                 numHL = numHH;
 745                 numHH.MakeUpperBound();
 746                 numHL.MakeLowerBound();
 747                 uHH = numHH.UMod1();
 748                 uHL = numHL.UMod1();
 749                 numLH = numLL;
 750                 numLH.MakeUpperBound();
 751                 numLL.MakeLowerBound();
 752                 uLH = numLH.UMod1();
 753                 uLL = numLL.UMod1();
 754                 Debug.Assert(uLL <= uLH && uLH <= uHL && uHL <= uHH);
 755
 756                 // Find the starting scale
 757                 uScale = 1;
 758                 if (uHH >= 100000000) {
 759                     uScale = 100000000;
 760                     wExp10 += 8;
 761                 } else {
 762                     if (uHH >= 10000) {
 763                         uScale = 10000;
 764                         wExp10 += 4;
 765                     }
 766                     if (uHH >= 100 * uScale) {
 767                         uScale *= 100;
 768                         wExp10 += 2;
 769                     }
 770                 }
 771                 if (uHH >= 10 * uScale) {
 772                     uScale *= 10;
 773                     wExp10++;
 774                 }
 775                 wExp10++;
 776                 Debug.Assert(uHH >= uScale && uHH / uScale < 10);
 777
 778                 for (ib = 0; ; ) {
 779                     Debug.Assert(uLL <= uHH);
 780                     bHH = (byte)(uHH / uScale);
 781                     uHH %= uScale;
 782                     bLL = (byte)(uLL / uScale);
 783                     uLL %= uScale;
 784
 785                     if (bHH != bLL) {
 786                         break;
 787                     }
 788
 789                     Debug.Assert(0 != uHH || !numHH.IsZero);
 790                     mantissa[ib++] = bHH;
 791
 792                     if (1 == uScale) {
 793                         // Multiply by 10^8.
 794                         uScale = 10000000;
 795
 796                         numHH.Mul(ref TenPowersPos[7]);
 797                         numHH.MakeUpperBound();
 798                         uHH = numHH.UMod1();
 799                         if (uHH >= 100000000) {
 800                             goto LFail;
 801                         }
 802                         numHL.Mul(ref TenPowersPos[7]);
 803                         numHL.MakeLowerBound();
 804                         uHL = numHL.UMod1();
 805
 806                         numLH.Mul(ref TenPowersPos[7]);
 807                         numLH.MakeUpperBound();
 808                         uLH = numLH.UMod1();
 809                         numLL.Mul(ref TenPowersPos[7]);
 810                         numLL.MakeLowerBound();
 811                         uLL = numLL.UMod1();
 812                     } else {
 813                         uScale /= 10;
 814                     }
 815                 }
 816
 817                 // LL and HH diverged. Get the digit values for LH and HL.
 818                 Debug.Assert(0 <= bLL && bLL < bHH && bHH <= 9);
 819                 bLH = (byte)((uLH / uScale) % 10);
 820                 uLH %= uScale;
 821                 bHL = (byte)((uHL / uScale) % 10);
 822                 uHL %= uScale;
 823
 824                 if (bLH >= bHL) {
 825                     goto LFail;
 826                 }
 827
 828                 // LH and HL also diverged.
 829
 830                 // We can get by with one fewer digit if: LL == LH and bLH is zero
 831                 // and the current value of LH is zero and the least significant bit of
 832                 // the double is zero. In this case, we have exactly the digit sequence
 833                 // for the original numLL and IEEE and will rounds numLL up to the double.
 834                 if (0 == bLH && 0 == uLH && numLH.IsZero && 0 == (dblLo & 1)) {
 835                 }
 836                 else if (bHL - bLH > 1) {
 837                     // HL and LH differ by at least two in this digit, so split
 838                     // the difference.
 839                     mantissa[ib++] = (byte)((bHL + bLH + 1) / 2);
 840                 } else if (0 != uHL || !numHL.IsZero || 0 == (dblLo & 1)) {
 841                     // We can just use bHL because this guarantees that we're bigger than
 842                     // LH and less than HL, so must convert to the double.
 843                     mantissa[ib++] = bHL;
 844                 } else {
 845                     goto LFail;
 846                 }
 847
 848                 exponent = wExp10;
 849                 mantissaSize = ib;
 850                 goto LSucceed;
 851
 852             LSmallInt:
 853                 // dbl should be an integer from 1 to (2^53 - 1).
 854                 dblInt = dbl;
 855                 Debug.Assert(dblInt == Math.Floor(dblInt) && 1 <= dblInt && dblInt <= 9007199254740991.0d);
 856
 857                 iT = 0;
 858                 if (dblInt >= C10toN[iT + 8]) {
 859                     iT += 8;
 860                 }
 861                 if (dblInt >= C10toN[iT + 4]) {
 862                     iT += 4;
 863                 }
 864                 if (dblInt >= C10toN[iT + 2]) {
 865                     iT += 2;
 866                 }
 867                 if (dblInt >= C10toN[iT + 1]) {
 868                     iT += 1;
 869                 }
 870                 Debug.Assert(iT >= 0 && iT <= 15);
 871                 Debug.Assert(dblInt >= C10toN[iT] && dblInt < C10toN[iT + 1]);
 872                 exponent = iT + 1;
 873
 874                 for (ib = 0; 0 != dblInt; iT--) {
 875                     Debug.Assert(iT >= 0);
 876                     bHH = (byte)(dblInt / C10toN[iT]);
 877                     dblInt -= bHH * C10toN[iT];
 878                     Debug.Assert(dblInt == Math.Floor(dblInt) && 0 <= dblInt && dblInt < C10toN[iT]);
 879                     mantissa[ib++] = bHH;
 880                 }
 881                 mantissaSize = ib;
 882
 883             LSucceed:
 884
 885             #if DEBUG
 886                 // Verify that precise is working and gives the same answer
 887                 if (mantissaSize > 0) {
 888                     byte[] mantissaPrec = new byte[20];
 889                     int    exponentPrec, mantissaSizePrec, idx;
 890
 891                     DblToRgbPrecise(dbl, mantissaPrec, out exponentPrec, out mantissaSizePrec);
 892                     Debug.Assert(exponent == exponentPrec && mantissaSize == mantissaSizePrec);
 893                     // Assert(!memcmp(mantissa, mantissaPrec, mantissaSizePrec - 1));
 894                     bool equal = true;
 895                     for (idx = 0; idx < mantissaSize; idx++) {
 896                         equal &= (
 897                             (mantissa[idx] == mantissaPrec[idx]) ||
 898                             (idx == mantissaSize - 1) && Math.Abs(mantissa[idx] - mantissaPrec[idx]) <= 1
 899                         );
 900                     }
 901                     Debug.Assert(equal, "DblToRgbFast and DblToRgbPrecise should give the same result");
 902                 }
 903             #endif
 904
 905                 return true;
 906
 907             LFail:
 908                 exponent = mantissaSize = 0;
 909                 return false;
 910             }
 911
 912             /*  ----------------------------------------------------------------------------
 913                 DblToRgbPrecise()
 914
 915                 Uses big integer arithmetic to get the sequence of digits.
 916             */
 917             public static void DblToRgbPrecise(double dbl, byte[] mantissa, out int exponent, out int mantissaSize) {
 918                 exponent = 0;
 919                 mantissaSize = 0;
 920                 BigInteger biNum = new BigInteger();
 921                 BigInteger biDen = new BigInteger();
 922                 BigInteger biHi  = new BigInteger();
 923                 BigInteger biLo  = new BigInteger();
 924                 BigInteger biT   = new BigInteger();
 925                 BigInteger biHiLo;
 926                 byte bT;
 927                 bool fPow2;
 928                 int ib, cu;
 929                 int wExp10, wExp2, w1, w2;
 930                 int c2Num, c2Den, c5Num, c5Den;
 931                 double dblT;
 932                 //uint *rgu = stackalloc uint[2];
 933                 uint rgu0, rgu1;
 934                 uint dblHi, dblLo;
 935
 936                 dblHi = DblHi(dbl);
 937                 dblLo = DblLo(dbl);
 938
 939                 // Caller should take care of 0, negative and non-finite values.
 940                 Debug.Assert(!IsSpecial(dbl));
 941                 Debug.Assert(0 < dbl);
 942
 943                 // Init the Denominator, Hi error and Lo error bigints.
 944                 biDen.InitFromDigits(1, 0, 1);
 945                 biHi.InitFromDigits(1, 0, 1);
 946
 947                 wExp2 = (int)(((dblHi & 0x7FF00000) >> 20) - 1075);
 948                 rgu1 = dblHi & 0x000FFFFF;
 949                 rgu0 = dblLo;
 950                 cu = 2;
 951                 fPow2 = false;
 952                 if (wExp2 == -1075) {
 953                     // dbl is denormalized.
 954                     Debug.Assert(0 == (dblHi & 0x7FF00000));
 955                     if (0 == rgu1) {
 956                         cu = 1;
 957                     }
 958
 959                     // Get dblT such that dbl / dblT is a power of 2 and 1 <= dblT < 2.
 960                     // First multiply by a power of 2 to get a normalized value.
 961                     dblT = BitConverter.Int64BitsToDouble(0x4FF00000L << 32);
 962                     dblT *= dbl;
 963                     Debug.Assert(0 != (DblHi(dblT) & 0x7FF00000));
 964
 965                     // This is the power of 2.
 966                     w1 = (int)((DblHi(dblT) & 0x7FF00000) >> 20) - (256 + 1023);
 967
 968                     dblHi = DblHi(dblT);
 969                     dblHi &= 0x000FFFFF;
 970                     dblHi |= 0x3FF00000;
 971                     dblT = BitConverter.Int64BitsToDouble((long)dblHi << 32 | DblLo(dblT));
 972
 973                     // Adjust wExp2 because we don't have the implicit bit.
 974                     wExp2++;
 975                 } else {
 976                     // Get dblT such that dbl / dblT is a power of 2 and 1 <= dblT < 2.
 977                     // First multiply by a power of 2 to get a normalized value.
 978                     dblHi &= 0x000FFFFF;
 979                     dblHi |= 0x3FF00000;
 980                     dblT = BitConverter.Int64BitsToDouble((long)dblHi << 32 | dblLo);
 981
 982                     // This is the power of 2.
 983                     w1 = wExp2 + 52;
 984
 985                     if (0 == rgu0 && 0 == rgu1 && wExp2 > -1074) {
 986                         // Power of 2 bigger than smallest normal. The next smaller
 987                         // representable value is closer than the next larger value.
 988                         rgu1 = 0x00200000;
 989                         wExp2--;
 990                         fPow2 = true;
 991                     } else {
 992                         // Normalized and not a power of 2 or the smallest normal. The
 993                         // representable values on either side are the same distance away.
 994                         rgu1 |= 0x00100000;
 995                     }
 996                 }
 997
 998                 // Compute an approximation to the base 10 log. This is borrowed from
 999                 // David ----'s paper.
1000                 Debug.Assert(1 <= dblT && dblT < 2);
1001                 dblT = (dblT - 1.5) * 0.289529654602168 + 0.1760912590558 +
1002                     w1 * 0.301029995663981;
1003                 wExp10 = (int)dblT;
1004                 if (dblT < 0 && dblT != wExp10) {
1005                     wExp10--;
1006                 }
1007
1008                 if (wExp2 >= 0) {
1009                     c2Num = wExp2;
1010                     c2Den = 0;
1011                 } else {
1012                     c2Num = 0;
1013                     c2Den = -wExp2;
1014                 }
1015
1016                 if (wExp10 >= 0) {
1017                     c5Num = 0;
1018                     c5Den = wExp10;
1019                     c2Den += wExp10;
1020                 } else {
1021                     c2Num -= wExp10;
1022                     c5Num = -wExp10;
1023                     c5Den = 0;
1024                 }
1025
1026                 if (c2Num > 0 && c2Den > 0) {
1027                     w1 = c2Num < c2Den ? c2Num : c2Den;
1028                     c2Num -= w1;
1029                     c2Den -= w1;
1030                 }
1031                 // We need a bit for the Hi and Lo values.
1032                 c2Num++;
1033                 c2Den++;
1034
1035                 // Initialize biNum and multiply by powers of 5.
1036                 if (c5Num > 0) {
1037                     Debug.Assert(0 == c5Den);
1038                     biHi.MulPow5(c5Num);
1039                     biNum.InitFromBigint(biHi);
1040                     if (1 == cu) {
1041                         biNum.MulAdd(rgu0, 0);
1042                     } else {
1043                         biNum.MulAdd(rgu1, 0);
1044                         biNum.ShiftLeft(32);
1045                         if (0 != rgu0) {
1046                             biT.InitFromBigint(biHi);
1047                             biT.MulAdd(rgu0, 0);
1048                             biNum.Add(biT);
1049                         }
1050                     }
1051                 } else {
1052                     Debug.Assert(cu <= 2);
1053                     biNum.InitFromDigits(rgu0, rgu1, cu);
1054                     if (c5Den > 0) {
1055                         biDen.MulPow5(c5Den);
1056                     }
1057                 }
1058
1059                 // BigInteger.DivRem only works if the 4 high bits of the divisor are 0.
1060                 // It works most efficiently if there are exactly 4 zero high bits.
1061                 // Adjust c2Den and c2Num to guarantee this.
1062                 w1 = CbitZeroLeft(biDen[biDen.Length - 1]);
1063                 w1 = (w1 + 28 - c2Den) & 0x1F;
1064                 c2Num += w1;
1065                 c2Den += w1;
1066
1067                 // Multiply by powers of 2.
1068                 Debug.Assert(c2Num > 0 && c2Den > 0);
1069                 biNum.ShiftLeft(c2Num);
1070                 if (c2Num > 1) {
1071                     biHi.ShiftLeft(c2Num - 1);
1072                 }
1073                 biDen.ShiftLeft(c2Den);
1074                 Debug.Assert(0 == (biDen[biDen.Length - 1] & 0xF0000000));
1075                 Debug.Assert(0 != (biDen[biDen.Length - 1] & 0x08000000));
1076
1077                 // Get biHiLo and handle the power of 2 case where biHi needs to be doubled.
1078                 if (fPow2) {
1079                     biHiLo = biLo;
1080                     biHiLo.InitFromBigint(biHi);
1081                     biHi.ShiftLeft(1);
1082                 } else {
1083                     biHiLo = biHi;
1084                 }
1085
1086                 for (ib = 0; ; ) {
1087                     bT = (byte)biNum.DivRem(biDen);
1088                     if (0 == ib && 0 == bT) {
1089                         // Our estimate of wExp10 was too big. Oh well.
1090                         wExp10--;
1091                         goto LSkip;
1092                     }
1093
1094                     // w1 = sign(biNum - biHiLo).
1095                     w1 = biNum.CompareTo(biHiLo);
1096
1097                     // w2 = sign(biNum + biHi - biDen).
1098                     if (biDen.CompareTo(biHi) < 0) {
1099                         w2 = 1;
1100                     } else {
1101                         //
1102                         biT.InitFromBigint(biDen);
1103                         biT.Subtract(biHi);
1104                         w2 = biNum.CompareTo(biT);
1105                     }
1106
1107                     // if (biNum + biHi == biDen && even)
1108                     if (0 == w2 && 0 == (dblLo & 1)) {
1109                         // Rounding up this digit produces exactly (biNum + biHi) which
1110                         // StrToDbl will round down to dbl.
1111                         if (bT == 9) {
1112                             goto LRoundUp9;
1113                         }
1114                         if (w1 > 0) {
1115                             bT++;
1116                         }
1117                         mantissa[ib++] = bT;
1118                         break;
1119                     }
1120
1121                     // if (biNum < biHiLo || biNum == biHiLo && even)
1122                     if (w1 < 0 || 0 == w1 && 0 == (dblLo & 1)) {
1123                         // if (biNum + biHi > biDen)
1124                         if (w2 > 0) {
1125                             // Decide whether to round up.
1126                             biNum.ShiftLeft(1);
1127                             w2 = biNum.CompareTo(biDen);
1128                             if ((w2 > 0 || 0 == w2 && (0 != (bT & 1))) && bT++ == 9) {
1129                                 goto LRoundUp9;
1130                             }
1131                         }
1132                         mantissa[ib++] = bT;
1133                         break;
1134                     }
1135
1136                     // if (biNum + biHi > biDen)
1137                     if (w2 > 0) {
1138                         // Round up and be done with it.
1139                         if (bT != 9) {
1140                             mantissa[ib++] = (byte)(bT + 1);
1141                             break;
1142                         }
1143                         goto LRoundUp9;
1144                     }
1145
1146                     // Save the digit.
1147                     mantissa[ib++] = bT;
1148
1149             LSkip:
1150                     biNum.MulAdd(10, 0);
1151                     biHi.MulAdd(10, 0);
1152                     if ((object) biHiLo != (object) biHi) {
1153                         biHiLo.MulAdd(10, 0);
1154                     }
1155                     continue;
1156
1157             LRoundUp9:
1158                     while (ib > 0) {
1159                         if (mantissa[--ib] != 9) {
1160                             mantissa[ib++]++;
1161                             goto LReturn;
1162                         }
1163                     }
1164                     wExp10++;
1165                     mantissa[ib++] = 1;
1166                     break;
1167                 }
1168
1169             LReturn:
1170                 exponent = wExp10 + 1;
1171                 mantissaSize = ib;
1172             }
1173
1174             #region Powers of ten
1175             private static readonly BigNumber[] TenPowersPos = new BigNumber[46] {
1176                 // Positive powers of 10 to 96 bits precision.
1177                 new BigNumber( 0x00000000, 0x00000000, 0xA0000000,     4, 0 ), // 10^1
1178                 new BigNumber( 0x00000000, 0x00000000, 0xC8000000,     7, 0 ), // 10^2
1179                 new BigNumber( 0x00000000, 0x00000000, 0xFA000000,    10, 0 ), // 10^3
1180                 new BigNumber( 0x00000000, 0x00000000, 0x9C400000,    14, 0 ), // 10^4
1181                 new BigNumber( 0x00000000, 0x00000000, 0xC3500000,    17, 0 ), // 10^5
1182                 new BigNumber( 0x00000000, 0x00000000, 0xF4240000,    20, 0 ), // 10^6
1183                 new BigNumber( 0x00000000, 0x00000000, 0x98968000,    24, 0 ), // 10^7
1184                 new BigNumber( 0x00000000, 0x00000000, 0xBEBC2000,    27, 0 ), // 10^8
1185                 new BigNumber( 0x00000000, 0x00000000, 0xEE6B2800,    30, 0 ), // 10^9
1186                 new BigNumber( 0x00000000, 0x00000000, 0x9502F900,    34, 0 ), // 10^10
1187                 new BigNumber( 0x00000000, 0x00000000, 0xBA43B740,    37, 0 ), // 10^11
1188                 new BigNumber( 0x00000000, 0x00000000, 0xE8D4A510,    40, 0 ), // 10^12
1189                 new BigNumber( 0x00000000, 0x00000000, 0x9184E72A,    44, 0 ), // 10^13
1190                 new BigNumber( 0x00000000, 0x80000000, 0xB5E620F4,    47, 0 ), // 10^14
1191                 new BigNumber( 0x00000000, 0xA0000000, 0xE35FA931,    50, 0 ), // 10^15
1192                 new BigNumber( 0x00000000, 0x04000000, 0x8E1BC9BF,    54, 0 ), // 10^16
1193                 new BigNumber( 0x00000000, 0xC5000000, 0xB1A2BC2E,    57, 0 ), // 10^17
1194                 new BigNumber( 0x00000000, 0x76400000, 0xDE0B6B3A,    60, 0 ), // 10^18
1195                 new BigNumber( 0x00000000, 0x89E80000, 0x8AC72304,    64, 0 ), // 10^19
1196                 new BigNumber( 0x00000000, 0xAC620000, 0xAD78EBC5,    67, 0 ), // 10^20
1197                 new BigNumber( 0x00000000, 0x177A8000, 0xD8D726B7,    70, 0 ), // 10^21
1198                 new BigNumber( 0x00000000, 0x6EAC9000, 0x87867832,    74, 0 ), // 10^22
1199                 new BigNumber( 0x00000000, 0x0A57B400, 0xA968163F,    77, 0 ), // 10^23
1200                 new BigNumber( 0x00000000, 0xCCEDA100, 0xD3C21BCE,    80, 0 ), // 10^24
1201                 new BigNumber( 0x00000000, 0x401484A0, 0x84595161,    84, 0 ), // 10^25
1202                 new BigNumber( 0x00000000, 0x9019A5C8, 0xA56FA5B9,    87, 0 ), // 10^26
1203                 new BigNumber( 0x00000000, 0xF4200F3A, 0xCECB8F27,    90, 0 ), // 10^27
1204                 new BigNumber( 0x40000000, 0xF8940984, 0x813F3978,    94, 0 ), // 10^28
1205                 new BigNumber( 0x50000000, 0x36B90BE5, 0xA18F07D7,    97, 0 ), // 10^29
1206                 new BigNumber( 0xA4000000, 0x04674EDE, 0xC9F2C9CD,   100, 0 ), // 10^30
1207                 new BigNumber( 0x4D000000, 0x45812296, 0xFC6F7C40,   103, 0 ), // 10^31
1208                 new BigNumber( 0xF0200000, 0x2B70B59D, 0x9DC5ADA8,   107, 0 ), // 10^32
1209                 new BigNumber( 0x3CBF6B72, 0xFFCFA6D5, 0xC2781F49,   213, 1 ), // 10^64   (rounded up)
1210                 new BigNumber( 0xC5CFE94F, 0xC59B14A2, 0xEFB3AB16,   319, 1 ), // 10^96   (rounded up)
1211                 new BigNumber( 0xC66F336C, 0x80E98CDF, 0x93BA47C9,   426, 1 ), // 10^128
1212                 new BigNumber( 0x577B986B, 0x7FE617AA, 0xB616A12B,   532, 1 ), // 10^160
1213                 new BigNumber( 0x85BBE254, 0x3927556A, 0xE070F78D,   638, 1 ), // 10^192  (rounded up)
1214                 new BigNumber( 0x82BD6B71, 0xE33CC92F, 0x8A5296FF,   745, 1 ), // 10^224  (rounded up)
1215                 new BigNumber( 0xDDBB901C, 0x9DF9DE8D, 0xAA7EEBFB,   851, 1 ), // 10^256  (rounded up)
1216                 new BigNumber( 0x73832EEC, 0x5C6A2F8C, 0xD226FC19,   957, 1 ), // 10^288
1217                 new BigNumber( 0xE6A11583, 0xF2CCE375, 0x81842F29,  1064, 1 ), // 10^320
1218                 new BigNumber( 0x5EBF18B7, 0xDB900AD2, 0x9FA42700,  1170, 1 ), // 10^352  (rounded up)
1219                 new BigNumber( 0x1027FFF5, 0xAEF8AA17, 0xC4C5E310,  1276, 1 ), // 10^384
1220                 new BigNumber( 0xB5E54F71, 0xE9B09C58, 0xF28A9C07,  1382, 1 ), // 10^416
1221                 new BigNumber( 0xA7EA9C88, 0xEBF7F3D3, 0x957A4AE1,  1489, 1 ), // 10^448
1222                 new BigNumber( 0x7DF40A74, 0x0795A262, 0xB83ED8DC,  1595, 1 ), // 10^480
1223             };
1224
1225             private static readonly BigNumber[] TenPowersNeg = new BigNumber[46] {
1226                 // Negative powers of 10 to 96 bits precision.
1227                 new BigNumber( 0xCCCCCCCD, 0xCCCCCCCC, 0xCCCCCCCC,    -3, 1 ), // 10^-1   (rounded up)
1228                 new BigNumber( 0x3D70A3D7, 0x70A3D70A, 0xA3D70A3D,    -6, 1 ), // 10^-2
1229                 new BigNumber( 0x645A1CAC, 0x8D4FDF3B, 0x83126E97,    -9, 1 ), // 10^-3
1230                 new BigNumber( 0xD3C36113, 0xE219652B, 0xD1B71758,   -13, 1 ), // 10^-4
1231                 new BigNumber( 0x0FCF80DC, 0x1B478423, 0xA7C5AC47,   -16, 1 ), // 10^-5
1232                 new BigNumber( 0xA63F9A4A, 0xAF6C69B5, 0x8637BD05,   -19, 1 ), // 10^-6   (rounded up)
1233                 new BigNumber( 0x3D329076, 0xE57A42BC, 0xD6BF94D5,   -23, 1 ), // 10^-7
1234                 new BigNumber( 0xFDC20D2B, 0x8461CEFC, 0xABCC7711,   -26, 1 ), // 10^-8
1235                 new BigNumber( 0x31680A89, 0x36B4A597, 0x89705F41,   -29, 1 ), // 10^-9   (rounded up)
1236                 new BigNumber( 0xB573440E, 0xBDEDD5BE, 0xDBE6FECE,   -33, 1 ), // 10^-10
1237                 new BigNumber( 0xF78F69A5, 0xCB24AAFE, 0xAFEBFF0B,   -36, 1 ), // 10^-11
1238                 new BigNumber( 0xF93F87B7, 0x6F5088CB, 0x8CBCCC09,   -39, 1 ), // 10^-12
1239                 new BigNumber( 0x2865A5F2, 0x4BB40E13, 0xE12E1342,   -43, 1 ), // 10^-13
1240                 new BigNumber( 0x538484C2, 0x095CD80F, 0xB424DC35,   -46, 1 ), // 10^-14  (rounded up)
1241                 new BigNumber( 0x0F9D3701, 0x3AB0ACD9, 0x901D7CF7,   -49, 1 ), // 10^-15
1242                 new BigNumber( 0x4C2EBE68, 0xC44DE15B, 0xE69594BE,   -53, 1 ), // 10^-16
1243                 new BigNumber( 0x09BEFEBA, 0x36A4B449, 0xB877AA32,   -56, 1 ), // 10^-17  (rounded up)
1244                 new BigNumber( 0x3AFF322E, 0x921D5D07, 0x9392EE8E,   -59, 1 ), // 10^-18
1245                 new BigNumber( 0x2B31E9E4, 0xB69561A5, 0xEC1E4A7D,   -63, 1 ), // 10^-19  (rounded up)
1246                 new BigNumber( 0x88F4BB1D, 0x92111AEA, 0xBCE50864,   -66, 1 ), // 10^-20  (rounded up)
1247                 new BigNumber( 0xD3F6FC17, 0x74DA7BEE, 0x971DA050,   -69, 1 ), // 10^-21  (rounded up)
1248                 new BigNumber( 0x5324C68B, 0xBAF72CB1, 0xF1C90080,   -73, 1 ), // 10^-22
1249                 new BigNumber( 0x75B7053C, 0x95928A27, 0xC16D9A00,   -76, 1 ), // 10^-23
1250                 new BigNumber( 0xC4926A96, 0x44753B52, 0x9ABE14CD,   -79, 1 ), // 10^-24
1251                 new BigNumber( 0x3A83DDBE, 0xD3EEC551, 0xF79687AE,   -83, 1 ), // 10^-25  (rounded up)
1252                 new BigNumber( 0x95364AFE, 0x76589DDA, 0xC6120625,   -86, 1 ), // 10^-26
1253                 new BigNumber( 0x775EA265, 0x91E07E48, 0x9E74D1B7,   -89, 1 ), // 10^-27  (rounded up)
1254                 new BigNumber( 0x8BCA9D6E, 0x8300CA0D, 0xFD87B5F2,   -93, 1 ), // 10^-28
1255                 new BigNumber( 0x096EE458, 0x359A3B3E, 0xCAD2F7F5,   -96, 1 ), // 10^-29
1256                 new BigNumber( 0xA125837A, 0x5E14FC31, 0xA2425FF7,   -99, 1 ), // 10^-30  (rounded up)
1257                 new BigNumber( 0x80EACF95, 0x4B43FCF4, 0x81CEB32C,  -102, 1 ), // 10^-31  (rounded up)
1258                 new BigNumber( 0x67DE18EE, 0x453994BA, 0xCFB11EAD,  -106, 1 ), // 10^-32  (rounded up)
1259                 new BigNumber( 0x3F2398D7, 0xA539E9A5, 0xA87FEA27,  -212, 1 ), // 10^-64
1260                 new BigNumber( 0x11DBCB02, 0xFD75539B, 0x88B402F7,  -318, 1 ), // 10^-96
1261                 new BigNumber( 0xAC7CB3F7, 0x64BCE4A0, 0xDDD0467C,  -425, 1 ), // 10^-128 (rounded up)
1262                 new BigNumber( 0x59ED2167, 0xDB73A093, 0xB3F4E093,  -531, 1 ), // 10^-160
1263                 new BigNumber( 0x7B6306A3, 0x5423CC06, 0x91FF8377,  -637, 1 ), // 10^-192
1264                 new BigNumber( 0xA4F8BF56, 0x4A314EBD, 0xECE53CEC,  -744, 1 ), // 10^-224
1265                 new BigNumber( 0xFA911156, 0x637A1939, 0xC0314325,  -850, 1 ), // 10^-256 (rounded up)
1266                 new BigNumber( 0x4EE367F9, 0x836AC577, 0x9BECCE62,  -956, 1 ), // 10^-288
1267                 new BigNumber( 0x8920B099, 0x478238D0, 0xFD00B897, -1063, 1 ), // 10^-320 (rounded up)
1268                 new BigNumber( 0x0092757C, 0x46F34F7D, 0xCD42A113, -1169, 1 ), // 10^-352 (rounded up)
1269                 new BigNumber( 0x88DBA000, 0xB11B0857, 0xA686E3E8, -1275, 1 ), // 10^-384 (rounded up)
1270                 new BigNumber( 0x1A4EB007, 0x3FFC68A6, 0x871A4981, -1381, 1 ), // 10^-416 (rounded up)
1271                 new BigNumber( 0x84C663CF, 0xB6074244, 0xDB377599, -1488, 1 ), // 10^-448 (rounded up)
1272                 new BigNumber( 0x61EB52E2, 0x79007736, 0xB1D983B4, -1594, 1 ), // 10^-480
1273             };
1274             #endregion
1275
1276 #if false
1277             /***************************************************************************
1278                 This is JScript code to compute the BigNumber values in the tables above.
1279             ***************************************************************************/
1280             var cbits = 100;
1281             var cbitsExact = 96;
1282             var arrPos = new Array;
1283             var arrNeg = new Array;
1284
1285             function main()
1286             {
1287                 Compute(0, 31);
1288                 Compute(5, 15);
1289
1290                 print();
1291                 Dump()
1292             }
1293
1294             function Dump()
1295             {
1296                 var i;
1297                 for (i = 0; i < arrPos.length; i++)
1298                     print(arrPos[i]);
1299                 print();
1300                 for (i = 0; i < arrNeg.length; i++)
1301                     print(arrNeg[i]);
1302             }
1303
1304             function Compute(p, n)
1305             {
1306                 var t;
1307                 var i;
1308                 var exp = 1;
1309                 var str = '101';
1310
1311                 for (i = 0; i < p; i++)
1312                     {
1313                     exp *= 2;
1314                     str = Mul(str, str);
1315                     }
1316
1317                 PrintNum(str, exp);
1318                 PrintInv(str, exp);
1319
1320                 t = str;
1321                 for (i = 2; i <= n; i++)
1322                     {
1323                     t = Mul(str, t);
1324                     PrintNum(t, i * exp);
1325                     PrintInv(t, i * exp);
1326                     }
1327             }
1328
1329             function Mul(a, b)
1330             {
1331                 var c;
1332                 var len;
1333                 var i;
1334                 var res;
1335                 var add;
1336
1337                 if (a.length > b.length)
1338                     {
1339                     c = a;
1340                     a = b;
1341                     b = c;
1342                     }
1343
1344                 res = '';
1345                 add = b;
1346                 len = a.length;
1347                 for (i = 1; i <= len; i++)
1348                     {
1349                     if (a.charAt(len - i) == '1')
1350                         res = Add(res, add);
1351                     add += '0';
1352                     }
1353
1354                 return res;
1355             }
1356
1357             function Add(a, b)
1358             {
1359                 var bit;
1360                 var i;
1361                 var c = 0;
1362                 var res = '';
1363                 var lena = a.length;
1364                 var lenb = b.length;
1365                 var lenm = Math.max(lena, lenb);
1366
1367                 for (i = 1; i <= lenm; i++)
1368                     {
1369                     bit = (a.charAt(lena - i) == '1') + (b.charAt(lenb - i) == '1') + c;
1370                     if (bit & 1)
1371                         res = '1' + res;
1372                     else
1373                         res = '0' + res;
1374                     c = bit >> 1;
1375                     }
1376                 if (c)
1377                     res = '1' + res;
1378
1379                 return res;
1380             }
1381
1382             function PrintNum(a, n)
1383             {
1384                 arrPos[arrPos.length] = PrintHex(a, a.length, n);
1385             }
1386
1387             function PrintHex(a, e, n)
1388             {
1389                 var arr;
1390                 var i;
1391                 var dig;
1392                 var fRoundUp = false;
1393                 var res = '';
1394
1395                 dig = 0;
1396                 for (i = 0; i < cbitsExact; )
1397                     {
1398                     dig *= 2;
1399                     if (a.charAt(i) == '1')
1400                         dig++;
1401                     if (0 == (++i & 0x1F) && i < cbitsExact)
1402                         {
1403                         strT = dig.toString(16).toUpperCase();
1404                         res += ' 0x' + '00000000'.substring(strT.length) + strT;
1405                         dig = 0;
1406                         }
1407                     }
1408
1409                 if (a.charAt(cbitsExact) == '1')
1410                     {
1411                     // Round up. Don't have to worry about overflowing.
1412                     fRoundUp = true;
1413                     dig++;
1414                     }
1415                 strT = dig.toString(16).toUpperCase();
1416                 res += ' 0x' + '00000000'.substring(strT.length) + strT;
1417
1418                 arr = res.split(' ');
1419                 res  = '\t{ ' + arr[3];
1420                 res += ', ' + arr[2];
1421                 res += ', ' + arr[1];
1422                 strT = '' + (e + n);
1423                 res += ', ' + '     '.substring(strT.length) + strT;
1424                 res += ', ' + (a.length <= cbitsExact ? 0 : 1);
1425                 strT = '' + n;
1426                 res += ' ), // 10^' + strT;
1427                 if (fRoundUp)
1428                     res +=  '     '.substring(strT.length) + '(rounded up)';
1429                 print(res);
1430
1431                 return res;
1432             }
1433
1434             function PrintInv(a, n)
1435             {
1436                 var exp;
1437                 var cdig;
1438                 var len = a.length;
1439                 var div = '1';
1440                 var res = '';
1441
1442                 for (exp = 0; div.length <= len; exp++)
1443                     div += '0';
1444
1445                 for (cdig = 0; cdig < cbits; cdig++)
1446                     {
1447                     if (div.length > len || div.length == len && div >= a)
1448                         {
1449                         res += '1';
1450                         div = Sub(div, a);
1451                         }
1452                     else
1453                         res += '0';
1454                     div += '0';
1455                     }
1456
1457                 arrNeg[arrNeg.length] = PrintHex(res, -exp + 1, -n);
1458             }
1459
1460             function Sub(a, b)
1461             {
1462                 var ad, bd;
1463                 var i;
1464                 var dig;
1465                 var cch;
1466                 var lena = a.length;
1467                 var lenb = b.length;
1468                 var c = false;
1469                 var res = '';
1470
1471                 for (i = 1; i <= lena; i++)
1472                     {
1473                     ad = (a.charAt(lena - i) == '1');
1474                     bd = (b.charAt(lenb - i) == '1');
1475                     if (ad == bd)
1476                         dig = c;
1477                     else
1478                         {
1479                         dig = !c;
1480                         c = bd;
1481                         }
1482                     if (dig)
1483                         {
1484                         cch = i;
1485                         res = '1' + res;
1486                         }
1487                     else
1488                         res = '0' + res;
1489                     }
1490                 return res.substring(res.length - cch, res.length);
1491             }
1492 #endif
1493         }
1494
1495         /**
1496         * Implementation of very large variable-precision non-negative integers.
1497         *
1498         * Hungarian: bi
1499         *
1500         */
1501         private class BigInteger : IComparable {
1502             // Make this big enough that we rarely have to reallocate.
1503             private const int InitCapacity = 30;
1504
1505             private int capacity;
1506             private int length;
1507             private uint[] digits;
1508
1509             public BigInteger() {
1510                 capacity = InitCapacity;
1511                 length = 0;
1512                 digits = new uint[InitCapacity];
1513                 AssertValid();
1514             }
1515
1516             public int Length {
1517                 get { return length; }
1518             }
1519
1520             public uint this[int idx] {
1521                 get {
1522                     AssertValid();
1523                     Debug.Assert(0 <= idx && idx < length);
1524                     return digits[idx];
1525                 }
1526             }
1527
1528             [Conditional("DEBUG")]
1529             private void AssertValidNoVal() {
1530                 Debug.Assert(capacity >= InitCapacity);
1531                 Debug.Assert(length >= 0 && length <= capacity);
1532             }
1533
1534             [Conditional("DEBUG")]
1535             private void AssertValid() {
1536                 AssertValidNoVal();
1537                 Debug.Assert(0 == length || 0 != digits[length - 1]);
1538             }
1539
1540             private void Ensure(int cu) {
1541                 AssertValidNoVal();
1542
1543                 if (cu <= capacity) {
1544                     return;
1545                 }
1546
1547                 cu += cu;
1548                 uint[] newDigits = new uint[cu];
1549                 digits.CopyTo(newDigits, 0);
1550                 digits = newDigits;
1551                 capacity = cu;
1552
1553                 AssertValidNoVal();
1554             }
1555
1556             /*  ----------------------------------------------------------------------------
1557                 InitFromRgu()
1558
1559                 Initialize this big integer from an array of uint values.
1560             */
1561             public void InitFromRgu(uint[] rgu, int cu) {
1562                 AssertValid();
1563                 Debug.Assert(cu >= 0);
1564
1565                 Ensure(cu);
1566                 length = cu;
1567                 //
1568                 for (int i = 0; i < cu; i++) {
1569                     digits[i] = rgu[i];
1570                 }
1571                 AssertValid();
1572             }
1573
1574             /*  ----------------------------------------------------------------------------
1575                 InitFromRgu()
1576
1577                 Initialize this big integer from 0, 1, or 2 uint values.
1578             */
1579             public void InitFromDigits(uint u0, uint u1, int cu) {
1580                 AssertValid();
1581                 Debug.Assert(2 <= capacity);
1582
1583                 length = cu;
1584                 digits[0] = u0;
1585                 digits[1] = u1;
1586                 AssertValid();
1587             }
1588
1589             /*  ----------------------------------------------------------------------------
1590                 InitFromBigint()
1591
1592                 Initialize this big integer from another BigInteger object.
1593             */
1594             public void InitFromBigint(BigInteger biSrc) {
1595                 AssertValid();
1596                 biSrc.AssertValid();
1597                 Debug.Assert((object)this != (object)biSrc);
1598
1599                 InitFromRgu(biSrc.digits, biSrc.length);
1600             }
1601
1602         #if !NOPARSE || DEBUG
1603             /*  ----------------------------------------------------------------------------
1604                 InitFromFloatingDecimal()
1605
1606                 Initialize this big integer from a FloatingDecimal object.
1607             */
1608             public void InitFromFloatingDecimal(FloatingDecimal dec) {
1609                 AssertValid();
1610                 Debug.Assert(dec.MantissaSize >= 0);
1611
1612                 uint uAdd, uMul;
1613                 int cu = (dec.MantissaSize + 8) / 9;
1614                 int mantissaSize = dec.MantissaSize;
1615
1616                 Ensure(cu);
1617                 length = 0;
1618
1619                 uAdd = 0;
1620                 uMul = 1;
1621                 for (int ib = 0; ib < mantissaSize; ib++) {
1622                     Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
1623                     if (1000000000 == uMul) {
1624                         MulAdd(uMul, uAdd);
1625                         uMul = 1;
1626                         uAdd = 0;
1627                     }
1628                     uMul *= 10;
1629                     uAdd = uAdd * 10 + dec[ib];
1630                 }
1631                 Debug.Assert(1 < uMul);
1632                 MulAdd(uMul, uAdd);
1633
1634                 AssertValid();
1635             }
1636         #endif
1637
1638             public void MulAdd(uint uMul, uint uAdd) {
1639                 AssertValid();
1640                 Debug.Assert(0 != uMul);
1641
1642                 for (int i = 0; i < length; i++) {
1643                     uint d, uT;
1644                     d = MulU(digits[i], uMul, out uT);
1645                     if (0 != uAdd) {
1646                         uT += AddU(ref d, uAdd);
1647                     }
1648                     digits[i] = d;
1649                     uAdd = uT;
1650                 }
1651                 if (0 != uAdd) {
1652                     Ensure(length + 1);
1653                     digits[length++] = uAdd;
1654                 }
1655                 AssertValid();
1656             }
1657
1658             public void MulPow5(int c5) {
1659                 AssertValid();
1660                 Debug.Assert(c5 >= 0);
1661
1662                 const uint C5to13 = 1220703125;
1663                 int cu = (c5 + 12) / 13;
1664
1665                 if (0 == length || 0 == c5) {
1666                     return;
1667                 }
1668
1669                 Ensure(length + cu);
1670
1671                 for ( ; c5 >= 13; c5 -= 13) {
1672                     MulAdd(C5to13, 0);
1673                 }
1674
1675                 if (c5 > 0) {
1676                     uint uT;
1677                     for (uT = 5; --c5 > 0; ) {
1678                         uT *= 5;
1679                     }
1680                     MulAdd(uT, 0);
1681                 }
1682                 AssertValid();
1683             }
1684
1685             public void ShiftLeft(int cbit) {
1686                 AssertValid();
1687                 Debug.Assert(cbit >= 0);
1688
1689                 int idx, cu;
1690                 uint uExtra;
1691
1692                 if (0 == cbit || 0 == length) {
1693                     return;
1694                 }
1695
1696                 cu = cbit >> 5;
1697                 cbit &= 0x001F;
1698
1699                 if (cbit > 0) {
1700                     idx = length - 1;
1701                     uExtra = digits[idx] >> (32 - cbit);
1702
1703                     for ( ; ; idx--) {
1704                         digits[idx] <<= cbit;
1705                         if (0 == idx) {
1706                             break;
1707                         }
1708                         digits[idx] |= digits[idx - 1] >> (32 - cbit);
1709                     }
1710                 } else {
1711                     uExtra = 0;
1712                 }
1713
1714                 if (cu > 0 || 0 != uExtra) {
1715                     // Make sure there's enough room.
1716                     idx = length + (0 != uExtra ? 1 : 0) + cu;
1717                     Ensure(idx);
1718
1719                     if (cu > 0) {
1720                         // Shift the digits.
1721                         //
1722                         for (int i = length; 0 != i--; ) {
1723                             digits[cu + i] = digits[i];
1724                         }
1725                         //
1726                         for (int i = 0; i < cu; i++) {
1727                             digits[i] = 0;
1728                         }
1729                         length += cu;
1730                     }
1731
1732                     // Throw on the extra one.
1733                     if (0 != uExtra) {
1734                         digits[length++] = uExtra;
1735                     }
1736                 }
1737                 AssertValid();
1738             }
1739
1740             public void ShiftUsRight(int cu) {
1741                 AssertValid();
1742                 Debug.Assert(cu >= 0);
1743
1744                 if (cu >= length) {
1745                     length = 0;
1746                 } else if (cu > 0) {
1747                     //
1748                     for (int i = 0; i < length - cu; i++) {
1749                         digits[i] = digits[cu + i];
1750                     }
1751                     length -= cu;
1752                 }
1753                 AssertValid();
1754             }
1755
1756             public void ShiftRight(int cbit) {
1757                 AssertValid();
1758                 Debug.Assert(cbit >= 0);
1759
1760                 int idx;
1761                 int cu = cbit >> 5;
1762                 cbit &= 0x001F;
1763
1764                 if (cu > 0) {
1765                     ShiftUsRight(cu);
1766                 }
1767
1768                 if (0 == cbit || 0 == length) {
1769                     AssertValid();
1770                     return;
1771                 }
1772
1773                 for (idx = 0; ; ) {
1774                     digits[idx] >>= cbit;
1775                     if (++idx >= length) {
1776                         // Last one.
1777                         if (0 == digits[idx - 1]) {
1778                             length--;
1779                         }
1780                         break;
1781                     }
1782                     digits[idx - 1] |= digits[idx] << (32 - cbit);
1783                 }
1784                 AssertValid();
1785             }
1786
1787             public int CompareTo(object obj) {
1788                 BigInteger bi = (BigInteger)obj;
1789                 AssertValid();
1790                 bi.AssertValid();
1791
1792                 if (length > bi.length) {
1793                     return 1;
1794                 } else if (length < bi.length) {
1795                     return -1;
1796                 } else if (0 == length) {
1797                     return 0;
1798                 }
1799
1800                 int idx;
1801
1802                 for (idx = length - 1; digits[idx] == bi.digits[idx]; idx--) {
1803                     if (0 == idx) {
1804                         return 0;
1805                     }
1806                 }
1807                 Debug.Assert(idx >= 0 && idx < length);
1808                 Debug.Assert(digits[idx] != bi.digits[idx]);
1809
1810                 return (digits[idx] > bi.digits[idx]) ? 1 : -1;
1811             }
1812
1813             public void Add(BigInteger bi) {
1814                 AssertValid();
1815                 bi.AssertValid();
1816                 Debug.Assert((object)this != (object)bi);
1817
1818                 int idx, cuMax, cuMin;
1819                 uint wCarry;
1820
1821                 if ((cuMax = length) < (cuMin = bi.length)) {
1822                     cuMax = bi.length;
1823                     cuMin = length;
1824                     Ensure(cuMax + 1);
1825                 }
1826
1827                 wCarry = 0;
1828                 for (idx = 0; idx < cuMin; idx++) {
1829                     if (0 != wCarry) {
1830                         wCarry = AddU(ref digits[idx], wCarry);
1831                     }
1832                     wCarry += AddU(ref digits[idx], bi.digits[idx]);
1833                 }
1834
1835                 if (length < bi.length) {
1836                     for ( ; idx < cuMax; idx++) {
1837                         digits[idx] = bi.digits[idx];
1838                         if (0 != wCarry) {
1839                             wCarry = AddU(ref digits[idx], wCarry);
1840                         }
1841                     }
1842                     length = cuMax;
1843                 } else {
1844                     for ( ; 0 != wCarry && idx < cuMax; idx++) {
1845                         wCarry = AddU(ref digits[idx], wCarry);
1846                     }
1847                 }
1848
1849                 if (0 != wCarry) {
1850                     Ensure(length + 1);
1851                     digits[length++] = wCarry;
1852                 }
1853                 AssertValid();
1854             }
1855
1856             public void Subtract(BigInteger bi) {
1857                 AssertValid();
1858                 bi.AssertValid();
1859                 Debug.Assert((object)this != (object)bi);
1860
1861                 int idx;
1862                 uint wCarry, uT;
1863
1864                 if (length < bi.length) {
1865                     goto LNegative;
1866                 }
1867
1868                 wCarry = 1;
1869                 for (idx = 0; idx < bi.length; idx++) {
1870                     Debug.Assert(0 == wCarry || 1 == wCarry);
1871                     uT = bi.digits[idx];
1872
1873                     // NOTE: We should really do:
1874                     //    wCarry = AddU(ref digits[idx], wCarry);
1875                     //    wCarry += AddU(ref digits[idx], ~uT);
1876                     // The only case where this is different than
1877                     //    wCarry = AddU(ref digits[idx], ~uT + wCarry);
1878                     // is when 0 == uT and 1 == wCarry, in which case we don't
1879                     // need to add anything and wCarry should still be 1, so we can
1880                     // just skip the operations.
1881
1882                     if (0 != uT || 0 == wCarry) {
1883                         wCarry = AddU(ref digits[idx], ~uT + wCarry);
1884                     }
1885                 }
1886                 while (0 == wCarry && idx < length) {
1887                     wCarry = AddU(ref digits[idx], 0xFFFFFFFF);
1888                 }
1889
1890                 if (0 != wCarry) {
1891                     if (idx == length) {
1892                         // Trim off zeros.
1893                         while (--idx >= 0 && 0 == digits[idx]) {
1894                         }
1895                         length = idx + 1;
1896                     }
1897                     AssertValid();
1898                     return;
1899                 }
1900
1901             LNegative:
1902                 // bi was bigger than this.
1903                 Debug.Assert(false, "Who's subtracting to negative?");
1904                 length = 0;
1905                 AssertValid();
1906             }
1907
1908             public uint DivRem(BigInteger bi) {
1909                 AssertValid();
1910                 bi.AssertValid();
1911                 Debug.Assert((object)this != (object)bi);
1912
1913                 int idx, cu;
1914                 uint uQuo, wCarry;
1915                 int wT;
1916                 uint uT, uHi, uLo;
1917
1918                 cu = bi.length;
1919                 Debug.Assert(length <= cu);
1920                 if (length < cu) {
1921                     return 0;
1922                 }
1923
1924                 // Get a lower bound on the quotient.
1925                 uQuo = (uint)(digits[cu - 1] / (bi.digits[cu - 1] + 1));
1926                 Debug.Assert(uQuo >= 0 && uQuo <= 9);
1927
1928                 // Handle 0 and 1 as special cases.
1929                 switch (uQuo) {
1930                 case 0:
1931                     break;
1932                 case 1:
1933                     Subtract(bi);
1934                     break;
1935                 default:
1936                     uHi = 0;
1937                     wCarry = 1;
1938                     for (idx = 0; idx < cu; idx++) {
1939                         Debug.Assert(0 == wCarry || 1 == wCarry);
1940
1941                         // Compute the product.
1942                         uLo = MulU(uQuo, bi.digits[idx], out uT);
1943                         uHi = uT + AddU(ref uLo, uHi);
1944
1945                         // Subtract the product. See note in BigInteger.Subtract.
1946                         if (0 != uLo || 0 == wCarry) {
1947                             wCarry = AddU(ref digits[idx], ~uLo + wCarry);
1948                         }
1949                     }
1950                     Debug.Assert(1 == wCarry);
1951                     Debug.Assert(idx == cu);
1952
1953                     // Trim off zeros.
1954                     while (--idx >= 0 && 0 == digits[idx]) {
1955                     }
1956                     length = idx + 1;
1957                     break;
1958                 }
1959
1960                 if (uQuo < 9 && (wT = CompareTo(bi)) >= 0) {
1961                     // Quotient was off too small (by one).
1962                     uQuo++;
1963                     if (0 == wT) {
1964                         length = 0;
1965                     } else {
1966                         Subtract(bi);
1967                     }
1968                 }
1969                 Debug.Assert(CompareTo(bi) < 0);
1970
1971                 return uQuo;
1972             }
1973
1974         #if NEVER
1975             public static explicit operator double(BigInteger bi) {
1976                 uint uHi, uLo;
1977                 uint u1, u2, u3;
1978                 int idx;
1979                 int cbit;
1980                 uint dblHi, dblLo;
1981
1982                 switch (bi.length) {
1983                 case 0:
1984                     return 0;
1985                 case 1:
1986                     return bi.digits[0];
1987                 case 2:
1988                     return (ulong)bi.digits[1] << 32 | bi.digits[0];
1989                 }
1990
1991                 Debug.Assert(3 <= bi.length);
1992                 if (bi.length > 32) {
1993                     // Result is infinite.
1994                     return BitConverter.Int64BitsToDouble(0x7FF00000L << 32);
1995                 }
1996
1997                 u1 = bi.digits[bi.length - 1];
1998                 u2 = bi.digits[bi.length - 2];
1999                 u3 = bi.digits[bi.length - 3];
2000                 Debug.Assert(0 != u1);
2001                 cbit = 31 - CbitZeroLeft(u1);
2002
2003                 if (0 == cbit) {
2004                     uHi = u2;
2005                     uLo = u3;
2006                 } else {
2007                     uHi = (u1 << (32 - cbit)) | (u2 >> cbit);
2008                     // Or 1 if there are any remaining nonzero bits in u3, so we take
2009                     // them into account when rounding.
2010                     uLo = (u2 << (32 - cbit)) | (u3 >> cbit) | NotZero(u3 << (32 - cbit));
2011                 }
2012
2013                 // Set the mantissa bits.
2014                 dblHi = uHi >> 12;
2015                 dblLo = (uHi << 20) | (uLo >> 12);
2016
2017                 // Set the exponent field.
2018                 dblHi |= (uint)(0x03FF + cbit + (bi.length - 1) * 0x0020) << 20;
2019
2020                 // Do IEEE rounding.
2021                 if (0 != (uLo & 0x0800)) {
2022                     if (0 != (uLo & 0x07FF) || 0 != (dblLo & 1)) {
2023                         if (0 == ++dblLo) {
2024                             ++dblHi;
2025                         }
2026                     } else {
2027                         // If there are any non-zero bits in digits from 0 to length - 4,
2028                         // round up.
2029                         for (idx = bi.length - 4; idx >= 0; idx--) {
2030                             if (0 != bi.digits[idx]) {
2031                                 if (0 == ++dblLo) {
2032                                     ++dblHi;
2033                                 }
2034                                 break;
2035                             }
2036                         }
2037                     }
2038                 }
2039                 return BitConverter.Int64BitsToDouble((long)dblHi << 32 | dblLo);
2040             }
2041         #endif
2042         };
2043
2044         /**
2045         * Floating point number represented in base-10.
2046         */
2047         private class FloatingDecimal {
2048             public  const int MaxDigits =   50;
2049             private const int MaxExp10  =  310;  // Upper bound on base 10 exponent
2050             private const int MinExp10  = -325;  // Lower bound on base 10 exponent
2051
2052             private int     exponent;             // Base-10 scaling factor (0 means decimal point immediately precedes first digit)
2053             private int     sign;                 // Sign is -1 or 1, depending on sign of number
2054             private int     mantissaSize;         // Size of mantissa
2055             private byte[]  mantissa = new byte[MaxDigits];    // Array of base-10 digits
2056
2057             public int      Exponent { get { return exponent; } set { exponent = value; } }
2058             public int      Sign     { get { return sign;     } set { sign = value;     } }
2059             public byte[]   Mantissa { get { return mantissa; } }
2060
2061             public int MantissaSize {
2062                 get {
2063                     return mantissaSize;
2064                 }
2065                 set {
2066                     Debug.Assert(value <= MaxDigits);
2067                     mantissaSize = value;
2068                 }
2069             }
2070
2071             public byte this[int ib] {
2072                 get {
2073                     Debug.Assert(0 <= ib && ib < mantissaSize);
2074                     return mantissa[ib];
2075                 }
2076             }
2077
2078             public FloatingDecimal() {
2079                 exponent = 0;
2080                 sign = 1;
2081                 mantissaSize = 0;
2082             }
2083
2084             public FloatingDecimal(double dbl) {
2085                 InitFromDouble(dbl);
2086
2087             #if DEBUG
2088                 if (0 != mantissaSize) {
2089                     Debug.Assert(dbl == (double)this);
2090
2091                     FloatingDecimal decAfter = new FloatingDecimal();
2092                     decAfter.InitFromDouble(Succ(dbl));
2093                     // Assert(memcmp(this, &decAfter, sizeof(*this) - MaxDigits + mantissaSize));
2094                     Debug.Assert(!this.Equals(decAfter));
2095
2096                     FloatingDecimal decBefore = new FloatingDecimal();
2097                     decBefore.InitFromDouble(Pred(dbl));
2098                     // Assert(memcmp(this, &decBefore, sizeof(*this) - MaxDigits + mantissaSize));
2099                     Debug.Assert(!this.Equals(decBefore));
2100                 }
2101             #endif
2102             }
2103
2104         #if DEBUG
2105             private bool Equals(FloatingDecimal other) {
2106                 if (exponent != other.exponent || sign != other.sign || mantissaSize != other.mantissaSize) {
2107                     return false;
2108                 }
2109                 for (int idx = 0; idx < mantissaSize; idx++) {
2110                     if (mantissa[idx] != other.mantissa[idx]) {
2111                         return false;
2112                     }
2113                 }
2114                 return true;
2115             }
2116         #endif
2117
2118         #if NEVER
2119             /*  ----------------------------------------------------------------------------
2120                 RoundTo()
2121
2122                 Rounds off the BCD representation of a number to a specified number of digits.
2123                 This may result in the exponent being incremented (e.g. if digits were 999).
2124             */
2125             public void RoundTo(int sizeMantissa) {
2126                 if (sizeMantissa >= mantissaSize) {
2127                     // No change required
2128                     return;
2129                 }
2130
2131                 if (sizeMantissa >= 0) {
2132                     bool fRoundUp = mantissa[sizeMantissa] >= 5;
2133                     mantissaSize = sizeMantissa;
2134
2135                     // Round up if necessary and trim trailing zeros
2136                     for (int idx = mantissaSize - 1; idx >= 0; idx--) {
2137                         if (fRoundUp) {
2138                             if (++(mantissa[idx]) <= 9) {
2139                                 // Trailing digit is non-zero, so break
2140                                 fRoundUp = false;
2141                                 break;
2142                             }
2143                         } else if (mantissa[idx] > 0) {
2144                             // Trailing digit is non-zero, so break
2145                             break;
2146                         }
2147
2148                         // Trim trailing zeros
2149                         mantissaSize--;
2150                     }
2151
2152                     if (fRoundUp) {
2153                         // Number consisted only of 9's
2154                         Debug.Assert(0 == mantissaSize);
2155                         mantissa[0] = 1;
2156                         mantissaSize = 1;
2157                         exponent++;
2158                     }
2159                 } else {
2160                     // Number was rounded past any significant digits (e.g. 0.001 rounded to 1 fractional place), so round to 0.0
2161                     mantissaSize = 0;
2162                 }
2163
2164                 if (0 == mantissaSize) {
2165                     // 0.0
2166                     sign = 1;
2167                     exponent = 0;
2168                 }
2169             }
2170         #endif
2171
2172         #if !NOPARSE || DEBUG
2173             /*  ----------------------------------------------------------------------------
2174                 explicit operator double()
2175
2176                 Returns the double value of this floating decimal.
2177             */
2178             public static explicit operator double(FloatingDecimal dec) {
2179                 BigNumber num, numHi, numLo;
2180                 uint   ul;
2181                 int    scale;
2182                 double dbl, dblLowPrec, dblLo;
2183                 int    mantissaSize = dec.mantissaSize;
2184
2185                 // Verify that there are no leading or trailing zeros in the mantissa
2186                 Debug.Assert(0 != mantissaSize && 0 != dec[0] && 0 != dec[mantissaSize - 1]);
2187
2188                 // See if we can just use IEEE double arithmetic.
2189                 scale = dec.exponent - mantissaSize;
2190                 if (mantissaSize <= 15 && scale >= -22 && dec.exponent <= 37) {
2191                     // These calculations are all exact since mantissaSize <= 15.
2192                     if (mantissaSize <= 9) {
2193                         // Can use the ALU to perform fast integer arithmetic
2194                         ul = 0;
2195                         for (int ib = 0; ib < mantissaSize; ib++) {
2196                             Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
2197                             ul = ul * 10 + dec[ib];
2198                         }
2199                         dbl = ul;
2200                     } else {
2201                         // Use floating point arithmetic
2202                         dbl = 0.0;
2203                         for (int ib = 0; ib < mantissaSize; ib++) {
2204                             Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
2205                             dbl = dbl * 10.0 + dec[ib];
2206                         }
2207                     }
2208
2209                     // This is the only (potential) rounding operation and we assume
2210                     // the compiler does the correct IEEE rounding.
2211                     if (scale > 0) {
2212                         // Need to scale upwards by powers of 10
2213                         if (scale > 22) {
2214                             // This one is exact. We're using the fact that mantissaSize < 15
2215                             // to handle exponents bigger than 22.
2216                             dbl *= C10toN[scale - 22];
2217                             dbl *= C10toN[22];
2218                         } else {
2219                             dbl *= C10toN[scale];
2220                         }
2221                     } else if (scale < 0) {
2222                         // Scale number by negative power of 10
2223                         dbl /= C10toN[-scale];
2224                     }
2225
2226             #if DEBUG
2227                     // In the debug version, execute the high precision code also and
2228                     // verify that the results are the same.
2229                     dblLowPrec = dbl;
2230                 }
2231                 else {
2232                     dblLowPrec = Double.NaN;
2233                 }
2234             #else
2235                     goto LDone;
2236                 }
2237             #endif
2238
2239                 if (dec.exponent >= MaxExp10) {
2240                     // Overflow to infinity.
2241                     dbl = Double.PositiveInfinity;
2242                     goto LDone;
2243                 }
2244
2245                 if (dec.exponent <= MinExp10) {
2246                     // Underflow to 0.
2247                     dbl = 0.0;
2248                     goto LDone;
2249                 }
2250
2251                 // Convert to a big number.
2252                 num = new BigNumber(dec);
2253
2254                 // If there is no error in the big number, just convert it to a double.
2255                 if (0 == num.Error) {
2256                     dbl = (double)num;
2257             #if DEBUG
2258                     dblLo = dec.AdjustDbl(dbl);
2259                     Debug.Assert(dbl == dblLo);
2260             #endif
2261                     goto LDone;
2262                 }
2263
2264                 // The big number has error in it, so see if the error matters.
2265                 // Get the upper bound and lower bound. If they convert to the same
2266                 // double we're done.
2267                 numHi = num;
2268                 numHi.MakeUpperBound();
2269                 numLo = num;
2270                 numLo.MakeLowerBound();
2271
2272                 dbl = (double)numHi;
2273                 dblLo = (double)numLo;
2274                 if (dbl == dblLo) {
2275             #if DEBUG
2276                     Debug.Assert(dbl == (double)num);
2277                     dblLo = dec.AdjustDbl(dbl);
2278                     Debug.Assert(dbl == dblLo || Double.IsNaN(dblLo));
2279             #endif
2280                     goto LDone;
2281                 }
2282
2283                 // Need to use big integer arithmetic. There's too much error in
2284                 // our result and it's close to a boundary value. This is rare,
2285                 // but does happen. Eg,
2286                 // x = 1.2345678901234568347913049445e+200;
2287                 //
2288                 dbl = dec.AdjustDbl((double)num);
2289
2290             LDone:
2291                 // This assert was removed because it would fire on VERY rare occasions. Not
2292                 // repro on all machines and very hard to repro even on machines that could repro it.
2293                 // The numbers (dblLowPrec and dbl) were different in their two least sig bits only
2294                 // which is _probably_ within expected errror. I did not take the time to fully
2295                 // investigate whether this really does meet the ECMA spec...
2296                 //
2297                 Debug.Assert(Double.IsNaN(dblLowPrec) || dblLowPrec == dbl);
2298                 return dec.sign < 0 ? -dbl : dbl;
2299             }
2300
2301             /*  ----------------------------------------------------------------------------
2302                 AdjustDbl()
2303
2304                 The double contains a binary value, M * 2^n, which is off by at most 1
2305                 in the least significant bit; this class' members represent a decimal
2306                 value, D * 10^e.
2307
2308                 The general scheme is to find an integer N (the smaller the better) such
2309                 that N * M * 2^n and N * D * 10^e are both integers. We then compare
2310                 N * M * 2^n to N * D * 10^e (at full precision). If the binary value is
2311                 greater, we adjust it to be exactly half way to the next value that can
2312                 come from a double. We then compare again to decided whether to bump the
2313                 double up to the next value. Similary if the binary value is smaller,
2314                 we adjust it to be exactly half way to the previous representable value
2315                 and recompare.
2316             */
2317             private double AdjustDbl(double dbl) {
2318                 BigInteger biDec = new BigInteger();
2319                 BigInteger biDbl = new BigInteger();
2320                 int c2Dec, c2Dbl;
2321                 int c5Dec, c5Dbl;
2322                 uint wAddHi, uT;
2323                 int wT, iT;
2324                 int lwExp, wExp2;
2325                 //uint *rgu = stackalloc uint[2];
2326                 uint rgu0, rgu1;
2327                 int cu;
2328
2329                 biDec.InitFromFloatingDecimal(this);
2330                 lwExp = exponent - mantissaSize;
2331
2332                 // lwExp is a base 10 exponent.
2333                 if (lwExp >= 0) {
2334                     c5Dec = c2Dec = lwExp;
2335                     c5Dbl = c2Dbl = 0;
2336                 } else {
2337                     c5Dec = c2Dec = 0;
2338                     c5Dbl = c2Dbl = -lwExp;
2339                 }
2340
2341                 rgu1 = DblHi(dbl);
2342                 rgu0 = DblLo(dbl);
2343                 wExp2 = (int)(rgu1 >> 20) & 0x07FF;
2344                 rgu1 &= 0x000FFFFF;
2345                 wAddHi = 1;
2346                 if (0 != wExp2) {
2347                     // Normal, so add implicit bit.
2348                     if (0 == rgu1 && 0 == rgu0 && 1 != wExp2) {
2349                         // Power of 2 (and not adjacent to the first denormal), so the
2350                         // adjacent low value is closer than the high value.
2351                         wAddHi = 2;
2352                         rgu1 = 0x00200000;
2353                         wExp2--;
2354                     } else {
2355                         rgu1 |= 0x00100000;
2356                     }
2357                     wExp2 -= 1076;
2358                 } else {
2359                     wExp2 = -1075;
2360                 }
2361
2362                 // Shift left by 1 bit : the adjustment values need the next lower bit.
2363                 rgu1 = (rgu1 << 1) | (rgu0 >> 31);
2364                 rgu0 <<= 1;
2365
2366                 // We must determine how many words of significant digits this requires.
2367                 if (0 == rgu0 && 0 == rgu1) {
2368                     cu = 0;
2369                 } else if (0 == rgu1) {
2370                     cu = 1;
2371                 } else {
2372                     cu = 2;
2373                 }
2374
2375                 biDbl.InitFromDigits(rgu0, rgu1, cu);
2376
2377                 if (wExp2 >= 0) {
2378                     c2Dbl += wExp2;
2379                 } else {
2380                     c2Dec += -wExp2;
2381                 }
2382
2383                 // Eliminate common powers of 2.
2384                 if (c2Dbl > c2Dec) {
2385                     c2Dbl -= c2Dec;
2386                     c2Dec = 0;
2387
2388                     // See if biDec has some powers of 2 that we can get rid of.
2389                     for (iT = 0; c2Dbl >= 32 && 0 == biDec[iT]; iT++) {
2390                         c2Dbl -= 32;
2391                     }
2392                     if (iT > 0) {
2393                         biDec.ShiftUsRight(iT);
2394                     }
2395                     Debug.Assert(c2Dbl < 32 || 0 != biDec[0]);
2396                     uT = biDec[0];
2397                     for (iT = 0; iT < c2Dbl && 0 == (uT & (1L << iT)); iT++) {
2398                     }
2399                     if (iT > 0) {
2400                         c2Dbl -= iT;
2401                         biDec.ShiftRight(iT);
2402                     }
2403                 } else {
2404                     c2Dec -= c2Dbl;
2405                     c2Dbl = 0;
2406                 }
2407
2408                 // There are no common powers of 2 or common powers of 5.
2409                 Debug.Assert(0 == c2Dbl || 0 == c2Dec);
2410                 Debug.Assert(0 == c5Dbl || 0 == c5Dec);
2411
2412                 // Fold in the powers of 5.
2413                 if (c5Dbl > 0) {
2414                     biDbl.MulPow5(c5Dbl);
2415                 } else if (c5Dec > 0) {
2416                     biDec.MulPow5(c5Dec);
2417                 }
2418
2419                 // Fold in the powers of 2.
2420                 if (c2Dbl > 0) {
2421                     biDbl.ShiftLeft(c2Dbl);
2422                 } else if (c2Dec > 0) {
2423                     biDec.ShiftLeft(c2Dec);
2424                 }
2425
2426                 // Now determine whether biDbl is above or below biDec.
2427                 wT = biDbl.CompareTo(biDec);
2428
2429                 if (0 == wT) {
2430                     return dbl;
2431                 } else if (wT > 0) {
2432                     // biDbl is greater. Recompute with the dbl minus half the distance
2433                     // to the next smaller double.
2434                     if (0 == AddU(ref rgu0, 0xFFFFFFFF)) {
2435                         AddU(ref rgu1, 0xFFFFFFFF);
2436                     }
2437                     biDbl.InitFromDigits(rgu0, rgu1, 1 + (0 != rgu1 ? 1 : 0));
2438                     if (c5Dbl > 0) {
2439                         biDbl.MulPow5(c5Dbl);
2440                     }
2441                     if (c2Dbl > 0) {
2442                         biDbl.ShiftLeft(c2Dbl);
2443                     }
2444
2445                     wT = biDbl.CompareTo(biDec);
2446                     if (wT > 0 || 0 == wT && 0 != (DblLo(dbl) & 1)) {
2447                         // Return the next lower value.
2448                         dbl = BitConverter.Int64BitsToDouble(BitConverter.DoubleToInt64Bits(dbl) - 1);
2449                     }
2450                 } else {
2451                     // biDbl is smaller. Recompute with the dbl plus half the distance
2452                     // to the next larger double.
2453                     if (0 != AddU(ref rgu0, wAddHi)) {
2454                         AddU(ref rgu1, 1);
2455                     }
2456                     biDbl.InitFromDigits(rgu0, rgu1, 1 + (0 != rgu1 ? 1 : 0));
2457                     if (c5Dbl > 0) {
2458                         biDbl.MulPow5(c5Dbl);
2459                     }
2460                     if (c2Dbl > 0) {
2461                         biDbl.ShiftLeft(c2Dbl);
2462                     }
2463
2464                     wT = biDbl.CompareTo(biDec);
2465                     if (wT < 0 || 0 == wT && 0 != (DblLo(dbl) & 1)) {
2466                         // Return the next higher value.
2467                         dbl = BitConverter.Int64BitsToDouble(BitConverter.DoubleToInt64Bits(dbl) + 1);
2468                     }
2469                 }
2470                 return dbl;
2471             }
2472         #endif
2473
2474             private void InitFromDouble(double dbl) {
2475                 if (0.0 == dbl || IsSpecial(dbl)) {
2476                     exponent = 0;
2477                     sign = 1;
2478                     mantissaSize = 0;
2479                 } else {
2480                     // Handle the sign.
2481                     if (dbl < 0) {
2482                         sign = -1;
2483                         dbl = -dbl;
2484                     } else {
2485                         sign = 1;
2486                     }
2487
2488                     if (!BigNumber.DblToRgbFast(dbl, mantissa, out exponent, out mantissaSize)) {
2489                         BigNumber.DblToRgbPrecise(dbl, mantissa, out exponent, out mantissaSize);
2490                     }
2491                 }
2492             }
2493         };
2494
2495         /*  ----------------------------------------------------------------------------
2496             IntToString()
2497
2498             Converts an integer to a string according to XPath rules.
2499         */
2500         private static unsafe string IntToString(int val) {
2501             // The maximum number of characters needed to represent any int value is 11
2502             const int BufSize = 12;
2503             char *pBuf = stackalloc char[BufSize];
2504             char *pch = pBuf += BufSize;
2505             uint u = (uint)(val < 0 ? -val : val);
2506
2507             while (u >= 10) {
2508                 // Fast division by 10
2509                 uint quot = (uint)((0x66666667L * u) >> 32) >> 2;
2510                 *(--pch) = (char)((u - quot * 10) + '0');
2511                 u = quot;
2512             }
2513
2514             *(--pch) = (char)(u + '0');
2515
2516             if (val < 0) {
2517                 *(--pch) = '-';
2518             }
2519
2520             return new string(pch, 0, (int)(pBuf - pch));
2521         }
2522
2523         /*  ----------------------------------------------------------------------------
2524             DoubleToString()
2525
2526             Converts a floating point number to a string according to XPath rules.
2527         */
2528         public static string DoubleToString(double dbl) {
2529             Debug.Assert(('0' & 0xF) == 0, "We use (char)(d |'0') to convert digit to char");
2530             int   maxSize, sizeInt, sizeFract, cntDigits, ib;
2531             int   iVal;
2532
2533             if (IsInteger(dbl, out iVal)) {
2534                 return IntToString(iVal);
2535             }
2536
2537             // Handle NaN and infinity
2538             if (IsSpecial(dbl)) {
2539                 if (Double.IsNaN(dbl)) {
2540                     return "NaN";
2541                 }
2542
2543                 Debug.Assert(Double.IsInfinity(dbl));
2544                 return dbl < 0 ? "-Infinity" : "Infinity";
2545             }
2546
2547             // Get decimal digits
2548             FloatingDecimal dec = new FloatingDecimal(dbl);
2549             Debug.Assert(0 != dec.MantissaSize);
2550
2551             // If exponent is negative, size of fraction increases
2552             sizeFract = dec.MantissaSize - dec.Exponent;
2553
2554             if (sizeFract > 0) {
2555                 // Decimal consists of a fraction part + possible integer part
2556                 sizeInt = (dec.Exponent > 0) ? dec.Exponent : 0;
2557             } else {
2558                 // Decimal consists of just an integer part
2559                 sizeInt = dec.Exponent;
2560                 sizeFract = 0;
2561             }
2562
2563             // Sign + integer + fraction + decimal point + leading zero + terminating null
2564             maxSize = sizeInt + sizeFract + 4;
2565
2566             unsafe {
2567                 // Allocate output memory
2568                 char *pBuf = stackalloc char[maxSize];
2569                 char *pch = pBuf;
2570
2571                 if (dec.Sign < 0) {
2572                     *pch++ = '-';
2573                 }
2574
2575                 cntDigits = dec.MantissaSize;
2576                 ib = 0;
2577
2578                 if (0 != sizeInt) {
2579                     do {
2580                         if (0 != cntDigits) {
2581                             // Still mantissa digits left to consume
2582                             Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
2583                             *pch++ = (char)(dec[ib++] | '0');
2584                             cntDigits--;
2585                         } else {
2586                             // Add trailing zeros
2587                             *pch++ = '0';
2588                         }
2589                     } while (0 != --sizeInt);
2590                 } else {
2591                     *pch++ = '0';
2592                 }
2593
2594                 if (0 != sizeFract) {
2595                     Debug.Assert(0 != cntDigits);
2596                     Debug.Assert(sizeFract == cntDigits || sizeFract > cntDigits && pch[-1] == '0');
2597                     *pch++ = '.';
2598
2599                     while (sizeFract > cntDigits) {
2600                         // Add leading zeros
2601                         *pch++ = '0';
2602                         sizeFract--;
2603                     }
2604
2605                     Debug.Assert(sizeFract == cntDigits);
2606                     while (0 != cntDigits) {
2607                         // Output remaining mantissa digits
2608                         Debug.Assert(dec[ib] >= 0 && dec[ib] <= 9);
2609                         *pch++ = (char)(dec[ib++] | '0');
2610                         cntDigits--;
2611                     }
2612                 }
2613
2614                 Debug.Assert(0 == sizeInt && 0 == cntDigits);
2615                 return new string(pBuf, 0, (int)(pch - pBuf));
2616             }
2617         }
2618
2619         private static bool IsAsciiDigit(char ch) {
2620             return (uint)(ch - '0') <= 9;
2621         }
2622
2623         private static bool IsWhitespace(char ch) {
2624             return ch == '\x20' || ch == '\x9' || ch == '\xA' || ch == '\xD';
2625         }
2626
2627         private static unsafe char *SkipWhitespace(char *pch) {
2628             while (IsWhitespace(*pch)) {
2629                 pch++;
2630             }
2631             return pch;
2632         }
2633
2634         /*  ----------------------------------------------------------------------------
2635             StringToDouble()
2636
2637             Converts a string to a floating point number according to XPath rules.
2638             NaN is returned if the entire string is not a valid number.
2639
2640             This code was stolen from MSXML6 DecimalFormat::parse(). The implementation
2641             depends on the fact that the String objects are always zero-terminated.
2642         */
2643         public static unsafe double StringToDouble(string s) {
2644             Debug.Assert(('0' & 0xF) == 0, "We use (ch & 0xF) to convert char to digit");
2645             // For the mantissa digits. After leaving the state machine, pchFirstDig
2646             // points to the first digit and pch points just past the last digit.
2647             // numDig is the number of digits. pch - pchFirstDig may be numDig + 1
2648             // (if there is a decimal point).
2649
2650             fixed (char *pchStart = s) {
2651                 int numDig = 0;
2652                 char *pch = pchStart;
2653                 char *pchFirstDig = null;
2654                 char ch;
2655
2656                 int sign = 1;       // sign of the mantissa
2657                 int expAdj = 0;     // exponent adjustment
2658
2659                 // Enter the state machine
2660
2661                 // Initialization
2662             LRestart:
2663                 ch = *pch++;
2664                 if (IsAsciiDigit(ch)) {
2665                     goto LGetLeft;
2666                 }
2667
2668                 switch (ch) {
2669                 case '-':
2670                     if (sign < 0) {
2671                         break;
2672                     }
2673                     sign = -1;
2674                     goto LRestart;
2675                 case '.':
2676                     if (IsAsciiDigit(*pch)) {
2677                         goto LGetRight;
2678                     }
2679                     break;
2680                 default:
2681                     // MSXML has a bug, we should not allow whitespace after a minus sign
2682                     if (IsWhitespace(ch) && sign > 0) {
2683                         pch = SkipWhitespace(pch);
2684                         goto LRestart;
2685                     }
2686                     break;
2687                 }
2688
2689                 // Nothing digested - set the result to NaN and exit.
2690                 return Double.NaN;
2691
2692             LGetLeft:
2693                 // Get digits to the left of the decimal point
2694                 if (ch == '0') {
2695                     do {
2696                         // Trim leading zeros
2697                         ch = *pch++;
2698                     } while (ch == '0');
2699
2700                     if (!IsAsciiDigit(ch)) {
2701                         goto LSkipNonZeroDigits;
2702                     }
2703                 }
2704
2705                 Debug.Assert(IsAsciiDigit(ch));
2706                 pchFirstDig = pch - 1;
2707                 do {
2708                     ch = *pch++;
2709                 } while (IsAsciiDigit(ch));
2710                 numDig = (int)(pch - pchFirstDig) - 1;
2711
2712             LSkipNonZeroDigits:
2713                 if (ch != '.') {
2714                     goto LEnd;
2715                 }
2716
2717             LGetRight:
2718                 Debug.Assert(ch == '.');
2719                 ch = *pch++;
2720                 if (pchFirstDig == null) {
2721                     // Count fractional leading zeros (e.g. six zeros in '0.0000005')
2722                     while (ch == '0') {
2723                         expAdj--;
2724                         ch = *pch++;
2725                     }
2726                     pchFirstDig = pch - 1;
2727                 }
2728
2729                 while (IsAsciiDigit(ch)) {
2730                     expAdj--;
2731                     numDig++;
2732                     ch = *pch++;
2733                 }
2734
2735             LEnd:
2736                 pch--;
2737                 char *pchEnd = pchStart + s.Length;
2738                 Debug.Assert(*pchEnd == '\0', "String objects must be zero-terminated");
2739
2740                 if (pch < pchEnd && SkipWhitespace(pch) < pchEnd) {
2741                     // If we're not at the end of the string, this is not a valid number
2742                     return Double.NaN;
2743                 }
2744
2745                 if (numDig == 0) {
2746                     return 0.0;
2747                 }
2748
2749                 Debug.Assert(pchFirstDig != null);
2750
2751                 if (expAdj == 0) {
2752                     // Assert(StrRChrW(pchFirstDig, &pchFirstDig[numDig], '.') == null);
2753
2754                     // Detect special case where number is an integer
2755                     if (numDig <= 9) {
2756                         Debug.Assert(pchFirstDig != pch);
2757                         int iNum = *pchFirstDig & 0xF; // - '0'
2758                         while (--numDig != 0) {
2759                             pchFirstDig++;
2760                             Debug.Assert(IsAsciiDigit(*pchFirstDig));
2761                             iNum = iNum * 10 + (*pchFirstDig & 0xF); // - '0'
2762                         }
2763                         return (double)(sign < 0 ? -iNum : iNum);
2764                     }
2765                 } else {
2766                     // The number has a fractional part
2767                     Debug.Assert(expAdj < 0);
2768                     // Assert(StrRChrW(pchStart, pch, '.') != null);
2769                 }
2770
2771                 // Limit to 50 digits (double is only precise up to 17 or so digits)
2772                 if (numDig > FloatingDecimal.MaxDigits) {
2773                     pch -= (numDig - FloatingDecimal.MaxDigits);
2774                     expAdj += (numDig - FloatingDecimal.MaxDigits);
2775                     numDig = FloatingDecimal.MaxDigits;
2776                 }
2777
2778                 // Remove trailing zero's from mantissa
2779                 Debug.Assert(IsAsciiDigit(*pchFirstDig) && *pchFirstDig != '0');
2780                 while (true) {
2781                     if (*--pch == '0') {
2782                         numDig--;
2783                         expAdj++;
2784                     } else if (*pch != '.') {
2785                         Debug.Assert(IsAsciiDigit(*pch) && *pch != '0');
2786                         pch++;
2787                         break;
2788                     }
2789                 }
2790                 Debug.Assert(pch - pchFirstDig == numDig || pch - pchFirstDig == numDig + 1);
2791
2792                 {
2793                     // Construct a floating decimal from the array of digits
2794                     Debug.Assert(numDig > 0 && numDig <= FloatingDecimal.MaxDigits);
2795
2796                     FloatingDecimal dec = new FloatingDecimal();
2797                     dec.Exponent = expAdj + numDig;
2798                     dec.Sign = sign;
2799                     dec.MantissaSize = numDig;
2800
2801                     fixed (byte *pin = dec.Mantissa) {
2802                         byte *mantissa = pin;
2803                         while (pchFirstDig < pch) {
2804                             if (*pchFirstDig != '.') {
2805                                 *mantissa = (byte)(*pchFirstDig & 0xF); // - '0'
2806                                 mantissa++;
2807                             }
2808                             pchFirstDig++;
2809                         }
2810                     }
2811
2812                     return (double)dec;
2813                 }
2814             }
2815         }
2816     }
2817 }