mcs/class/System.Numerics/System.Numerics/BigInteger.cs

   1 //
   2 // System.Numerics.BigInteger
   3 //
   4 // Rodrigo Kumpera (rkumpera@novell.com)
   5
   6 //
   7 // Copyright (C) 2010 Novell, Inc (http://www.novell.com)
   8 //
   9 // Permission is hereby granted, free of charge, to any person obtaining
  10 // a copy of this software and associated documentation files (the
  11 // "Software"), to deal in the Software without restriction, including
  12 // without limitation the rights to use, copy, modify, merge, publish,
  13 // distribute, sublicense, and/or sell copies of the Software, and to
  14 // permit persons to whom the Software is furnished to do so, subject to
  15 // the following conditions:
  16 //
  17 // The above copyright notice and this permission notice shall be
  18 // included in all copies or substantial portions of the Software.
  19 //
  20 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  21 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  22 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  23 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  24 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  25 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  26 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  27 //
  28 // The Multiply and Div code comes the DLR (as hosted in http://ironpython.codeplex.com),
  29 // which has the following License:
  30 //
  31 /* ****************************************************************************
  32  *
  33  * Copyright (c) Microsoft Corporation.
  34  *
  35  * This source code is subject to terms and conditions of the Microsoft Public License. A
  36  * copy of the license can be found in the License.html file at the root of this distribution. If
  37  * you cannot locate the  Microsoft Public License, please send an email to
  38  * dlr@microsoft.com. By using this source code in any fashion, you are agreeing to be bound
  39  * by the terms of the Microsoft Public License.
  40  *
  41  * You must not remove this notice, or any other, from this software.
  42  *
  43  *
  44  * ***************************************************************************/
  45
  46 using System;\r
  47 using System.Collections.Generic;\r
  48 using System.Diagnostics;\r
  49 using System.Diagnostics.CodeAnalysis;\r
  50 using System.Globalization;\r
  51 using System.Text;\r
  52
  53 /*
  54 Optimization
  55         Have proper popcount function for IsPowerOfTwo
  56         Use unsafe ops to avoid bounds check
  57         CoreAdd could avoid some resizes by checking for equal sized array that top overflow
  58         For bitwise operators, hoist the conditionals out of their main loop
  59         Optimize BitScanBackward
  60         Use a carry variable to make shift opts do half the number of array ops.
  61         Schoolbook multiply is O(n^2), use Karatsuba /Toom-3 for large numbers
  62 */\r
  63 namespace System.Numerics {\r
  64         public struct BigInteger : IComparable, IComparable<BigInteger>, IEquatable<BigInteger>
  65         {
  66                 //LSB on [0]
  67                 readonly uint[] data;
  68                 readonly short sign;
  69
  70                 static readonly uint[] ZERO = new uint [1];
  71                 static readonly uint[] ONE = new uint [1] { 1 };
  72
  73                 BigInteger (short sign, uint[] data)
  74                 {
  75                         this.sign = sign;
  76                         this.data = data;
  77                 }
  78
  79                 public BigInteger (int value)
  80                 {
  81                         if (value == 0) {
  82                                 sign = 0;
  83                                 data = ZERO;
  84                         } else if (value > 0) {
  85                                 sign = 1;
  86                                 data = new uint[] { (uint) value };
  87                         } else {
  88                                 sign = -1;
  89                                 data = new uint[1] { (uint)-value };
  90                         }
  91                 }
  92
  93                 [CLSCompliantAttribute (false)]
  94                 public BigInteger (uint value)
  95                 {
  96                         if (value == 0) {
  97                                 sign = 0;
  98                                 data = ZERO;
  99                         } else {
 100                                 sign = 1;
 101                                 data = new uint [1] { value };
 102                         }
 103                 }
 104
 105                 public BigInteger (long value)
 106                 {
 107                         if (value == 0) {
 108                                 sign = 0;
 109                                 data = ZERO;
 110                         } else if (value > 0) {
 111                                 sign = 1;
 112                                 uint low = (uint)value;
 113                                 uint high = (uint)(value >> 32);
 114
 115                                 data = new uint [high != 0 ? 2 : 1];
 116                                 data [0] = low;
 117                                 if (high != 0)
 118                                         data [1] = high;
 119                         } else {
 120                                 sign = -1;
 121                                 value = -value;
 122                                 uint low = (uint)value;
 123                                 uint high = (uint)((ulong)value >> 32);
 124
 125                                 data = new uint [high != 0 ? 2 : 1];
 126                                 data [0] = low;
 127                                 if (high != 0)
 128                                         data [1] = high;
 129                         }
 130                 }
 131
 132                 [CLSCompliantAttribute (false)]
 133                 public BigInteger (ulong value)
 134                 {
 135                         if (value == 0) {
 136                                 sign = 0;
 137                                 data = ZERO;
 138                         } else {
 139                                 sign = 1;
 140                                 uint low = (uint)value;
 141                                 uint high = (uint)(value >> 32);
 142
 143                                 data = new uint [high != 0 ? 2 : 1];
 144                                 data [0] = low;
 145                                 if (high != 0)
 146                                         data [1] = high;
 147                         }
 148                 }
 149
 150                 [CLSCompliantAttribute (false)]
 151                 public BigInteger (byte[] value)
 152                 {
 153                         if (value == null)
 154                                 throw new ArgumentNullException ("value");
 155
 156                         int len = value.Length;
 157
 158                         if (len == 0 || (len == 1 && value [0] == 0)) {
 159                                 sign = 0;
 160                                 data = ZERO;
 161                                 return;
 162                         }
 163
 164                         if ((value [len - 1] & 0x80) != 0)
 165                                 sign = -1;
 166                         else
 167                                 sign = 1;
 168
 169                         if (sign == 1) {
 170                                 while (value [len - 1] == 0)
 171                                         --len;
 172
 173                                 int full_words, size;
 174                                 full_words = size = len / 4;
 175                                 if ((len & 0x3) != 0)
 176                                         ++size;
 177
 178                                 data = new uint [size];
 179                                 int j = 0;
 180                                 for (int i = 0; i < full_words; ++i) {
 181                                         data [i] =      (uint)value [j++] |
 182                                                                 (uint)(value [j++] << 8) |
 183                                                                 (uint)(value [j++] << 16) |
 184                                                                 (uint)(value [j++] << 24);
 185                                 }
 186                                 size = len & 0x3;
 187                                 if (size > 0) {
 188                                         int idx = data.Length - 1;
 189                                         for (int i = 0; i < size; ++i)
 190                                                 data [idx] |= (uint)(value [j++] << (i * 8));
 191                                 }
 192                         } else {
 193                                 int full_words, size;
 194                                 full_words = size = len / 4;
 195                                 if ((len & 0x3) != 0)
 196                                         ++size;
 197
 198                                 data = new uint [size];
 199
 200                                 uint word, borrow = 1;
 201                                 ulong sub = 0;
 202                                 int j = 0;
 203
 204                                 for (int i = 0; i < full_words; ++i) {
 205                                         word =  (uint)value [j++] |
 206                                                         (uint)(value [j++] << 8) |
 207                                                         (uint)(value [j++] << 16) |
 208                                                         (uint)(value [j++] << 24);
 209
 210                                         sub = (ulong)word - borrow;
 211                                         word = (uint)sub;
 212                                         borrow = (uint)(sub >> 32) & 0x1u;
 213                                         data [i] = ~word;
 214                                 }
 215                                 size = len & 0x3;
 216
 217                                 if (size > 0) {
 218                                         word = 0;
 219                                         uint store_mask = 0;
 220                                         for (int i = 0; i < size; ++i) {
 221                                                 word |= (uint)(value [j++] << (i * 8));
 222                                                 store_mask = (store_mask << 8) | 0xFF;
 223                                         }
 224
 225                                         sub = word - borrow;
 226                                         word = (uint)sub;
 227                                         borrow = (uint)(sub >> 32) & 0x1u;
 228
 229                                         data [data.Length - 1] = ~word & store_mask;
 230                                 }
 231                                 if (borrow != 0) //FIXME I believe this can't happen, can someone write a test for it?
 232                                         throw new Exception ("non zero final carry");
 233                         }
 234
 235                 }
 236
 237                 public bool IsEven {
 238                         get { return (data [0] & 0x1) == 0; }
 239                 }
 240
 241                 public bool IsOne {
 242                         get { return sign == 1 && data.Length == 1 && data [0] == 1; }
 243                 }
 244
 245
 246                 //Gem from Hacker's Delight
 247                 //Returns the number of bits set in @x
 248                 static int PopulationCount (uint x)
 249                 {
 250                         x = x - ((x >> 1) & 0x55555555);
 251                         x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
 252                         x = (x + (x >> 4)) & 0x0F0F0F0F;
 253                         x = x + (x >> 8);
 254                         x = x + (x >> 16);
 255                         return (int)(x & 0x0000003F);
 256                 }
 257
 258                 public bool IsPowerOfTwo {
 259                         get {
 260                                 bool foundBit = false;
 261                                 if (sign != 1)
 262                                         return false;
 263                                 //This function is pop count == 1 for positive numbers
 264                                 for (int i = 0; i < data.Length; ++i) {
 265                                         int p = PopulationCount (data [i]);
 266                                         if (p > 0) {
 267                                                 if (p > 1 || foundBit)
 268                                                         return false;
 269                                                 foundBit = true;
 270                                         }
 271                                 }
 272                                 return foundBit;
 273                         }
 274                 }
 275
 276                 public bool IsZero {
 277                         get { return sign == 0; }
 278                 }
 279
 280                 public int Sign {
 281                         get { return sign; }
 282                 }
 283
 284                 public static BigInteger MinusOne {
 285                         get { return new BigInteger (-1, ONE); }
 286                 }
 287
 288                 public static BigInteger One {
 289                         get { return new BigInteger (1, ONE); }
 290                 }
 291
 292                 public static BigInteger Zero {
 293                         get { return new BigInteger (0, ZERO); }
 294                 }
 295
 296                 public static explicit operator int (BigInteger value)
 297                 {
 298                         if (value.data.Length > 1)
 299                                 throw new OverflowException ();
 300                         uint data = value.data [0];
 301
 302                         if (value.sign == 1) {
 303                                 if (data > (uint)int.MaxValue)
 304                                         throw new OverflowException ();
 305                                 return (int)data;
 306                         } else if (value.sign == -1) {
 307                                 if (data > 0x80000000u)
 308                                         throw new OverflowException ();
 309                                 return -(int)data;
 310                         }
 311
 312                         return 0;
 313                 }
 314
 315                 [CLSCompliantAttribute (false)]
 316                 public static explicit operator uint (BigInteger value)
 317                 {
 318                         if (value.data.Length > 1 || value.sign == -1)
 319                                 throw new OverflowException ();
 320                         return value.data [0];
 321                 }
 322
 323                 public static explicit operator short (BigInteger value)
 324                 {
 325                         int val = (int)value;
 326                         if (val < short.MinValue || val > short.MaxValue)
 327                                 throw new OverflowException ();
 328                         return (short)val;
 329                 }
 330
 331                 [CLSCompliantAttribute (false)]
 332                 public static explicit operator ushort (BigInteger value)
 333                 {
 334                         uint val = (uint)value;
 335                         if (val > ushort.MaxValue)
 336                                 throw new OverflowException ();
 337                         return (ushort)val;
 338                 }
 339
 340                 public static explicit operator byte (BigInteger value)
 341                 {
 342                         uint val = (uint)value;
 343                         if (val > byte.MaxValue)
 344                                 throw new OverflowException ();
 345                         return (byte)val;
 346                 }
 347
 348                 [CLSCompliantAttribute (false)]
 349                 public static explicit operator sbyte (BigInteger value)
 350                 {
 351                         int val = (int)value;
 352                         if (val < sbyte.MinValue || val > sbyte.MaxValue)
 353                                 throw new OverflowException ();
 354                         return (sbyte)val;
 355                 }
 356
 357
 358                 public static explicit operator long (BigInteger value)
 359                 {
 360                         if (value.sign == 0)
 361                                 return 0;
 362
 363                         if (value.data.Length > 2)
 364                                 throw new OverflowException ();
 365
 366                         uint low = value.data [0];
 367
 368                         if (value.data.Length == 1) {
 369                                 if (value.sign == 1)
 370                                         return (long)low;
 371                                 long res = (long)low;
 372                                 return -res;
 373                         }
 374
 375                         uint high = value.data [1];
 376
 377                         if (value.sign == 1) {
 378                                 if (high >= 0x80000000u)
 379                                         throw new OverflowException ();
 380                                 return (((long)high) << 32) | low;
 381                         }
 382
 383                         if (high > 0x80000000u)
 384                                 throw new OverflowException ();
 385
 386                         return - ((((long)high) << 32) | (long)low);
 387                 }
 388
 389                 [CLSCompliantAttribute (false)]
 390                 public static explicit operator ulong (BigInteger value)
 391                 {
 392                         if (value.data.Length > 2 || value.sign == -1)
 393                                 throw new OverflowException ();
 394
 395                         uint low = value.data [0];
 396                         if (value.data.Length == 1)
 397                                 return low;
 398
 399                         uint high = value.data [1];
 400                         return (((ulong)high) << 32) | low;
 401                 }
 402
 403
 404                 public static implicit operator BigInteger (int value)
 405                 {
 406                         return new BigInteger (value);
 407                 }
 408
 409                 [CLSCompliantAttribute (false)]
 410                 public static implicit operator BigInteger (uint value)
 411                 {
 412                         return new BigInteger (value);
 413                 }
 414
 415                 public static implicit operator BigInteger (short value)
 416                 {
 417                         return new BigInteger (value);
 418                 }
 419
 420                 [CLSCompliantAttribute (false)]
 421                 public static implicit operator BigInteger (ushort value)
 422                 {
 423                         return new BigInteger (value);
 424                 }
 425
 426                 public static implicit operator BigInteger (byte value)
 427                 {
 428                         return new BigInteger (value);
 429                 }
 430
 431                 [CLSCompliantAttribute (false)]
 432                 public static implicit operator BigInteger (sbyte value)
 433                 {
 434                         return new BigInteger (value);
 435                 }
 436
 437                 public static implicit operator BigInteger (long value)
 438                 {
 439                         return new BigInteger (value);
 440                 }
 441
 442                 [CLSCompliantAttribute (false)]
 443                 public static implicit operator BigInteger (ulong value)
 444                 {
 445                         return new BigInteger (value);
 446                 }
 447
 448                 public static BigInteger operator+ (BigInteger left, BigInteger right)
 449                 {
 450                         if (left.sign == 0)
 451                                 return right;
 452                         if (right.sign == 0)
 453                                 return left;
 454
 455                         if (left.sign == right.sign)
 456                                 return new BigInteger (left.sign, CoreAdd (left.data, right.data));
 457
 458                         int r = CoreCompare (left.data, right.data);
 459
 460                         if (r == 0)
 461                                 return new BigInteger (0, ZERO);
 462
 463                         if (r > 0) //left > right
 464                                 return new BigInteger (left.sign, CoreSub (left.data, right.data));
 465
 466                         return new BigInteger (right.sign, CoreSub (right.data, left.data));
 467                 }
 468
 469                 public static BigInteger operator- (BigInteger left, BigInteger right)
 470                 {
 471                         if (right.sign == 0)
 472                                 return left;
 473                         if (left.sign == 0)
 474                                 return new BigInteger ((short)-right.sign, right.data);
 475
 476                         if (left.sign == right.sign) {
 477                                 int r = CoreCompare (left.data, right.data);
 478
 479                                 if (r == 0)
 480                                         return new BigInteger (0, ZERO);
 481
 482                                 if (r > 0) //left > right
 483                                         return new BigInteger (left.sign, CoreSub (left.data, right.data));
 484
 485                                 return new BigInteger ((short)-right.sign, CoreSub (right.data, left.data));
 486                         }
 487
 488                         return new BigInteger (left.sign, CoreAdd (left.data, right.data));
 489                 }
 490
 491                 public static BigInteger operator* (BigInteger left, BigInteger right)
 492                 {
 493                         if (left.sign == 0 || right.sign == 0)
 494                                 return new BigInteger (0, ZERO);
 495
 496                         if (left.data [0] == 1 && left.data.Length == 1) {
 497                                 if (left.sign == 1)
 498                                         return right;
 499                                 return new BigInteger ((short)-right.sign, right.data);
 500                         }
 501
 502                         if (right.data [0] == 1 && right.data.Length == 1) {
 503                                 if (right.sign == 1)
 504                                         return left;
 505                                 return new BigInteger ((short)-left.sign, left.data);
 506                         }
 507
 508                         uint[] a = left.data;
 509                         uint[] b = right.data;
 510
 511                         uint[] res = new uint [a.Length + b.Length];
 512
 513             for (int i = 0; i < a.Length; ++i) {
 514                 uint ai = a [i];
 515                 int k = i;
 516
 517                 ulong carry = 0;
 518                 for (int j = 0; j < b.Length; ++j) {
 519                     carry = carry + ((ulong)ai) * b [j] + res [k];
 520                     res[k++] = (uint)carry;
 521                     carry >>= 32;
 522                 }
 523
 524                 while (carry != 0) {
 525                     carry += res [k];
 526                     res[k++] = (uint)carry;
 527                     carry >>= 32;
 528                 }
 529             }
 530
 531                         int m;
 532                         for (m = res.Length - 1; m >= 0 && res [m] == 0; --m) ;
 533                         if (m < res.Length - 1)
 534                                 res = Resize (res, m + 1);
 535
 536                         return new BigInteger ((short)(left.sign * right.sign), res);
 537                 }
 538
 539                 public static BigInteger operator/ (BigInteger dividend, BigInteger divisor)
 540                 {
 541                         if (divisor.sign == 0)
 542                                 throw new DivideByZeroException ();
 543
 544                         if (dividend.sign == 0)
 545                                 return dividend;
 546
 547                         uint[] quotient;
 548                         uint[] remainder_value;
 549
 550                         DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
 551
 552                         int i;
 553                         for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
 554                         if (i == -1)
 555                                 return new BigInteger (0, ZERO);
 556                         if (i < quotient.Length - 1)
 557                                 quotient = Resize (quotient, i + 1);
 558
 559                         return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
 560                 }
 561
 562                 public static BigInteger operator% (BigInteger dividend, BigInteger divisor)
 563                 {
 564                         if (divisor.sign == 0)
 565                                 throw new DivideByZeroException ();
 566
 567                         if (dividend.sign == 0)
 568                                 return dividend;
 569
 570                         uint[] quotient;
 571                         uint[] remainder_value;
 572
 573                         DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
 574
 575                         int i;
 576                         for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
 577                         if (i == -1)
 578                                 return new BigInteger (0, ZERO);
 579
 580                         if (i < remainder_value.Length - 1)
 581                                 remainder_value = Resize (remainder_value, i + 1);
 582                         return new BigInteger (dividend.sign, remainder_value);
 583                 }
 584
 585                 public static BigInteger operator- (BigInteger value)
 586                 {
 587                         if (value.sign == 0)
 588                                 return value;
 589                         return new BigInteger ((short)-value.sign, value.data);
 590                 }
 591
 592                 public static BigInteger operator+ (BigInteger value)
 593                 {
 594                         return value;
 595                 }
 596
 597                 public static BigInteger operator++ (BigInteger value)
 598                 {
 599                         short sign = value.sign;
 600                         uint[] data = value.data;
 601                         if (data.Length == 1) {
 602                                 if (sign == -1 && data [0] == 1)
 603                                         return new BigInteger (0, ZERO);
 604                                 if (sign == 0)
 605                                         return new BigInteger (1, ONE);
 606                         }
 607
 608                         if (sign == -1)
 609                                 data = CoreSub (data, 1);
 610                         else
 611                                 data = CoreAdd (data, 1);
 612
 613                         return new BigInteger (sign, data);
 614                 }
 615
 616                 public static BigInteger operator-- (BigInteger value)
 617                 {
 618                         short sign = value.sign;
 619                         uint[] data = value.data;
 620                         if (data.Length == 1) {
 621                                 if (sign == 1 && data [0] == 1)
 622                                         return new BigInteger (0, ZERO);
 623                                 if (sign == 0)
 624                                         return new BigInteger (-1, ONE);
 625                         }
 626
 627                         if (sign == -1)
 628                                 data = CoreAdd (data, 1);
 629                         else
 630                                 data = CoreSub (data, 1);
 631
 632                         return new BigInteger (sign, data);
 633                 }
 634
 635                 public static BigInteger operator& (BigInteger left, BigInteger right)
 636                 {
 637                         if (left.sign == 0)
 638                                 return left;
 639
 640                         if (right.sign == 0)
 641                                 return right;
 642
 643                         uint[] a = left.data;
 644                         uint[] b = right.data;
 645                         int ls = left.sign;
 646                         int rs = right.sign;
 647
 648                         bool neg_res = (ls == rs) && (ls == -1);
 649
 650                         uint[] result = new uint [Math.Max (a.Length, b.Length)];
 651
 652                         ulong ac = 1, bc = 1, borrow = 1;
 653
 654                         int i;
 655                         for (i = 0; i < result.Length; ++i) {
 656                                 uint va = 0;
 657                                 if (i < a.Length)
 658                                         va = a [i];
 659                                 if (ls == -1) {
 660                                         ac = ~va + ac;
 661                                         va = (uint)ac;
 662                                         ac = (uint)(ac >> 32);
 663                                 }
 664
 665                                 uint vb = 0;
 666                                 if (i < b.Length)
 667                                         vb = b [i];
 668                                 if (rs == -1) {
 669                                         bc = ~vb + bc;
 670                                         vb = (uint)bc;
 671                                         bc = (uint)(bc >> 32);
 672                                 }
 673
 674                                 uint word = va & vb;
 675
 676                                 if (neg_res) {
 677                                         borrow = word - borrow;
 678                                         word = ~(uint)borrow;
 679                                         borrow = (uint)(borrow >> 32) & 0x1u;
 680                                 }
 681
 682                                 result [i] = word;
 683                         }
 684
 685                         for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
 686                         if (i == -1)
 687                                 return new BigInteger (0, ZERO);
 688
 689                         if (i < result.Length - 1)
 690                                 result = Resize (result, i + 1);
 691
 692                         return new BigInteger (neg_res ? (short)-1 : (short)1, result);
 693                 }
 694
 695                 public static BigInteger operator| (BigInteger left, BigInteger right)
 696                 {
 697                         if (left.sign == 0)
 698                                 return right;
 699
 700                         if (right.sign == 0)
 701                                 return left;
 702
 703                         uint[] a = left.data;
 704                         uint[] b = right.data;
 705                         int ls = left.sign;
 706                         int rs = right.sign;
 707
 708                         bool neg_res = (ls == -1) || (rs == -1);
 709
 710                         uint[] result = new uint [Math.Max (a.Length, b.Length)];
 711
 712                         ulong ac = 1, bc = 1, borrow = 1;
 713
 714                         int i;
 715                         for (i = 0; i < result.Length; ++i) {
 716                                 uint va = 0;
 717                                 if (i < a.Length)
 718                                         va = a [i];
 719                                 if (ls == -1) {
 720                                         ac = ~va + ac;
 721                                         va = (uint)ac;
 722                                         ac = (uint)(ac >> 32);
 723                                 }
 724
 725                                 uint vb = 0;
 726                                 if (i < b.Length)
 727                                         vb = b [i];
 728                                 if (rs == -1) {
 729                                         bc = ~vb + bc;
 730                                         vb = (uint)bc;
 731                                         bc = (uint)(bc >> 32);
 732                                 }
 733
 734                                 uint word = va | vb;
 735
 736                                 if (neg_res) {
 737                                         borrow = word - borrow;
 738                                         word = ~(uint)borrow;
 739                                         borrow = (uint)(borrow >> 32) & 0x1u;
 740                                 }
 741
 742                                 result [i] = word;
 743                         }
 744
 745                         for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
 746                         if (i == -1)
 747                                 return new BigInteger (0, ZERO);
 748
 749                         if (i < result.Length - 1)
 750                                 result = Resize (result, i + 1);
 751
 752                         return new BigInteger (neg_res ? (short)-1 : (short)1, result);
 753                 }
 754
 755                 public static BigInteger operator^ (BigInteger left, BigInteger right)
 756                 {
 757                         if (left.sign == 0)
 758                                 return right;
 759
 760                         if (right.sign == 0)
 761                                 return left;
 762
 763                         uint[] a = left.data;
 764                         uint[] b = right.data;
 765                         int ls = left.sign;
 766                         int rs = right.sign;
 767
 768                         bool neg_res = (ls == -1) ^ (rs == -1);
 769
 770                         uint[] result = new uint [Math.Max (a.Length, b.Length)];
 771
 772                         ulong ac = 1, bc = 1, borrow = 1;
 773
 774                         int i;
 775                         for (i = 0; i < result.Length; ++i) {
 776                                 uint va = 0;
 777                                 if (i < a.Length)
 778                                         va = a [i];
 779                                 if (ls == -1) {
 780                                         ac = ~va + ac;
 781                                         va = (uint)ac;
 782                                         ac = (uint)(ac >> 32);
 783                                 }
 784
 785                                 uint vb = 0;
 786                                 if (i < b.Length)
 787                                         vb = b [i];
 788                                 if (rs == -1) {
 789                                         bc = ~vb + bc;
 790                                         vb = (uint)bc;
 791                                         bc = (uint)(bc >> 32);
 792                                 }
 793
 794                                 uint word = va ^ vb;
 795
 796                                 if (neg_res) {
 797                                         borrow = word - borrow;
 798                                         word = ~(uint)borrow;
 799                                         borrow = (uint)(borrow >> 32) & 0x1u;
 800                                 }
 801
 802                                 result [i] = word;
 803                         }
 804
 805                         for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
 806                         if (i == -1)
 807                                 return new BigInteger (0, ZERO);
 808
 809                         if (i < result.Length - 1)
 810                                 result = Resize (result, i + 1);
 811
 812                         return new BigInteger (neg_res ? (short)-1 : (short)1, result);
 813                 }
 814
 815                 public static BigInteger operator~ (BigInteger value)
 816                 {
 817                         if (value.sign == 0)
 818                                 return new BigInteger (-1, ONE);
 819
 820                         uint[] data = value.data;
 821                         int sign = value.sign;
 822
 823                         bool neg_res = sign == 1;
 824
 825                         uint[] result = new uint [data.Length];
 826
 827                         ulong carry = 1, borrow = 1;
 828
 829                         int i;
 830                         for (i = 0; i < result.Length; ++i) {
 831                                 uint word = data [i];
 832                                 if (sign == -1) {
 833                                         carry = ~word + carry;
 834                                         word = (uint)carry;
 835                                         carry = (uint)(carry >> 32);
 836                                 }
 837
 838                                 word = ~word;
 839
 840                                 if (neg_res) {
 841                                         borrow = word - borrow;
 842                                         word = ~(uint)borrow;
 843                                         borrow = (uint)(borrow >> 32) & 0x1u;
 844                                 }
 845
 846                                 result [i] = word;
 847                         }
 848
 849                         for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
 850                         if (i == -1)
 851                                 return new BigInteger (0, ZERO);
 852
 853                         if (i < result.Length - 1)
 854                                 result = Resize (result, i + 1);
 855
 856                         return new BigInteger (neg_res ? (short)-1 : (short)1, result);
 857                 }
 858
 859                 //returns the 0-based index of the most significant set bit
 860                 //returns 0 if no bit is set, so extra care when using it
 861                 static int BitScanBackward (uint word)
 862                 {
 863                         for (int i = 31; i >= 0; --i) {
 864                                 uint mask = 1u << i;
 865                                 if ((word & mask) == mask)
 866                                         return i;
 867                         }
 868                         return 0;
 869                 }
 870
 871                 public static BigInteger operator<< (BigInteger value, int shift)
 872                 {
 873                         if (shift == 0 || value.sign == 0)
 874                                 return value;
 875                         if (shift < 0)
 876                                 return value >> -shift;
 877
 878                         uint[] data = value.data;
 879                         int sign = value.sign;
 880
 881                         int topMostIdx = BitScanBackward (data [data.Length - 1]);
 882                         int bits = shift - (31 - topMostIdx);
 883                         int extra_words = (bits >> 5) + ((bits & 0x1F) != 0 ? 1 : 0);
 884
 885                         uint[] res = new uint [data.Length + extra_words];
 886
 887                         int idx_shift = shift >> 5;
 888                         int bit_shift = shift & 0x1F;
 889                         int carry_shift = 32 - bit_shift;
 890
 891                         for (int i = 0; i < data.Length; ++i) {
 892                                 uint word = data [i];
 893                                 res [i + idx_shift] |= word << bit_shift;
 894                                 if (i + idx_shift + 1 < res.Length)
 895                                         res [i + idx_shift + 1] = word >> carry_shift;
 896                         }
 897
 898                         return new BigInteger ((short)sign, res);
 899                 }
 900
 901                 public static BigInteger operator>> (BigInteger value, int shift)
 902                 {
 903                         if (shift == 0 || value.sign == 0)
 904                                 return value;
 905                         if (shift < 0)
 906                                 return value << -shift;
 907
 908                         uint[] data = value.data;
 909                         int sign = value.sign;
 910
 911                         int topMostIdx = BitScanBackward (data [data.Length - 1]);
 912                         int idx_shift = shift >> 5;
 913                         int bit_shift = shift & 0x1F;
 914
 915                         int extra_words = idx_shift;
 916                         if (bit_shift > topMostIdx)
 917                                 ++extra_words;
 918                         int size = data.Length - extra_words;
 919
 920                         if (size <= 0) {
 921                                 if (sign == 1)
 922                                         return new BigInteger (0, ZERO);
 923                                 return new BigInteger (-1, ONE);
 924                         }
 925
 926                         uint[] res = new uint [size];
 927                         int carry_shift = 32 - bit_shift;
 928
 929                         for (int i = data.Length - 1; i >= idx_shift; --i) {
 930                                 uint word = data [i];
 931
 932                                 if (i - idx_shift < res.Length)
 933                                         res [i - idx_shift] |= word >> bit_shift;
 934                                 if (i - idx_shift - 1 >= 0)
 935                                         res [i - idx_shift - 1] = word << carry_shift;
 936                         }
 937
 938                         //Round down instead of toward zero
 939                         if (sign == -1) {
 940                                 for (int i = 0; i < idx_shift; i++) {
 941                                         if (data [i] != 0u) {
 942                                                 var tmp = new BigInteger ((short)sign, res);
 943                                                 --tmp;
 944                                                 return tmp;
 945                                         }
 946                                 }
 947                                 if (bit_shift > 0 && (data [idx_shift] << carry_shift) != 0u) {
 948                                         var tmp = new BigInteger ((short)sign, res);
 949                                         --tmp;
 950                                         return tmp;
 951                                 }
 952                         }
 953                         return new BigInteger ((short)sign, res);
 954                 }
 955
 956                 public static bool operator< (BigInteger left, BigInteger right)
 957                 {
 958                         return Compare (left, right) < 0;
 959                 }
 960
 961                 public static bool operator< (BigInteger left, long right)
 962                 {
 963                         return left.CompareTo (right) < 0;
 964                 }
 965
 966
 967                 public static bool operator< (long left, BigInteger right)
 968                 {
 969                         return right.CompareTo (left) > 0;
 970                 }
 971
 972
 973                 [CLSCompliantAttribute (false)]
 974                 public static bool operator< (BigInteger left, ulong right)
 975                 {
 976                         return left.CompareTo (right) < 0;
 977                 }
 978
 979                 [CLSCompliantAttribute (false)]
 980                 public static bool operator< (ulong left, BigInteger right)
 981                 {
 982                         return right.CompareTo (left) > 0;
 983                 }
 984
 985                 public static bool operator<= (BigInteger left, BigInteger right)
 986                 {
 987                         return Compare (left, right) <= 0;
 988                 }
 989
 990                 public static bool operator<= (BigInteger left, long right)
 991                 {
 992                         return left.CompareTo (right) <= 0;
 993                 }
 994
 995                 public static bool operator<= (long left, BigInteger right)
 996                 {
 997                         return right.CompareTo (left) >= 0;
 998                 }
 999
1000                 [CLSCompliantAttribute (false)]
1001                 public static bool operator<= (BigInteger left, ulong right)
1002                 {
1003                         return left.CompareTo (right) <= 0;
1004                 }
1005
1006                 [CLSCompliantAttribute (false)]
1007                 public static bool operator<= (ulong left, BigInteger right)
1008                 {
1009                         return right.CompareTo (left) >= 0;
1010                 }
1011
1012                 public static bool operator> (BigInteger left, BigInteger right)
1013                 {
1014                         return Compare (left, right) > 0;
1015                 }
1016
1017                 public static bool operator> (BigInteger left, long right)
1018                 {
1019                         return left.CompareTo (right) > 0;
1020                 }
1021
1022                 public static bool operator> (long left, BigInteger right)
1023                 {
1024                         return right.CompareTo (left) < 0;
1025                 }
1026
1027                 [CLSCompliantAttribute (false)]
1028                 public static bool operator> (BigInteger left, ulong right)
1029                 {
1030                         return left.CompareTo (right) > 0;
1031                 }
1032
1033                 [CLSCompliantAttribute (false)]
1034                 public static bool operator> (ulong left, BigInteger right)
1035                 {
1036                         return right.CompareTo (left) < 0;
1037                 }
1038
1039                 public static bool operator>= (BigInteger left, BigInteger right)
1040                 {
1041                         return Compare (left, right) >= 0;
1042                 }
1043
1044                 public static bool operator>= (BigInteger left, long right)
1045                 {
1046                         return left.CompareTo (right) >= 0;
1047                 }
1048
1049                 public static bool operator>= (long left, BigInteger right)
1050                 {
1051                         return right.CompareTo (left) <= 0;
1052                 }
1053
1054                 [CLSCompliantAttribute (false)]
1055                 public static bool operator>= (BigInteger left, ulong right)
1056                 {
1057                         return left.CompareTo (right) >= 0;
1058                 }
1059
1060                 [CLSCompliantAttribute (false)]
1061                 public static bool operator>= (ulong left, BigInteger right)
1062                 {
1063                         return right.CompareTo (left) <= 0;
1064                 }
1065
1066                 public static bool operator== (BigInteger left, BigInteger right)
1067                 {
1068                         return Compare (left, right) == 0;
1069                 }
1070
1071                 public static bool operator== (BigInteger left, long right)
1072                 {
1073                         return left.CompareTo (right) == 0;
1074                 }
1075
1076                 public static bool operator== (long left, BigInteger right)
1077                 {
1078                         return right.CompareTo (left) == 0;
1079                 }
1080
1081                 [CLSCompliantAttribute (false)]
1082                 public static bool operator== (BigInteger left, ulong right)
1083                 {
1084                         return left.CompareTo (right) == 0;
1085                 }
1086
1087                 [CLSCompliantAttribute (false)]
1088                 public static bool operator== (ulong left, BigInteger right)
1089                 {
1090                         return right.CompareTo (left) == 0;
1091                 }
1092
1093                 public static bool operator!= (BigInteger left, BigInteger right)
1094                 {
1095                         return Compare (left, right) != 0;
1096                 }
1097
1098                 public static bool operator!= (BigInteger left, long right)
1099                 {
1100                         return left.CompareTo (right) != 0;
1101                 }
1102
1103                 public static bool operator!= (long left, BigInteger right)
1104                 {
1105                         return right.CompareTo (left) != 0;
1106                 }
1107
1108                 [CLSCompliantAttribute (false)]
1109                 public static bool operator!= (BigInteger left, ulong right)
1110                 {
1111                         return left.CompareTo (right) != 0;
1112                 }
1113
1114                 [CLSCompliantAttribute (false)]
1115                 public static bool operator!= (ulong left, BigInteger right)
1116                 {
1117                         return right.CompareTo (left) != 0;
1118                 }
1119
1120                 public override bool Equals (object obj)
1121                 {
1122                         if (!(obj is BigInteger))
1123                                 return false;
1124                         return Equals ((BigInteger)obj);
1125                 }
1126
1127                 public bool Equals (BigInteger other)
1128                 {
1129                         if (sign != other.sign)
1130                                 return false;
1131                         if (data.Length != other.data.Length)
1132                                 return false;
1133                         for (int i = 0; i < data.Length; ++i) {
1134                                 if (data [i] != other.data [i])
1135                                         return false;
1136                         }
1137                         return true;
1138                 }
1139
1140                 public bool Equals (long other)
1141                 {
1142                         return CompareTo (other) == 0;
1143                 }
1144
1145                 public static BigInteger Min (BigInteger left, BigInteger right)
1146                 {
1147                         int ls = left.sign;
1148                         int rs = right.sign;
1149
1150                         if (ls < rs)
1151                                 return left;
1152                         if (rs < ls)
1153                                 return right;
1154
1155                         int r = CoreCompare (left.data, right.data);
1156                         if (ls == -1)
1157                                 r = -r;
1158
1159                         if (r <= 0)
1160                                 return left;
1161                         return right;
1162                 }
1163
1164
1165                 public static BigInteger Max (BigInteger left, BigInteger right)
1166                 {
1167                         int ls = left.sign;
1168                         int rs = right.sign;
1169
1170                         if (ls > rs)
1171                                 return left;
1172                         if (rs > ls)
1173                                 return right;
1174
1175                         int r = CoreCompare (left.data, right.data);
1176                         if (ls == -1)
1177                                 r = -r;
1178
1179                         if (r >= 0)
1180                                 return left;
1181                         return right;
1182                 }
1183
1184                 public static BigInteger Abs (BigInteger value)
1185                 {
1186                         return new BigInteger ((short)Math.Abs (value.sign), value.data);
1187                 }
1188
1189
1190                 public static BigInteger DivRem (BigInteger dividend, BigInteger divisor, out BigInteger remainder)
1191                 {
1192                         if (divisor.sign == 0)
1193                                 throw new DivideByZeroException ();
1194
1195                         if (dividend.sign == 0) {
1196                                 remainder = dividend;
1197                                 return dividend;
1198                         }
1199
1200                         uint[] quotient;
1201                         uint[] remainder_value;
1202
1203                         DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
1204
1205                         int i;
1206                         for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
1207                         if (i == -1) {
1208                                 remainder = new BigInteger (0, ZERO);
1209                         } else {
1210                                 if (i < remainder_value.Length - 1)
1211                                         remainder_value = Resize (remainder_value, i + 1);
1212                                 remainder = new BigInteger (dividend.sign, remainder_value);
1213                         }
1214
1215                         for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
1216                         if (i == -1)
1217                                 return new BigInteger (0, ZERO);
1218                         if (i < quotient.Length - 1)
1219                                 quotient = Resize (quotient, i + 1);
1220
1221                         return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
1222                 }
1223
1224         public static BigInteger Pow (BigInteger value, int exponent)
1225                 {
1226                         if (exponent < 0)
1227                                 throw new ArgumentOutOfRangeException("exp", "exp must be >= 0");
1228                         if (exponent == 0)
1229                                 return One;
1230                         if (exponent == 1)
1231                                 return value;
1232
1233                         BigInteger result = One;
1234                         while (exponent != 0) {
1235                                 if ((exponent & 1) != 0)
1236                                         result = result * value;
1237                                 if (exponent == 1)
1238                                         break;
1239
1240                                 value = value * value;
1241                                 exponent >>= 1;
1242                         }
1243                         return result;
1244         }
1245
1246
1247                 [CLSCompliantAttribute (false)]
1248                 public bool Equals (ulong other)
1249                 {
1250                         return CompareTo (other) == 0;
1251                 }
1252
1253                 public override int GetHashCode ()
1254                 {
1255                         uint hash = (uint)(sign * 0x01010101u);
1256
1257                         for (int i = 0; i < data.Length; ++i)
1258                                 hash ^= data [i];
1259                         return (int)hash;
1260                 }
1261
1262                 public static BigInteger Add (BigInteger left, BigInteger right)
1263                 {
1264                         return left + right;
1265                 }
1266
1267                 public static BigInteger Subtract (BigInteger left, BigInteger right)
1268                 {
1269                         return left - right;
1270                 }
1271
1272                 public static BigInteger Multiply (BigInteger left, BigInteger right)
1273                 {
1274                         return left * right;
1275                 }
1276
1277                 public static BigInteger Divide (BigInteger dividend, BigInteger divisor)
1278                 {
1279                         return dividend / divisor;
1280                 }
1281
1282                 public static BigInteger Remainder (BigInteger dividend, BigInteger divisor)
1283                 {
1284                         return dividend % divisor;
1285                 }
1286
1287                 public static BigInteger Negate (BigInteger value)
1288                 {
1289                         return - value;
1290                 }
1291
1292                 public int CompareTo (object obj)
1293                 {
1294                         if (obj == null)
1295                                 return 1;
1296
1297                         if (!(obj is BigInteger))
1298                                 return -1;
1299
1300                         return Compare (this, (BigInteger)obj);
1301                 }
1302
1303                 public int CompareTo (BigInteger other)
1304                 {
1305                         return Compare (this, other);
1306                 }
1307
1308                 [CLSCompliantAttribute (false)]
1309                 public int CompareTo (ulong other)
1310                 {
1311                         if (sign < 0)
1312                                 return -1;
1313                         if (sign == 0)
1314                                 return other == 0 ? 0 : -1;
1315
1316                         if (data.Length > 2)
1317                                 return 1;
1318
1319                         uint high = (uint)(other >> 32);
1320                         uint low = (uint)other;
1321
1322                         return LongCompare (low, high);
1323                 }
1324
1325                 int LongCompare (uint low, uint high)
1326                 {
1327                         uint h = 0;
1328                         if (data.Length > 1)
1329                                 h = data [1];
1330
1331                         if (h > high)
1332                                 return 1;
1333                         if (h < high)
1334                                 return -1;
1335
1336                         uint l = data [0];
1337
1338                         if (l > low)
1339                                 return 1;
1340                         if (l < low)
1341                                 return -1;
1342
1343                         return 0;
1344                 }
1345
1346                 public int CompareTo (long other)
1347                 {
1348                         int ls = sign;
1349                         int rs = Math.Sign (other);
1350
1351                         if (ls != rs)
1352                                 return ls > rs ? 1 : -1;
1353
1354                         if (ls == 0)
1355                                 return 0;
1356
1357                         if (data.Length > 2)
1358                                 return -sign;
1359
1360                         if (other < 0)
1361                                 other = -other;
1362                         uint low = (uint)other;
1363                         uint high = (uint)((ulong)other >> 32);
1364
1365                         int r = LongCompare (low, high);
1366                         if (ls == -1)
1367                                 r = -r;
1368
1369                         return r;
1370                 }
1371
1372                 public static int Compare (BigInteger left, BigInteger right)
1373                 {
1374                         int ls = left.sign;
1375                         int rs = right.sign;
1376
1377                         if (ls != rs)
1378                                 return ls > rs ? 1 : -1;
1379
1380                         int r = CoreCompare (left.data, right.data);
1381                         if (ls < 0)
1382                                 r = -r;
1383                         return r;
1384                 }
1385
1386
1387                 static int TopByte (uint x)
1388                 {
1389                         if ((x & 0xFFFF0000u) != 0) {
1390                                 if ((x & 0xFF000000u) != 0)
1391                                         return 4;
1392                                 return 3;
1393                         }
1394                         if ((x & 0xFF00u) != 0)
1395                                 return 2;
1396                         return 1;
1397                 }
1398
1399                 static int FirstNonFFByte (uint word)
1400                 {
1401                         if ((word & 0xFF000000u) != 0xFF000000u)
1402                                 return 4;
1403                         else if ((word & 0xFF0000u) != 0xFF0000u)
1404                                 return 3;
1405                         else if ((word & 0xFF00u) != 0xFF00u)
1406                                 return 2;
1407                         return 1;
1408                 }
1409
1410                 public byte[] ToByteArray ()
1411                 {
1412                         if (sign == 0)
1413                                 return new byte [1];
1414
1415                         //number of bytes not counting upper word
1416                         int bytes = (data.Length - 1) * 4;
1417                         bool needExtraZero = false;
1418
1419                         uint topWord = data [data.Length - 1];
1420                         int extra;
1421
1422                         //if the topmost bit is set we need an extra
1423                         if (sign == 1) {
1424                                 extra = TopByte (topWord);
1425                                 uint mask = 0x80u << ((extra - 1) * 8);
1426                                 if ((topWord & mask) != 0) {
1427                                         needExtraZero = true;
1428                                 }
1429                         } else {
1430                                 extra = TopByte (topWord);
1431                         }
1432
1433                         byte[] res = new byte [bytes + extra + (needExtraZero ? 1 : 0) ];
1434                         if (sign == 1) {
1435                                 int j = 0;
1436                                 int end = data.Length - 1;
1437                                 for (int i = 0; i < end; ++i) {
1438                                         uint word = data [i];
1439
1440                                         res [j++] = (byte)word;
1441                                         res [j++] = (byte)(word >> 8);
1442                                         res [j++] = (byte)(word >> 16);
1443                                         res [j++] = (byte)(word >> 24);
1444                                 }
1445                                 while (extra-- > 0) {
1446                                         res [j++] = (byte)topWord;
1447                                         topWord >>= 8;
1448                                 }
1449                         } else {
1450                                 int j = 0;
1451                                 int end = data.Length - 1;
1452
1453                                 uint carry = 1, word;
1454                                 ulong add;
1455                                 for (int i = 0; i < end; ++i) {
1456                                         word = data [i];
1457                                         add = (ulong)~word + carry;
1458                                         word = (uint)add;
1459                                         carry = (uint)(add >> 32);
1460
1461                                         res [j++] = (byte)word;
1462                                         res [j++] = (byte)(word >> 8);
1463                                         res [j++] = (byte)(word >> 16);
1464                                         res [j++] = (byte)(word >> 24);
1465                                 }
1466
1467                                 add = (ulong)~topWord + (carry);
1468                                 word = (uint)add;
1469                                 carry = (uint)(add >> 32);
1470                                 if (carry == 0) {
1471                                         int ex = FirstNonFFByte (word);
1472                                         bool needExtra = (word & (1 << (ex * 8 - 1))) == 0;
1473                                         int to = ex + (needExtra ? 1 : 0);
1474
1475                                         if (to != extra)
1476                                                 res = Resize (res, bytes + to);
1477
1478                                         while (ex-- > 0) {
1479                                                 res [j++] = (byte)word;
1480                                                 word >>= 8;
1481                                         }
1482                                         if (needExtra)
1483                                                 res [j++] = 0xFF;
1484                                 } else {
1485                                         res = Resize (res, bytes + 5);
1486                                         res [j++] = (byte)word;
1487                                         res [j++] = (byte)(word >> 8);
1488                                         res [j++] = (byte)(word >> 16);
1489                                         res [j++] = (byte)(word >> 24);
1490                                         res [j++] = 0xFF;
1491                                 }
1492                         }
1493
1494                         return res;
1495                 }
1496
1497                 static byte[] Resize (byte[] v, int len)
1498                 {\r
1499                         byte[] res = new byte [len];\r
1500                         Array.Copy (v, res, Math.Min (v.Length, len));\r
1501                         return res;\r
1502                 }
1503
1504                 static uint[] Resize (uint[] v, int len)
1505                 {
1506                         uint[] res = new uint [len];\r
1507                         Array.Copy (v, res, Math.Min (v.Length, len));\r
1508                         return res;\r
1509                 }
1510
1511                 static uint[] CoreAdd (uint[] a, uint[] b)
1512                 {
1513                         if (a.Length < b.Length) {
1514                                 uint[] tmp = a;
1515                                 a = b;
1516                                 b = tmp;
1517                         }
1518
1519                         int bl = a.Length;
1520                         int sl = b.Length;
1521
1522                         uint[] res = new uint [bl];
1523
1524                         ulong sum = 0;
1525
1526                         int i = 0;
1527                         for (; i < sl; i++) {
1528                                 sum = sum + a [i] + b [i];
1529                                 res [i] = (uint)sum;
1530                                 sum >>= 32;
1531                         }
1532
1533                         for (; i < bl; i++) {
1534                                 sum = sum + a [i];
1535                                 res [i] = (uint)sum;
1536                                 sum >>= 32;
1537                         }
1538
1539                         if (sum != 0) {
1540                                 res = Resize (res, bl + 1);
1541                                 res [i] = (uint)sum;
1542                         }
1543
1544                         return res;
1545                 }
1546
1547                 /*invariant a > b*/
1548                 static uint[] CoreSub (uint[] a, uint[] b)
1549                 {
1550                         int bl = a.Length;
1551                         int sl = b.Length;
1552
1553                         uint[] res = new uint [bl];
1554
1555                         ulong borrow = 0;
1556                         int i;
1557                         for (i = 0; i < sl; ++i) {
1558                                 borrow = (ulong)a [i] - b [i] - borrow;
1559
1560                                 res [i] = (uint)borrow;
1561                                 borrow = (borrow >> 32) & 0x1;
1562                         }
1563
1564                         for (; i < bl; i++) {
1565                                 borrow = (ulong)a [i] - borrow;
1566                                 res [i] = (uint)borrow;
1567                                 borrow = (borrow >> 32) & 0x1;
1568                         }
1569
1570                         //remove extra zeroes
1571                         for (i = bl - 1; i >= 0 && res [i] == 0; --i) ;
1572                         if (i < bl - 1)
1573                                 res = Resize (res, i + 1);
1574
1575             return res;
1576                 }
1577
1578
1579                 static uint[] CoreAdd (uint[] a, uint b)
1580                 {
1581                         int len = a.Length;
1582                         uint[] res = new uint [len];
1583
1584                         ulong sum = b;
1585                         int i;
1586                         for (i = 0; i < len; i++) {
1587                                 sum = sum + a [i];
1588                                 res [i] = (uint)sum;
1589                                 sum >>= 32;
1590                         }
1591
1592                         if (sum != 0) {
1593                                 res = Resize (res, len + 1);
1594                                 res [i] = (uint)sum;
1595                         }
1596
1597                         return res;
1598                 }
1599
1600                 static uint[] CoreSub (uint[] a, uint b)
1601                 {
1602                         int len = a.Length;
1603                         uint[] res = new uint [len];
1604
1605                         ulong borrow = b;
1606                         int i;
1607                         for (i = 0; i < len; i++) {
1608                                 borrow = (ulong)a [i] - borrow;
1609                                 res [i] = (uint)borrow;
1610                                 borrow = (borrow >> 32) & 0x1;
1611                         }
1612
1613                         //remove extra zeroes
1614                         for (i = len - 1; i >= 0 && res [i] == 0; --i) ;
1615                         if (i < len - 1)
1616                                 res = Resize (res, i + 1);
1617
1618             return res;
1619                 }
1620
1621                 static int CoreCompare (uint[] a, uint[] b)
1622                 {
1623                         int     al = a.Length;
1624                         int bl = b.Length;
1625
1626                         if (al > bl)
1627                                 return 1;
1628                         if (bl > al)
1629                                 return -1;
1630
1631                         for (int i = al - 1; i >= 0; --i) {
1632                                 uint ai = a [i];
1633                                 uint bi = b [i];
1634                                 if (ai > bi)
1635                                         return 1;
1636                                 if (ai < bi)
1637                                         return -1;
1638                         }
1639                         return 0;
1640                 }
1641
1642                 static int GetNormalizeShift(uint value) {
1643                         int shift = 0;
1644
1645                         if ((value & 0xFFFF0000) == 0) { value <<= 16; shift += 16; }
1646                         if ((value & 0xFF000000) == 0) { value <<= 8; shift += 8; }
1647                         if ((value & 0xF0000000) == 0) { value <<= 4; shift += 4; }
1648                         if ((value & 0xC0000000) == 0) { value <<= 2; shift += 2; }
1649                         if ((value & 0x80000000) == 0) { value <<= 1; shift += 1; }
1650
1651                         return shift;
1652                 }
1653
1654                 static void Normalize (uint[] u, int l, uint[] un, int shift)
1655                 {
1656                         uint carry = 0;
1657                         int i;
1658                         if (shift > 0) {
1659                                 int rshift = 32 - shift;
1660                                 for (i = 0; i < l; i++) {
1661                                         uint ui = u [i];
1662                                         un [i] = (ui << shift) | carry;
1663                                         carry = ui >> rshift;
1664                                 }
1665                         } else {
1666                                 for (i = 0; i < l; i++) {
1667                                         un [i] = u [i];
1668                                 }
1669                         }
1670
1671                         while (i < un.Length) {
1672                                 un [i++] = 0;
1673                         }
1674
1675                         if (carry != 0) {
1676                                 un [l] = carry;
1677                         }
1678                 }
1679
1680                 static void Unnormalize (uint[] un, out uint[] r, int shift)
1681                 {
1682                         int length = un.Length;
1683                         r = new uint [length];
1684
1685                         if (shift > 0) {
1686                                 int lshift = 32 - shift;
1687                                 uint carry = 0;
1688                                 for (int i = length - 1; i >= 0; i--) {
1689                                         uint uni = un [i];
1690                                         r [i] = (uni >> shift) | carry;
1691                                         carry = (uni << lshift);
1692                                 }
1693                         } else {
1694                                 for (int i = 0; i < length; i++) {
1695                                         r [i] = un [i];
1696                                 }
1697                         }
1698                 }
1699
1700                 const ulong Base = 0x100000000;
1701                 static void DivModUnsigned (uint[] u, uint[] v, out uint[] q, out uint[] r)
1702                 {
1703                         int m = u.Length;
1704                         int n = v.Length;
1705
1706                         if (n <= 1) {
1707                                 //  Divide by single digit
1708                                 //
1709                                 ulong rem = 0;
1710                                 uint v0 = v [0];
1711                                 q = new uint[m];
1712                                 r = new uint [1];
1713
1714                                 for (int j = m - 1; j >= 0; j--) {
1715                                         rem *= Base;
1716                                         rem += u[j];
1717
1718                                         ulong div = rem / v0;
1719                                         rem -= div * v0;
1720                                         q[j] = (uint)div;
1721                                 }
1722                                 r [0] = (uint)rem;
1723                         } else if (m >= n) {
1724                                 int shift = GetNormalizeShift (v [n - 1]);
1725
1726                                 uint[] un = new uint [m + 1];
1727                                 uint[] vn = new uint [n];
1728
1729                                 Normalize (u, m, un, shift);
1730                                 Normalize (v, n, vn, shift);
1731
1732                                 q = new uint [m - n + 1];
1733                                 r = null;
1734
1735                                 //  Main division loop
1736                                 //
1737                                 for (int j = m - n; j >= 0; j--) {
1738                                         ulong rr, qq;
1739                                         int i;
1740
1741                                         rr = Base * un [j + n] + un [j + n - 1];
1742                                         qq = rr / vn [n - 1];
1743                                         rr -= qq * vn [n - 1];
1744
1745                                         for (; ; ) {
1746                                                 // Estimate too big ?
1747                                                 //
1748                                                 if ((qq >= Base) || (qq * vn [n - 2] > (rr * Base + un [j + n - 2]))) {
1749                                                         qq--;
1750                                                         rr += (ulong)vn [n - 1];
1751                                                         if (rr < Base)
1752                                                                 continue;
1753                                                 }
1754                                                 break;
1755                                         }
1756
1757
1758                                         //  Multiply and subtract
1759                                         //
1760                                         long b = 0;
1761                                         long t = 0;
1762                                         for (i = 0; i < n; i++) {
1763                                                 ulong p = vn [i] * qq;
1764                                                 t = (long)un [i + j] - (long)(uint)p - b;
1765                                                 un [i + j] = (uint)t;
1766                                                 p >>= 32;
1767                                                 t >>= 32;
1768                                                 b = (long)p - t;
1769                                         }
1770                                         t = (long)un [j + n] - b;
1771                                         un [j + n] = (uint)t;
1772
1773                                         //  Store the calculated value
1774                                         //
1775                                         q [j] = (uint)qq;
1776
1777                                         //  Add back vn[0..n] to un[j..j+n]
1778                                         //
1779                                         if (t < 0) {
1780                                                 q [j]--;
1781                                                 ulong c = 0;
1782                                                 for (i = 0; i < n; i++) {
1783                                                         c = (ulong)vn [i] + un [j + i] + c;
1784                                                         un [j + i] = (uint)c;
1785                                                         c >>= 32;
1786                                                 }
1787                                                 c += (ulong)un [j + n];
1788                                                 un [j + n] = (uint)c;
1789                                         }
1790                                 }
1791
1792                                 Unnormalize (un, out r, shift);
1793                         } else {
1794                                 q = new uint [] { 0 };
1795                                 r = u;
1796                         }
1797                 }
1798         }
1799 }\r