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