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("exponent", "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                 public static BigInteger ModPow (BigInteger value, BigInteger exponent, BigInteger modulus) {
1247                         if (exponent.sign == -1)
1248                                 throw new ArgumentOutOfRangeException("exponent", "power must be >= 0");
1249                         if (modulus.sign == 0)
1250                                 throw new DivideByZeroException ();
1251
1252                         BigInteger result = One % modulus;
1253                         while (exponent.sign != 0) {
1254                                 if (!exponent.IsEven) {
1255                                         result = result * value;
1256                                         result = result % modulus;
1257                                 }
1258                                 if (exponent.IsOne)
1259                                         break;
1260                                 value = value * value;
1261                                 value = value % modulus;
1262                                 exponent >>= 1;
1263                         }
1264                         return result;
1265                 }
1266
1267                 public static BigInteger GreatestCommonDivisor (BigInteger left, BigInteger right)
1268                 {
1269                         if (left.data.Length == 1 && left.data [0] == 1)
1270                                 return new BigInteger (1, ONE);
1271                         if (right.data.Length == 1 && right.data [0] == 1)
1272                                 return new BigInteger (1, ONE);
1273                         if (left.IsZero)
1274                                 return right;
1275                         if (right.IsZero)
1276                                 return left;
1277
1278                         BigInteger x = new BigInteger (1, left.data);
1279                         BigInteger y = new BigInteger (1, right.data);
1280
1281                         BigInteger g = y;
1282
1283                         while (x.data.Length > 1) {
1284                                 g = x;
1285                                 x = y % x;
1286                                 y = g;
1287
1288                         }
1289                         if (x.IsZero) return g;
1290
1291                         // TODO: should we have something here if we can convert to long?
1292
1293                         //
1294                         // Now we can just do it with single precision. I am using the binary gcd method,
1295                         // as it should be faster.
1296                         //
1297
1298                         uint yy = x.data [0];
1299                         uint xx = (uint)(y % yy);
1300
1301                         int t = 0;
1302
1303                         while (((xx | yy) & 1) == 0) {
1304                                 xx >>= 1; yy >>= 1; t++;
1305                         }
1306                         while (xx != 0) {
1307                                 while ((xx & 1) == 0) xx >>= 1;
1308                                 while ((yy & 1) == 0) yy >>= 1;
1309                                 if (xx >= yy)
1310                                         xx = (xx - yy) >> 1;
1311                                 else
1312                                         yy = (yy - xx) >> 1;
1313                         }
1314
1315                         return yy << t;
1316                 }
1317
1318                 [CLSCompliantAttribute (false)]
1319                 public bool Equals (ulong other)
1320                 {
1321                         return CompareTo (other) == 0;
1322                 }
1323
1324                 public override int GetHashCode ()
1325                 {
1326                         uint hash = (uint)(sign * 0x01010101u);
1327
1328                         for (int i = 0; i < data.Length; ++i)
1329                                 hash ^= data [i];
1330                         return (int)hash;
1331                 }
1332
1333                 public static BigInteger Add (BigInteger left, BigInteger right)
1334                 {
1335                         return left + right;
1336                 }
1337
1338                 public static BigInteger Subtract (BigInteger left, BigInteger right)
1339                 {
1340                         return left - right;
1341                 }
1342
1343                 public static BigInteger Multiply (BigInteger left, BigInteger right)
1344                 {
1345                         return left * right;
1346                 }
1347
1348                 public static BigInteger Divide (BigInteger dividend, BigInteger divisor)
1349                 {
1350                         return dividend / divisor;
1351                 }
1352
1353                 public static BigInteger Remainder (BigInteger dividend, BigInteger divisor)
1354                 {
1355                         return dividend % divisor;
1356                 }
1357
1358                 public static BigInteger Negate (BigInteger value)
1359                 {
1360                         return - value;
1361                 }
1362
1363                 public int CompareTo (object obj)
1364                 {
1365                         if (obj == null)
1366                                 return 1;
1367
1368                         if (!(obj is BigInteger))
1369                                 return -1;
1370
1371                         return Compare (this, (BigInteger)obj);
1372                 }
1373
1374                 public int CompareTo (BigInteger other)
1375                 {
1376                         return Compare (this, other);
1377                 }
1378
1379                 [CLSCompliantAttribute (false)]
1380                 public int CompareTo (ulong other)
1381                 {
1382                         if (sign < 0)
1383                                 return -1;
1384                         if (sign == 0)
1385                                 return other == 0 ? 0 : -1;
1386
1387                         if (data.Length > 2)
1388                                 return 1;
1389
1390                         uint high = (uint)(other >> 32);
1391                         uint low = (uint)other;
1392
1393                         return LongCompare (low, high);
1394                 }
1395
1396                 int LongCompare (uint low, uint high)
1397                 {
1398                         uint h = 0;
1399                         if (data.Length > 1)
1400                                 h = data [1];
1401
1402                         if (h > high)
1403                                 return 1;
1404                         if (h < high)
1405                                 return -1;
1406
1407                         uint l = data [0];
1408
1409                         if (l > low)
1410                                 return 1;
1411                         if (l < low)
1412                                 return -1;
1413
1414                         return 0;
1415                 }
1416
1417                 public int CompareTo (long other)
1418                 {
1419                         int ls = sign;
1420                         int rs = Math.Sign (other);
1421
1422                         if (ls != rs)
1423                                 return ls > rs ? 1 : -1;
1424
1425                         if (ls == 0)
1426                                 return 0;
1427
1428                         if (data.Length > 2)
1429                                 return -sign;
1430
1431                         if (other < 0)
1432                                 other = -other;
1433                         uint low = (uint)other;
1434                         uint high = (uint)((ulong)other >> 32);
1435
1436                         int r = LongCompare (low, high);
1437                         if (ls == -1)
1438                                 r = -r;
1439
1440                         return r;
1441                 }
1442
1443                 public static int Compare (BigInteger left, BigInteger right)
1444                 {
1445                         int ls = left.sign;
1446                         int rs = right.sign;
1447
1448                         if (ls != rs)
1449                                 return ls > rs ? 1 : -1;
1450
1451                         int r = CoreCompare (left.data, right.data);
1452                         if (ls < 0)
1453                                 r = -r;
1454                         return r;
1455                 }
1456
1457
1458                 static int TopByte (uint x)
1459                 {
1460                         if ((x & 0xFFFF0000u) != 0) {
1461                                 if ((x & 0xFF000000u) != 0)
1462                                         return 4;
1463                                 return 3;
1464                         }
1465                         if ((x & 0xFF00u) != 0)
1466                                 return 2;
1467                         return 1;
1468                 }
1469
1470                 static int FirstNonFFByte (uint word)
1471                 {
1472                         if ((word & 0xFF000000u) != 0xFF000000u)
1473                                 return 4;
1474                         else if ((word & 0xFF0000u) != 0xFF0000u)
1475                                 return 3;
1476                         else if ((word & 0xFF00u) != 0xFF00u)
1477                                 return 2;
1478                         return 1;
1479                 }
1480
1481                 public byte[] ToByteArray ()
1482                 {
1483                         if (sign == 0)
1484                                 return new byte [1];
1485
1486                         //number of bytes not counting upper word
1487                         int bytes = (data.Length - 1) * 4;
1488                         bool needExtraZero = false;
1489
1490                         uint topWord = data [data.Length - 1];
1491                         int extra;
1492
1493                         //if the topmost bit is set we need an extra
1494                         if (sign == 1) {
1495                                 extra = TopByte (topWord);
1496                                 uint mask = 0x80u << ((extra - 1) * 8);
1497                                 if ((topWord & mask) != 0) {
1498                                         needExtraZero = true;
1499                                 }
1500                         } else {
1501                                 extra = TopByte (topWord);
1502                         }
1503
1504                         byte[] res = new byte [bytes + extra + (needExtraZero ? 1 : 0) ];
1505                         if (sign == 1) {
1506                                 int j = 0;
1507                                 int end = data.Length - 1;
1508                                 for (int i = 0; i < end; ++i) {
1509                                         uint word = data [i];
1510
1511                                         res [j++] = (byte)word;
1512                                         res [j++] = (byte)(word >> 8);
1513                                         res [j++] = (byte)(word >> 16);
1514                                         res [j++] = (byte)(word >> 24);
1515                                 }
1516                                 while (extra-- > 0) {
1517                                         res [j++] = (byte)topWord;
1518                                         topWord >>= 8;
1519                                 }
1520                         } else {
1521                                 int j = 0;
1522                                 int end = data.Length - 1;
1523
1524                                 uint carry = 1, word;
1525                                 ulong add;
1526                                 for (int i = 0; i < end; ++i) {
1527                                         word = data [i];
1528                                         add = (ulong)~word + carry;
1529                                         word = (uint)add;
1530                                         carry = (uint)(add >> 32);
1531
1532                                         res [j++] = (byte)word;
1533                                         res [j++] = (byte)(word >> 8);
1534                                         res [j++] = (byte)(word >> 16);
1535                                         res [j++] = (byte)(word >> 24);
1536                                 }
1537
1538                                 add = (ulong)~topWord + (carry);
1539                                 word = (uint)add;
1540                                 carry = (uint)(add >> 32);
1541                                 if (carry == 0) {
1542                                         int ex = FirstNonFFByte (word);
1543                                         bool needExtra = (word & (1 << (ex * 8 - 1))) == 0;
1544                                         int to = ex + (needExtra ? 1 : 0);
1545
1546                                         if (to != extra)
1547                                                 res = Resize (res, bytes + to);
1548
1549                                         while (ex-- > 0) {
1550                                                 res [j++] = (byte)word;
1551                                                 word >>= 8;
1552                                         }
1553                                         if (needExtra)
1554                                                 res [j++] = 0xFF;
1555                                 } else {
1556                                         res = Resize (res, bytes + 5);
1557                                         res [j++] = (byte)word;
1558                                         res [j++] = (byte)(word >> 8);
1559                                         res [j++] = (byte)(word >> 16);
1560                                         res [j++] = (byte)(word >> 24);
1561                                         res [j++] = 0xFF;
1562                                 }
1563                         }
1564
1565                         return res;
1566                 }
1567
1568                 static byte[] Resize (byte[] v, int len)
1569                 {\r
1570                         byte[] res = new byte [len];\r
1571                         Array.Copy (v, res, Math.Min (v.Length, len));\r
1572                         return res;\r
1573                 }
1574
1575                 static uint[] Resize (uint[] v, int len)
1576                 {
1577                         uint[] res = new uint [len];\r
1578                         Array.Copy (v, res, Math.Min (v.Length, len));\r
1579                         return res;\r
1580                 }
1581
1582                 static uint[] CoreAdd (uint[] a, uint[] b)
1583                 {
1584                         if (a.Length < b.Length) {
1585                                 uint[] tmp = a;
1586                                 a = b;
1587                                 b = tmp;
1588                         }
1589
1590                         int bl = a.Length;
1591                         int sl = b.Length;
1592
1593                         uint[] res = new uint [bl];
1594
1595                         ulong sum = 0;
1596
1597                         int i = 0;
1598                         for (; i < sl; i++) {
1599                                 sum = sum + a [i] + b [i];
1600                                 res [i] = (uint)sum;
1601                                 sum >>= 32;
1602                         }
1603
1604                         for (; i < bl; i++) {
1605                                 sum = sum + a [i];
1606                                 res [i] = (uint)sum;
1607                                 sum >>= 32;
1608                         }
1609
1610                         if (sum != 0) {
1611                                 res = Resize (res, bl + 1);
1612                                 res [i] = (uint)sum;
1613                         }
1614
1615                         return res;
1616                 }
1617
1618                 /*invariant a > b*/
1619                 static uint[] CoreSub (uint[] a, uint[] b)
1620                 {
1621                         int bl = a.Length;
1622                         int sl = b.Length;
1623
1624                         uint[] res = new uint [bl];
1625
1626                         ulong borrow = 0;
1627                         int i;
1628                         for (i = 0; i < sl; ++i) {
1629                                 borrow = (ulong)a [i] - b [i] - borrow;
1630
1631                                 res [i] = (uint)borrow;
1632                                 borrow = (borrow >> 32) & 0x1;
1633                         }
1634
1635                         for (; i < bl; i++) {
1636                                 borrow = (ulong)a [i] - borrow;
1637                                 res [i] = (uint)borrow;
1638                                 borrow = (borrow >> 32) & 0x1;
1639                         }
1640
1641                         //remove extra zeroes
1642                         for (i = bl - 1; i >= 0 && res [i] == 0; --i) ;
1643                         if (i < bl - 1)
1644                                 res = Resize (res, i + 1);
1645
1646             return res;
1647                 }
1648
1649
1650                 static uint[] CoreAdd (uint[] a, uint b)
1651                 {
1652                         int len = a.Length;
1653                         uint[] res = new uint [len];
1654
1655                         ulong sum = b;
1656                         int i;
1657                         for (i = 0; i < len; i++) {
1658                                 sum = sum + a [i];
1659                                 res [i] = (uint)sum;
1660                                 sum >>= 32;
1661                         }
1662
1663                         if (sum != 0) {
1664                                 res = Resize (res, len + 1);
1665                                 res [i] = (uint)sum;
1666                         }
1667
1668                         return res;
1669                 }
1670
1671                 static uint[] CoreSub (uint[] a, uint b)
1672                 {
1673                         int len = a.Length;
1674                         uint[] res = new uint [len];
1675
1676                         ulong borrow = b;
1677                         int i;
1678                         for (i = 0; i < len; i++) {
1679                                 borrow = (ulong)a [i] - borrow;
1680                                 res [i] = (uint)borrow;
1681                                 borrow = (borrow >> 32) & 0x1;
1682                         }
1683
1684                         //remove extra zeroes
1685                         for (i = len - 1; i >= 0 && res [i] == 0; --i) ;
1686                         if (i < len - 1)
1687                                 res = Resize (res, i + 1);
1688
1689             return res;
1690                 }
1691
1692                 static int CoreCompare (uint[] a, uint[] b)
1693                 {
1694                         int     al = a.Length;
1695                         int bl = b.Length;
1696
1697                         if (al > bl)
1698                                 return 1;
1699                         if (bl > al)
1700                                 return -1;
1701
1702                         for (int i = al - 1; i >= 0; --i) {
1703                                 uint ai = a [i];
1704                                 uint bi = b [i];
1705                                 if (ai > bi)
1706                                         return 1;
1707                                 if (ai < bi)
1708                                         return -1;
1709                         }
1710                         return 0;
1711                 }
1712
1713                 static int GetNormalizeShift(uint value) {
1714                         int shift = 0;
1715
1716                         if ((value & 0xFFFF0000) == 0) { value <<= 16; shift += 16; }
1717                         if ((value & 0xFF000000) == 0) { value <<= 8; shift += 8; }
1718                         if ((value & 0xF0000000) == 0) { value <<= 4; shift += 4; }
1719                         if ((value & 0xC0000000) == 0) { value <<= 2; shift += 2; }
1720                         if ((value & 0x80000000) == 0) { value <<= 1; shift += 1; }
1721
1722                         return shift;
1723                 }
1724
1725                 static void Normalize (uint[] u, int l, uint[] un, int shift)
1726                 {
1727                         uint carry = 0;
1728                         int i;
1729                         if (shift > 0) {
1730                                 int rshift = 32 - shift;
1731                                 for (i = 0; i < l; i++) {
1732                                         uint ui = u [i];
1733                                         un [i] = (ui << shift) | carry;
1734                                         carry = ui >> rshift;
1735                                 }
1736                         } else {
1737                                 for (i = 0; i < l; i++) {
1738                                         un [i] = u [i];
1739                                 }
1740                         }
1741
1742                         while (i < un.Length) {
1743                                 un [i++] = 0;
1744                         }
1745
1746                         if (carry != 0) {
1747                                 un [l] = carry;
1748                         }
1749                 }
1750
1751                 static void Unnormalize (uint[] un, out uint[] r, int shift)
1752                 {
1753                         int length = un.Length;
1754                         r = new uint [length];
1755
1756                         if (shift > 0) {
1757                                 int lshift = 32 - shift;
1758                                 uint carry = 0;
1759                                 for (int i = length - 1; i >= 0; i--) {
1760                                         uint uni = un [i];
1761                                         r [i] = (uni >> shift) | carry;
1762                                         carry = (uni << lshift);
1763                                 }
1764                         } else {
1765                                 for (int i = 0; i < length; i++) {
1766                                         r [i] = un [i];
1767                                 }
1768                         }
1769                 }
1770
1771                 const ulong Base = 0x100000000;
1772                 static void DivModUnsigned (uint[] u, uint[] v, out uint[] q, out uint[] r)
1773                 {
1774                         int m = u.Length;
1775                         int n = v.Length;
1776
1777                         if (n <= 1) {
1778                                 //  Divide by single digit
1779                                 //
1780                                 ulong rem = 0;
1781                                 uint v0 = v [0];
1782                                 q = new uint[m];
1783                                 r = new uint [1];
1784
1785                                 for (int j = m - 1; j >= 0; j--) {
1786                                         rem *= Base;
1787                                         rem += u[j];
1788
1789                                         ulong div = rem / v0;
1790                                         rem -= div * v0;
1791                                         q[j] = (uint)div;
1792                                 }
1793                                 r [0] = (uint)rem;
1794                         } else if (m >= n) {
1795                                 int shift = GetNormalizeShift (v [n - 1]);
1796
1797                                 uint[] un = new uint [m + 1];
1798                                 uint[] vn = new uint [n];
1799
1800                                 Normalize (u, m, un, shift);
1801                                 Normalize (v, n, vn, shift);
1802
1803                                 q = new uint [m - n + 1];
1804                                 r = null;
1805
1806                                 //  Main division loop
1807                                 //
1808                                 for (int j = m - n; j >= 0; j--) {
1809                                         ulong rr, qq;
1810                                         int i;
1811
1812                                         rr = Base * un [j + n] + un [j + n - 1];
1813                                         qq = rr / vn [n - 1];
1814                                         rr -= qq * vn [n - 1];
1815
1816                                         for (; ; ) {
1817                                                 // Estimate too big ?
1818                                                 //
1819                                                 if ((qq >= Base) || (qq * vn [n - 2] > (rr * Base + un [j + n - 2]))) {
1820                                                         qq--;
1821                                                         rr += (ulong)vn [n - 1];
1822                                                         if (rr < Base)
1823                                                                 continue;
1824                                                 }
1825                                                 break;
1826                                         }
1827
1828
1829                                         //  Multiply and subtract
1830                                         //
1831                                         long b = 0;
1832                                         long t = 0;
1833                                         for (i = 0; i < n; i++) {
1834                                                 ulong p = vn [i] * qq;
1835                                                 t = (long)un [i + j] - (long)(uint)p - b;
1836                                                 un [i + j] = (uint)t;
1837                                                 p >>= 32;
1838                                                 t >>= 32;
1839                                                 b = (long)p - t;
1840                                         }
1841                                         t = (long)un [j + n] - b;
1842                                         un [j + n] = (uint)t;
1843
1844                                         //  Store the calculated value
1845                                         //
1846                                         q [j] = (uint)qq;
1847
1848                                         //  Add back vn[0..n] to un[j..j+n]
1849                                         //
1850                                         if (t < 0) {
1851                                                 q [j]--;
1852                                                 ulong c = 0;
1853                                                 for (i = 0; i < n; i++) {
1854                                                         c = (ulong)vn [i] + un [j + i] + c;
1855                                                         un [j + i] = (uint)c;
1856                                                         c >>= 32;
1857                                                 }
1858                                                 c += (ulong)un [j + n];
1859                                                 un [j + n] = (uint)c;
1860                                         }
1861                                 }
1862
1863                                 Unnormalize (un, out r, shift);
1864                         } else {
1865                                 q = new uint [] { 0 };
1866                                 r = u;
1867                         }
1868                 }
1869         }
1870 }\r