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) {
262 int full_words, size;
263 full_words = size = len / 4;
264 if ((len & 0x3) != 0)
267 data = new uint [size];
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);
277 int idx = data.Length - 1;
278 for (int i = 0; i < size; ++i)
279 data [idx] |= (uint)(value [j++] << (i * 8));
282 int full_words, size;
283 full_words = size = len / 4;
284 if ((len & 0x3) != 0)
287 data = new uint [size];
289 uint word, borrow = 1;
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);
299 sub = (ulong)word - borrow;
301 borrow = (uint)(sub >> 32) & 0x1u;
309 for (int i = 0; i < size; ++i) {
310 word |= (uint)(value [j++] << (i * 8));
311 store_mask = (store_mask << 8) | 0xFF;
316 borrow = (uint)(sub >> 32) & 0x1u;
318 data [data.Length - 1] = ~word & store_mask;
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");
327 get { return sign == 0 || (data [0] & 0x1) == 0; }
331 get { return sign == 1 && data.Length == 1 && data [0] == 1; }
335 //Gem from Hacker's Delight
336 //Returns the number of bits set in @x
337 static int PopulationCount (uint x)
339 x = x - ((x >> 1) & 0x55555555);
340 x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
341 x = (x + (x >> 4)) & 0x0F0F0F0F;
344 return (int)(x & 0x0000003F);
347 public bool IsPowerOfTwo {
349 bool foundBit = 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]);
356 if (p > 1 || foundBit)
366 get { return sign == 0; }
373 public static BigInteger MinusOne {
374 get { return new BigInteger (-1, ONE); }
377 public static BigInteger One {
378 get { return new BigInteger (1, ONE); }
381 public static BigInteger Zero {
382 get { return new BigInteger (0, ZERO); }
385 public static explicit operator int (BigInteger value)
389 if (value.data.Length > 1)
390 throw new OverflowException ();
391 uint data = value.data [0];
393 if (value.sign == 1) {
394 if (data > (uint)int.MaxValue)
395 throw new OverflowException ();
397 } else if (value.sign == -1) {
398 if (data > 0x80000000u)
399 throw new OverflowException ();
406 [CLSCompliantAttribute (false)]
407 public static explicit operator uint (BigInteger value)
411 if (value.data.Length > 1 || value.sign == -1)
412 throw new OverflowException ();
413 return value.data [0];
416 public static explicit operator short (BigInteger value)
418 int val = (int)value;
419 if (val < short.MinValue || val > short.MaxValue)
420 throw new OverflowException ();
424 [CLSCompliantAttribute (false)]
425 public static explicit operator ushort (BigInteger value)
427 uint val = (uint)value;
428 if (val > ushort.MaxValue)
429 throw new OverflowException ();
433 public static explicit operator byte (BigInteger value)
435 uint val = (uint)value;
436 if (val > byte.MaxValue)
437 throw new OverflowException ();
441 [CLSCompliantAttribute (false)]
442 public static explicit operator sbyte (BigInteger value)
444 int val = (int)value;
445 if (val < sbyte.MinValue || val > sbyte.MaxValue)
446 throw new OverflowException ();
451 public static explicit operator long (BigInteger value)
456 if (value.data.Length > 2)
457 throw new OverflowException ();
459 uint low = value.data [0];
461 if (value.data.Length == 1) {
464 long res = (long)low;
468 uint high = value.data [1];
470 if (value.sign == 1) {
471 if (high >= 0x80000000u)
472 throw new OverflowException ();
473 return (((long)high) << 32) | low;
476 if (high > 0x80000000u)
477 throw new OverflowException ();
479 return - ((((long)high) << 32) | (long)low);
482 [CLSCompliantAttribute (false)]
483 public static explicit operator ulong (BigInteger value)
487 if (value.data.Length > 2 || value.sign == -1)
488 throw new OverflowException ();
490 uint low = value.data [0];
491 if (value.data.Length == 1)
494 uint high = value.data [1];
495 return (((ulong)high) << 32) | low;
498 public static explicit operator double (BigInteger value)
502 return double.Parse (value.ToString (),
503 System.Globalization.CultureInfo.InvariantCulture.NumberFormat);
504 } catch (OverflowException) {
505 return value.sign == -1 ? double.NegativeInfinity : double.PositiveInfinity;
509 public static explicit operator float (BigInteger value)
513 return float.Parse (value.ToString (),
514 System.Globalization.CultureInfo.InvariantCulture.NumberFormat);
515 } catch (OverflowException) {
516 return value.sign == -1 ? float.NegativeInfinity : float.PositiveInfinity;
520 public static explicit operator decimal (BigInteger value)
525 uint[] data = value.data;
527 throw new OverflowException ();
529 int lo = 0, mi = 0, hi = 0;
531 hi = (Int32)data [2];
533 mi = (Int32)data [1];
535 lo = (Int32)data [0];
537 return new Decimal(lo, mi, hi, value.sign < 0, 0);
540 public static implicit operator BigInteger (int value)
542 return new BigInteger (value);
545 [CLSCompliantAttribute (false)]
546 public static implicit operator BigInteger (uint value)
548 return new BigInteger (value);
551 public static implicit operator BigInteger (short value)
553 return new BigInteger (value);
556 [CLSCompliantAttribute (false)]
557 public static implicit operator BigInteger (ushort value)
559 return new BigInteger (value);
562 public static implicit operator BigInteger (byte value)
564 return new BigInteger (value);
567 [CLSCompliantAttribute (false)]
568 public static implicit operator BigInteger (sbyte value)
570 return new BigInteger (value);
573 public static implicit operator BigInteger (long value)
575 return new BigInteger (value);
578 [CLSCompliantAttribute (false)]
579 public static implicit operator BigInteger (ulong value)
581 return new BigInteger (value);
584 public static explicit operator BigInteger (double value)
586 return new BigInteger (value);
589 public static explicit operator BigInteger (float value)
591 return new BigInteger (value);
594 public static explicit operator BigInteger (decimal value)
596 return new BigInteger (value);
599 public static BigInteger operator+ (BigInteger left, BigInteger right)
606 if (left.sign == right.sign)
607 return new BigInteger (left.sign, CoreAdd (left.data, right.data));
609 int r = CoreCompare (left.data, right.data);
612 return new BigInteger (0, ZERO);
614 if (r > 0) //left > right
615 return new BigInteger (left.sign, CoreSub (left.data, right.data));
617 return new BigInteger (right.sign, CoreSub (right.data, left.data));
620 public static BigInteger operator- (BigInteger left, BigInteger right)
625 return new BigInteger ((short)-right.sign, right.data);
627 if (left.sign == right.sign) {
628 int r = CoreCompare (left.data, right.data);
631 return new BigInteger (0, ZERO);
633 if (r > 0) //left > right
634 return new BigInteger (left.sign, CoreSub (left.data, right.data));
636 return new BigInteger ((short)-right.sign, CoreSub (right.data, left.data));
639 return new BigInteger (left.sign, CoreAdd (left.data, right.data));
642 public static BigInteger operator* (BigInteger left, BigInteger right)
644 if (left.sign == 0 || right.sign == 0)
645 return new BigInteger (0, ZERO);
647 if (left.data [0] == 1 && left.data.Length == 1) {
650 return new BigInteger ((short)-right.sign, right.data);
653 if (right.data [0] == 1 && right.data.Length == 1) {
656 return new BigInteger ((short)-left.sign, left.data);
659 uint[] a = left.data;
660 uint[] b = right.data;
662 uint[] res = new uint [a.Length + b.Length];
664 for (int i = 0; i < a.Length; ++i) {
669 for (int j = 0; j < b.Length; ++j) {
670 carry = carry + ((ulong)ai) * b [j] + res [k];
671 res[k++] = (uint)carry;
677 res[k++] = (uint)carry;
683 for (m = res.Length - 1; m >= 0 && res [m] == 0; --m) ;
684 if (m < res.Length - 1)
685 res = Resize (res, m + 1);
687 return new BigInteger ((short)(left.sign * right.sign), res);
690 public static BigInteger operator/ (BigInteger dividend, BigInteger divisor)
692 if (divisor.sign == 0)
693 throw new DivideByZeroException ();
695 if (dividend.sign == 0)
699 uint[] remainder_value;
701 DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
704 for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
706 return new BigInteger (0, ZERO);
707 if (i < quotient.Length - 1)
708 quotient = Resize (quotient, i + 1);
710 return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
713 public static BigInteger operator% (BigInteger dividend, BigInteger divisor)
715 if (divisor.sign == 0)
716 throw new DivideByZeroException ();
718 if (dividend.sign == 0)
722 uint[] remainder_value;
724 DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
727 for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
729 return new BigInteger (0, ZERO);
731 if (i < remainder_value.Length - 1)
732 remainder_value = Resize (remainder_value, i + 1);
733 return new BigInteger (dividend.sign, remainder_value);
736 public static BigInteger operator- (BigInteger value)
740 return new BigInteger ((short)-value.sign, value.data);
743 public static BigInteger operator+ (BigInteger value)
748 public static BigInteger operator++ (BigInteger value)
753 short sign = value.sign;
754 uint[] data = value.data;
755 if (data.Length == 1) {
756 if (sign == -1 && data [0] == 1)
757 return new BigInteger (0, ZERO);
759 return new BigInteger (1, ONE);
763 data = CoreSub (data, 1);
765 data = CoreAdd (data, 1);
767 return new BigInteger (sign, data);
770 public static BigInteger operator-- (BigInteger value)
775 short sign = value.sign;
776 uint[] data = value.data;
777 if (data.Length == 1) {
778 if (sign == 1 && data [0] == 1)
779 return new BigInteger (0, ZERO);
781 return new BigInteger (-1, ONE);
785 data = CoreAdd (data, 1);
787 data = CoreSub (data, 1);
789 return new BigInteger (sign, data);
792 public static BigInteger operator& (BigInteger left, BigInteger right)
800 uint[] a = left.data;
801 uint[] b = right.data;
805 bool neg_res = (ls == rs) && (ls == -1);
807 uint[] result = new uint [Math.Max (a.Length, b.Length)];
809 ulong ac = 1, bc = 1, borrow = 1;
812 for (i = 0; i < result.Length; ++i) {
819 ac = (uint)(ac >> 32);
828 bc = (uint)(bc >> 32);
834 borrow = word - borrow;
835 word = ~(uint)borrow;
836 borrow = (uint)(borrow >> 32) & 0x1u;
842 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
844 return new BigInteger (0, ZERO);
846 if (i < result.Length - 1)
847 result = Resize (result, i + 1);
849 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
852 public static BigInteger operator| (BigInteger left, BigInteger right)
860 uint[] a = left.data;
861 uint[] b = right.data;
865 bool neg_res = (ls == -1) || (rs == -1);
867 uint[] result = new uint [Math.Max (a.Length, b.Length)];
869 ulong ac = 1, bc = 1, borrow = 1;
872 for (i = 0; i < result.Length; ++i) {
879 ac = (uint)(ac >> 32);
888 bc = (uint)(bc >> 32);
894 borrow = word - borrow;
895 word = ~(uint)borrow;
896 borrow = (uint)(borrow >> 32) & 0x1u;
902 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
904 return new BigInteger (0, ZERO);
906 if (i < result.Length - 1)
907 result = Resize (result, i + 1);
909 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
912 public static BigInteger operator^ (BigInteger left, BigInteger right)
920 uint[] a = left.data;
921 uint[] b = right.data;
925 bool neg_res = (ls == -1) ^ (rs == -1);
927 uint[] result = new uint [Math.Max (a.Length, b.Length)];
929 ulong ac = 1, bc = 1, borrow = 1;
932 for (i = 0; i < result.Length; ++i) {
939 ac = (uint)(ac >> 32);
948 bc = (uint)(bc >> 32);
954 borrow = word - borrow;
955 word = ~(uint)borrow;
956 borrow = (uint)(borrow >> 32) & 0x1u;
962 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
964 return new BigInteger (0, ZERO);
966 if (i < result.Length - 1)
967 result = Resize (result, i + 1);
969 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
972 public static BigInteger operator~ (BigInteger value)
975 return new BigInteger (-1, ONE);
977 uint[] data = value.data;
978 int sign = value.sign;
980 bool neg_res = sign == 1;
982 uint[] result = new uint [data.Length];
984 ulong carry = 1, borrow = 1;
987 for (i = 0; i < result.Length; ++i) {
988 uint word = data [i];
990 carry = ~word + carry;
992 carry = (uint)(carry >> 32);
998 borrow = word - borrow;
999 word = ~(uint)borrow;
1000 borrow = (uint)(borrow >> 32) & 0x1u;
1006 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
1008 return new BigInteger (0, ZERO);
1010 if (i < result.Length - 1)
1011 result = Resize (result, i + 1);
1013 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
1016 //returns the 0-based index of the most significant set bit
1017 //returns 0 if no bit is set, so extra care when using it
1018 static int BitScanBackward (uint word)
1020 for (int i = 31; i >= 0; --i) {
1021 uint mask = 1u << i;
1022 if ((word & mask) == mask)
1028 public static BigInteger operator<< (BigInteger value, int shift)
1030 if (shift == 0 || value.sign == 0)
1033 return value >> -shift;
1035 uint[] data = value.data;
1036 int sign = value.sign;
1038 int topMostIdx = BitScanBackward (data [data.Length - 1]);
1039 int bits = shift - (31 - topMostIdx);
1040 int extra_words = (bits >> 5) + ((bits & 0x1F) != 0 ? 1 : 0);
1042 uint[] res = new uint [data.Length + extra_words];
1044 int idx_shift = shift >> 5;
1045 int bit_shift = shift & 0x1F;
1046 int carry_shift = 32 - bit_shift;
1048 if (carry_shift == 32) {
1049 for (int i = 0; i < data.Length; ++i) {
1050 uint word = data [i];
1051 res [i + idx_shift] |= word << bit_shift;
1054 for (int i = 0; i < data.Length; ++i) {
1055 uint word = data [i];
1056 res [i + idx_shift] |= word << bit_shift;
1057 if (i + idx_shift + 1 < res.Length)
1058 res [i + idx_shift + 1] = word >> carry_shift;
1062 return new BigInteger ((short)sign, res);
1065 public static BigInteger operator>> (BigInteger value, int shift)
1067 if (shift == 0 || value.sign == 0)
1070 return value << -shift;
1072 uint[] data = value.data;
1073 int sign = value.sign;
1075 int topMostIdx = BitScanBackward (data [data.Length - 1]);
1076 int idx_shift = shift >> 5;
1077 int bit_shift = shift & 0x1F;
1079 int extra_words = idx_shift;
1080 if (bit_shift > topMostIdx)
1082 int size = data.Length - extra_words;
1086 return new BigInteger (0, ZERO);
1087 return new BigInteger (-1, ONE);
1090 uint[] res = new uint [size];
1091 int carry_shift = 32 - bit_shift;
1093 if (carry_shift == 32) {
1094 for (int i = data.Length - 1; i >= idx_shift; --i) {
1095 uint word = data [i];
1097 if (i - idx_shift < res.Length)
1098 res [i - idx_shift] |= word >> bit_shift;
1101 for (int i = data.Length - 1; i >= idx_shift; --i) {
1102 uint word = data [i];
1104 if (i - idx_shift < res.Length)
1105 res [i - idx_shift] |= word >> bit_shift;
1106 if (i - idx_shift - 1 >= 0)
1107 res [i - idx_shift - 1] = word << carry_shift;
1112 //Round down instead of toward zero
1114 for (int i = 0; i < idx_shift; i++) {
1115 if (data [i] != 0u) {
1116 var tmp = new BigInteger ((short)sign, res);
1121 if (bit_shift > 0 && (data [idx_shift] << carry_shift) != 0u) {
1122 var tmp = new BigInteger ((short)sign, res);
1127 return new BigInteger ((short)sign, res);
1130 public static bool operator< (BigInteger left, BigInteger right)
1132 return Compare (left, right) < 0;
1135 public static bool operator< (BigInteger left, long right)
1137 return left.CompareTo (right) < 0;
1141 public static bool operator< (long left, BigInteger right)
1143 return right.CompareTo (left) > 0;
1147 [CLSCompliantAttribute (false)]
1148 public static bool operator< (BigInteger left, ulong right)
1150 return left.CompareTo (right) < 0;
1153 [CLSCompliantAttribute (false)]
1154 public static bool operator< (ulong left, BigInteger right)
1156 return right.CompareTo (left) > 0;
1159 public static bool operator<= (BigInteger left, BigInteger right)
1161 return Compare (left, right) <= 0;
1164 public static bool operator<= (BigInteger left, long right)
1166 return left.CompareTo (right) <= 0;
1169 public static bool operator<= (long left, BigInteger right)
1171 return right.CompareTo (left) >= 0;
1174 [CLSCompliantAttribute (false)]
1175 public static bool operator<= (BigInteger left, ulong right)
1177 return left.CompareTo (right) <= 0;
1180 [CLSCompliantAttribute (false)]
1181 public static bool operator<= (ulong left, BigInteger right)
1183 return right.CompareTo (left) >= 0;
1186 public static bool operator> (BigInteger left, BigInteger right)
1188 return Compare (left, right) > 0;
1191 public static bool operator> (BigInteger left, long right)
1193 return left.CompareTo (right) > 0;
1196 public static bool operator> (long left, BigInteger right)
1198 return right.CompareTo (left) < 0;
1201 [CLSCompliantAttribute (false)]
1202 public static bool operator> (BigInteger left, ulong right)
1204 return left.CompareTo (right) > 0;
1207 [CLSCompliantAttribute (false)]
1208 public static bool operator> (ulong left, BigInteger right)
1210 return right.CompareTo (left) < 0;
1213 public static bool operator>= (BigInteger left, BigInteger right)
1215 return Compare (left, right) >= 0;
1218 public static bool operator>= (BigInteger left, long right)
1220 return left.CompareTo (right) >= 0;
1223 public static bool operator>= (long left, BigInteger right)
1225 return right.CompareTo (left) <= 0;
1228 [CLSCompliantAttribute (false)]
1229 public static bool operator>= (BigInteger left, ulong right)
1231 return left.CompareTo (right) >= 0;
1234 [CLSCompliantAttribute (false)]
1235 public static bool operator>= (ulong left, BigInteger right)
1237 return right.CompareTo (left) <= 0;
1240 public static bool operator== (BigInteger left, BigInteger right)
1242 return Compare (left, right) == 0;
1245 public static bool operator== (BigInteger left, long right)
1247 return left.CompareTo (right) == 0;
1250 public static bool operator== (long left, BigInteger right)
1252 return right.CompareTo (left) == 0;
1255 [CLSCompliantAttribute (false)]
1256 public static bool operator== (BigInteger left, ulong right)
1258 return left.CompareTo (right) == 0;
1261 [CLSCompliantAttribute (false)]
1262 public static bool operator== (ulong left, BigInteger right)
1264 return right.CompareTo (left) == 0;
1267 public static bool operator!= (BigInteger left, BigInteger right)
1269 return Compare (left, right) != 0;
1272 public static bool operator!= (BigInteger left, long right)
1274 return left.CompareTo (right) != 0;
1277 public static bool operator!= (long left, BigInteger right)
1279 return right.CompareTo (left) != 0;
1282 [CLSCompliantAttribute (false)]
1283 public static bool operator!= (BigInteger left, ulong right)
1285 return left.CompareTo (right) != 0;
1288 [CLSCompliantAttribute (false)]
1289 public static bool operator!= (ulong left, BigInteger right)
1291 return right.CompareTo (left) != 0;
1294 public override bool Equals (object obj)
1296 if (!(obj is BigInteger))
1298 return Equals ((BigInteger)obj);
1301 public bool Equals (BigInteger other)
1303 if (sign != other.sign)
1306 int alen = data != null ? data.Length : 0;
1307 int blen = other.data != null ? other.data.Length : 0;
1311 for (int i = 0; i < alen; ++i) {
1312 if (data [i] != other.data [i])
1318 public bool Equals (long other)
1320 return CompareTo (other) == 0;
1323 public override string ToString ()
1325 return ToString (10, null);
1328 string ToStringWithPadding (string format, uint radix, IFormatProvider provider)
1330 if (format.Length > 1) {
1331 int precision = Convert.ToInt32(format.Substring (1), CultureInfo.InvariantCulture.NumberFormat);
1332 string baseStr = ToString (radix, provider);
1333 if (baseStr.Length < precision) {
1334 string additional = new String ('0', precision - baseStr.Length);
1335 if (baseStr[0] != '-') {
1336 return additional + baseStr;
1338 return "-" + additional + baseStr.Substring (1);
1343 return ToString (radix, provider);
1346 public string ToString (string format)
1348 return ToString (format, null);
1351 public string ToString (IFormatProvider provider)
1353 return ToString (null, provider);
1357 public string ToString (string format, IFormatProvider provider)
1359 if (format == null || format == "")
1360 return ToString (10, provider);
1362 switch (format[0]) {
1369 return ToStringWithPadding (format, 10, provider);
1372 return ToStringWithPadding (format, 16, null);
1374 throw new FormatException (string.Format ("format '{0}' not implemented", format));
1378 static uint[] MakeTwoComplement (uint[] v)
1380 uint[] res = new uint [v.Length];
1383 for (int i = 0; i < v.Length; ++i) {
1385 carry = (ulong)~word + carry;
1387 carry = (uint)(carry >> 32);
1391 uint last = res [res.Length - 1];
1392 int idx = FirstNonFFByte (last);
1394 for (int i = 1; i < idx; ++i)
1395 mask = (mask << 8) | 0xFF;
1397 res [res.Length - 1] = last & mask;
1401 string ToString (uint radix, IFormatProvider provider)
1403 const string characterSet = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
1405 if (characterSet.Length < radix)
1406 throw new ArgumentException ("charSet length less than radix", "characterSet");
1408 throw new ArgumentException ("There is no such thing as radix one notation", "radix");
1412 if (data.Length == 1 && data [0] == 1)
1413 return sign == 1 ? "1" : "-1";
1415 List<char> digits = new List<char> (1 + data.Length * 3 / 10);
1423 dt = MakeTwoComplement (dt);
1424 a = new BigInteger (1, dt);
1429 a = DivRem (a, radix, out rem);
1430 digits.Add (characterSet [(int) rem]);
1433 if (sign == -1 && radix == 10) {
1434 NumberFormatInfo info = null;
1435 if (provider != null)
1436 info = provider.GetFormat (typeof (NumberFormatInfo)) as NumberFormatInfo;
1438 string str = info.NegativeSign;
1439 for (int i = str.Length - 1; i >= 0; --i)
1440 digits.Add (str [i]);
1446 char last = digits [digits.Count - 1];
1447 if (sign == 1 && radix > 10 && (last < '0' || last > '9'))
1452 return new String (digits.ToArray ());
1455 public static BigInteger Parse (string value)
1460 if (!Parse (value, false, out result, out ex))
1465 public static bool TryParse (string value, out BigInteger result)
1468 return Parse (value, true, out result, out ex);
1473 public static BigInteger Parse (string value, NumberStyles style)
1475 throw new NotImplementedException ();
1479 public static BigInteger Parse (string value, IFormatProvider provider)
1481 throw new NotImplementedException ();
1485 public static BigInteger Parse (
1486 string value, NumberStyles style, IFormatProvider provider)
1488 throw new InvalidOperationException ();
1492 public static bool TryParse (
1493 string value, NumberStyles style, IFormatProvider provider,
1494 out BigInteger result)
1496 throw new NotImplementedException ();
1501 static Exception GetFormatException ()
1503 return new FormatException ("Input string was not in the correct format");
1506 static bool ProcessTrailingWhitespace (bool tryParse, string s, int position, ref Exception exc)
1510 for (int i = position; i < len; i++){
1513 if (c != 0 && !Char.IsWhiteSpace (c)){
1515 exc = GetFormatException ();
1522 static bool Parse (string s, bool tryParse, out BigInteger result, out Exception exc)
1526 bool digits_seen = false;
1533 exc = new ArgumentNullException ("value");
1540 for (i = 0; i < len; i++){
1542 if (!Char.IsWhiteSpace (c))
1548 exc = GetFormatException ();
1552 var info = Thread.CurrentThread.CurrentCulture.NumberFormat;
1554 string negative = info.NegativeSign;
1555 string positive = info.PositiveSign;
1557 if (string.CompareOrdinal (s, i, positive, 0, positive.Length) == 0)
1558 i += positive.Length;
1559 else if (string.CompareOrdinal (s, i, negative, 0, negative.Length) == 0) {
1561 i += negative.Length;
1564 BigInteger val = Zero;
1565 for (; i < len; i++){
1573 if (c >= '0' && c <= '9'){
1574 byte d = (byte) (c - '0');
1579 } else if (!ProcessTrailingWhitespace (tryParse, s, i, ref exc))
1585 exc = GetFormatException ();
1591 else if (sign == -1)
1592 result = new BigInteger (-1, val.data);
1594 result = new BigInteger (1, val.data);
1599 public static BigInteger Min (BigInteger left, BigInteger right)
1602 int rs = right.sign;
1609 int r = CoreCompare (left.data, right.data);
1619 public static BigInteger Max (BigInteger left, BigInteger right)
1622 int rs = right.sign;
1629 int r = CoreCompare (left.data, right.data);
1638 public static BigInteger Abs (BigInteger value)
1640 return new BigInteger ((short)Math.Abs (value.sign), value.data);
1644 public static BigInteger DivRem (BigInteger dividend, BigInteger divisor, out BigInteger remainder)
1646 if (divisor.sign == 0)
1647 throw new DivideByZeroException ();
1649 if (dividend.sign == 0) {
1650 remainder = dividend;
1655 uint[] remainder_value;
1657 DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
1660 for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
1662 remainder = new BigInteger (0, ZERO);
1664 if (i < remainder_value.Length - 1)
1665 remainder_value = Resize (remainder_value, i + 1);
1666 remainder = new BigInteger (dividend.sign, remainder_value);
1669 for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
1671 return new BigInteger (0, ZERO);
1672 if (i < quotient.Length - 1)
1673 quotient = Resize (quotient, i + 1);
1675 return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
1678 public static BigInteger Pow (BigInteger value, int exponent)
1681 throw new ArgumentOutOfRangeException("exponent", "exp must be >= 0");
1687 BigInteger result = One;
1688 while (exponent != 0) {
1689 if ((exponent & 1) != 0)
1690 result = result * value;
1694 value = value * value;
1700 public static BigInteger ModPow (BigInteger value, BigInteger exponent, BigInteger modulus) {
1701 if (exponent.sign == -1)
1702 throw new ArgumentOutOfRangeException("exponent", "power must be >= 0");
1703 if (modulus.sign == 0)
1704 throw new DivideByZeroException ();
1706 BigInteger result = One % modulus;
1707 while (exponent.sign != 0) {
1708 if (!exponent.IsEven) {
1709 result = result * value;
1710 result = result % modulus;
1714 value = value * value;
1715 value = value % modulus;
1721 public static BigInteger GreatestCommonDivisor (BigInteger left, BigInteger right)
1723 if (left.sign != 0 && left.data.Length == 1 && left.data [0] == 1)
1724 return new BigInteger (1, ONE);
1725 if (right.sign != 0 && right.data.Length == 1 && right.data [0] == 1)
1726 return new BigInteger (1, ONE);
1732 BigInteger x = new BigInteger (1, left.data);
1733 BigInteger y = new BigInteger (1, right.data);
1737 while (x.data.Length > 1) {
1743 if (x.IsZero) return g;
1745 // TODO: should we have something here if we can convert to long?
1748 // Now we can just do it with single precision. I am using the binary gcd method,
1749 // as it should be faster.
1752 uint yy = x.data [0];
1753 uint xx = (uint)(y % yy);
1757 while (((xx | yy) & 1) == 0) {
1758 xx >>= 1; yy >>= 1; t++;
1761 while ((xx & 1) == 0) xx >>= 1;
1762 while ((yy & 1) == 0) yy >>= 1;
1764 xx = (xx - yy) >> 1;
1766 yy = (yy - xx) >> 1;
1772 /*LAMESPEC Log doesn't specify to how many ulp is has to be precise
1773 We are equilavent to MS with about 2 ULP
1775 public static double Log (BigInteger value, Double baseValue)
1777 if (value.sign == -1 || baseValue == 1.0d || baseValue == -1.0d ||
1778 baseValue == Double.NegativeInfinity || double.IsNaN (baseValue))
1781 if (baseValue == 0.0d || baseValue == Double.PositiveInfinity)
1782 return value.IsOne ? 0 : double.NaN;
1784 if (value.sign == 0)
1785 return double.NegativeInfinity;
1787 int length = value.data.Length - 1;
1789 for (int curBit = 31; curBit >= 0; curBit--) {
1790 if ((value.data [length] & (1 << curBit)) != 0) {
1791 bitCount = curBit + length * 32;
1796 long bitlen = bitCount;
1797 Double c = 0, d = 1;
1799 BigInteger testBit = One;
1800 long tempBitlen = bitlen;
1801 while (tempBitlen > Int32.MaxValue) {
1802 testBit = testBit << Int32.MaxValue;
1803 tempBitlen -= Int32.MaxValue;
1805 testBit = testBit << (int)tempBitlen;
1807 for (long curbit = bitlen; curbit >= 0; --curbit) {
1808 if ((value & testBit).sign != 0)
1811 testBit = testBit >> 1;
1813 return (System.Math.Log (c) + System.Math.Log (2) * bitlen) / System.Math.Log (baseValue);
1817 public static double Log (BigInteger value)
1819 return Log (value, Math.E);
1823 public static double Log10 (BigInteger value)
1825 return Log (value, 10);
1828 [CLSCompliantAttribute (false)]
1829 public bool Equals (ulong other)
1831 return CompareTo (other) == 0;
1834 public override int GetHashCode ()
1836 uint hash = (uint)(sign * 0x01010101u);
1838 for (int i = 0; i < data.Length; ++i)
1843 public static BigInteger Add (BigInteger left, BigInteger right)
1845 return left + right;
1848 public static BigInteger Subtract (BigInteger left, BigInteger right)
1850 return left - right;
1853 public static BigInteger Multiply (BigInteger left, BigInteger right)
1855 return left * right;
1858 public static BigInteger Divide (BigInteger dividend, BigInteger divisor)
1860 return dividend / divisor;
1863 public static BigInteger Remainder (BigInteger dividend, BigInteger divisor)
1865 return dividend % divisor;
1868 public static BigInteger Negate (BigInteger value)
1873 public int CompareTo (object obj)
1878 if (!(obj is BigInteger))
1881 return Compare (this, (BigInteger)obj);
1884 public int CompareTo (BigInteger other)
1886 return Compare (this, other);
1889 [CLSCompliantAttribute (false)]
1890 public int CompareTo (ulong other)
1895 return other == 0 ? 0 : -1;
1897 if (data.Length > 2)
1900 uint high = (uint)(other >> 32);
1901 uint low = (uint)other;
1903 return LongCompare (low, high);
1906 int LongCompare (uint low, uint high)
1909 if (data.Length > 1)
1927 public int CompareTo (long other)
1930 int rs = Math.Sign (other);
1933 return ls > rs ? 1 : -1;
1938 if (data.Length > 2)
1943 uint low = (uint)other;
1944 uint high = (uint)((ulong)other >> 32);
1946 int r = LongCompare (low, high);
1953 public static int Compare (BigInteger left, BigInteger right)
1956 int rs = right.sign;
1959 return ls > rs ? 1 : -1;
1961 int r = CoreCompare (left.data, right.data);
1968 static int TopByte (uint x)
1970 if ((x & 0xFFFF0000u) != 0) {
1971 if ((x & 0xFF000000u) != 0)
1975 if ((x & 0xFF00u) != 0)
1980 static int FirstNonFFByte (uint word)
1982 if ((word & 0xFF000000u) != 0xFF000000u)
1984 else if ((word & 0xFF0000u) != 0xFF0000u)
1986 else if ((word & 0xFF00u) != 0xFF00u)
1991 public byte[] ToByteArray ()
1994 return new byte [1];
1996 //number of bytes not counting upper word
1997 int bytes = (data.Length - 1) * 4;
1998 bool needExtraZero = false;
2000 uint topWord = data [data.Length - 1];
2003 //if the topmost bit is set we need an extra
2005 extra = TopByte (topWord);
2006 uint mask = 0x80u << ((extra - 1) * 8);
2007 if ((topWord & mask) != 0) {
2008 needExtraZero = true;
2011 extra = TopByte (topWord);
2014 byte[] res = new byte [bytes + extra + (needExtraZero ? 1 : 0) ];
2017 int end = data.Length - 1;
2018 for (int i = 0; i < end; ++i) {
2019 uint word = data [i];
2021 res [j++] = (byte)word;
2022 res [j++] = (byte)(word >> 8);
2023 res [j++] = (byte)(word >> 16);
2024 res [j++] = (byte)(word >> 24);
2026 while (extra-- > 0) {
2027 res [j++] = (byte)topWord;
2032 int end = data.Length - 1;
2034 uint carry = 1, word;
2036 for (int i = 0; i < end; ++i) {
2038 add = (ulong)~word + carry;
2040 carry = (uint)(add >> 32);
2042 res [j++] = (byte)word;
2043 res [j++] = (byte)(word >> 8);
2044 res [j++] = (byte)(word >> 16);
2045 res [j++] = (byte)(word >> 24);
2048 add = (ulong)~topWord + (carry);
2050 carry = (uint)(add >> 32);
2052 int ex = FirstNonFFByte (word);
2053 bool needExtra = (word & (1 << (ex * 8 - 1))) == 0;
2054 int to = ex + (needExtra ? 1 : 0);
2057 res = Resize (res, bytes + to);
2060 res [j++] = (byte)word;
2066 res = Resize (res, bytes + 5);
2067 res [j++] = (byte)word;
2068 res [j++] = (byte)(word >> 8);
2069 res [j++] = (byte)(word >> 16);
2070 res [j++] = (byte)(word >> 24);
2078 static byte[] Resize (byte[] v, int len)
2080 byte[] res = new byte [len];
2081 Array.Copy (v, res, Math.Min (v.Length, len));
2085 static uint[] Resize (uint[] v, int len)
2087 uint[] res = new uint [len];
2088 Array.Copy (v, res, Math.Min (v.Length, len));
2092 static uint[] CoreAdd (uint[] a, uint[] b)
2094 if (a.Length < b.Length) {
2103 uint[] res = new uint [bl];
2108 for (; i < sl; i++) {
2109 sum = sum + a [i] + b [i];
2110 res [i] = (uint)sum;
2114 for (; i < bl; i++) {
2116 res [i] = (uint)sum;
2121 res = Resize (res, bl + 1);
2122 res [i] = (uint)sum;
2129 static uint[] CoreSub (uint[] a, uint[] b)
2134 uint[] res = new uint [bl];
2138 for (i = 0; i < sl; ++i) {
2139 borrow = (ulong)a [i] - b [i] - borrow;
2141 res [i] = (uint)borrow;
2142 borrow = (borrow >> 32) & 0x1;
2145 for (; i < bl; i++) {
2146 borrow = (ulong)a [i] - borrow;
2147 res [i] = (uint)borrow;
2148 borrow = (borrow >> 32) & 0x1;
2151 //remove extra zeroes
2152 for (i = bl - 1; i >= 0 && res [i] == 0; --i) ;
2154 res = Resize (res, i + 1);
2160 static uint[] CoreAdd (uint[] a, uint b)
2163 uint[] res = new uint [len];
2167 for (i = 0; i < len; i++) {
2169 res [i] = (uint)sum;
2174 res = Resize (res, len + 1);
2175 res [i] = (uint)sum;
2181 static uint[] CoreSub (uint[] a, uint b)
2184 uint[] res = new uint [len];
2188 for (i = 0; i < len; i++) {
2189 borrow = (ulong)a [i] - borrow;
2190 res [i] = (uint)borrow;
2191 borrow = (borrow >> 32) & 0x1;
2194 //remove extra zeroes
2195 for (i = len - 1; i >= 0 && res [i] == 0; --i) ;
2197 res = Resize (res, i + 1);
2202 static int CoreCompare (uint[] a, uint[] b)
2204 int al = a != null ? a.Length : 0;
2205 int bl = b != null ? b.Length : 0;
2212 for (int i = al - 1; i >= 0; --i) {
2223 static int GetNormalizeShift(uint value) {
2226 if ((value & 0xFFFF0000) == 0) { value <<= 16; shift += 16; }
2227 if ((value & 0xFF000000) == 0) { value <<= 8; shift += 8; }
2228 if ((value & 0xF0000000) == 0) { value <<= 4; shift += 4; }
2229 if ((value & 0xC0000000) == 0) { value <<= 2; shift += 2; }
2230 if ((value & 0x80000000) == 0) { value <<= 1; shift += 1; }
2235 static void Normalize (uint[] u, int l, uint[] un, int shift)
2240 int rshift = 32 - shift;
2241 for (i = 0; i < l; i++) {
2243 un [i] = (ui << shift) | carry;
2244 carry = ui >> rshift;
2247 for (i = 0; i < l; i++) {
2252 while (i < un.Length) {
2261 static void Unnormalize (uint[] un, out uint[] r, int shift)
2263 int length = un.Length;
2264 r = new uint [length];
2267 int lshift = 32 - shift;
2269 for (int i = length - 1; i >= 0; i--) {
2271 r [i] = (uni >> shift) | carry;
2272 carry = (uni << lshift);
2275 for (int i = 0; i < length; i++) {
2281 const ulong Base = 0x100000000;
2282 static void DivModUnsigned (uint[] u, uint[] v, out uint[] q, out uint[] r)
2288 // Divide by single digit
2295 for (int j = m - 1; j >= 0; j--) {
2299 ulong div = rem / v0;
2304 } else if (m >= n) {
2305 int shift = GetNormalizeShift (v [n - 1]);
2307 uint[] un = new uint [m + 1];
2308 uint[] vn = new uint [n];
2310 Normalize (u, m, un, shift);
2311 Normalize (v, n, vn, shift);
2313 q = new uint [m - n + 1];
2316 // Main division loop
2318 for (int j = m - n; j >= 0; j--) {
2322 rr = Base * un [j + n] + un [j + n - 1];
2323 qq = rr / vn [n - 1];
2324 rr -= qq * vn [n - 1];
2327 // Estimate too big ?
2329 if ((qq >= Base) || (qq * vn [n - 2] > (rr * Base + un [j + n - 2]))) {
2331 rr += (ulong)vn [n - 1];
2339 // Multiply and subtract
2343 for (i = 0; i < n; i++) {
2344 ulong p = vn [i] * qq;
2345 t = (long)un [i + j] - (long)(uint)p - b;
2346 un [i + j] = (uint)t;
2351 t = (long)un [j + n] - b;
2352 un [j + n] = (uint)t;
2354 // Store the calculated value
2358 // Add back vn[0..n] to un[j..j+n]
2363 for (i = 0; i < n; i++) {
2364 c = (ulong)vn [i] + un [j + i] + c;
2365 un [j + i] = (uint)c;
2368 c += (ulong)un [j + n];
2369 un [j + n] = (uint)c;
2373 Unnormalize (un, out r, shift);
2375 q = new uint [] { 0 };