mcs/class/dlr/Runtime/Microsoft.Dynamic/Utils/MathUtils.cs

   1 /* ****************************************************************************
   2  *
   3  * Copyright (c) Microsoft Corporation.
   4  *
   5  * This source code is subject to terms and conditions of the Apache License, Version 2.0. A
   6  * copy of the license can be found in the License.html file at the root of this distribution. If
   7  * you cannot locate the  Apache License, Version 2.0, please send an email to
   8  * dlr@microsoft.com. By using this source code in any fashion, you are agreeing to be bound
   9  * by the terms of the Apache License, Version 2.0.
  10  *
  11  * You must not remove this notice, or any other, from this software.
  12  *
  13  *
  14  * ***************************************************************************/
  15
  16 #if FEATURE_NUMERICS
  17 using BigInt = System.Numerics.BigInteger;
  18 using Complex = System.Numerics.Complex;
  19 #endif
  20
  21 using System;
  22 using System.Text;
  23 using System.Collections.Generic;
  24 using Microsoft.Scripting.Math;
  25 using Microsoft.Scripting.Runtime;
  26
  27 namespace Microsoft.Scripting.Utils {
  28     using Math = System.Math;
  29
  30     public static class MathUtils {
  31         /// <summary>
  32         /// Calculates the quotient of two 32-bit signed integers rounded towards negative infinity.
  33         /// </summary>
  34         /// <param name="x">Dividend.</param>
  35         /// <param name="y">Divisor.</param>
  36         /// <returns>The quotient of the specified numbers rounded towards negative infinity, or <code>(int)Floor((double)x/(double)y)</code>.</returns>
  37         /// <exception cref="DivideByZeroException"><paramref name="y"/> is 0.</exception>
  38         /// <remarks>The caller must check for overflow (x = Int32.MinValue, y = -1)</remarks>
  39         public static int FloorDivideUnchecked(int x, int y) {
  40             int q = x / y;
  41
  42             if (x >= 0) {
  43                 if (y > 0) {
  44                     return q;
  45                 } else if (x % y == 0) {
  46                     return q;
  47                 } else {
  48                     return q - 1;
  49                 }
  50             } else {
  51                 if (y > 0) {
  52                     if (x % y == 0) {
  53                         return q;
  54                     } else {
  55                         return q - 1;
  56                     }
  57                 } else {
  58                     return q;
  59                 }
  60             }
  61         }
  62
  63         /// <summary>
  64         /// Calculates the quotient of two 32-bit signed integers rounded towards negative infinity.
  65         /// </summary>
  66         /// <param name="x">Dividend.</param>
  67         /// <param name="y">Divisor.</param>
  68         /// <returns>The quotient of the specified numbers rounded towards negative infinity, or <code>(int)Floor((double)x/(double)y)</code>.</returns>
  69         /// <exception cref="DivideByZeroException"><paramref name="y"/> is 0.</exception>
  70         /// <remarks>The caller must check for overflow (x = Int64.MinValue, y = -1)</remarks>
  71         public static long FloorDivideUnchecked(long x, long y) {
  72             long q = x / y;
  73
  74             if (x >= 0) {
  75                 if (y > 0) {
  76                     return q;
  77                 } else if (x % y == 0) {
  78                     return q;
  79                 } else {
  80                     return q - 1;
  81                 }
  82             } else {
  83                 if (y > 0) {
  84                     if (x % y == 0) {
  85                         return q;
  86                     } else {
  87                         return q - 1;
  88                     }
  89                 } else {
  90                     return q;
  91                 }
  92             }
  93         }
  94
  95         /// <summary>
  96         /// Calculates the remainder of floor division of two 32-bit signed integers.
  97         /// </summary>
  98         /// <param name="x">Dividend.</param>
  99         /// <param name="y">Divisor.</param>
 100         /// <returns>The remainder of of floor division of the specified numbers, or <code>x - (int)Floor((double)x/(double)y) * y</code>.</returns>
 101         /// <exception cref="DivideByZeroException"><paramref name="y"/> is 0.</exception>
 102         public static int FloorRemainder(int x, int y) {
 103             if (y == -1) return 0;
 104             int r = x % y;
 105
 106             if (x >= 0) {
 107                 if (y > 0) {
 108                     return r;
 109                 } else if (r == 0) {
 110                     return 0;
 111                 } else {
 112                     return r + y;
 113                 }
 114             } else {
 115                 if (y > 0) {
 116                     if (r == 0) {
 117                         return 0;
 118                     } else {
 119                         return r + y;
 120                     }
 121                 } else {
 122                     return r;
 123                 }
 124             }
 125         }
 126
 127         /// <summary>
 128         /// Calculates the remainder of floor division of two 32-bit signed integers.
 129         /// </summary>
 130         /// <param name="x">Dividend.</param>
 131         /// <param name="y">Divisor.</param>
 132         /// <returns>The remainder of of floor division of the specified numbers, or <code>x - (int)Floor((double)x/(double)y) * y</code>.</returns>
 133         /// <exception cref="DivideByZeroException"><paramref name="y"/> is 0.</exception>
 134         public static long FloorRemainder(long x, long y) {
 135             if (y == -1) return 0;
 136             long r = x % y;
 137
 138             if (x >= 0) {
 139                 if (y > 0) {
 140                     return r;
 141                 } else if (r == 0) {
 142                     return 0;
 143                 } else {
 144                     return r + y;
 145                 }
 146             } else {
 147                 if (y > 0) {
 148                     if (r == 0) {
 149                         return 0;
 150                     } else {
 151                         return r + y;
 152                     }
 153                 } else {
 154                     return r;
 155                 }
 156             }
 157         }
 158
 159         /// <summary>
 160         /// Behaves like Math.Round(value, MidpointRounding.AwayFromZero)
 161         /// Needed because CoreCLR doesn't support this particular overload of Math.Round
 162         /// </summary>
 163         public static double RoundAwayFromZero(double value) {
 164 #if !SILVERLIGHT && !WP75
 165             return Math.Round(value, MidpointRounding.AwayFromZero);
 166 #else
 167             if (value < 0) {
 168                 return -RoundAwayFromZero(-value);
 169             }
 170
 171             // we can assume positive value
 172             double result = Math.Floor(value);
 173             if (value - result >= 0.5) {
 174                 result += 1.0;
 175             }
 176             return result;
 177 #endif
 178         }
 179
 180         private static readonly double[] _RoundPowersOfTens = new double[] { 1E0, 1E1, 1E2, 1E3, 1E4, 1E5, 1E6, 1E7, 1E8, 1E9, 1E10, 1E11, 1E12, 1E13, 1E14, 1E15 };
 181
 182         private static double GetPowerOf10(int precision) {
 183             return (precision < 16) ? _RoundPowersOfTens[precision] : Math.Pow(10, precision);
 184         }
 185
 186         /// <summary>
 187         /// Behaves like Math.Round(value, precision, MidpointRounding.AwayFromZero)
 188         /// However, it works correctly on negative precisions and cases where precision is
 189         /// outside of the [-15, 15] range.
 190         ///
 191         /// (This function is also needed because CoreCLR lacks this overload.)
 192         /// </summary>
 193         public static double RoundAwayFromZero(double value, int precision) {
 194             if (double.IsInfinity(value) || double.IsNaN(value)) {
 195                 return value;
 196             }
 197
 198             if (precision >= 0) {
 199                 if (precision > 308) {
 200                     return value;
 201                 }
 202
 203                 double num = GetPowerOf10(precision);
 204                 return RoundAwayFromZero(value * num) / num;
 205             } else if (precision >= -308) {
 206                 // Note: this code path could be merged with the precision >= 0 path,
 207                 // (by extending the cache to negative powers of 10)
 208                 // but the results seem to be more precise if we do it this way
 209                 double num = GetPowerOf10(-precision);
 210                 return RoundAwayFromZero(value / num) * num;
 211             } else {
 212                 // Preserve the sign of the input, including +/-0.0
 213                 return value < 0.0 || 1.0 / value < 0.0 ? -0.0 : 0.0;
 214             }
 215         }
 216
 217         public static bool IsNegativeZero(double self) {
 218 #if SILVERLIGHT // BitConverter.DoubleToInt64Bits
 219             if ( self != 0.0 ) {
 220               return false;
 221             }
 222             byte[] bits = BitConverter.GetBytes(self);
 223             return (bits[7] == 0x80 && bits[6] == 0x00 && bits[5] == 0x00 && bits[4] == 0x00
 224                 && bits[3] == 0x00 && bits[2] == 0x00 && bits[1] == 0x00 && bits[0] == 0x00);
 225 #else
 226             return (self == 0.0 && 1.0 / self < 0);
 227 #endif
 228         }
 229
 230         #region Special Functions
 231
 232         public static double Erf(double v0) {
 233             // Calculate the error function using the approximation method outlined in
 234             // W. J. Cody's "Rational Chebyshev Approximations for the Error Function"
 235
 236             if (v0 >= 10.0) {
 237                 return 1.0;
 238             } else if (v0 <= -10.0) {
 239                 return -1.0;
 240             }
 241
 242             if (v0 > 0.47 || v0 < -0.47) {
 243                 return 1.0 - ErfComplement(v0);
 244             }
 245
 246             double sq = v0 * v0;
 247             double numer = EvalPolynomial(sq, ErfNumerCoeffs);
 248             double denom = EvalPolynomial(sq, ErfDenomCoeffs);
 249
 250             return v0 * numer / denom;
 251         }
 252
 253         public static double ErfComplement(double v0) {
 254             if (v0 >= 30.0) {
 255                 return 0.0;
 256             } else if (v0 <= -10.0) {
 257                 return 2.0;
 258             }
 259
 260             double a = Math.Abs(v0);
 261             if (a < 0.47) {
 262                 return 1.0 - Erf(v0);
 263             }
 264
 265             // Different approximations are required for different ranges of v0
 266             double res;
 267             if (a <= 4.0) {
 268                 // Use the approximation method outlined in W. J. Cody's "Rational Chebyshev
 269                 // Approximations for the Error Function"
 270                 double numer = EvalPolynomial(a, ErfcNumerCoeffs);
 271                 double denom = EvalPolynomial(a, ErfcDenomCoeffs);
 272
 273                 res = Math.Exp(-a * a) * numer / denom;
 274             } else {
 275                 // Use the approximation method introduced by C. Tellambura and A. Annamalai
 276                 // in "Efficient Computation of erfc(x) for Large Arguments"
 277                 const double h = 0.5;
 278                 const double hSquared = 0.25;
 279                 const int nTerms = 10;
 280                 double sq = a * a;
 281                 res = 0.0;
 282                 for (int i = nTerms; i > 0; i--) {
 283                     double term = i * i * hSquared;
 284                     res += Math.Exp(-term) / (term + sq);
 285                 }
 286
 287                 res = h * a * Math.Exp(-sq) / Math.PI * (res * 2 + 1.0 / sq);
 288             }
 289
 290             if (v0 < 0.0) {
 291                 res = 2.0 - res;
 292             }
 293             return res;
 294         }
 295
 296         public static double Gamma(double v0) {
 297             // Calculate the Gamma function using the Lanczos approximation
 298
 299             if (double.IsNegativeInfinity(v0)) {
 300                 return double.NaN;
 301             }
 302             double a = Math.Abs(v0);
 303
 304             // Special-case integers
 305             if (a % 1.0 == 0.0) {
 306                 // Gamma is undefined on non-positive integers
 307                 if (v0 <= 0.0) {
 308                     return double.NaN;
 309                 }
 310
 311                 // factorial(v0 - 1)
 312                 if (a <= 25.0) {
 313                     if (a <= 2.0) {
 314                         return 1.0;
 315                     }
 316                     a -= 1.0;
 317                     v0 -= 1.0;
 318                     while (--v0 > 1.0) {
 319                         a *= v0;
 320                     }
 321                     return a;
 322                 }
 323             }
 324
 325             // lim(Gamma(v0)) = 1.0 / v0 as v0 approaches 0.0
 326             if (a < 1e-50) {
 327                 return 1.0 / v0;
 328             }
 329
 330             double res;
 331             if (v0 < -150.0) {
 332                 // If Gamma(1 - v0) could overflow for large v0, use the duplication formula to
 333                 // compute Gamma(1 - v0):
 334                 //     Gamma(x) * Gamma(x + 0,5) = sqrt(pi) * 2**(1 - 2x) * Gamma(2x)
 335                 // ==> Gamma(1 - x) = Gamma((1-x)/2) * Gamma((2-x)/2) / (2**x * sqrt(pi))
 336                 // Then apply the reflection formula:
 337                 //     Gamma(x) = pi / sin(pi * x) / Gamma(1 - x)
 338                 double halfV0 = v0 / 2.0;
 339                 res = Math.Pow(Math.PI, 1.5) / SinPi(v0);
 340                 res *= Math.Pow(2.0, v0);
 341                 res /= PositiveGamma(0.5 - halfV0);
 342                 res /= PositiveGamma(1.0 - halfV0);
 343             } else if (v0 < 0.001) {
 344                 // For values less than or close to zero, just use the reflection formula
 345                 res = Math.PI / SinPi(v0);
 346                 double v1 = 1.0 - v0;
 347                 if (v0 == 1.0 - v1) {
 348                     res /= PositiveGamma(v1);
 349                 } else {
 350                     // Computing v1 has resulted in a loss of precision. To avoid this, use the
 351                     // recurrence relation Gamma(x + 1) = x * Gamma(x).
 352                     res /= -v0 * PositiveGamma(-v0);
 353                 }
 354             } else {
 355                 res = PositiveGamma(v0);
 356             }
 357
 358             return res;
 359         }
 360
 361         public static double LogGamma(double v0) {
 362             // Calculate the log of the Gamma function using the Lanczos approximation
 363
 364             if (double.IsInfinity(v0)) {
 365                 return double.PositiveInfinity;
 366             }
 367             double a = Math.Abs(v0);
 368
 369             // Gamma is undefined on non-positive integers
 370             if (v0 <= 0.0 && a % 1.0 == 0.0) {
 371                 return double.NaN;
 372             }
 373
 374             // lim(LGamma(v0)) = -log|v0| as v0 approaches 0.0
 375             if (a < 1e-50) {
 376                 return -Math.Log(a);
 377             }
 378
 379             double res;
 380             if (v0 < 0.0) {
 381                 // For negative values, use the reflection formula:
 382                 //     Gamma(x) = pi / sin(pi * x) / Gamma(1 - x)
 383                 // ==> LGamma(x) = log(pi / |sin(pi * x)|) - LGamma(1 - x)
 384                 res = Math.Log(Math.PI / AbsSinPi(v0));
 385                 res -= PositiveLGamma(1.0 - v0);
 386             } else {
 387                 res = PositiveLGamma(v0);
 388             }
 389
 390             return res;
 391         }
 392
 393         public static double Hypot(double x, double y) {
 394             //
 395             // sqrt(x*x + y*y) == sqrt(x*x * (1 + (y*y)/(x*x))) ==
 396             // sqrt(x*x) * sqrt(1 + (y/x)*(y/x)) ==
 397             // abs(x) * sqrt(1 + (y/x)*(y/x))
 398             //
 399
 400             // Handle infinities
 401             if (double.IsInfinity(x) || double.IsInfinity(y)) {
 402                 return double.PositiveInfinity;
 403             }
 404
 405             //  First, get abs
 406             if (x < 0.0) x = -x;
 407             if (y < 0.0) y = -y;
 408
 409             // Obvious cases
 410             if (x == 0.0) return y;
 411             if (y == 0.0) return x;
 412
 413             // Divide smaller number by bigger number to safeguard the (y/x)*(y/x)
 414             if (x < y) {
 415                 double temp = y; y = x; x = temp;
 416             }
 417
 418             y /= x;
 419
 420             // calculate abs(x) * sqrt(1 + (y/x)*(y/x))
 421             return x * System.Math.Sqrt(1 + y * y);
 422         }
 423
 424         /// <summary>
 425         /// Evaluates a polynomial in v0 where the coefficients are ordered in increasing degree
 426         /// </summary>
 427         private static double EvalPolynomial(double v0, double[] coeffs) {
 428             double res = 0.0;
 429             for (int i = coeffs.Length - 1; i >= 0; i--) {
 430                 res = checked(res * v0 + coeffs[i]);
 431             }
 432
 433             return res;
 434         }
 435
 436         /// <summary>
 437         /// Evaluates a polynomial in v0 where the coefficients are ordered in increasing degree
 438         /// if reverse is false, and increasing degree if reverse is true.
 439         /// </summary>
 440         private static double EvalPolynomial(double v0, double[] coeffs, bool reverse) {
 441             if (!reverse) {
 442                 return EvalPolynomial(v0, coeffs);
 443             }
 444
 445             double res = 0.0;
 446             for (int i = 0; i < coeffs.Length; i++) {
 447                 res = checked(res * v0 + coeffs[i]);
 448             }
 449
 450             return res;
 451         }
 452
 453         /// <summary>
 454         /// A numerically precise version of sin(v0 * pi)
 455         /// </summary>
 456         private static double SinPi(double v0) {
 457             double res = Math.Abs(v0) % 2.0;
 458
 459             if (res < 0.25) {
 460                 res = Math.Sin(res * Math.PI);
 461             } else if (res < 0.75) {
 462                 res = Math.Cos((res - 0.5) * Math.PI);
 463             } else if (res < 1.25) {
 464                 res = -Math.Sin((res - 1.0) * Math.PI);
 465             } else if (res < 1.75) {
 466                 res = -Math.Cos((res - 1.5) * Math.PI);
 467             } else {
 468                 res = Math.Sin((res - 2.0) * Math.PI);
 469             }
 470
 471             return v0 < 0 ? -res : res;
 472         }
 473
 474         /// <summary>
 475         /// A numerically precise version of |sin(v0 * pi)|
 476         /// </summary>
 477         private static double AbsSinPi(double v0) {
 478             double res = Math.Abs(v0) % 1.0;
 479
 480             if (res < 0.25) {
 481                 res = Math.Sin(res * Math.PI);
 482             } else if (res < 0.75) {
 483                 res = Math.Cos((res - 0.5) * Math.PI);
 484             } else {
 485                 res = Math.Sin((res - 1.0) * Math.PI);
 486             }
 487
 488             return Math.Abs(res);
 489         }
 490
 491         // polynomial coefficients ordered by increasing degree
 492         private static double[] ErfNumerCoeffs = {
 493             2.4266795523053175e02, 2.1979261618294152e01,
 494             6.9963834886191355, -3.5609843701815385e-02
 495         };
 496         private static double[] ErfDenomCoeffs = {
 497             2.1505887586986120e02, 9.1164905404514901e01,
 498             1.5082797630407787e01, 1.0
 499         };
 500         private static double[] ErfcNumerCoeffs = {
 501             3.004592610201616005e02, 4.519189537118729422e02,
 502             3.393208167343436870e02, 1.529892850469404039e02,
 503             4.316222722205673530e01, 7.211758250883093659,
 504             5.641955174789739711e-01, -1.368648573827167067e-07
 505         };
 506         private static double[] ErfcDenomCoeffs = {
 507             3.004592609569832933e02, 7.909509253278980272e02,
 508             9.313540948506096211e02, 6.389802644656311665e02,
 509             2.775854447439876434e02, 7.700015293522947295e01,
 510             1.278272731962942351e01, 1.0
 511         };
 512         private static double[] GammaNumerCoeffs = {
 513             4.401213842800460895436e13, 4.159045335859320051581e13,
 514             1.801384278711799677796e13, 4.728736263475388896889e12,
 515             8.379100836284046470415e11, 1.055837072734299344907e11,
 516             9.701363618494999493386e09, 6.549143975482052641016e08,
 517             3.223832294213356530668e07, 1.128514219497091438040e06,
 518             2.666579378459858944762e04, 3.818801248632926870394e02,
 519             2.506628274631000502415
 520         };
 521         private static double[] GammaDenomCoeffs = {
 522             0.0, 39916800.0, 120543840.0, 150917976.0,
 523             105258076.0, 45995730.0, 13339535.0, 2637558.0,
 524             357423.0, 32670.0, 1925.0, 66.0, 1.0
 525         };
 526
 527         /// <summary>
 528         /// Take the quotient of the 2 polynomials forming the Lanczos approximation
 529         /// with N=13 and G=13.144565
 530         /// </summary>
 531         private static double GammaRationalFunc(double v0) {
 532             double numer = 0.0;
 533             double denom = 0.0;
 534
 535             if (v0 < 1e15) {
 536                 numer = EvalPolynomial(v0, GammaNumerCoeffs);
 537                 denom = EvalPolynomial(v0, GammaDenomCoeffs);
 538             } else {
 539                 double vRecip = 1.0 / v0;
 540                 numer = EvalPolynomial(vRecip, GammaNumerCoeffs, true);
 541                 denom = EvalPolynomial(vRecip, GammaDenomCoeffs, true);
 542             }
 543
 544             return numer / denom;
 545         }
 546
 547         /// <summary>
 548         /// Computes the Gamma function on positive values, using the Lanczos approximation.
 549         /// Lanczos parameters are N=13 and G=13.144565.
 550         /// </summary>
 551         private static double PositiveGamma(double v0) {
 552             if (v0 > 200.0) {
 553                 return Double.PositiveInfinity;
 554             }
 555
 556             double vg = v0 + 12.644565; // v0 + g - 0.5
 557             double res = GammaRationalFunc(v0);
 558             res /= Math.Exp(vg);
 559             if (v0 < 120.0) {
 560                 res *= Math.Pow(vg, v0 - 0.5);
 561             } else {
 562                 // Use a smaller exponent if we're in danger of overflowing Math.Pow
 563                 double sqrt = Math.Pow(vg, v0 / 2.0 - 0.25);
 564                 res *= sqrt;
 565                 res *= sqrt;
 566             }
 567
 568             return res;
 569         }
 570
 571         /// <summary>
 572         /// Computes the Log-Gamma function on positive values, using the Lanczos approximation.
 573         /// Lanczos parameters are N=13 and G=13.144565.
 574         /// </summary>
 575         private static double PositiveLGamma(double v0) {
 576             double vg = v0 + 12.644565; // v0 + g - 0.5
 577             double res = Math.Log(GammaRationalFunc(v0)) - vg;
 578             res += (v0 - 0.5) * Math.Log(vg);
 579
 580             return res;
 581         }
 582
 583         #endregion
 584
 585         #region BigInteger
 586
 587         // generated by scripts/radix_generator.py
 588         private static readonly uint[] maxCharsPerDigit = { 0, 0, 31, 20, 15, 13, 12, 11, 10, 10, 9, 9, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 };
 589         private static readonly uint[] groupRadixValues = { 0, 0, 2147483648, 3486784401, 1073741824, 1220703125, 2176782336, 1977326743, 1073741824, 3486784401, 1000000000, 2357947691, 429981696, 815730721, 1475789056, 2562890625, 268435456, 410338673, 612220032, 893871739, 1280000000, 1801088541, 2494357888, 3404825447, 191102976, 244140625, 308915776, 387420489, 481890304, 594823321, 729000000, 887503681, 1073741824, 1291467969, 1544804416, 1838265625, 2176782336 };
 590
 591         internal static string BigIntegerToString(uint[] d, int sign, int radix, bool lowerCase) {
 592             if (radix < 2) {
 593                 throw ExceptionUtils.MakeArgumentOutOfRangeException("radix", radix, "radix must be >= 2");
 594             }
 595             if (radix > 36) {
 596                 throw ExceptionUtils.MakeArgumentOutOfRangeException("radix", radix, "radix must be <= 36");
 597             }
 598
 599             int dl = d.Length;
 600             if (dl == 0) {
 601                 return "0";
 602             }
 603
 604             List<uint> digitGroups = new List<uint>();
 605
 606             uint groupRadix = groupRadixValues[radix];
 607             while (dl > 0) {
 608                 uint rem = div(d, ref dl, groupRadix);
 609                 digitGroups.Add(rem);
 610             }
 611
 612             StringBuilder ret = new StringBuilder();
 613             if (sign == -1) {
 614                 ret.Append("-");
 615             }
 616
 617             int digitIndex = digitGroups.Count - 1;
 618
 619             char[] tmpDigits = new char[maxCharsPerDigit[radix]];
 620
 621             AppendRadix((uint)digitGroups[digitIndex--], (uint)radix, tmpDigits, ret, false, lowerCase);
 622             while (digitIndex >= 0) {
 623                 AppendRadix((uint)digitGroups[digitIndex--], (uint)radix, tmpDigits, ret, true, lowerCase);
 624             }
 625             return ret.Length == 0 ? "0" : ret.ToString();
 626         }
 627
 628         private const int BitsPerDigit = 32;
 629
 630         private static uint div(uint[] n, ref int nl, uint d) {
 631             ulong rem = 0;
 632             int i = nl;
 633             bool seenNonZero = false;
 634             while (--i >= 0) {
 635                 rem <<= BitsPerDigit;
 636                 rem |= n[i];
 637                 uint v = (uint)(rem / d);
 638                 n[i] = v;
 639                 if (v == 0) {
 640                     if (!seenNonZero) nl--;
 641                 } else {
 642                     seenNonZero = true;
 643                 }
 644                 rem %= d;
 645             }
 646             return (uint)rem;
 647         }
 648
 649         private static void AppendRadix(uint rem, uint radix, char[] tmp, StringBuilder buf, bool leadingZeros, bool lowerCase) {
 650             string symbols = lowerCase ? "0123456789abcdefghijklmnopqrstuvwxyz" : "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
 651
 652             int digits = tmp.Length;
 653             int i = digits;
 654             while (i > 0 && (leadingZeros || rem != 0)) {
 655                 uint digit = rem % radix;
 656                 rem /= radix;
 657                 tmp[--i] = symbols[(int)digit];
 658             }
 659             if (leadingZeros) buf.Append(tmp);
 660             else buf.Append(tmp, i, digits - i);
 661         }
 662
 663         // Helper for GetRandBits
 664         private static uint GetWord(byte[] bytes, int start, int end) {
 665             uint four = 0;
 666             int bits = end - start;
 667             int shift = 0;
 668             if (bits > 32) {
 669                 bits = 32;
 670             }
 671             start /= 8;
 672             while (bits > 0) {
 673                 uint value = bytes[start];
 674                 if (bits < 8) {
 675                     value &= (1u << bits) - 1u;
 676                 }
 677                 value <<= shift;
 678                 four |= value;
 679                 bits -= 8;
 680                 shift += 8;
 681                 start++;
 682             }
 683
 684             return four;
 685         }
 686 #if !MONO_INTERPRETER
 687 #if !FEATURE_NUMERICS
 688         public static BigInteger GetRandBits(this Random generator, int bits) {
 689             ContractUtils.Requires(bits > 0);
 690
 691             // equivalent to (bits + 7) / 8 without possibility of overflow
 692             int count = bits % 8 == 0 ? bits / 8 : bits / 8 + 1;
 693
 694             // Pad the end (most significant) with zeros if we align to the byte
 695             // to ensure that we end up with a positive value
 696             byte[] bytes = new byte[bits % 8 == 0 ? count + 1 : count];
 697             generator.NextBytes(bytes);
 698             if (bits % 8 == 0) {
 699                 bytes[bytes.Length - 1] = 0;
 700             } else {
 701                 bytes[bytes.Length - 1] = (byte)(bytes[bytes.Length - 1] & ((1 << (bits % 8)) - 1));
 702             }
 703
 704             if (bits <= 32) {
 705                 return (BigInteger)GetWord(bytes, 0, bits);
 706             } else if (bits <= 64) {
 707                 ulong a = GetWord(bytes, 0, bits);
 708                 ulong b = GetWord(bytes, 32, bits);
 709                 return (BigInteger)(a | (b << 32));
 710             } else {
 711                 count = (count + 3) / 4;
 712                 uint[] data = new uint[count];
 713                 for (int i = 0; i < count; i++) {
 714                     data[i] = GetWord(bytes, i * 32, bits);
 715                 }
 716                 return new BigInteger(1, data);
 717             }
 718         }
 719
 720         public static BigInteger Random(this Random generator, BigInteger limit) {
 721             ContractUtils.Requires(limit.Sign > 0, "limit");
 722             ContractUtils.RequiresNotNull(generator, "generator");
 723
 724             // TODO: this doesn't yield a uniform distribution (small numbers will be picked more frequently):
 725             uint[] result = new uint[limit.GetWordCount() + 1];
 726             for (int i = 0; i < result.Length; i++) {
 727                 result[i] = unchecked((uint)generator.Next());
 728             }
 729             return new BigInteger(1, result) % limit;
 730         }
 731 #else
 732         public static BigInt GetRandBits(this Random generator, int bits) {
 733             ContractUtils.Requires(bits > 0);
 734
 735             // equivalent to (bits + 7) / 8 without possibility of overflow
 736             int count = bits % 8 == 0 ? bits / 8 : bits / 8 + 1;
 737
 738             // Pad the end (most significant) with zeros if we align to the byte
 739             // to ensure that we end up with a positive value
 740             byte[] bytes = new byte[bits % 8 == 0 ? count + 1 : count];
 741             generator.NextBytes(bytes);
 742             if (bits % 8 == 0) {
 743                 bytes[bytes.Length - 1] = 0;
 744             } else {
 745                 bytes[bytes.Length - 1] = (byte)(bytes[bytes.Length - 1] & ((1 << (bits % 8)) - 1));
 746             }
 747
 748             if (bits <= 32) {
 749                 return (BigInt)GetWord(bytes, 0, bits);
 750             } else if (bits <= 64) {
 751                 ulong a = GetWord(bytes, 0, bits);
 752                 ulong b = GetWord(bytes, 32, bits);
 753                 return (BigInt)(a | (b << 32));
 754             }
 755
 756             return new BigInt(bytes);
 757         }
 758
 759         public static BigInteger Random(this Random generator, BigInteger limit) {
 760             return new BigInteger(generator.Random(limit.Value));
 761         }
 762
 763         public static BigInt Random(this Random generator, BigInt limit) {
 764             ContractUtils.Requires(limit.Sign > 0, "limit");
 765             ContractUtils.RequiresNotNull(generator, "generator");
 766
 767             BigInt res = BigInt.Zero;
 768
 769             while (true) {
 770                 // if we've run out of significant digits, we can return the total
 771                 if (limit == BigInt.Zero) {
 772                     return res;
 773                 }
 774
 775                 // if we're small enough to fit in an int, do so
 776                 int iLimit;
 777                 if (limit.AsInt32(out iLimit)) {
 778                     return res + generator.Next(iLimit);
 779                 }
 780
 781                 // get the 3 or 4 uppermost bytes that fit into an int
 782                 int hiData;
 783                 byte[] data = limit.ToByteArray();
 784                 int index = data.Length;
 785                 while (data[--index] == 0) ;
 786                 if (data[index] < 0x80) {
 787                     hiData = data[index] << 24;
 788                     data[index--] = (byte)0;
 789                 } else {
 790                     hiData = 0;
 791                 }
 792                 hiData |= data[index] << 16;
 793                 data[index--] = (byte)0;
 794                 hiData |= data[index] << 8;
 795                 data[index--] = (byte)0;
 796                 hiData |= data[index];
 797                 data[index--] = (byte)0;
 798
 799                 // get a uniform random number for the uppermost portion of the bigint
 800                 byte[] randomData = new byte[index + 2];
 801                 generator.NextBytes(randomData);
 802                 randomData[index + 1] = (byte)0;
 803                 res += new BigInt(randomData);
 804                 res += (BigInt)generator.Next(hiData) << ((index + 1) * 8);
 805
 806                 // sum it with a uniform random number for the remainder of the bigint
 807                 limit = new BigInt(data);
 808             }
 809         }
 810
 811         public static bool TryToFloat64(this BigInt self, out double result) {
 812             return StringUtils.TryParseDouble(
 813                 self.ToString(),
 814                 System.Globalization.NumberStyles.Number,
 815                 System.Globalization.CultureInfo.InvariantCulture.NumberFormat,
 816                 out result
 817             );
 818         }
 819
 820         public static double ToFloat64(this BigInt self) {
 821             return double.Parse(
 822                 self.ToString(),
 823                 System.Globalization.NumberStyles.Number,
 824                 System.Globalization.CultureInfo.InvariantCulture.NumberFormat
 825             );
 826         }
 827 #endif
 828         public static bool TryToFloat64(this BigInteger self, out double result) {
 829             return StringUtils.TryParseDouble(
 830                 self.ToString(10),
 831                 System.Globalization.NumberStyles.Number,
 832                 System.Globalization.CultureInfo.InvariantCulture.NumberFormat,
 833                 out result
 834             );
 835         }
 836
 837         public static double ToFloat64(this BigInteger self) {
 838             return double.Parse(
 839                 self.ToString(10),
 840                 System.Globalization.NumberStyles.Number,
 841                 System.Globalization.CultureInfo.InvariantCulture.NumberFormat
 842             );
 843         }
 844
 845         // Like GetBitCount(Abs(x)), except 0 maps to 0
 846         public static int BitLength(BigInteger x) {
 847             if (x.IsZero()) {
 848                 return 0;
 849             }
 850
 851             return x.Abs().GetBitCount();
 852         }
 853
 854 #if FEATURE_NUMERICS
 855         public static int BitLength(BigInt x) {
 856             if (x.IsZero) {
 857                 return 0;
 858             }
 859
 860             byte[] bytes = BigInt.Abs(x).ToByteArray();
 861             int index = bytes.Length;
 862             while (bytes[--index] == 0) ;
 863
 864             return index * 8 + BitLength((int)bytes[index]);
 865         }
 866 #endif
 867 #endif
 868         // Like GetBitCount(Abs(x)), except 0 maps to 0
 869         public static int BitLength(long x) {
 870             if (x == 0) {
 871                 return 0;
 872             }
 873             if (x == Int64.MinValue) {
 874                 return 64;
 875             }
 876
 877             x = Math.Abs(x);
 878             int res = 1;
 879             if (x >= 1L << 32) {
 880                 x >>= 32;
 881                 res += 32;
 882             }
 883             if (x >= 1L << 16) {
 884                 x >>= 16;
 885                 res += 16;
 886             }
 887             if (x >= 1L << 8) {
 888                 x >>= 8;
 889                 res += 8;
 890             }
 891             if (x >= 1L << 4) {
 892                 x >>= 4;
 893                 res += 4;
 894             }
 895             if (x >= 1L << 2) {
 896                 x >>= 2;
 897                 res += 2;
 898             }
 899             if (x >= 1L << 1) {
 900                 res += 1;
 901             }
 902
 903             return res;
 904         }
 905
 906         // Like GetBitCount(Abs(x)), except 0 maps to 0
 907         [CLSCompliant(false)]
 908         public static int BitLengthUnsigned(ulong x) {
 909             if (x >= 1uL << 63) {
 910                 return 64;
 911             }
 912             return BitLength((long)x);
 913         }
 914
 915         // Like GetBitCount(Abs(x)), except 0 maps to 0
 916         public static int BitLength(int x) {
 917             if (x == 0) {
 918                 return 0;
 919             }
 920             if (x == Int32.MinValue) {
 921                 return 32;
 922             }
 923
 924             x = Math.Abs(x);
 925             int res = 1;
 926             if (x >= 1 << 16) {
 927                 x >>= 16;
 928                 res += 16;
 929             }
 930             if (x >= 1 << 8) {
 931                 x >>= 8;
 932                 res += 8;
 933             }
 934             if (x >= 1 << 4) {
 935                 x >>= 4;
 936                 res += 4;
 937             }
 938             if (x >= 1 << 2) {
 939                 x >>= 2;
 940                 res += 2;
 941             }
 942             if (x >= 1 << 1) {
 943                 res += 1;
 944             }
 945
 946             return res;
 947         }
 948
 949         // Like GetBitCount(Abs(x)), except 0 maps to 0
 950         [CLSCompliant(false)]
 951         public static int BitLengthUnsigned(uint x) {
 952             if (x >= 1u << 31) {
 953                 return 32;
 954             }
 955             return BitLength((int)x);
 956         }
 957
 958         #region Extending BigInt with BigInteger API
 959 #if FEATURE_NUMERICS
 960
 961         public static bool AsInt32(this BigInt self, out int ret) {
 962             if (self >= Int32.MinValue && self <= Int32.MaxValue) {
 963                 ret = (Int32)self;
 964                 return true;
 965             }
 966             ret = 0;
 967             return false;
 968         }
 969
 970         public static bool AsInt64(this BigInt self, out long ret) {
 971             if (self >= Int64.MinValue && self <= Int64.MaxValue) {
 972                 ret = (long)self;
 973                 return true;
 974             }
 975             ret = 0;
 976             return false;
 977         }
 978
 979         [CLSCompliant(false)]
 980         public static bool AsUInt32(this BigInt self, out uint ret) {
 981             if (self >= UInt32.MinValue && self <= UInt32.MaxValue) {
 982                 ret = (UInt32)self;
 983                 return true;
 984             }
 985             ret = 0;
 986             return false;
 987         }
 988
 989         [CLSCompliant(false)]
 990         public static bool AsUInt64(this BigInt self, out ulong ret) {
 991             if (self >= UInt64.MinValue && self <= UInt64.MaxValue) {
 992                 ret = (UInt64)self;
 993                 return true;
 994             }
 995             ret = 0;
 996             return false;
 997         }
 998
 999         public static BigInt Abs(this BigInt self) {
1000             return BigInt.Abs(self);
1001         }
1002
1003         public static bool IsZero(this BigInt self) {
1004             return self.IsZero;
1005         }
1006
1007         public static bool IsPositive(this BigInt self) {
1008             return self.Sign > 0;
1009         }
1010
1011         public static bool IsNegative(this BigInt self) {
1012             return self.Sign < 0;
1013         }
1014
1015         public static double Log(this BigInt self) {
1016             return BigInt.Log(self);
1017         }
1018
1019         public static double Log(this BigInt self, double baseValue) {
1020             return BigInt.Log(self, baseValue);
1021         }
1022
1023         public static double Log10(this BigInt self) {
1024             return BigInt.Log10(self);
1025         }
1026
1027         public static BigInt Power(this BigInt self, int exp) {
1028             return BigInt.Pow(self, exp);
1029         }
1030
1031         public static BigInt ModPow(this BigInt self, int power, BigInt mod) {
1032             return BigInt.ModPow(self, power, mod);
1033         }
1034
1035         public static BigInt ModPow(this BigInt self, BigInt power, BigInt mod) {
1036             return BigInt.ModPow(self, power, mod);
1037         }
1038
1039         public static string ToString(this BigInt self, int radix) {
1040             const string symbols = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
1041
1042             if (radix < 2) {
1043                 throw ExceptionUtils.MakeArgumentOutOfRangeException("radix", radix, "radix must be >= 2");
1044             }
1045             if (radix > 36) {
1046                 throw ExceptionUtils.MakeArgumentOutOfRangeException("radix", radix, "radix must be <= 36");
1047             }
1048
1049             bool isNegative = false;
1050             if (self < BigInt.Zero) {
1051                 self = -self;
1052                 isNegative = true;
1053             } else if (self == BigInt.Zero) {
1054                 return "0";
1055             }
1056
1057             List<char> digits = new List<char>();
1058             while (self > 0) {
1059                 digits.Add(symbols[(int)(self % radix)]);
1060                 self /= radix;
1061             }
1062
1063             StringBuilder ret = new StringBuilder();
1064             if (isNegative) {
1065                 ret.Append('-');
1066             }
1067             for (int digitIndex = digits.Count - 1; digitIndex >= 0; digitIndex--) {
1068                 ret.Append(digits[digitIndex]);
1069             }
1070             return ret.ToString();
1071         }
1072 #endif
1073         #endregion
1074
1075         #region Exposing underlying data
1076 #if FEATURE_NUMERICS
1077
1078         [CLSCompliant(false)]
1079         public static uint[] GetWords(this BigInt self) {
1080             if (self.IsZero) {
1081                 return new uint[] { 0 };
1082             }
1083
1084             int hi;
1085             byte[] bytes;
1086             GetHighestByte(self, out hi, out bytes);
1087
1088             uint[] result = new uint[(hi + 1 + 3) / 4];
1089             int i = 0;
1090             int j = 0;
1091             uint u = 0;
1092             int shift = 0;
1093             while (i < bytes.Length) {
1094                 u |= (uint)bytes[i++] << shift;
1095                 if (i % 4 == 0) {
1096                     result[j++] = u;
1097                     u = 0;
1098                 }
1099                 shift += 8;
1100             }
1101             if (u != 0) {
1102                 result[j] = u;
1103             }
1104             return result;
1105         }
1106
1107         [CLSCompliant(false)]
1108         public static uint GetWord(this BigInt self, int index) {
1109             return GetWords(self)[index];
1110         }
1111
1112         public static int GetWordCount(this BigInt self) {
1113             int index;
1114             byte[] bytes;
1115             GetHighestByte(self, out index, out bytes);
1116             return index / 4 + 1; // return (index + 1 + 3) / 4;
1117         }
1118
1119         public static int GetByteCount(this BigInt self) {
1120             int index;
1121             byte[] bytes;
1122             GetHighestByte(self, out index, out bytes);
1123             return index + 1;
1124         }
1125
1126         public static int GetBitCount(this BigInt self) {
1127             if (self.IsZero) {
1128                 return 1;
1129             }
1130             byte[] bytes = BigInt.Abs(self).ToByteArray();
1131
1132             int index = bytes.Length;
1133             while (bytes[--index] == 0) ;
1134
1135             int count = index * 8;
1136             for (int hiByte = bytes[index]; hiByte > 0; hiByte >>= 1) {
1137                 count++;
1138             }
1139             return count;
1140         }
1141
1142         private static byte GetHighestByte(BigInt self, out int index, out byte[] byteArray) {
1143             byte[] bytes = BigInt.Abs(self).ToByteArray();
1144             if (self.IsZero) {
1145                 byteArray = bytes;
1146                 index = 0;
1147                 return 1;
1148             }
1149
1150             int hi = bytes.Length;
1151             byte b;
1152             do {
1153                 b = bytes[--hi];
1154             } while (b == 0);
1155             index = hi;
1156             byteArray = bytes;
1157             return b;
1158         }
1159
1160 #endif
1161         #endregion
1162
1163         #endregion
1164
1165 #if !MONO_INTERPRETER
1166         #region Complex
1167
1168 #if !FEATURE_NUMERICS
1169         public static Complex64 MakeReal(double real) {
1170             return new Complex64(real, 0.0);
1171         }
1172
1173         public static Complex64 MakeImaginary(double imag) {
1174             return new Complex64(0.0, imag);
1175         }
1176
1177         public static Complex64 MakeComplex(double real, double imag) {
1178             return new Complex64(real, imag);
1179         }
1180
1181         public static double Imaginary(this Complex64 self) {
1182             return self.Imag;
1183         }
1184
1185         public static bool IsZero(this Complex64 self) {
1186             return self.IsZero;
1187         }
1188
1189         public static Complex64 Pow(this Complex64 self, Complex64 power) {
1190             return self.Power(power);
1191         }
1192 #else
1193         public static Complex MakeReal(double real) {
1194             return new Complex(real, 0.0);
1195         }
1196
1197         public static Complex MakeImaginary(double imag) {
1198             return new Complex(0.0, imag);
1199         }
1200
1201         public static Complex MakeComplex(double real, double imag) {
1202             return new Complex(real, imag);
1203         }
1204
1205         public static double Imaginary(this Complex self) {
1206             return self.Imaginary;
1207         }
1208
1209         public static bool IsZero(this Complex self) {
1210             return self.Equals(Complex.Zero);
1211         }
1212
1213         public static Complex Conjugate(this Complex self) {
1214             return new Complex(self.Real, -self.Imaginary);
1215         }
1216
1217         public static double Abs(this Complex self) {
1218             return Complex.Abs(self);
1219         }
1220
1221         public static Complex Pow(this Complex self, Complex power) {
1222             return Complex.Pow(self, power);
1223         }
1224 #endif
1225
1226         #endregion
1227 #endif
1228     }
1229
1230 }