2 // System.Numerics.BigInteger
4 // Rodrigo Kumpera (rkumpera@novell.com)
7 // Copyright (C) 2010 Novell, Inc (http://www.novell.com)
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:
17 // The above copyright notice and this permission notice shall be
18 // included in all copies or substantial portions of the Software.
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.
28 // A big chuck of code comes the DLR (as hosted in http://ironpython.codeplex.com),
29 // which has the following License:
31 /* ****************************************************************************
33 * Copyright (c) Microsoft Corporation.
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.
41 * You must not remove this notice, or any other, from this software.
44 * ***************************************************************************/
47 using System.Collections.Generic;
48 using System.Diagnostics;
49 using System.Diagnostics.CodeAnalysis;
50 using System.Globalization;
52 using System.Threading;
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
64 namespace System.Numerics {
65 public struct BigInteger : IComparable, IFormattable, IComparable<BigInteger>, IEquatable<BigInteger>
71 static readonly uint[] ZERO = new uint [1];
72 static readonly uint[] ONE = new uint [1] { 1 };
74 BigInteger (short sign, uint[] data)
80 public BigInteger (int value)
85 } else if (value > 0) {
87 data = new uint[] { (uint) value };
90 data = new uint[1] { (uint)-value };
94 [CLSCompliantAttribute (false)]
95 public BigInteger (uint value)
102 data = new uint [1] { value };
106 public BigInteger (long value)
111 } else if (value > 0) {
113 uint low = (uint)value;
114 uint high = (uint)(value >> 32);
116 data = new uint [high != 0 ? 2 : 1];
123 uint low = (uint)value;
124 uint high = (uint)((ulong)value >> 32);
126 data = new uint [high != 0 ? 2 : 1];
133 [CLSCompliantAttribute (false)]
134 public BigInteger (ulong value)
141 uint low = (uint)value;
142 uint high = (uint)(value >> 32);
144 data = new uint [high != 0 ? 2 : 1];
152 static bool Negative (byte[] v)
154 return ((v[7] & 0x80) != 0);
157 static ushort Exponent (byte[] v)
159 return (ushort)((((ushort)(v[7] & 0x7F)) << (ushort)4) | (((ushort)(v[6] & 0xF0)) >> 4));
162 static ulong Mantissa(byte[] v)
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));
167 return (ulong)((ulong)i1 | ((ulong)i2 << 32));
170 const int bias = 1075;
171 public BigInteger (double value)
173 if (double.IsNaN (value) || Double.IsInfinity (value))
174 throw new OverflowException ();
176 byte[] bytes = BitConverter.GetBytes (value);
177 ulong mantissa = Mantissa (bytes);
179 // 1.0 * 2**exp, we have a power of 2
180 int exponent = Exponent (bytes);
187 BigInteger res = Negative (bytes) ? MinusOne : One;
188 res = res << (exponent - 0x3ff);
189 this.sign = res.sign;
190 this.data = res.data;
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);
198 this.sign = (short) (Negative (bytes) ? -1 : 1);
199 this.data = res.data;
203 public BigInteger (float value) : this ((double)value)
207 const Int32 DecimalScaleFactorMask = 0x00FF0000;
208 const Int32 DecimalSignMask = unchecked((Int32)0x80000000);
210 public BigInteger (decimal value)
212 // First truncate to get scale to 0 and extract bits
213 int[] bits = Decimal.GetBits(Decimal.Truncate(value));
216 while (size > 0 && bits[size - 1] == 0) size--;
224 sign = (short) ((bits [3] & DecimalSignMask) != 0 ? -1 : 1);
226 data = new uint [size];
227 data [0] = (uint)bits [0];
229 data [1] = (uint)bits [1];
231 data [2] = (uint)bits [2];
234 [CLSCompliantAttribute (false)]
235 public BigInteger (byte[] value)
238 throw new ArgumentNullException ("value");
240 int len = value.Length;
242 if (len == 0 || (len == 1 && value [0] == 0)) {
248 if ((value [len - 1] & 0x80) != 0)
254 while (value [len - 1] == 0)
257 int full_words, size;
258 full_words = size = len / 4;
259 if ((len & 0x3) != 0)
262 data = new uint [size];
264 for (int i = 0; i < full_words; ++i) {
265 data [i] = (uint)value [j++] |
266 (uint)(value [j++] << 8) |
267 (uint)(value [j++] << 16) |
268 (uint)(value [j++] << 24);
272 int idx = data.Length - 1;
273 for (int i = 0; i < size; ++i)
274 data [idx] |= (uint)(value [j++] << (i * 8));
277 int full_words, size;
278 full_words = size = len / 4;
279 if ((len & 0x3) != 0)
282 data = new uint [size];
284 uint word, borrow = 1;
288 for (int i = 0; i < full_words; ++i) {
289 word = (uint)value [j++] |
290 (uint)(value [j++] << 8) |
291 (uint)(value [j++] << 16) |
292 (uint)(value [j++] << 24);
294 sub = (ulong)word - borrow;
296 borrow = (uint)(sub >> 32) & 0x1u;
304 for (int i = 0; i < size; ++i) {
305 word |= (uint)(value [j++] << (i * 8));
306 store_mask = (store_mask << 8) | 0xFF;
311 borrow = (uint)(sub >> 32) & 0x1u;
313 data [data.Length - 1] = ~word & store_mask;
315 if (borrow != 0) //FIXME I believe this can't happen, can someone write a test for it?
316 throw new Exception ("non zero final carry");
322 get { return (data [0] & 0x1) == 0; }
326 get { return sign == 1 && data.Length == 1 && data [0] == 1; }
330 //Gem from Hacker's Delight
331 //Returns the number of bits set in @x
332 static int PopulationCount (uint x)
334 x = x - ((x >> 1) & 0x55555555);
335 x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
336 x = (x + (x >> 4)) & 0x0F0F0F0F;
339 return (int)(x & 0x0000003F);
342 public bool IsPowerOfTwo {
344 bool foundBit = false;
347 //This function is pop count == 1 for positive numbers
348 for (int i = 0; i < data.Length; ++i) {
349 int p = PopulationCount (data [i]);
351 if (p > 1 || foundBit)
361 get { return sign == 0; }
368 public static BigInteger MinusOne {
369 get { return new BigInteger (-1, ONE); }
372 public static BigInteger One {
373 get { return new BigInteger (1, ONE); }
376 public static BigInteger Zero {
377 get { return new BigInteger (0, ZERO); }
380 public static explicit operator int (BigInteger value)
382 if (value.data.Length > 1)
383 throw new OverflowException ();
384 uint data = value.data [0];
386 if (value.sign == 1) {
387 if (data > (uint)int.MaxValue)
388 throw new OverflowException ();
390 } else if (value.sign == -1) {
391 if (data > 0x80000000u)
392 throw new OverflowException ();
399 [CLSCompliantAttribute (false)]
400 public static explicit operator uint (BigInteger value)
402 if (value.data.Length > 1 || value.sign == -1)
403 throw new OverflowException ();
404 return value.data [0];
407 public static explicit operator short (BigInteger value)
409 int val = (int)value;
410 if (val < short.MinValue || val > short.MaxValue)
411 throw new OverflowException ();
415 [CLSCompliantAttribute (false)]
416 public static explicit operator ushort (BigInteger value)
418 uint val = (uint)value;
419 if (val > ushort.MaxValue)
420 throw new OverflowException ();
424 public static explicit operator byte (BigInteger value)
426 uint val = (uint)value;
427 if (val > byte.MaxValue)
428 throw new OverflowException ();
432 [CLSCompliantAttribute (false)]
433 public static explicit operator sbyte (BigInteger value)
435 int val = (int)value;
436 if (val < sbyte.MinValue || val > sbyte.MaxValue)
437 throw new OverflowException ();
442 public static explicit operator long (BigInteger value)
447 if (value.data.Length > 2)
448 throw new OverflowException ();
450 uint low = value.data [0];
452 if (value.data.Length == 1) {
455 long res = (long)low;
459 uint high = value.data [1];
461 if (value.sign == 1) {
462 if (high >= 0x80000000u)
463 throw new OverflowException ();
464 return (((long)high) << 32) | low;
467 if (high > 0x80000000u)
468 throw new OverflowException ();
470 return - ((((long)high) << 32) | (long)low);
473 [CLSCompliantAttribute (false)]
474 public static explicit operator ulong (BigInteger value)
476 if (value.data.Length > 2 || value.sign == -1)
477 throw new OverflowException ();
479 uint low = value.data [0];
480 if (value.data.Length == 1)
483 uint high = value.data [1];
484 return (((ulong)high) << 32) | low;
487 public static explicit operator double (BigInteger value)
491 return double.Parse (value.ToString (),
492 System.Globalization.CultureInfo.InvariantCulture.NumberFormat);
493 } catch (OverflowException) {
494 return value.sign == -1 ? double.NegativeInfinity : double.PositiveInfinity;
498 public static explicit operator float (BigInteger value)
502 return float.Parse (value.ToString (),
503 System.Globalization.CultureInfo.InvariantCulture.NumberFormat);
504 } catch (OverflowException) {
505 return value.sign == -1 ? float.NegativeInfinity : float.PositiveInfinity;
509 public static explicit operator decimal (BigInteger value)
514 uint[] data = value.data;
516 throw new OverflowException ();
518 int lo = 0, mi = 0, hi = 0;
520 hi = (Int32)data [2];
522 mi = (Int32)data [1];
524 lo = (Int32)data [0];
526 return new Decimal(lo, mi, hi, value.sign < 0, 0);
529 public static implicit operator BigInteger (int value)
531 return new BigInteger (value);
534 [CLSCompliantAttribute (false)]
535 public static implicit operator BigInteger (uint value)
537 return new BigInteger (value);
540 public static implicit operator BigInteger (short value)
542 return new BigInteger (value);
545 [CLSCompliantAttribute (false)]
546 public static implicit operator BigInteger (ushort value)
548 return new BigInteger (value);
551 public static implicit operator BigInteger (byte value)
553 return new BigInteger (value);
556 [CLSCompliantAttribute (false)]
557 public static implicit operator BigInteger (sbyte value)
559 return new BigInteger (value);
562 public static implicit operator BigInteger (long value)
564 return new BigInteger (value);
567 [CLSCompliantAttribute (false)]
568 public static implicit operator BigInteger (ulong value)
570 return new BigInteger (value);
573 public static explicit operator BigInteger (double value)
575 return new BigInteger (value);
578 public static explicit operator BigInteger (float value)
580 return new BigInteger (value);
583 public static explicit operator BigInteger (decimal value)
585 return new BigInteger (value);
588 public static BigInteger operator+ (BigInteger left, BigInteger right)
595 if (left.sign == right.sign)
596 return new BigInteger (left.sign, CoreAdd (left.data, right.data));
598 int r = CoreCompare (left.data, right.data);
601 return new BigInteger (0, ZERO);
603 if (r > 0) //left > right
604 return new BigInteger (left.sign, CoreSub (left.data, right.data));
606 return new BigInteger (right.sign, CoreSub (right.data, left.data));
609 public static BigInteger operator- (BigInteger left, BigInteger right)
614 return new BigInteger ((short)-right.sign, right.data);
616 if (left.sign == right.sign) {
617 int r = CoreCompare (left.data, right.data);
620 return new BigInteger (0, ZERO);
622 if (r > 0) //left > right
623 return new BigInteger (left.sign, CoreSub (left.data, right.data));
625 return new BigInteger ((short)-right.sign, CoreSub (right.data, left.data));
628 return new BigInteger (left.sign, CoreAdd (left.data, right.data));
631 public static BigInteger operator* (BigInteger left, BigInteger right)
633 if (left.sign == 0 || right.sign == 0)
634 return new BigInteger (0, ZERO);
636 if (left.data [0] == 1 && left.data.Length == 1) {
639 return new BigInteger ((short)-right.sign, right.data);
642 if (right.data [0] == 1 && right.data.Length == 1) {
645 return new BigInteger ((short)-left.sign, left.data);
648 uint[] a = left.data;
649 uint[] b = right.data;
651 uint[] res = new uint [a.Length + b.Length];
653 for (int i = 0; i < a.Length; ++i) {
658 for (int j = 0; j < b.Length; ++j) {
659 carry = carry + ((ulong)ai) * b [j] + res [k];
660 res[k++] = (uint)carry;
666 res[k++] = (uint)carry;
672 for (m = res.Length - 1; m >= 0 && res [m] == 0; --m) ;
673 if (m < res.Length - 1)
674 res = Resize (res, m + 1);
676 return new BigInteger ((short)(left.sign * right.sign), res);
679 public static BigInteger operator/ (BigInteger dividend, BigInteger divisor)
681 if (divisor.sign == 0)
682 throw new DivideByZeroException ();
684 if (dividend.sign == 0)
688 uint[] remainder_value;
690 DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
693 for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
695 return new BigInteger (0, ZERO);
696 if (i < quotient.Length - 1)
697 quotient = Resize (quotient, i + 1);
699 return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
702 public static BigInteger operator% (BigInteger dividend, BigInteger divisor)
704 if (divisor.sign == 0)
705 throw new DivideByZeroException ();
707 if (dividend.sign == 0)
711 uint[] remainder_value;
713 DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
716 for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
718 return new BigInteger (0, ZERO);
720 if (i < remainder_value.Length - 1)
721 remainder_value = Resize (remainder_value, i + 1);
722 return new BigInteger (dividend.sign, remainder_value);
725 public static BigInteger operator- (BigInteger value)
729 return new BigInteger ((short)-value.sign, value.data);
732 public static BigInteger operator+ (BigInteger value)
737 public static BigInteger operator++ (BigInteger value)
739 short sign = value.sign;
740 uint[] data = value.data;
741 if (data.Length == 1) {
742 if (sign == -1 && data [0] == 1)
743 return new BigInteger (0, ZERO);
745 return new BigInteger (1, ONE);
749 data = CoreSub (data, 1);
751 data = CoreAdd (data, 1);
753 return new BigInteger (sign, data);
756 public static BigInteger operator-- (BigInteger value)
758 short sign = value.sign;
759 uint[] data = value.data;
760 if (data.Length == 1) {
761 if (sign == 1 && data [0] == 1)
762 return new BigInteger (0, ZERO);
764 return new BigInteger (-1, ONE);
768 data = CoreAdd (data, 1);
770 data = CoreSub (data, 1);
772 return new BigInteger (sign, data);
775 public static BigInteger operator& (BigInteger left, BigInteger right)
783 uint[] a = left.data;
784 uint[] b = right.data;
788 bool neg_res = (ls == rs) && (ls == -1);
790 uint[] result = new uint [Math.Max (a.Length, b.Length)];
792 ulong ac = 1, bc = 1, borrow = 1;
795 for (i = 0; i < result.Length; ++i) {
802 ac = (uint)(ac >> 32);
811 bc = (uint)(bc >> 32);
817 borrow = word - borrow;
818 word = ~(uint)borrow;
819 borrow = (uint)(borrow >> 32) & 0x1u;
825 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
827 return new BigInteger (0, ZERO);
829 if (i < result.Length - 1)
830 result = Resize (result, i + 1);
832 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
835 public static BigInteger operator| (BigInteger left, BigInteger right)
843 uint[] a = left.data;
844 uint[] b = right.data;
848 bool neg_res = (ls == -1) || (rs == -1);
850 uint[] result = new uint [Math.Max (a.Length, b.Length)];
852 ulong ac = 1, bc = 1, borrow = 1;
855 for (i = 0; i < result.Length; ++i) {
862 ac = (uint)(ac >> 32);
871 bc = (uint)(bc >> 32);
877 borrow = word - borrow;
878 word = ~(uint)borrow;
879 borrow = (uint)(borrow >> 32) & 0x1u;
885 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
887 return new BigInteger (0, ZERO);
889 if (i < result.Length - 1)
890 result = Resize (result, i + 1);
892 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
895 public static BigInteger operator^ (BigInteger left, BigInteger right)
903 uint[] a = left.data;
904 uint[] b = right.data;
908 bool neg_res = (ls == -1) ^ (rs == -1);
910 uint[] result = new uint [Math.Max (a.Length, b.Length)];
912 ulong ac = 1, bc = 1, borrow = 1;
915 for (i = 0; i < result.Length; ++i) {
922 ac = (uint)(ac >> 32);
931 bc = (uint)(bc >> 32);
937 borrow = word - borrow;
938 word = ~(uint)borrow;
939 borrow = (uint)(borrow >> 32) & 0x1u;
945 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
947 return new BigInteger (0, ZERO);
949 if (i < result.Length - 1)
950 result = Resize (result, i + 1);
952 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
955 public static BigInteger operator~ (BigInteger value)
958 return new BigInteger (-1, ONE);
960 uint[] data = value.data;
961 int sign = value.sign;
963 bool neg_res = sign == 1;
965 uint[] result = new uint [data.Length];
967 ulong carry = 1, borrow = 1;
970 for (i = 0; i < result.Length; ++i) {
971 uint word = data [i];
973 carry = ~word + carry;
975 carry = (uint)(carry >> 32);
981 borrow = word - borrow;
982 word = ~(uint)borrow;
983 borrow = (uint)(borrow >> 32) & 0x1u;
989 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
991 return new BigInteger (0, ZERO);
993 if (i < result.Length - 1)
994 result = Resize (result, i + 1);
996 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
999 //returns the 0-based index of the most significant set bit
1000 //returns 0 if no bit is set, so extra care when using it
1001 static int BitScanBackward (uint word)
1003 for (int i = 31; i >= 0; --i) {
1004 uint mask = 1u << i;
1005 if ((word & mask) == mask)
1011 public static BigInteger operator<< (BigInteger value, int shift)
1013 if (shift == 0 || value.sign == 0)
1016 return value >> -shift;
1018 uint[] data = value.data;
1019 int sign = value.sign;
1021 int topMostIdx = BitScanBackward (data [data.Length - 1]);
1022 int bits = shift - (31 - topMostIdx);
1023 int extra_words = (bits >> 5) + ((bits & 0x1F) != 0 ? 1 : 0);
1025 uint[] res = new uint [data.Length + extra_words];
1027 int idx_shift = shift >> 5;
1028 int bit_shift = shift & 0x1F;
1029 int carry_shift = 32 - bit_shift;
1031 for (int i = 0; i < data.Length; ++i) {
1032 uint word = data [i];
1033 res [i + idx_shift] |= word << bit_shift;
1034 if (i + idx_shift + 1 < res.Length)
1035 res [i + idx_shift + 1] = word >> carry_shift;
1038 return new BigInteger ((short)sign, res);
1041 public static BigInteger operator>> (BigInteger value, int shift)
1043 if (shift == 0 || value.sign == 0)
1046 return value << -shift;
1048 uint[] data = value.data;
1049 int sign = value.sign;
1051 int topMostIdx = BitScanBackward (data [data.Length - 1]);
1052 int idx_shift = shift >> 5;
1053 int bit_shift = shift & 0x1F;
1055 int extra_words = idx_shift;
1056 if (bit_shift > topMostIdx)
1058 int size = data.Length - extra_words;
1062 return new BigInteger (0, ZERO);
1063 return new BigInteger (-1, ONE);
1066 uint[] res = new uint [size];
1067 int carry_shift = 32 - bit_shift;
1069 for (int i = data.Length - 1; i >= idx_shift; --i) {
1070 uint word = data [i];
1072 if (i - idx_shift < res.Length)
1073 res [i - idx_shift] |= word >> bit_shift;
1074 if (i - idx_shift - 1 >= 0)
1075 res [i - idx_shift - 1] = word << carry_shift;
1078 //Round down instead of toward zero
1080 for (int i = 0; i < idx_shift; i++) {
1081 if (data [i] != 0u) {
1082 var tmp = new BigInteger ((short)sign, res);
1087 if (bit_shift > 0 && (data [idx_shift] << carry_shift) != 0u) {
1088 var tmp = new BigInteger ((short)sign, res);
1093 return new BigInteger ((short)sign, res);
1096 public static bool operator< (BigInteger left, BigInteger right)
1098 return Compare (left, right) < 0;
1101 public static bool operator< (BigInteger left, long right)
1103 return left.CompareTo (right) < 0;
1107 public static bool operator< (long left, BigInteger right)
1109 return right.CompareTo (left) > 0;
1113 [CLSCompliantAttribute (false)]
1114 public static bool operator< (BigInteger left, ulong right)
1116 return left.CompareTo (right) < 0;
1119 [CLSCompliantAttribute (false)]
1120 public static bool operator< (ulong left, BigInteger right)
1122 return right.CompareTo (left) > 0;
1125 public static bool operator<= (BigInteger left, BigInteger right)
1127 return Compare (left, right) <= 0;
1130 public static bool operator<= (BigInteger left, long right)
1132 return left.CompareTo (right) <= 0;
1135 public static bool operator<= (long left, BigInteger right)
1137 return right.CompareTo (left) >= 0;
1140 [CLSCompliantAttribute (false)]
1141 public static bool operator<= (BigInteger left, ulong right)
1143 return left.CompareTo (right) <= 0;
1146 [CLSCompliantAttribute (false)]
1147 public static bool operator<= (ulong left, BigInteger right)
1149 return right.CompareTo (left) >= 0;
1152 public static bool operator> (BigInteger left, BigInteger right)
1154 return Compare (left, right) > 0;
1157 public static bool operator> (BigInteger left, long right)
1159 return left.CompareTo (right) > 0;
1162 public static bool operator> (long left, BigInteger right)
1164 return right.CompareTo (left) < 0;
1167 [CLSCompliantAttribute (false)]
1168 public static bool operator> (BigInteger left, ulong right)
1170 return left.CompareTo (right) > 0;
1173 [CLSCompliantAttribute (false)]
1174 public static bool operator> (ulong left, BigInteger right)
1176 return right.CompareTo (left) < 0;
1179 public static bool operator>= (BigInteger left, BigInteger right)
1181 return Compare (left, right) >= 0;
1184 public static bool operator>= (BigInteger left, long right)
1186 return left.CompareTo (right) >= 0;
1189 public static bool operator>= (long left, BigInteger right)
1191 return right.CompareTo (left) <= 0;
1194 [CLSCompliantAttribute (false)]
1195 public static bool operator>= (BigInteger left, ulong right)
1197 return left.CompareTo (right) >= 0;
1200 [CLSCompliantAttribute (false)]
1201 public static bool operator>= (ulong left, BigInteger right)
1203 return right.CompareTo (left) <= 0;
1206 public static bool operator== (BigInteger left, BigInteger right)
1208 return Compare (left, right) == 0;
1211 public static bool operator== (BigInteger left, long right)
1213 return left.CompareTo (right) == 0;
1216 public static bool operator== (long left, BigInteger right)
1218 return right.CompareTo (left) == 0;
1221 [CLSCompliantAttribute (false)]
1222 public static bool operator== (BigInteger left, ulong right)
1224 return left.CompareTo (right) == 0;
1227 [CLSCompliantAttribute (false)]
1228 public static bool operator== (ulong left, BigInteger right)
1230 return right.CompareTo (left) == 0;
1233 public static bool operator!= (BigInteger left, BigInteger right)
1235 return Compare (left, right) != 0;
1238 public static bool operator!= (BigInteger left, long right)
1240 return left.CompareTo (right) != 0;
1243 public static bool operator!= (long left, BigInteger right)
1245 return right.CompareTo (left) != 0;
1248 [CLSCompliantAttribute (false)]
1249 public static bool operator!= (BigInteger left, ulong right)
1251 return left.CompareTo (right) != 0;
1254 [CLSCompliantAttribute (false)]
1255 public static bool operator!= (ulong left, BigInteger right)
1257 return right.CompareTo (left) != 0;
1260 public override bool Equals (object obj)
1262 if (!(obj is BigInteger))
1264 return Equals ((BigInteger)obj);
1267 public bool Equals (BigInteger other)
1269 if (sign != other.sign)
1271 if (data.Length != other.data.Length)
1273 for (int i = 0; i < data.Length; ++i) {
1274 if (data [i] != other.data [i])
1280 public bool Equals (long other)
1282 return CompareTo (other) == 0;
1285 public override string ToString ()
1287 return ToString (10, null);
1290 string ToStringWithPadding (string format, uint radix, IFormatProvider provider)
1292 if (format.Length > 1) {
1293 int precision = Convert.ToInt32(format.Substring (1), CultureInfo.InvariantCulture.NumberFormat);
1294 string baseStr = ToString (radix, provider);
1295 if (baseStr.Length < precision) {
1296 string additional = new String ('0', precision - baseStr.Length);
1297 if (baseStr[0] != '-') {
1298 return additional + baseStr;
1300 return "-" + additional + baseStr.Substring (1);
1305 return ToString (radix, provider);
1308 public string ToString (string format)
1310 return ToString (format, null);
1313 public string ToString (IFormatProvider provider)
1315 return ToString (null, provider);
1319 public string ToString (string format, IFormatProvider provider)
1321 if (format == null || format == "")
1322 return ToString (10, provider);
1324 switch (format[0]) {
1331 return ToStringWithPadding (format, 10, provider);
1334 return ToStringWithPadding (format, 16, null);
1336 throw new FormatException (string.Format ("format '{0}' not implemented", format));
1340 static uint[] MakeTwoComplement (uint[] v)
1342 uint[] res = new uint [v.Length];
1345 for (int i = 0; i < v.Length; ++i) {
1347 carry = (ulong)~word + carry;
1349 carry = (uint)(carry >> 32);
1353 uint last = res [res.Length - 1];
1354 int idx = FirstNonFFByte (last);
1356 for (int i = 1; i < idx; ++i)
1357 mask = (mask << 8) | 0xFF;
1359 res [res.Length - 1] = last & mask;
1363 string ToString (uint radix, IFormatProvider provider)
1365 const string characterSet = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
1367 if (characterSet.Length < radix)
1368 throw new ArgumentException ("charSet length less than radix", "characterSet");
1370 throw new ArgumentException ("There is no such thing as radix one notation", "radix");
1374 if (data.Length == 1 && data [0] == 1)
1375 return sign == 1 ? "1" : "-1";
1377 List<char> digits = new List<char> (1 + data.Length * 3 / 10);
1385 dt = MakeTwoComplement (dt);
1386 a = new BigInteger (1, dt);
1391 a = DivRem (a, radix, out rem);
1392 digits.Add (characterSet [(int) rem]);
1395 if (sign == -1 && radix == 10) {
1396 NumberFormatInfo info = null;
1397 if (provider != null)
1398 info = provider.GetFormat (typeof (NumberFormatInfo)) as NumberFormatInfo;
1400 string str = info.NegativeSign;
1401 for (int i = str.Length - 1; i >= 0; --i)
1402 digits.Add (str [i]);
1408 char last = digits [digits.Count - 1];
1409 if (sign == 1 && radix > 10 && (last < '0' || last > '9'))
1414 return new String (digits.ToArray ());
1417 public static BigInteger Parse (string value)
1422 if (!Parse (value, false, out result, out ex))
1428 public static bool TryParse (string value, out BigInteger result)
1431 return Parse (value, true, out result, out ex);
1434 static Exception GetFormatException ()
1436 return new FormatException ("Input string was not in the correct format");
1439 static bool ProcessTrailingWhitespace (bool tryParse, string s, int position, ref Exception exc)
1443 for (int i = position; i < len; i++){
1446 if (c != 0 && !Char.IsWhiteSpace (c)){
1448 exc = GetFormatException ();
1455 static bool Parse (string s, bool tryParse, out BigInteger result, out Exception exc)
1459 bool digits_seen = false;
1466 exc = new ArgumentNullException ("value");
1473 for (i = 0; i < len; i++){
1475 if (!Char.IsWhiteSpace (c))
1481 exc = GetFormatException ();
1485 var info = Thread.CurrentThread.CurrentCulture.NumberFormat;
1487 string negative = info.NegativeSign;
1488 string positive = info.PositiveSign;
1490 if (string.CompareOrdinal (s, i, positive, 0, positive.Length) == 0)
1491 i += positive.Length;
1492 else if (string.CompareOrdinal (s, i, negative, 0, negative.Length) == 0) {
1494 i += negative.Length;
1497 BigInteger val = Zero;
1498 for (; i < len; i++){
1506 if (c >= '0' && c <= '9'){
1507 byte d = (byte) (c - '0');
1512 } else if (!ProcessTrailingWhitespace (tryParse, s, i, ref exc))
1518 exc = GetFormatException ();
1524 else if (sign == -1)
1525 result = new BigInteger (-1, val.data);
1527 result = new BigInteger (1, val.data);
1532 public static BigInteger Min (BigInteger left, BigInteger right)
1535 int rs = right.sign;
1542 int r = CoreCompare (left.data, right.data);
1552 public static BigInteger Max (BigInteger left, BigInteger right)
1555 int rs = right.sign;
1562 int r = CoreCompare (left.data, right.data);
1571 public static BigInteger Abs (BigInteger value)
1573 return new BigInteger ((short)Math.Abs (value.sign), value.data);
1577 public static BigInteger DivRem (BigInteger dividend, BigInteger divisor, out BigInteger remainder)
1579 if (divisor.sign == 0)
1580 throw new DivideByZeroException ();
1582 if (dividend.sign == 0) {
1583 remainder = dividend;
1588 uint[] remainder_value;
1590 DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
1593 for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
1595 remainder = new BigInteger (0, ZERO);
1597 if (i < remainder_value.Length - 1)
1598 remainder_value = Resize (remainder_value, i + 1);
1599 remainder = new BigInteger (dividend.sign, remainder_value);
1602 for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
1604 return new BigInteger (0, ZERO);
1605 if (i < quotient.Length - 1)
1606 quotient = Resize (quotient, i + 1);
1608 return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
1611 public static BigInteger Pow (BigInteger value, int exponent)
1614 throw new ArgumentOutOfRangeException("exponent", "exp must be >= 0");
1620 BigInteger result = One;
1621 while (exponent != 0) {
1622 if ((exponent & 1) != 0)
1623 result = result * value;
1627 value = value * value;
1633 public static BigInteger ModPow (BigInteger value, BigInteger exponent, BigInteger modulus) {
1634 if (exponent.sign == -1)
1635 throw new ArgumentOutOfRangeException("exponent", "power must be >= 0");
1636 if (modulus.sign == 0)
1637 throw new DivideByZeroException ();
1639 BigInteger result = One % modulus;
1640 while (exponent.sign != 0) {
1641 if (!exponent.IsEven) {
1642 result = result * value;
1643 result = result % modulus;
1647 value = value * value;
1648 value = value % modulus;
1654 public static BigInteger GreatestCommonDivisor (BigInteger left, BigInteger right)
1656 if (left.data.Length == 1 && left.data [0] == 1)
1657 return new BigInteger (1, ONE);
1658 if (right.data.Length == 1 && right.data [0] == 1)
1659 return new BigInteger (1, ONE);
1665 BigInteger x = new BigInteger (1, left.data);
1666 BigInteger y = new BigInteger (1, right.data);
1670 while (x.data.Length > 1) {
1676 if (x.IsZero) return g;
1678 // TODO: should we have something here if we can convert to long?
1681 // Now we can just do it with single precision. I am using the binary gcd method,
1682 // as it should be faster.
1685 uint yy = x.data [0];
1686 uint xx = (uint)(y % yy);
1690 while (((xx | yy) & 1) == 0) {
1691 xx >>= 1; yy >>= 1; t++;
1694 while ((xx & 1) == 0) xx >>= 1;
1695 while ((yy & 1) == 0) yy >>= 1;
1697 xx = (xx - yy) >> 1;
1699 yy = (yy - xx) >> 1;
1705 /*LAMESPEC Log doesn't specify to how many ulp is has to be precise
1706 We are equilavent to MS with about 2 ULP
1708 public static double Log (BigInteger value, Double baseValue)
1710 if (value.sign == -1 || baseValue == 1.0d || baseValue == -1.0d ||
1711 baseValue == Double.NegativeInfinity || double.IsNaN (baseValue))
1714 if (baseValue == 0.0d || baseValue == Double.PositiveInfinity)
1715 return value.IsOne ? 0 : double.NaN;
1717 if (value.sign == 0)
1718 return double.NegativeInfinity;
1720 int length = value.data.Length - 1;
1722 for (int curBit = 31; curBit >= 0; curBit--) {
1723 if ((value.data [length] & (1 << curBit)) != 0) {
1724 bitCount = curBit + length * 32;
1729 long bitlen = bitCount;
1730 Double c = 0, d = 1;
1732 BigInteger testBit = One;
1733 long tempBitlen = bitlen;
1734 while (tempBitlen > Int32.MaxValue) {
1735 testBit = testBit << Int32.MaxValue;
1736 tempBitlen -= Int32.MaxValue;
1738 testBit = testBit << (int)tempBitlen;
1740 for (long curbit = bitlen; curbit >= 0; --curbit) {
1741 if ((value & testBit).sign != 0)
1744 testBit = testBit >> 1;
1746 return (System.Math.Log (c) + System.Math.Log (2) * bitlen) / System.Math.Log (baseValue);
1750 public static double Log (BigInteger value)
1752 return Log (value, Math.E);
1756 public static double Log10 (BigInteger value)
1758 return Log (value, 10);
1761 [CLSCompliantAttribute (false)]
1762 public bool Equals (ulong other)
1764 return CompareTo (other) == 0;
1767 public override int GetHashCode ()
1769 uint hash = (uint)(sign * 0x01010101u);
1771 for (int i = 0; i < data.Length; ++i)
1776 public static BigInteger Add (BigInteger left, BigInteger right)
1778 return left + right;
1781 public static BigInteger Subtract (BigInteger left, BigInteger right)
1783 return left - right;
1786 public static BigInteger Multiply (BigInteger left, BigInteger right)
1788 return left * right;
1791 public static BigInteger Divide (BigInteger dividend, BigInteger divisor)
1793 return dividend / divisor;
1796 public static BigInteger Remainder (BigInteger dividend, BigInteger divisor)
1798 return dividend % divisor;
1801 public static BigInteger Negate (BigInteger value)
1806 public int CompareTo (object obj)
1811 if (!(obj is BigInteger))
1814 return Compare (this, (BigInteger)obj);
1817 public int CompareTo (BigInteger other)
1819 return Compare (this, other);
1822 [CLSCompliantAttribute (false)]
1823 public int CompareTo (ulong other)
1828 return other == 0 ? 0 : -1;
1830 if (data.Length > 2)
1833 uint high = (uint)(other >> 32);
1834 uint low = (uint)other;
1836 return LongCompare (low, high);
1839 int LongCompare (uint low, uint high)
1842 if (data.Length > 1)
1860 public int CompareTo (long other)
1863 int rs = Math.Sign (other);
1866 return ls > rs ? 1 : -1;
1871 if (data.Length > 2)
1876 uint low = (uint)other;
1877 uint high = (uint)((ulong)other >> 32);
1879 int r = LongCompare (low, high);
1886 public static int Compare (BigInteger left, BigInteger right)
1889 int rs = right.sign;
1892 return ls > rs ? 1 : -1;
1894 int r = CoreCompare (left.data, right.data);
1901 static int TopByte (uint x)
1903 if ((x & 0xFFFF0000u) != 0) {
1904 if ((x & 0xFF000000u) != 0)
1908 if ((x & 0xFF00u) != 0)
1913 static int FirstNonFFByte (uint word)
1915 if ((word & 0xFF000000u) != 0xFF000000u)
1917 else if ((word & 0xFF0000u) != 0xFF0000u)
1919 else if ((word & 0xFF00u) != 0xFF00u)
1924 public byte[] ToByteArray ()
1927 return new byte [1];
1929 //number of bytes not counting upper word
1930 int bytes = (data.Length - 1) * 4;
1931 bool needExtraZero = false;
1933 uint topWord = data [data.Length - 1];
1936 //if the topmost bit is set we need an extra
1938 extra = TopByte (topWord);
1939 uint mask = 0x80u << ((extra - 1) * 8);
1940 if ((topWord & mask) != 0) {
1941 needExtraZero = true;
1944 extra = TopByte (topWord);
1947 byte[] res = new byte [bytes + extra + (needExtraZero ? 1 : 0) ];
1950 int end = data.Length - 1;
1951 for (int i = 0; i < end; ++i) {
1952 uint word = data [i];
1954 res [j++] = (byte)word;
1955 res [j++] = (byte)(word >> 8);
1956 res [j++] = (byte)(word >> 16);
1957 res [j++] = (byte)(word >> 24);
1959 while (extra-- > 0) {
1960 res [j++] = (byte)topWord;
1965 int end = data.Length - 1;
1967 uint carry = 1, word;
1969 for (int i = 0; i < end; ++i) {
1971 add = (ulong)~word + carry;
1973 carry = (uint)(add >> 32);
1975 res [j++] = (byte)word;
1976 res [j++] = (byte)(word >> 8);
1977 res [j++] = (byte)(word >> 16);
1978 res [j++] = (byte)(word >> 24);
1981 add = (ulong)~topWord + (carry);
1983 carry = (uint)(add >> 32);
1985 int ex = FirstNonFFByte (word);
1986 bool needExtra = (word & (1 << (ex * 8 - 1))) == 0;
1987 int to = ex + (needExtra ? 1 : 0);
1990 res = Resize (res, bytes + to);
1993 res [j++] = (byte)word;
1999 res = Resize (res, bytes + 5);
2000 res [j++] = (byte)word;
2001 res [j++] = (byte)(word >> 8);
2002 res [j++] = (byte)(word >> 16);
2003 res [j++] = (byte)(word >> 24);
2011 static byte[] Resize (byte[] v, int len)
2013 byte[] res = new byte [len];
2014 Array.Copy (v, res, Math.Min (v.Length, len));
2018 static uint[] Resize (uint[] v, int len)
2020 uint[] res = new uint [len];
2021 Array.Copy (v, res, Math.Min (v.Length, len));
2025 static uint[] CoreAdd (uint[] a, uint[] b)
2027 if (a.Length < b.Length) {
2036 uint[] res = new uint [bl];
2041 for (; i < sl; i++) {
2042 sum = sum + a [i] + b [i];
2043 res [i] = (uint)sum;
2047 for (; i < bl; i++) {
2049 res [i] = (uint)sum;
2054 res = Resize (res, bl + 1);
2055 res [i] = (uint)sum;
2062 static uint[] CoreSub (uint[] a, uint[] b)
2067 uint[] res = new uint [bl];
2071 for (i = 0; i < sl; ++i) {
2072 borrow = (ulong)a [i] - b [i] - borrow;
2074 res [i] = (uint)borrow;
2075 borrow = (borrow >> 32) & 0x1;
2078 for (; i < bl; i++) {
2079 borrow = (ulong)a [i] - borrow;
2080 res [i] = (uint)borrow;
2081 borrow = (borrow >> 32) & 0x1;
2084 //remove extra zeroes
2085 for (i = bl - 1; i >= 0 && res [i] == 0; --i) ;
2087 res = Resize (res, i + 1);
2093 static uint[] CoreAdd (uint[] a, uint b)
2096 uint[] res = new uint [len];
2100 for (i = 0; i < len; i++) {
2102 res [i] = (uint)sum;
2107 res = Resize (res, len + 1);
2108 res [i] = (uint)sum;
2114 static uint[] CoreSub (uint[] a, uint b)
2117 uint[] res = new uint [len];
2121 for (i = 0; i < len; i++) {
2122 borrow = (ulong)a [i] - borrow;
2123 res [i] = (uint)borrow;
2124 borrow = (borrow >> 32) & 0x1;
2127 //remove extra zeroes
2128 for (i = len - 1; i >= 0 && res [i] == 0; --i) ;
2130 res = Resize (res, i + 1);
2135 static int CoreCompare (uint[] a, uint[] b)
2145 for (int i = al - 1; i >= 0; --i) {
2156 static int GetNormalizeShift(uint value) {
2159 if ((value & 0xFFFF0000) == 0) { value <<= 16; shift += 16; }
2160 if ((value & 0xFF000000) == 0) { value <<= 8; shift += 8; }
2161 if ((value & 0xF0000000) == 0) { value <<= 4; shift += 4; }
2162 if ((value & 0xC0000000) == 0) { value <<= 2; shift += 2; }
2163 if ((value & 0x80000000) == 0) { value <<= 1; shift += 1; }
2168 static void Normalize (uint[] u, int l, uint[] un, int shift)
2173 int rshift = 32 - shift;
2174 for (i = 0; i < l; i++) {
2176 un [i] = (ui << shift) | carry;
2177 carry = ui >> rshift;
2180 for (i = 0; i < l; i++) {
2185 while (i < un.Length) {
2194 static void Unnormalize (uint[] un, out uint[] r, int shift)
2196 int length = un.Length;
2197 r = new uint [length];
2200 int lshift = 32 - shift;
2202 for (int i = length - 1; i >= 0; i--) {
2204 r [i] = (uni >> shift) | carry;
2205 carry = (uni << lshift);
2208 for (int i = 0; i < length; i++) {
2214 const ulong Base = 0x100000000;
2215 static void DivModUnsigned (uint[] u, uint[] v, out uint[] q, out uint[] r)
2221 // Divide by single digit
2228 for (int j = m - 1; j >= 0; j--) {
2232 ulong div = rem / v0;
2237 } else if (m >= n) {
2238 int shift = GetNormalizeShift (v [n - 1]);
2240 uint[] un = new uint [m + 1];
2241 uint[] vn = new uint [n];
2243 Normalize (u, m, un, shift);
2244 Normalize (v, n, vn, shift);
2246 q = new uint [m - n + 1];
2249 // Main division loop
2251 for (int j = m - n; j >= 0; j--) {
2255 rr = Base * un [j + n] + un [j + n - 1];
2256 qq = rr / vn [n - 1];
2257 rr -= qq * vn [n - 1];
2260 // Estimate too big ?
2262 if ((qq >= Base) || (qq * vn [n - 2] > (rr * Base + un [j + n - 2]))) {
2264 rr += (ulong)vn [n - 1];
2272 // Multiply and subtract
2276 for (i = 0; i < n; i++) {
2277 ulong p = vn [i] * qq;
2278 t = (long)un [i + j] - (long)(uint)p - b;
2279 un [i + j] = (uint)t;
2284 t = (long)un [j + n] - b;
2285 un [j + n] = (uint)t;
2287 // Store the calculated value
2291 // Add back vn[0..n] to un[j..j+n]
2296 for (i = 0; i < n; i++) {
2297 c = (ulong)vn [i] + un [j + i] + c;
2298 un [j + i] = (uint)c;
2301 c += (ulong)un [j + n];
2302 un [j + n] = (uint)c;
2306 Unnormalize (un, out r, shift);
2308 q = new uint [] { 0 };