Handle CompareTo(long) and huge numbers.
[mono.git] / 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 // A big chuck of 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;
47 using System.Collections.Generic;
48 using System.Diagnostics;
49 using System.Diagnostics.CodeAnalysis;
50 using System.Globalization;
51 using System.Text;
52 using System.Threading;
53
54 /*
55 Optimization
56         Have proper popcount function for IsPowerOfTwo
57         Use unsafe ops to avoid bounds check
58         CoreAdd could avoid some resizes by checking for equal sized array that top overflow
59         For bitwise operators, hoist the conditionals out of their main loop
60         Optimize BitScanBackward
61         Use a carry variable to make shift opts do half the number of array ops.
62         Schoolbook multiply is O(n^2), use Karatsuba /Toom-3 for large numbers
63 */
64 namespace System.Numerics {
65         public struct BigInteger : IComparable, IFormattable, IComparable<BigInteger>, IEquatable<BigInteger>
66         {
67                 //LSB on [0]
68                 readonly uint[] data;
69                 readonly short sign;
70
71                 static readonly uint[] ZERO = new uint [1];
72                 static readonly uint[] ONE = new uint [1] { 1 };
73
74                 BigInteger (short sign, uint[] data)
75                 {
76                         this.sign = sign;
77                         this.data = data;
78                 }
79
80                 public BigInteger (int value)
81                 {
82                         if (value == 0) {
83                                 sign = 0;
84                                 data = ZERO;
85                         } else if (value > 0) {
86                                 sign = 1;
87                                 data = new uint[] { (uint) value };
88                         } else {
89                                 sign = -1;
90                                 data = new uint[1] { (uint)-value };
91                         }
92                 }
93
94                 [CLSCompliantAttribute (false)]
95                 public BigInteger (uint value)
96                 {
97                         if (value == 0) {
98                                 sign = 0;
99                                 data = ZERO;
100                         } else {
101                                 sign = 1;
102                                 data = new uint [1] { value };
103                         }
104                 }
105
106                 public BigInteger (long value)
107                 {
108                         if (value == 0) {
109                                 sign = 0;
110                                 data = ZERO;
111                         } else if (value > 0) {
112                                 sign = 1;
113                                 uint low = (uint)value;
114                                 uint high = (uint)(value >> 32);
115
116                                 data = new uint [high != 0 ? 2 : 1];
117                                 data [0] = low;
118                                 if (high != 0)
119                                         data [1] = high;
120                         } else {
121                                 sign = -1;
122                                 value = -value;
123                                 uint low = (uint)value;
124                                 uint high = (uint)((ulong)value >> 32);
125
126                                 data = new uint [high != 0 ? 2 : 1];
127                                 data [0] = low;
128                                 if (high != 0)
129                                         data [1] = high;
130                         }                       
131                 }
132
133                 [CLSCompliantAttribute (false)]
134                 public BigInteger (ulong value)
135                 {
136                         if (value == 0) {
137                                 sign = 0;
138                                 data = ZERO;
139                         } else {
140                                 sign = 1;
141                                 uint low = (uint)value;
142                                 uint high = (uint)(value >> 32);
143
144                                 data = new uint [high != 0 ? 2 : 1];
145                                 data [0] = low;
146                                 if (high != 0)
147                                         data [1] = high;
148                         }
149                 }
150
151
152                 static bool Negative (byte[] v)
153                 {
154                         return ((v[7] & 0x80) != 0);
155                 }
156
157                 static ushort Exponent (byte[] v)
158                 {
159                         return (ushort)((((ushort)(v[7] & 0x7F)) << (ushort)4) | (((ushort)(v[6] & 0xF0)) >> 4));
160                 }
161
162                 static ulong Mantissa(byte[] v)
163                 {
164                         uint i1 = ((uint)v[0] | ((uint)v[1] << 8) | ((uint)v[2] << 16) | ((uint)v[3] << 24));
165                         uint i2 = ((uint)v[4] | ((uint)v[5] << 8) | ((uint)(v[6] & 0xF) << 16));
166
167                         return (ulong)((ulong)i1 | ((ulong)i2 << 32));
168                 }
169
170                 const int bias = 1075;
171                 public BigInteger (double value)
172                 {
173                         if (double.IsNaN (value) || Double.IsInfinity (value))
174                                 throw new OverflowException ();
175
176                         byte[] bytes = BitConverter.GetBytes (value);
177                         ulong mantissa = Mantissa (bytes);
178                         if (mantissa == 0) {
179                                 // 1.0 * 2**exp, we have a power of 2
180                                 int exponent = Exponent (bytes);
181                                 if (exponent == 0) {
182                                         sign = 0;
183                                         data = ZERO;
184                                         return;
185                                 }
186
187                                 BigInteger res = Negative (bytes) ? MinusOne : One;
188                                 res = res << (exponent - 0x3ff);
189                                 this.sign = res.sign;
190                                 this.data = res.data;
191                         } else {
192                                 // 1.mantissa * 2**exp
193                                 int exponent = Exponent(bytes);
194                                 mantissa |= 0x10000000000000ul;
195                                 BigInteger res = mantissa;
196                                 res = exponent > bias ? res << (exponent - bias) : res >> (bias - exponent);
197
198                                 this.sign = (short) (Negative (bytes) ? -1 : 1);
199                                 this.data = res.data;
200                         }
201                 }
202
203                 public BigInteger (float value) : this ((double)value)
204                 {
205                 }
206
207                 const Int32 DecimalScaleFactorMask = 0x00FF0000;
208                 const Int32 DecimalSignMask = unchecked((Int32)0x80000000);
209
210                 public BigInteger (decimal value)
211                 {
212                         // First truncate to get scale to 0 and extract bits
213                         int[] bits = Decimal.GetBits(Decimal.Truncate(value));
214
215                         int size = 3;
216                         while (size > 0 && bits[size - 1] == 0) size--;
217
218                         if (size == 0) {
219                                 sign = 0;
220                                 data = ZERO;
221                                 return;
222                         }
223
224                         sign = (short) ((bits [3] & DecimalSignMask) != 0 ? -1 : 1);
225
226                         data = new uint [size];
227                         data [0] = (uint)bits [0];
228                         if (size > 1)
229                                 data [1] = (uint)bits [1];
230                         if (size > 2)
231                                 data [2] = (uint)bits [2];
232                 }
233
234                 [CLSCompliantAttribute (false)]
235                 public BigInteger (byte[] value)
236                 {
237                         if (value == null)
238                                 throw new ArgumentNullException ("value");
239
240                         int len = value.Length;
241
242                         if (len == 0 || (len == 1 && value [0] == 0)) {
243                                 sign = 0;
244                                 data = ZERO;
245                                 return;
246                         }
247
248                         if ((value [len - 1] & 0x80) != 0)
249                                 sign = -1;
250                         else
251                                 sign = 1;
252
253                         if (sign == 1) {
254                                 while (value [len - 1] == 0)
255                                         --len;
256
257                                 int full_words, size;
258                                 full_words = size = len / 4;
259                                 if ((len & 0x3) != 0)
260                                         ++size;
261
262                                 data = new uint [size];
263                                 int j = 0;
264                                 for (int i = 0; i < full_words; ++i) {
265                                         data [i] =      (uint)value [j++] |
266                                                                 (uint)(value [j++] << 8) |
267                                                                 (uint)(value [j++] << 16) |
268                                                                 (uint)(value [j++] << 24);
269                                 }
270                                 size = len & 0x3;
271                                 if (size > 0) {
272                                         int idx = data.Length - 1;
273                                         for (int i = 0; i < size; ++i)
274                                                 data [idx] |= (uint)(value [j++] << (i * 8));
275                                 }
276                         } else {
277                                 int full_words, size;
278                                 full_words = size = len / 4;
279                                 if ((len & 0x3) != 0)
280                                         ++size;
281
282                                 data = new uint [size];
283
284                                 uint word, borrow = 1;
285                                 ulong sub = 0;
286                                 int j = 0;
287
288                                 for (int i = 0; i < full_words; ++i) {
289                                         word =  (uint)value [j++] |
290                                                         (uint)(value [j++] << 8) |
291                                                         (uint)(value [j++] << 16) |
292                                                         (uint)(value [j++] << 24);
293
294                                         sub = (ulong)word - borrow;
295                                         word = (uint)sub;
296                                         borrow = (uint)(sub >> 32) & 0x1u;
297                                         data [i] = ~word;
298                                 }
299                                 size = len & 0x3;
300
301                                 if (size > 0) {
302                                         word = 0;
303                                         uint store_mask = 0;
304                                         for (int i = 0; i < size; ++i) {
305                                                 word |= (uint)(value [j++] << (i * 8));
306                                                 store_mask = (store_mask << 8) | 0xFF;
307                                         }
308
309                                         sub = word - borrow;
310                                         word = (uint)sub;
311                                         borrow = (uint)(sub >> 32) & 0x1u;
312
313                                         data [data.Length - 1] = ~word & store_mask;
314                                 }
315                                 if (borrow != 0) //FIXME I believe this can't happen, can someone write a test for it?
316                                         throw new Exception ("non zero final carry");
317                         }
318
319                 }
320
321                 public bool IsEven {
322                         get { return (data [0] & 0x1) == 0; }
323                 }               
324
325                 public bool IsOne {
326                         get { return sign == 1 && data.Length == 1 && data [0] == 1; }
327                 }               
328
329
330                 //Gem from Hacker's Delight
331                 //Returns the number of bits set in @x
332                 static int PopulationCount (uint x)
333                 {
334                         x = x - ((x >> 1) & 0x55555555);
335                         x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
336                         x = (x + (x >> 4)) & 0x0F0F0F0F;
337                         x = x + (x >> 8);
338                         x = x + (x >> 16);
339                         return (int)(x & 0x0000003F);
340                 }
341
342                 public bool IsPowerOfTwo {
343                         get {
344                                 bool foundBit = false;
345                                 if (sign != 1)
346                                         return false;
347                                 //This function is pop count == 1 for positive numbers
348                                 for (int i = 0; i < data.Length; ++i) {
349                                         int p = PopulationCount (data [i]);
350                                         if (p > 0) {
351                                                 if (p > 1 || foundBit)
352                                                         return false;
353                                                 foundBit = true;
354                                         }
355                                 }
356                                 return foundBit;
357                         }
358                 }               
359
360                 public bool IsZero {
361                         get { return sign == 0; }
362                 }               
363
364                 public int Sign {
365                         get { return sign; }
366                 }
367
368                 public static BigInteger MinusOne {
369                         get { return new BigInteger (-1, ONE); }
370                 }
371
372                 public static BigInteger One {
373                         get { return new BigInteger (1, ONE); }
374                 }
375
376                 public static BigInteger Zero {
377                         get { return new BigInteger (0, ZERO); }
378                 }
379
380                 public static explicit operator int (BigInteger value)
381                 {
382                         if (value.data.Length > 1)
383                                 throw new OverflowException ();
384                         uint data = value.data [0];
385
386                         if (value.sign == 1) {
387                                 if (data > (uint)int.MaxValue)
388                                         throw new OverflowException ();
389                                 return (int)data;
390                         } else if (value.sign == -1) {
391                                 if (data > 0x80000000u)
392                                         throw new OverflowException ();
393                                 return -(int)data;
394                         }
395
396                         return 0;
397                 }
398
399                 [CLSCompliantAttribute (false)]
400                 public static explicit operator uint (BigInteger value)
401                 {
402                         if (value.data.Length > 1 || value.sign == -1)
403                                 throw new OverflowException ();
404                         return value.data [0];
405                 }
406
407                 public static explicit operator short (BigInteger value)
408                 {
409                         int val = (int)value;
410                         if (val < short.MinValue || val > short.MaxValue)
411                                 throw new OverflowException ();
412                         return (short)val;
413                 }
414
415                 [CLSCompliantAttribute (false)]
416                 public static explicit operator ushort (BigInteger value)
417                 {
418                         uint val = (uint)value;
419                         if (val > ushort.MaxValue)
420                                 throw new OverflowException ();
421                         return (ushort)val;
422                 }
423
424                 public static explicit operator byte (BigInteger value)
425                 {
426                         uint val = (uint)value;
427                         if (val > byte.MaxValue)
428                                 throw new OverflowException ();
429                         return (byte)val;
430                 }
431
432                 [CLSCompliantAttribute (false)]
433                 public static explicit operator sbyte (BigInteger value)
434                 {
435                         int val = (int)value;
436                         if (val < sbyte.MinValue || val > sbyte.MaxValue)
437                                 throw new OverflowException ();
438                         return (sbyte)val;
439                 }
440
441
442                 public static explicit operator long (BigInteger value)
443                 {
444                         if (value.sign == 0)
445                                 return 0;
446
447                         if (value.data.Length > 2)
448                                 throw new OverflowException ();
449
450                         uint low = value.data [0];
451
452                         if (value.data.Length == 1) {
453                                 if (value.sign == 1)
454                                         return (long)low;
455                                 long res = (long)low;
456                                 return -res;
457                         }
458
459                         uint high = value.data [1];
460
461                         if (value.sign == 1) {
462                                 if (high >= 0x80000000u)
463                                         throw new OverflowException ();
464                                 return (((long)high) << 32) | low;
465                         }
466
467                         if (high > 0x80000000u)
468                                 throw new OverflowException ();
469
470                         return - ((((long)high) << 32) | (long)low);
471                 }
472
473                 [CLSCompliantAttribute (false)]
474                 public static explicit operator ulong (BigInteger value)
475                 {
476                         if (value.data.Length > 2 || value.sign == -1)
477                                 throw new OverflowException ();
478
479                         uint low = value.data [0];
480                         if (value.data.Length == 1)
481                                 return low;
482
483                         uint high = value.data [1];
484                         return (((ulong)high) << 32) | low;
485                 }
486
487                 public static explicit operator double (BigInteger value)
488                 {
489                         //FIXME
490                         try {
491                     return double.Parse (value.ToString (),
492                     System.Globalization.CultureInfo.InvariantCulture.NumberFormat);
493                         } catch (OverflowException) {
494                                 return value.sign == -1 ? double.NegativeInfinity : double.PositiveInfinity;
495                         }
496         }
497
498                 public static explicit operator float (BigInteger value)
499                 {
500                         //FIXME
501                         try {
502                                 return float.Parse (value.ToString (),
503                                 System.Globalization.CultureInfo.InvariantCulture.NumberFormat);
504                         } catch (OverflowException) {
505                                 return value.sign == -1 ? float.NegativeInfinity : float.PositiveInfinity;
506                         }
507                 }
508
509                 public static explicit operator decimal (BigInteger value)
510                 {
511                         if (value.sign == 0)
512                         return Decimal.Zero;
513
514                         uint[] data = value.data;
515                         if (data.Length > 3) 
516                                 throw new OverflowException ();
517
518                         int lo = 0, mi = 0, hi = 0;
519                         if (data.Length > 2)
520                                 hi = (Int32)data [2];
521                         if (data.Length > 1)
522                                 mi = (Int32)data [1];
523                         if (data.Length > 0)
524                                 lo = (Int32)data [0];
525
526                         return new Decimal(lo, mi, hi, value.sign < 0, 0);
527                 }
528
529                 public static implicit operator BigInteger (int value)
530                 {
531                         return new BigInteger (value);
532                 }
533
534                 [CLSCompliantAttribute (false)]
535                 public static implicit operator BigInteger (uint value)
536                 {
537                         return new BigInteger (value);
538                 }
539
540                 public static implicit operator BigInteger (short value)
541                 {
542                         return new BigInteger (value);
543                 }
544
545                 [CLSCompliantAttribute (false)]
546                 public static implicit operator BigInteger (ushort value)
547                 {
548                         return new BigInteger (value);
549                 }
550
551                 public static implicit operator BigInteger (byte value)
552                 {
553                         return new BigInteger (value);
554                 }
555
556                 [CLSCompliantAttribute (false)]
557                 public static implicit operator BigInteger (sbyte value)
558                 {
559                         return new BigInteger (value);
560                 }
561
562                 public static implicit operator BigInteger (long value)
563                 {
564                         return new BigInteger (value);
565                 }
566
567                 [CLSCompliantAttribute (false)]
568                 public static implicit operator BigInteger (ulong value)
569                 {
570                         return new BigInteger (value);
571                 }
572
573                 public static explicit operator BigInteger (double value)
574                 {
575                         return new BigInteger (value);
576                 }
577
578                 public static explicit operator BigInteger (float value)
579                 {
580                         return new BigInteger (value);
581                 }
582
583                 public static explicit operator BigInteger (decimal value)
584                 {
585                         return new BigInteger (value);
586                 }
587
588                 public static BigInteger operator+ (BigInteger left, BigInteger right)
589                 {
590                         if (left.sign == 0)
591                                 return right;
592                         if (right.sign == 0)
593                                 return left;
594
595                         if (left.sign == right.sign)
596                                 return new BigInteger (left.sign, CoreAdd (left.data, right.data));
597
598                         int r = CoreCompare (left.data, right.data);
599
600                         if (r == 0)     
601                                 return new BigInteger (0, ZERO);
602
603                         if (r > 0) //left > right
604                                 return new BigInteger (left.sign, CoreSub (left.data, right.data));
605
606                         return new BigInteger (right.sign, CoreSub (right.data, left.data));
607                 }
608
609                 public static BigInteger operator- (BigInteger left, BigInteger right)
610                 {
611                         if (right.sign == 0)
612                                 return left;
613                         if (left.sign == 0)
614                                 return new BigInteger ((short)-right.sign, right.data);
615
616                         if (left.sign == right.sign) {
617                                 int r = CoreCompare (left.data, right.data);
618
619                                 if (r == 0)     
620                                         return new BigInteger (0, ZERO);
621
622                                 if (r > 0) //left > right
623                                         return new BigInteger (left.sign, CoreSub (left.data, right.data));
624
625                                 return new BigInteger ((short)-right.sign, CoreSub (right.data, left.data));
626                         }
627
628                         return new BigInteger (left.sign, CoreAdd (left.data, right.data));
629                 }
630
631                 public static BigInteger operator* (BigInteger left, BigInteger right)
632                 {
633                         if (left.sign == 0 || right.sign == 0)
634                                 return new BigInteger (0, ZERO);
635
636                         if (left.data [0] == 1 && left.data.Length == 1) {
637                                 if (left.sign == 1)
638                                         return right;
639                                 return new BigInteger ((short)-right.sign, right.data);
640                         }
641
642                         if (right.data [0] == 1 && right.data.Length == 1) {
643                                 if (right.sign == 1)
644                                         return left;
645                                 return new BigInteger ((short)-left.sign, left.data);
646                         }
647
648                         uint[] a = left.data;
649                         uint[] b = right.data;
650
651                         uint[] res = new uint [a.Length + b.Length];
652
653             for (int i = 0; i < a.Length; ++i) {
654                 uint ai = a [i];
655                 int k = i;
656
657                 ulong carry = 0;
658                 for (int j = 0; j < b.Length; ++j) {
659                     carry = carry + ((ulong)ai) * b [j] + res [k];
660                     res[k++] = (uint)carry;
661                     carry >>= 32;
662                 }
663
664                 while (carry != 0) {
665                     carry += res [k];
666                     res[k++] = (uint)carry;
667                     carry >>= 32;
668                 }
669             }
670
671                         int m;
672                         for (m = res.Length - 1; m >= 0 && res [m] == 0; --m) ;
673                         if (m < res.Length - 1)
674                                 res = Resize (res, m + 1);
675
676                         return new BigInteger ((short)(left.sign * right.sign), res);
677                 }
678
679                 public static BigInteger operator/ (BigInteger dividend, BigInteger divisor)
680                 {
681                         if (divisor.sign == 0)
682                                 throw new DivideByZeroException ();
683
684                         if (dividend.sign == 0) 
685                                 return dividend;
686
687                         uint[] quotient;
688                         uint[] remainder_value;
689
690                         DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
691
692                         int i;
693                         for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
694                         if (i == -1)
695                                 return new BigInteger (0, ZERO);
696                         if (i < quotient.Length - 1)
697                                 quotient = Resize (quotient, i + 1);
698
699                         return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
700                 }
701
702                 public static BigInteger operator% (BigInteger dividend, BigInteger divisor)
703                 {
704                         if (divisor.sign == 0)
705                                 throw new DivideByZeroException ();
706
707                         if (dividend.sign == 0)
708                                 return dividend;
709
710                         uint[] quotient;
711                         uint[] remainder_value;
712
713                         DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
714
715                         int i;
716                         for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
717                         if (i == -1)
718                                 return new BigInteger (0, ZERO);
719
720                         if (i < remainder_value.Length - 1)
721                                 remainder_value = Resize (remainder_value, i + 1);
722                         return new BigInteger (dividend.sign, remainder_value);
723                 }
724
725                 public static BigInteger operator- (BigInteger value)
726                 {
727                         if (value.sign == 0)
728                                 return value;
729                         return new BigInteger ((short)-value.sign, value.data);
730                 }
731
732                 public static BigInteger operator+ (BigInteger value)
733                 {
734                         return value;
735                 }
736
737                 public static BigInteger operator++ (BigInteger value)
738                 {
739                         short sign = value.sign;
740                         uint[] data = value.data;
741                         if (data.Length == 1) {
742                                 if (sign == -1 && data [0] == 1)
743                                         return new BigInteger (0, ZERO);
744                                 if (sign == 0)
745                                         return new BigInteger (1, ONE);
746                         }
747
748                         if (sign == -1)
749                                 data = CoreSub (data, 1);
750                         else
751                                 data = CoreAdd (data, 1);
752                 
753                         return new BigInteger (sign, data);
754                 }
755
756                 public static BigInteger operator-- (BigInteger value)
757                 {
758                         short sign = value.sign;
759                         uint[] data = value.data;
760                         if (data.Length == 1) {
761                                 if (sign == 1 && data [0] == 1)
762                                         return new BigInteger (0, ZERO);
763                                 if (sign == 0)
764                                         return new BigInteger (-1, ONE);
765                         }
766
767                         if (sign == -1)
768                                 data = CoreAdd (data, 1);
769                         else
770                                 data = CoreSub (data, 1);
771                 
772                         return new BigInteger (sign, data);
773                 }
774
775                 public static BigInteger operator& (BigInteger left, BigInteger right)
776                 {
777                         if (left.sign == 0)
778                                 return left;
779
780                         if (right.sign == 0)
781                                 return right;
782
783                         uint[] a = left.data;
784                         uint[] b = right.data;
785                         int ls = left.sign;
786                         int rs = right.sign;
787
788                         bool neg_res = (ls == rs) && (ls == -1);
789
790                         uint[] result = new uint [Math.Max (a.Length, b.Length)];
791
792                         ulong ac = 1, bc = 1, borrow = 1;
793
794                         int i;
795                         for (i = 0; i < result.Length; ++i) {
796                                 uint va = 0;
797                                 if (i < a.Length)
798                                         va = a [i];
799                                 if (ls == -1) {
800                                         ac = ~va + ac;
801                                         va = (uint)ac;
802                                         ac = (uint)(ac >> 32);
803                                 }
804
805                                 uint vb = 0;
806                                 if (i < b.Length)
807                                         vb = b [i];
808                                 if (rs == -1) {
809                                         bc = ~vb + bc;
810                                         vb = (uint)bc;
811                                         bc = (uint)(bc >> 32);
812                                 }
813
814                                 uint word = va & vb;
815
816                                 if (neg_res) {
817                                         borrow = word - borrow;
818                                         word = ~(uint)borrow;
819                                         borrow = (uint)(borrow >> 32) & 0x1u;
820                                 }
821
822                                 result [i] = word;
823                         }
824
825                         for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
826                         if (i == -1)
827                                 return new BigInteger (0, ZERO);
828         
829                         if (i < result.Length - 1)
830                                 result = Resize (result, i + 1);
831
832                         return new BigInteger (neg_res ? (short)-1 : (short)1, result);
833                 }
834
835                 public static BigInteger operator| (BigInteger left, BigInteger right)
836                 {
837                         if (left.sign == 0)
838                                 return right;
839
840                         if (right.sign == 0)
841                                 return left;
842
843                         uint[] a = left.data;
844                         uint[] b = right.data;
845                         int ls = left.sign;
846                         int rs = right.sign;
847
848                         bool neg_res = (ls == -1) || (rs == -1);
849
850                         uint[] result = new uint [Math.Max (a.Length, b.Length)];
851
852                         ulong ac = 1, bc = 1, borrow = 1;
853
854                         int i;
855                         for (i = 0; i < result.Length; ++i) {
856                                 uint va = 0;
857                                 if (i < a.Length)
858                                         va = a [i];
859                                 if (ls == -1) {
860                                         ac = ~va + ac;
861                                         va = (uint)ac;
862                                         ac = (uint)(ac >> 32);
863                                 }
864
865                                 uint vb = 0;
866                                 if (i < b.Length)
867                                         vb = b [i];
868                                 if (rs == -1) {
869                                         bc = ~vb + bc;
870                                         vb = (uint)bc;
871                                         bc = (uint)(bc >> 32);
872                                 }
873
874                                 uint word = va | vb;
875
876                                 if (neg_res) {
877                                         borrow = word - borrow;
878                                         word = ~(uint)borrow;
879                                         borrow = (uint)(borrow >> 32) & 0x1u;
880                                 }
881
882                                 result [i] = word;
883                         }
884
885                         for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
886                         if (i == -1)
887                                 return new BigInteger (0, ZERO);
888         
889                         if (i < result.Length - 1)
890                                 result = Resize (result, i + 1);
891
892                         return new BigInteger (neg_res ? (short)-1 : (short)1, result);
893                 }
894
895                 public static BigInteger operator^ (BigInteger left, BigInteger right)
896                 {
897                         if (left.sign == 0)
898                                 return right;
899
900                         if (right.sign == 0)
901                                 return left;
902
903                         uint[] a = left.data;
904                         uint[] b = right.data;
905                         int ls = left.sign;
906                         int rs = right.sign;
907
908                         bool neg_res = (ls == -1) ^ (rs == -1);
909
910                         uint[] result = new uint [Math.Max (a.Length, b.Length)];
911
912                         ulong ac = 1, bc = 1, borrow = 1;
913
914                         int i;
915                         for (i = 0; i < result.Length; ++i) {
916                                 uint va = 0;
917                                 if (i < a.Length)
918                                         va = a [i];
919                                 if (ls == -1) {
920                                         ac = ~va + ac;
921                                         va = (uint)ac;
922                                         ac = (uint)(ac >> 32);
923                                 }
924
925                                 uint vb = 0;
926                                 if (i < b.Length)
927                                         vb = b [i];
928                                 if (rs == -1) {
929                                         bc = ~vb + bc;
930                                         vb = (uint)bc;
931                                         bc = (uint)(bc >> 32);
932                                 }
933
934                                 uint word = va ^ vb;
935
936                                 if (neg_res) {
937                                         borrow = word - borrow;
938                                         word = ~(uint)borrow;
939                                         borrow = (uint)(borrow >> 32) & 0x1u;
940                                 }
941
942                                 result [i] = word;
943                         }
944
945                         for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
946                         if (i == -1)
947                                 return new BigInteger (0, ZERO);
948         
949                         if (i < result.Length - 1)
950                                 result = Resize (result, i + 1);
951
952                         return new BigInteger (neg_res ? (short)-1 : (short)1, result);
953                 }
954
955                 public static BigInteger operator~ (BigInteger value)
956                 {
957                         if (value.sign == 0)
958                                 return new BigInteger (-1, ONE);
959
960                         uint[] data = value.data;
961                         int sign = value.sign;
962
963                         bool neg_res = sign == 1;
964
965                         uint[] result = new uint [data.Length];
966
967                         ulong carry = 1, borrow = 1;
968
969                         int i;
970                         for (i = 0; i < result.Length; ++i) {
971                                 uint word = data [i];
972                                 if (sign == -1) {
973                                         carry = ~word + carry;
974                                         word = (uint)carry;
975                                         carry = (uint)(carry >> 32);
976                                 }
977
978                                 word = ~word;
979
980                                 if (neg_res) {
981                                         borrow = word - borrow;
982                                         word = ~(uint)borrow;
983                                         borrow = (uint)(borrow >> 32) & 0x1u;
984                                 }
985
986                                 result [i] = word;
987                         }
988
989                         for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
990                         if (i == -1)
991                                 return new BigInteger (0, ZERO);
992         
993                         if (i < result.Length - 1)
994                                 result = Resize (result, i + 1);
995
996                         return new BigInteger (neg_res ? (short)-1 : (short)1, result);
997                 }
998
999                 //returns the 0-based index of the most significant set bit
1000                 //returns 0 if no bit is set, so extra care when using it
1001                 static int BitScanBackward (uint word)
1002                 {
1003                         for (int i = 31; i >= 0; --i) {
1004                                 uint mask = 1u << i;
1005                                 if ((word & mask) == mask)
1006                                         return i;
1007                         }
1008                         return 0;
1009                 }
1010
1011                 public static BigInteger operator<< (BigInteger value, int shift)
1012                 {
1013                         if (shift == 0 || value.sign == 0)
1014                                 return value;
1015                         if (shift < 0)
1016                                 return value >> -shift;
1017
1018                         uint[] data = value.data;
1019                         int sign = value.sign;
1020
1021                         int topMostIdx = BitScanBackward (data [data.Length - 1]);
1022                         int bits = shift - (31 - topMostIdx);
1023                         int extra_words = (bits >> 5) + ((bits & 0x1F) != 0 ? 1 : 0);
1024
1025                         uint[] res = new uint [data.Length + extra_words];
1026
1027                         int idx_shift = shift >> 5;
1028                         int bit_shift = shift & 0x1F;
1029                         int carry_shift = 32 - bit_shift;
1030
1031                         for (int i = 0; i < data.Length; ++i) {
1032                                 uint word = data [i];
1033                                 res [i + idx_shift] |= word << bit_shift;
1034                                 if (i + idx_shift + 1 < res.Length)
1035                                         res [i + idx_shift + 1] = word >> carry_shift;
1036                         }
1037
1038                         return new BigInteger ((short)sign, res);
1039                 }
1040
1041                 public static BigInteger operator>> (BigInteger value, int shift)
1042                 {
1043                         if (shift == 0 || value.sign == 0)
1044                                 return value;
1045                         if (shift < 0)
1046                                 return value << -shift;
1047
1048                         uint[] data = value.data;
1049                         int sign = value.sign;
1050
1051                         int topMostIdx = BitScanBackward (data [data.Length - 1]);
1052                         int idx_shift = shift >> 5;
1053                         int bit_shift = shift & 0x1F;
1054
1055                         int extra_words = idx_shift;
1056                         if (bit_shift > topMostIdx)
1057                                 ++extra_words;
1058                         int size = data.Length - extra_words;
1059
1060                         if (size <= 0) {
1061                                 if (sign == 1)
1062                                         return new BigInteger (0, ZERO);
1063                                 return new BigInteger (-1, ONE);
1064                         }
1065
1066                         uint[] res = new uint [size];
1067                         int carry_shift = 32 - bit_shift;
1068
1069                         for (int i = data.Length - 1; i >= idx_shift; --i) {
1070                                 uint word = data [i];
1071
1072                                 if (i - idx_shift < res.Length)
1073                                         res [i - idx_shift] |= word >> bit_shift;
1074                                 if (i - idx_shift - 1 >= 0)
1075                                         res [i - idx_shift - 1] = word << carry_shift;
1076                         }
1077
1078                         //Round down instead of toward zero
1079                         if (sign == -1) {
1080                                 for (int i = 0; i < idx_shift; i++) {
1081                                         if (data [i] != 0u) {
1082                                                 var tmp = new BigInteger ((short)sign, res);
1083                                                 --tmp;
1084                                                 return tmp;
1085                                         }
1086                                 }
1087                                 if (bit_shift > 0 && (data [idx_shift] << carry_shift) != 0u) {
1088                                         var tmp = new BigInteger ((short)sign, res);
1089                                         --tmp;
1090                                         return tmp;
1091                                 }
1092                         }
1093                         return new BigInteger ((short)sign, res);
1094                 }
1095
1096                 public static bool operator< (BigInteger left, BigInteger right)
1097                 {
1098                         return Compare (left, right) < 0;
1099                 }
1100
1101                 public static bool operator< (BigInteger left, long right)
1102                 {
1103                         return left.CompareTo (right) < 0;
1104                 }
1105
1106
1107                 public static bool operator< (long left, BigInteger right)
1108                 {
1109                         return right.CompareTo (left) > 0;
1110                 }
1111
1112
1113                 [CLSCompliantAttribute (false)]
1114                 public static bool operator< (BigInteger left, ulong right)
1115                 {
1116                         return left.CompareTo (right) < 0;
1117                 }
1118
1119                 [CLSCompliantAttribute (false)]
1120                 public static bool operator< (ulong left, BigInteger right)
1121                 {
1122                         return right.CompareTo (left) > 0;
1123                 }
1124
1125                 public static bool operator<= (BigInteger left, BigInteger right)
1126                 {
1127                         return Compare (left, right) <= 0;
1128                 }
1129
1130                 public static bool operator<= (BigInteger left, long right)
1131                 {
1132                         return left.CompareTo (right) <= 0;
1133                 }
1134
1135                 public static bool operator<= (long left, BigInteger right)
1136                 {
1137                         return right.CompareTo (left) >= 0;
1138                 }
1139
1140                 [CLSCompliantAttribute (false)]
1141                 public static bool operator<= (BigInteger left, ulong right)
1142                 {
1143                         return left.CompareTo (right) <= 0;
1144                 }
1145
1146                 [CLSCompliantAttribute (false)]
1147                 public static bool operator<= (ulong left, BigInteger right)
1148                 {
1149                         return right.CompareTo (left) >= 0;
1150                 }
1151
1152                 public static bool operator> (BigInteger left, BigInteger right)
1153                 {
1154                         return Compare (left, right) > 0;
1155                 }
1156
1157                 public static bool operator> (BigInteger left, long right)
1158                 {
1159                         return left.CompareTo (right) > 0;
1160                 }
1161
1162                 public static bool operator> (long left, BigInteger right)
1163                 {
1164                         return right.CompareTo (left) < 0;
1165                 }
1166
1167                 [CLSCompliantAttribute (false)]
1168                 public static bool operator> (BigInteger left, ulong right)
1169                 {
1170                         return left.CompareTo (right) > 0;
1171                 }
1172
1173                 [CLSCompliantAttribute (false)]
1174                 public static bool operator> (ulong left, BigInteger right)
1175                 {
1176                         return right.CompareTo (left) < 0;
1177                 }
1178
1179                 public static bool operator>= (BigInteger left, BigInteger right)
1180                 {
1181                         return Compare (left, right) >= 0;
1182                 }
1183
1184                 public static bool operator>= (BigInteger left, long right)
1185                 {
1186                         return left.CompareTo (right) >= 0;
1187                 }
1188
1189                 public static bool operator>= (long left, BigInteger right)
1190                 {
1191                         return right.CompareTo (left) <= 0;
1192                 }
1193
1194                 [CLSCompliantAttribute (false)]
1195                 public static bool operator>= (BigInteger left, ulong right)
1196                 {
1197                         return left.CompareTo (right) >= 0;
1198                 }
1199
1200                 [CLSCompliantAttribute (false)]
1201                 public static bool operator>= (ulong left, BigInteger right)
1202                 {
1203                         return right.CompareTo (left) <= 0;
1204                 }
1205
1206                 public static bool operator== (BigInteger left, BigInteger right)
1207                 {
1208                         return Compare (left, right) == 0;
1209                 }
1210
1211                 public static bool operator== (BigInteger left, long right)
1212                 {
1213                         return left.CompareTo (right) == 0;
1214                 }
1215
1216                 public static bool operator== (long left, BigInteger right)
1217                 {
1218                         return right.CompareTo (left) == 0;
1219                 }
1220
1221                 [CLSCompliantAttribute (false)]
1222                 public static bool operator== (BigInteger left, ulong right)
1223                 {
1224                         return left.CompareTo (right) == 0;
1225                 }
1226
1227                 [CLSCompliantAttribute (false)]
1228                 public static bool operator== (ulong left, BigInteger right)
1229                 {
1230                         return right.CompareTo (left) == 0;
1231                 }
1232
1233                 public static bool operator!= (BigInteger left, BigInteger right)
1234                 {
1235                         return Compare (left, right) != 0;
1236                 }
1237
1238                 public static bool operator!= (BigInteger left, long right)
1239                 {
1240                         return left.CompareTo (right) != 0;
1241                 }
1242
1243                 public static bool operator!= (long left, BigInteger right)
1244                 {
1245                         return right.CompareTo (left) != 0;
1246                 }
1247
1248                 [CLSCompliantAttribute (false)]
1249                 public static bool operator!= (BigInteger left, ulong right)
1250                 {
1251                         return left.CompareTo (right) != 0;
1252                 }
1253
1254                 [CLSCompliantAttribute (false)]
1255                 public static bool operator!= (ulong left, BigInteger right)
1256                 {
1257                         return right.CompareTo (left) != 0;
1258                 }
1259
1260                 public override bool Equals (object obj)
1261                 {
1262                         if (!(obj is BigInteger))
1263                                 return false;
1264                         return Equals ((BigInteger)obj);
1265                 }
1266
1267                 public bool Equals (BigInteger other)
1268                 {
1269                         if (sign != other.sign)
1270                                 return false;
1271                         if (data.Length != other.data.Length)
1272                                 return false;
1273                         for (int i = 0; i < data.Length; ++i) {
1274                                 if (data [i] != other.data [i])
1275                                         return false;
1276                         }
1277                         return true;
1278                 }
1279
1280                 public bool Equals (long other)
1281                 {
1282                         return CompareTo (other) == 0;
1283                 }
1284
1285                 public override string ToString ()
1286                 {
1287                         return ToString (10, null);
1288                 }
1289
1290                 string ToStringWithPadding (string format, uint radix, IFormatProvider provider)
1291                 {
1292                         if (format.Length > 1) {
1293                                 int precision = Convert.ToInt32(format.Substring (1), CultureInfo.InvariantCulture.NumberFormat);
1294                                 string baseStr = ToString (radix, provider);
1295                                 if (baseStr.Length < precision) {
1296                                         string additional = new String ('0', precision - baseStr.Length);
1297                                         if (baseStr[0] != '-') {
1298                                                 return additional + baseStr;
1299                                         } else {
1300                                                         return "-" + additional + baseStr.Substring (1);
1301                                         }
1302                                 }
1303                                 return baseStr;
1304                         }
1305                         return ToString (radix, provider);
1306                 }
1307
1308                 public string ToString (string format)
1309                 {
1310                         return ToString (format, null);
1311                 }
1312
1313                 public string ToString (IFormatProvider provider)
1314                 {
1315                         return ToString (null, provider);
1316                 }
1317
1318
1319                 public string ToString (string format, IFormatProvider provider)
1320                 {
1321                         if (format == null || format == "")
1322                                 return ToString (10, provider);
1323
1324                         switch (format[0]) {
1325                         case 'd':
1326                         case 'D':
1327                         case 'g':
1328                         case 'G':
1329                         case 'r':
1330                         case 'R':
1331                                 return ToStringWithPadding (format, 10, provider);
1332                         case 'x':
1333                         case 'X':
1334                                 return ToStringWithPadding (format, 16, null);
1335                         default:
1336                                 throw new FormatException (string.Format ("format '{0}' not implemented", format));
1337                         }
1338                 }
1339
1340                 static uint[] MakeTwoComplement (uint[] v)
1341                 {
1342                         uint[] res = new uint [v.Length];
1343
1344                         ulong carry = 1;
1345                         for (int i = 0; i < v.Length; ++i) {
1346                                 uint word = v [i];
1347                                 carry = (ulong)~word + carry;
1348                                 word = (uint)carry;
1349                                 carry = (uint)(carry >> 32);
1350                                 res [i] = word;
1351                         }
1352
1353                         uint last = res [res.Length - 1];
1354                         int idx = FirstNonFFByte (last);
1355                         uint mask = 0xFF;
1356                         for (int i = 1; i < idx; ++i)
1357                                 mask = (mask << 8) | 0xFF;
1358
1359                         res [res.Length - 1] = last & mask;
1360                         return res;
1361                 }
1362
1363                 string ToString (uint radix, IFormatProvider provider)
1364                 {
1365                         const string characterSet = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
1366
1367                         if (characterSet.Length < radix)
1368                                 throw new ArgumentException ("charSet length less than radix", "characterSet");
1369                         if (radix == 1)
1370                                 throw new ArgumentException ("There is no such thing as radix one notation", "radix");
1371
1372                         if (sign == 0)
1373                                 return "0";
1374                         if (data.Length == 1 && data [0] == 1)
1375                                 return sign == 1 ? "1" : "-1";
1376
1377                         List<char> digits = new List<char> (1 + data.Length * 3 / 10);
1378
1379                         BigInteger a;
1380                         if (sign == 1)
1381                                 a = this;
1382                         else {
1383                                 uint[] dt = data;
1384                                 if (radix > 10)
1385                                         dt = MakeTwoComplement (dt);
1386                                 a = new BigInteger (1, dt);
1387                         }               
1388
1389                         while (a != 0) {
1390                                 BigInteger rem;
1391                                 a = DivRem (a, radix, out rem);
1392                                 digits.Add (characterSet [(int) rem]);
1393                         }
1394
1395                         if (sign == -1 && radix == 10) {
1396                                 NumberFormatInfo info = null;
1397                                 if (provider != null)
1398                                         info = provider.GetFormat (typeof (NumberFormatInfo)) as NumberFormatInfo;
1399                                 if (info != null) {
1400                                         string str = info.NegativeSign;
1401                                         for (int i = str.Length - 1; i >= 0; --i)
1402                                                 digits.Add (str [i]);
1403                                 } else {
1404                                         digits.Add ('-');
1405                                 }
1406                         }
1407
1408                         char last = digits [digits.Count - 1];
1409                         if (sign == 1 && radix > 10 && (last < '0' || last > '9'))
1410                                 digits.Add ('0');
1411                 
1412                         digits.Reverse ();
1413
1414                         return new String (digits.ToArray ());
1415                 }
1416
1417                 public static BigInteger Parse (string value)
1418                 {
1419                         Exception ex;
1420                         BigInteger result;
1421
1422                         if (!Parse (value, false, out result, out ex))
1423                                 throw ex;
1424                         return result;
1425                 }
1426
1427
1428                 public static bool TryParse (string value, out BigInteger result)
1429                 {
1430                         Exception ex;
1431                         return Parse (value, true, out result, out ex);
1432                 }
1433
1434                 static Exception GetFormatException ()
1435                 {
1436                         return new FormatException ("Input string was not in the correct format");
1437                 }
1438
1439                 static bool ProcessTrailingWhitespace (bool tryParse, string s, int position, ref Exception exc)
1440                 {
1441                         int len = s.Length;
1442                         
1443                         for (int i = position; i < len; i++){
1444                                 char c = s [i];
1445                                 
1446                                 if (c != 0 && !Char.IsWhiteSpace (c)){
1447                                         if (!tryParse)
1448                                                 exc = GetFormatException ();
1449                                         return false;
1450                                 }
1451                         }
1452                         return true;
1453                 }
1454
1455                 static bool Parse (string s, bool tryParse, out BigInteger result, out Exception exc)
1456                 {
1457                         int len;
1458                         int i, sign = 1;
1459                         bool digits_seen = false;
1460
1461                         result = Zero;
1462                         exc = null;
1463
1464                         if (s == null) {
1465                                 if (!tryParse)
1466                                         exc = new ArgumentNullException ("value");
1467                                 return false;
1468                         }
1469
1470                         len = s.Length;
1471
1472                         char c;
1473                         for (i = 0; i < len; i++){
1474                                 c = s [i];
1475                                 if (!Char.IsWhiteSpace (c))
1476                                         break;
1477                         }
1478                         
1479                         if (i == len) {
1480                                 if (!tryParse)
1481                                         exc = GetFormatException ();
1482                                 return false;
1483                         }
1484
1485                         var info = Thread.CurrentThread.CurrentCulture.NumberFormat;
1486                         
1487                         string negative = info.NegativeSign;
1488                         string positive = info.PositiveSign;
1489
1490                         if (string.CompareOrdinal (s, i, positive, 0, positive.Length) == 0)
1491                                 i += positive.Length;
1492                         else if (string.CompareOrdinal (s, i, negative, 0, negative.Length) == 0) {
1493                                 sign = -1;
1494                                 i += negative.Length;
1495                         }
1496
1497                         BigInteger val = Zero;
1498                         for (; i < len; i++){
1499                                 c = s [i];
1500
1501                                 if (c == '\0') {
1502                                         i = len;
1503                                         continue;
1504                                 }
1505
1506                                 if (c >= '0' && c <= '9'){
1507                                         byte d = (byte) (c - '0');
1508
1509                                         val = val * 10 + d;
1510
1511                                         digits_seen = true;
1512                                 } else if (!ProcessTrailingWhitespace (tryParse, s, i, ref exc))
1513                                         return false;
1514                         }
1515
1516                         if (!digits_seen) {
1517                                 if (!tryParse)
1518                                         exc = GetFormatException ();
1519                                 return false;
1520                         }
1521
1522                         if (val.sign == 0)
1523                                 result = val;
1524                         else if (sign == -1)
1525                                 result = new BigInteger (-1, val.data);
1526                         else
1527                                 result = new BigInteger (1, val.data);
1528
1529                         return true;
1530                 }
1531
1532                 public static BigInteger Min (BigInteger left, BigInteger right)
1533                 {
1534                         int ls = left.sign;
1535                         int rs = right.sign;
1536
1537                         if (ls < rs)
1538                                 return left;
1539                         if (rs < ls)
1540                                 return right;
1541
1542                         int r = CoreCompare (left.data, right.data);
1543                         if (ls == -1)
1544                                 r = -r;
1545
1546                         if (r <= 0)
1547                                 return left;
1548                         return right;
1549                 }
1550
1551
1552                 public static BigInteger Max (BigInteger left, BigInteger right)
1553                 {
1554                         int ls = left.sign;
1555                         int rs = right.sign;
1556
1557                         if (ls > rs)
1558                                 return left;
1559                         if (rs > ls)
1560                                 return right;
1561
1562                         int r = CoreCompare (left.data, right.data);
1563                         if (ls == -1)
1564                                 r = -r;
1565
1566                         if (r >= 0)
1567                                 return left;
1568                         return right;
1569                 }
1570
1571                 public static BigInteger Abs (BigInteger value)
1572                 {
1573                         return new BigInteger ((short)Math.Abs (value.sign), value.data);
1574                 }
1575
1576
1577                 public static BigInteger DivRem (BigInteger dividend, BigInteger divisor, out BigInteger remainder)
1578                 {
1579                         if (divisor.sign == 0)
1580                                 throw new DivideByZeroException ();
1581
1582                         if (dividend.sign == 0) {
1583                                 remainder = dividend;
1584                                 return dividend;
1585                         }
1586
1587                         uint[] quotient;
1588                         uint[] remainder_value;
1589
1590                         DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
1591
1592                         int i;
1593                         for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
1594                         if (i == -1) {
1595                                 remainder = new BigInteger (0, ZERO);
1596                         } else {
1597                                 if (i < remainder_value.Length - 1)
1598                                         remainder_value = Resize (remainder_value, i + 1);
1599                                 remainder = new BigInteger (dividend.sign, remainder_value);
1600                         }
1601
1602                         for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
1603                         if (i == -1)
1604                                 return new BigInteger (0, ZERO);
1605                         if (i < quotient.Length - 1)
1606                                 quotient = Resize (quotient, i + 1);
1607
1608                         return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
1609                 }
1610
1611         public static BigInteger Pow (BigInteger value, int exponent)
1612                 {
1613                         if (exponent < 0)
1614                                 throw new ArgumentOutOfRangeException("exponent", "exp must be >= 0");
1615                         if (exponent == 0)
1616                                 return One;
1617                         if (exponent == 1)
1618                                 return value;
1619
1620                         BigInteger result = One;
1621                         while (exponent != 0) {
1622                                 if ((exponent & 1) != 0)
1623                                         result = result * value;
1624                                 if (exponent == 1)
1625                                         break;
1626
1627                                 value = value * value;
1628                                 exponent >>= 1;
1629                         }
1630                         return result;
1631         }
1632
1633                 public static BigInteger ModPow (BigInteger value, BigInteger exponent, BigInteger modulus) {
1634                         if (exponent.sign == -1)
1635                                 throw new ArgumentOutOfRangeException("exponent", "power must be >= 0");
1636                         if (modulus.sign == 0)
1637                                 throw new DivideByZeroException ();
1638
1639                         BigInteger result = One % modulus;
1640                         while (exponent.sign != 0) {
1641                                 if (!exponent.IsEven) {
1642                                         result = result * value;
1643                                         result = result % modulus;
1644                                 }
1645                                 if (exponent.IsOne)
1646                                         break;
1647                                 value = value * value;
1648                                 value = value % modulus;
1649                                 exponent >>= 1;
1650                         }
1651                         return result;
1652                 }
1653
1654                 public static BigInteger GreatestCommonDivisor (BigInteger left, BigInteger right)
1655                 {
1656                         if (left.data.Length == 1 && left.data [0] == 1)
1657                                 return new BigInteger (1, ONE);
1658                         if (right.data.Length == 1 && right.data [0] == 1)
1659                                 return new BigInteger (1, ONE);
1660                         if (left.IsZero)
1661                                 return right;
1662                         if (right.IsZero)
1663                                 return left;
1664
1665                         BigInteger x = new BigInteger (1, left.data);
1666                         BigInteger y = new BigInteger (1, right.data);
1667
1668                         BigInteger g = y;
1669
1670                         while (x.data.Length > 1) {
1671                                 g = x;
1672                                 x = y % x;
1673                                 y = g;
1674
1675                         }
1676                         if (x.IsZero) return g;
1677
1678                         // TODO: should we have something here if we can convert to long?
1679
1680                         //
1681                         // Now we can just do it with single precision. I am using the binary gcd method,
1682                         // as it should be faster.
1683                         //
1684
1685                         uint yy = x.data [0];
1686                         uint xx = (uint)(y % yy);
1687
1688                         int t = 0;
1689
1690                         while (((xx | yy) & 1) == 0) {
1691                                 xx >>= 1; yy >>= 1; t++;
1692                         }
1693                         while (xx != 0) {
1694                                 while ((xx & 1) == 0) xx >>= 1;
1695                                 while ((yy & 1) == 0) yy >>= 1;
1696                                 if (xx >= yy)
1697                                         xx = (xx - yy) >> 1;
1698                                 else
1699                                         yy = (yy - xx) >> 1;
1700                         }
1701
1702                         return yy << t;
1703                 }
1704
1705                 /*LAMESPEC Log doesn't specify to how many ulp is has to be precise
1706                 We are equilavent to MS with about 2 ULP
1707                 */
1708                 public static double Log (BigInteger value, Double baseValue)
1709                 {
1710                         if (value.sign == -1 || baseValue == 1.0d || baseValue == -1.0d ||
1711                                         baseValue == Double.NegativeInfinity || double.IsNaN (baseValue))
1712                                 return double.NaN;
1713
1714                         if (baseValue == 0.0d || baseValue == Double.PositiveInfinity)
1715                                 return value.IsOne ? 0 : double.NaN;
1716         
1717                         if (value.sign == 0)
1718                                 return double.NegativeInfinity;
1719
1720                         int length = value.data.Length - 1;
1721                         int bitCount = -1;
1722                         for (int curBit = 31; curBit >= 0; curBit--) {
1723                                 if ((value.data [length] & (1 << curBit)) != 0) {
1724                                         bitCount = curBit + length * 32;
1725                                         break;
1726                                 }
1727                         }
1728
1729                         long bitlen = bitCount;
1730                         Double c = 0, d = 1;
1731
1732                         BigInteger testBit = One;
1733                         long tempBitlen = bitlen;
1734                         while (tempBitlen > Int32.MaxValue) {
1735                                 testBit = testBit << Int32.MaxValue;
1736                                 tempBitlen -= Int32.MaxValue;
1737                         }
1738                         testBit = testBit << (int)tempBitlen;
1739
1740                         for (long curbit = bitlen; curbit >= 0; --curbit) {
1741                                 if ((value & testBit).sign != 0)
1742                                         c += d;
1743                                 d *= 0.5;
1744                                 testBit = testBit >> 1;
1745                         }
1746                         return (System.Math.Log (c) + System.Math.Log (2) * bitlen) / System.Math.Log (baseValue);
1747                 }
1748
1749
1750         public static double Log (BigInteger value)
1751                 {
1752             return Log (value, Math.E);
1753         }
1754
1755
1756         public static double Log10 (BigInteger value)
1757                 {
1758             return Log (value, 10);
1759         }
1760
1761                 [CLSCompliantAttribute (false)]
1762                 public bool Equals (ulong other)
1763                 {
1764                         return CompareTo (other) == 0;
1765                 }
1766
1767                 public override int GetHashCode ()
1768                 {
1769                         uint hash = (uint)(sign * 0x01010101u);
1770
1771                         for (int i = 0; i < data.Length; ++i)
1772                                 hash ^= data [i];
1773                         return (int)hash;
1774                 }
1775
1776                 public static BigInteger Add (BigInteger left, BigInteger right)
1777                 {
1778                         return left + right;
1779                 }
1780
1781                 public static BigInteger Subtract (BigInteger left, BigInteger right)
1782                 {
1783                         return left - right;
1784                 }
1785
1786                 public static BigInteger Multiply (BigInteger left, BigInteger right)
1787                 {
1788                         return left * right;
1789                 }
1790
1791                 public static BigInteger Divide (BigInteger dividend, BigInteger divisor)
1792                 {
1793                         return dividend / divisor;
1794                 }
1795
1796                 public static BigInteger Remainder (BigInteger dividend, BigInteger divisor)
1797                 {
1798                         return dividend % divisor;
1799                 }
1800
1801                 public static BigInteger Negate (BigInteger value)
1802                 {
1803                         return - value;
1804                 }
1805
1806                 public int CompareTo (object obj)
1807                 {
1808                         if (obj == null)
1809                                 return 1;
1810                         
1811                         if (!(obj is BigInteger))
1812                                 return -1;
1813
1814                         return Compare (this, (BigInteger)obj);
1815                 }
1816
1817                 public int CompareTo (BigInteger other)
1818                 {
1819                         return Compare (this, other);
1820                 }
1821
1822                 [CLSCompliantAttribute (false)]
1823                 public int CompareTo (ulong other)
1824                 {
1825                         if (sign < 0)
1826                                 return -1;
1827                         if (sign == 0)
1828                                 return other == 0 ? 0 : -1;
1829
1830                         if (data.Length > 2)
1831                                 return 1;
1832
1833                         uint high = (uint)(other >> 32);
1834                         uint low = (uint)other;
1835
1836                         return LongCompare (low, high);
1837                 }
1838
1839                 int LongCompare (uint low, uint high)
1840                 {
1841                         uint h = 0;
1842                         if (data.Length > 1)
1843                                 h = data [1];
1844
1845                         if (h > high)
1846                                 return 1;
1847                         if (h < high)
1848                                 return -1;
1849
1850                         uint l = data [0];
1851
1852                         if (l > low)
1853                                 return 1;
1854                         if (l < low)
1855                                 return -1;
1856
1857                         return 0;
1858                 }
1859
1860                 public int CompareTo (long other)
1861                 {
1862                         int ls = sign;
1863                         int rs = Math.Sign (other);
1864
1865                         if (ls != rs)
1866                                 return ls > rs ? 1 : -1;
1867
1868                         if (ls == 0)
1869                                 return 0;
1870
1871                         if (data.Length > 2)
1872                                 return sign;
1873
1874                         if (other < 0)
1875                                 other = -other;
1876                         uint low = (uint)other;
1877                         uint high = (uint)((ulong)other >> 32);
1878
1879                         int r = LongCompare (low, high);
1880                         if (ls == -1)
1881                                 r = -r;
1882
1883                         return r;
1884                 }
1885
1886                 public static int Compare (BigInteger left, BigInteger right)
1887                 {
1888                         int ls = left.sign;
1889                         int rs = right.sign;
1890
1891                         if (ls != rs)
1892                                 return ls > rs ? 1 : -1;
1893
1894                         int r = CoreCompare (left.data, right.data);
1895                         if (ls < 0)
1896                                 r = -r;
1897                         return r;
1898                 }
1899
1900
1901                 static int TopByte (uint x)
1902                 {
1903                         if ((x & 0xFFFF0000u) != 0) {
1904                                 if ((x & 0xFF000000u) != 0)
1905                                         return 4;
1906                                 return 3;
1907                         }
1908                         if ((x & 0xFF00u) != 0)
1909                                 return 2;
1910                         return 1;       
1911                 }
1912
1913                 static int FirstNonFFByte (uint word)
1914                 {
1915                         if ((word & 0xFF000000u) != 0xFF000000u)
1916                                 return 4;
1917                         else if ((word & 0xFF0000u) != 0xFF0000u)
1918                                 return 3;
1919                         else if ((word & 0xFF00u) != 0xFF00u)
1920                                 return 2;
1921                         return 1;
1922                 }
1923
1924                 public byte[] ToByteArray ()
1925                 {
1926                         if (sign == 0)
1927                                 return new byte [1];
1928
1929                         //number of bytes not counting upper word
1930                         int bytes = (data.Length - 1) * 4;
1931                         bool needExtraZero = false;
1932
1933                         uint topWord = data [data.Length - 1];
1934                         int extra;
1935
1936                         //if the topmost bit is set we need an extra 
1937                         if (sign == 1) {
1938                                 extra = TopByte (topWord);
1939                                 uint mask = 0x80u << ((extra - 1) * 8);
1940                                 if ((topWord & mask) != 0) {
1941                                         needExtraZero = true;
1942                                 }
1943                         } else {
1944                                 extra = TopByte (topWord);
1945                         }
1946
1947                         byte[] res = new byte [bytes + extra + (needExtraZero ? 1 : 0) ];
1948                         if (sign == 1) {
1949                                 int j = 0;
1950                                 int end = data.Length - 1;
1951                                 for (int i = 0; i < end; ++i) {
1952                                         uint word = data [i];
1953
1954                                         res [j++] = (byte)word;
1955                                         res [j++] = (byte)(word >> 8);
1956                                         res [j++] = (byte)(word >> 16);
1957                                         res [j++] = (byte)(word >> 24);
1958                                 }
1959                                 while (extra-- > 0) {
1960                                         res [j++] = (byte)topWord;
1961                                         topWord >>= 8;
1962                                 }
1963                         } else {
1964                                 int j = 0;
1965                                 int end = data.Length - 1;
1966
1967                                 uint carry = 1, word;
1968                                 ulong add;
1969                                 for (int i = 0; i < end; ++i) {
1970                                         word = data [i];
1971                                         add = (ulong)~word + carry;
1972                                         word = (uint)add;
1973                                         carry = (uint)(add >> 32);
1974
1975                                         res [j++] = (byte)word;
1976                                         res [j++] = (byte)(word >> 8);
1977                                         res [j++] = (byte)(word >> 16);
1978                                         res [j++] = (byte)(word >> 24);
1979                                 }
1980
1981                                 add = (ulong)~topWord + (carry);
1982                                 word = (uint)add;
1983                                 carry = (uint)(add >> 32);
1984                                 if (carry == 0) {
1985                                         int ex = FirstNonFFByte (word);
1986                                         bool needExtra = (word & (1 << (ex * 8 - 1))) == 0;
1987                                         int to = ex + (needExtra ? 1 : 0);
1988
1989                                         if (to != extra)
1990                                                 res = Resize (res, bytes + to);
1991
1992                                         while (ex-- > 0) {
1993                                                 res [j++] = (byte)word;
1994                                                 word >>= 8;
1995                                         }
1996                                         if (needExtra)
1997                                                 res [j++] = 0xFF;
1998                                 } else {
1999                                         res = Resize (res, bytes + 5);
2000                                         res [j++] = (byte)word;
2001                                         res [j++] = (byte)(word >> 8);
2002                                         res [j++] = (byte)(word >> 16);
2003                                         res [j++] = (byte)(word >> 24);
2004                                         res [j++] = 0xFF;
2005                                 }
2006                         }
2007
2008                         return res;
2009                 }
2010
2011                 static byte[] Resize (byte[] v, int len)
2012                 {
2013                         byte[] res = new byte [len];
2014                         Array.Copy (v, res, Math.Min (v.Length, len));
2015                         return res;
2016                 }
2017
2018                 static uint[] Resize (uint[] v, int len)
2019                 {
2020                         uint[] res = new uint [len];
2021                         Array.Copy (v, res, Math.Min (v.Length, len));
2022                         return res;
2023                 }
2024
2025                 static uint[] CoreAdd (uint[] a, uint[] b)
2026                 {
2027                         if (a.Length < b.Length) {
2028                                 uint[] tmp = a;
2029                                 a = b;
2030                                 b = tmp;
2031                         }
2032
2033                         int bl = a.Length;
2034                         int sl = b.Length;
2035
2036                         uint[] res = new uint [bl];
2037
2038                         ulong sum = 0;
2039
2040                         int i = 0;
2041                         for (; i < sl; i++) {
2042                                 sum = sum + a [i] + b [i];
2043                                 res [i] = (uint)sum;
2044                                 sum >>= 32;
2045                         }
2046
2047                         for (; i < bl; i++) {
2048                                 sum = sum + a [i];
2049                                 res [i] = (uint)sum;
2050                                 sum >>= 32;
2051                         }
2052
2053                         if (sum != 0) {
2054                                 res = Resize (res, bl + 1);
2055                                 res [i] = (uint)sum;
2056                         }
2057
2058                         return res;
2059                 }
2060
2061                 /*invariant a > b*/
2062                 static uint[] CoreSub (uint[] a, uint[] b)
2063                 {
2064                         int bl = a.Length;
2065                         int sl = b.Length;
2066
2067                         uint[] res = new uint [bl];
2068
2069                         ulong borrow = 0;
2070                         int i;
2071                         for (i = 0; i < sl; ++i) {
2072                                 borrow = (ulong)a [i] - b [i] - borrow;
2073
2074                                 res [i] = (uint)borrow;
2075                                 borrow = (borrow >> 32) & 0x1;
2076                         }
2077
2078                         for (; i < bl; i++) {
2079                                 borrow = (ulong)a [i] - borrow;
2080                                 res [i] = (uint)borrow;
2081                                 borrow = (borrow >> 32) & 0x1;
2082                         }
2083
2084                         //remove extra zeroes
2085                         for (i = bl - 1; i >= 0 && res [i] == 0; --i) ;
2086                         if (i < bl - 1)
2087                                 res = Resize (res, i + 1);
2088
2089             return res;
2090                 }
2091
2092
2093                 static uint[] CoreAdd (uint[] a, uint b)
2094                 {
2095                         int len = a.Length;
2096                         uint[] res = new uint [len];
2097
2098                         ulong sum = b;
2099                         int i;
2100                         for (i = 0; i < len; i++) {
2101                                 sum = sum + a [i];
2102                                 res [i] = (uint)sum;
2103                                 sum >>= 32;
2104                         }
2105
2106                         if (sum != 0) {
2107                                 res = Resize (res, len + 1);
2108                                 res [i] = (uint)sum;
2109                         }
2110
2111                         return res;
2112                 }
2113
2114                 static uint[] CoreSub (uint[] a, uint b)
2115                 {
2116                         int len = a.Length;
2117                         uint[] res = new uint [len];
2118
2119                         ulong borrow = b;
2120                         int i;
2121                         for (i = 0; i < len; i++) {
2122                                 borrow = (ulong)a [i] - borrow;
2123                                 res [i] = (uint)borrow;
2124                                 borrow = (borrow >> 32) & 0x1;
2125                         }
2126
2127                         //remove extra zeroes
2128                         for (i = len - 1; i >= 0 && res [i] == 0; --i) ;
2129                         if (i < len - 1)
2130                                 res = Resize (res, i + 1);
2131
2132             return res;
2133                 }
2134
2135                 static int CoreCompare (uint[] a, uint[] b)
2136                 {
2137                         int     al = a.Length;
2138                         int bl = b.Length;
2139
2140                         if (al > bl)
2141                                 return 1;
2142                         if (bl > al)
2143                                 return -1;
2144
2145                         for (int i = al - 1; i >= 0; --i) {
2146                                 uint ai = a [i];
2147                                 uint bi = b [i];
2148                                 if (ai > bi)    
2149                                         return 1;
2150                                 if (ai < bi)    
2151                                         return -1;
2152                         }
2153                         return 0;
2154                 }
2155
2156                 static int GetNormalizeShift(uint value) {
2157                         int shift = 0;
2158
2159                         if ((value & 0xFFFF0000) == 0) { value <<= 16; shift += 16; }
2160                         if ((value & 0xFF000000) == 0) { value <<= 8; shift += 8; }
2161                         if ((value & 0xF0000000) == 0) { value <<= 4; shift += 4; }
2162                         if ((value & 0xC0000000) == 0) { value <<= 2; shift += 2; }
2163                         if ((value & 0x80000000) == 0) { value <<= 1; shift += 1; }
2164
2165                         return shift;
2166                 }
2167
2168                 static void Normalize (uint[] u, int l, uint[] un, int shift)
2169                 {
2170                         uint carry = 0;
2171                         int i;
2172                         if (shift > 0) {
2173                                 int rshift = 32 - shift;
2174                                 for (i = 0; i < l; i++) {
2175                                         uint ui = u [i];
2176                                         un [i] = (ui << shift) | carry;
2177                                         carry = ui >> rshift;
2178                                 }
2179                         } else {
2180                                 for (i = 0; i < l; i++) {
2181                                         un [i] = u [i];
2182                                 }
2183                         }
2184
2185                         while (i < un.Length) {
2186                                 un [i++] = 0;
2187                         }
2188
2189                         if (carry != 0) {
2190                                 un [l] = carry;
2191                         }
2192                 }
2193
2194                 static void Unnormalize (uint[] un, out uint[] r, int shift)
2195                 {
2196                         int length = un.Length;
2197                         r = new uint [length];
2198
2199                         if (shift > 0) {
2200                                 int lshift = 32 - shift;
2201                                 uint carry = 0;
2202                                 for (int i = length - 1; i >= 0; i--) {
2203                                         uint uni = un [i];
2204                                         r [i] = (uni >> shift) | carry;
2205                                         carry = (uni << lshift);
2206                                 }
2207                         } else {
2208                                 for (int i = 0; i < length; i++) {
2209                                         r [i] = un [i];
2210                                 }
2211                         }
2212                 }
2213
2214                 const ulong Base = 0x100000000;
2215                 static void DivModUnsigned (uint[] u, uint[] v, out uint[] q, out uint[] r)
2216                 {
2217                         int m = u.Length;
2218                         int n = v.Length;
2219
2220                         if (n <= 1) {
2221                                 //  Divide by single digit
2222                                 //
2223                                 ulong rem = 0;
2224                                 uint v0 = v [0];
2225                                 q = new uint[m];
2226                                 r = new uint [1];
2227
2228                                 for (int j = m - 1; j >= 0; j--) {
2229                                         rem *= Base;
2230                                         rem += u[j];
2231
2232                                         ulong div = rem / v0;
2233                                         rem -= div * v0;
2234                                         q[j] = (uint)div;
2235                                 }
2236                                 r [0] = (uint)rem;
2237                         } else if (m >= n) {
2238                                 int shift = GetNormalizeShift (v [n - 1]);
2239
2240                                 uint[] un = new uint [m + 1];
2241                                 uint[] vn = new uint [n];
2242
2243                                 Normalize (u, m, un, shift);
2244                                 Normalize (v, n, vn, shift);
2245
2246                                 q = new uint [m - n + 1];
2247                                 r = null;
2248
2249                                 //  Main division loop
2250                                 //
2251                                 for (int j = m - n; j >= 0; j--) {
2252                                         ulong rr, qq;
2253                                         int i;
2254
2255                                         rr = Base * un [j + n] + un [j + n - 1];
2256                                         qq = rr / vn [n - 1];
2257                                         rr -= qq * vn [n - 1];
2258
2259                                         for (; ; ) {
2260                                                 // Estimate too big ?
2261                                                 //
2262                                                 if ((qq >= Base) || (qq * vn [n - 2] > (rr * Base + un [j + n - 2]))) {
2263                                                         qq--;
2264                                                         rr += (ulong)vn [n - 1];
2265                                                         if (rr < Base)
2266                                                                 continue;
2267                                                 }
2268                                                 break;
2269                                         }
2270
2271
2272                                         //  Multiply and subtract
2273                                         //
2274                                         long b = 0;
2275                                         long t = 0;
2276                                         for (i = 0; i < n; i++) {
2277                                                 ulong p = vn [i] * qq;
2278                                                 t = (long)un [i + j] - (long)(uint)p - b;
2279                                                 un [i + j] = (uint)t;
2280                                                 p >>= 32;
2281                                                 t >>= 32;
2282                                                 b = (long)p - t;
2283                                         }
2284                                         t = (long)un [j + n] - b;
2285                                         un [j + n] = (uint)t;
2286
2287                                         //  Store the calculated value
2288                                         //
2289                                         q [j] = (uint)qq;
2290
2291                                         //  Add back vn[0..n] to un[j..j+n]
2292                                         //
2293                                         if (t < 0) {
2294                                                 q [j]--;
2295                                                 ulong c = 0;
2296                                                 for (i = 0; i < n; i++) {
2297                                                         c = (ulong)vn [i] + un [j + i] + c;
2298                                                         un [j + i] = (uint)c;
2299                                                         c >>= 32;
2300                                                 }
2301                                                 c += (ulong)un [j + n];
2302                                                 un [j + n] = (uint)c;
2303                                         }
2304                                 }
2305
2306                                 Unnormalize (un, out r, shift);
2307                         } else {
2308                                 q = new uint [] { 0 };
2309                                 r = u;
2310                         }
2311                 }
2312         }
2313 }