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