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 (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)
387 if (value.data.Length > 1)
388 throw new OverflowException ();
389 uint data = value.data [0];
391 if (value.sign == 1) {
392 if (data > (uint)int.MaxValue)
393 throw new OverflowException ();
395 } else if (value.sign == -1) {
396 if (data > 0x80000000u)
397 throw new OverflowException ();
404 [CLSCompliantAttribute (false)]
405 public static explicit operator uint (BigInteger value)
407 if (value.data.Length > 1 || value.sign == -1)
408 throw new OverflowException ();
409 return value.data [0];
412 public static explicit operator short (BigInteger value)
414 int val = (int)value;
415 if (val < short.MinValue || val > short.MaxValue)
416 throw new OverflowException ();
420 [CLSCompliantAttribute (false)]
421 public static explicit operator ushort (BigInteger value)
423 uint val = (uint)value;
424 if (val > ushort.MaxValue)
425 throw new OverflowException ();
429 public static explicit operator byte (BigInteger value)
431 uint val = (uint)value;
432 if (val > byte.MaxValue)
433 throw new OverflowException ();
437 [CLSCompliantAttribute (false)]
438 public static explicit operator sbyte (BigInteger value)
440 int val = (int)value;
441 if (val < sbyte.MinValue || val > sbyte.MaxValue)
442 throw new OverflowException ();
447 public static explicit operator long (BigInteger value)
452 if (value.data.Length > 2)
453 throw new OverflowException ();
455 uint low = value.data [0];
457 if (value.data.Length == 1) {
460 long res = (long)low;
464 uint high = value.data [1];
466 if (value.sign == 1) {
467 if (high >= 0x80000000u)
468 throw new OverflowException ();
469 return (((long)high) << 32) | low;
472 if (high > 0x80000000u)
473 throw new OverflowException ();
475 return - ((((long)high) << 32) | (long)low);
478 [CLSCompliantAttribute (false)]
479 public static explicit operator ulong (BigInteger value)
481 if (value.data.Length > 2 || value.sign == -1)
482 throw new OverflowException ();
484 uint low = value.data [0];
485 if (value.data.Length == 1)
488 uint high = value.data [1];
489 return (((ulong)high) << 32) | low;
492 public static explicit operator double (BigInteger value)
496 return double.Parse (value.ToString (),
497 System.Globalization.CultureInfo.InvariantCulture.NumberFormat);
498 } catch (OverflowException) {
499 return value.sign == -1 ? double.NegativeInfinity : double.PositiveInfinity;
503 public static explicit operator float (BigInteger value)
507 return float.Parse (value.ToString (),
508 System.Globalization.CultureInfo.InvariantCulture.NumberFormat);
509 } catch (OverflowException) {
510 return value.sign == -1 ? float.NegativeInfinity : float.PositiveInfinity;
514 public static explicit operator decimal (BigInteger value)
519 uint[] data = value.data;
521 throw new OverflowException ();
523 int lo = 0, mi = 0, hi = 0;
525 hi = (Int32)data [2];
527 mi = (Int32)data [1];
529 lo = (Int32)data [0];
531 return new Decimal(lo, mi, hi, value.sign < 0, 0);
534 public static implicit operator BigInteger (int value)
536 return new BigInteger (value);
539 [CLSCompliantAttribute (false)]
540 public static implicit operator BigInteger (uint value)
542 return new BigInteger (value);
545 public static implicit operator BigInteger (short value)
547 return new BigInteger (value);
550 [CLSCompliantAttribute (false)]
551 public static implicit operator BigInteger (ushort value)
553 return new BigInteger (value);
556 public static implicit operator BigInteger (byte value)
558 return new BigInteger (value);
561 [CLSCompliantAttribute (false)]
562 public static implicit operator BigInteger (sbyte value)
564 return new BigInteger (value);
567 public static implicit operator BigInteger (long value)
569 return new BigInteger (value);
572 [CLSCompliantAttribute (false)]
573 public static implicit operator BigInteger (ulong value)
575 return new BigInteger (value);
578 public static explicit operator BigInteger (double value)
580 return new BigInteger (value);
583 public static explicit operator BigInteger (float value)
585 return new BigInteger (value);
588 public static explicit operator BigInteger (decimal value)
590 return new BigInteger (value);
593 public static BigInteger operator+ (BigInteger left, BigInteger right)
600 if (left.sign == right.sign)
601 return new BigInteger (left.sign, CoreAdd (left.data, right.data));
603 int r = CoreCompare (left.data, right.data);
606 return new BigInteger (0, ZERO);
608 if (r > 0) //left > right
609 return new BigInteger (left.sign, CoreSub (left.data, right.data));
611 return new BigInteger (right.sign, CoreSub (right.data, left.data));
614 public static BigInteger operator- (BigInteger left, BigInteger right)
619 return new BigInteger ((short)-right.sign, right.data);
621 if (left.sign == right.sign) {
622 int r = CoreCompare (left.data, right.data);
625 return new BigInteger (0, ZERO);
627 if (r > 0) //left > right
628 return new BigInteger (left.sign, CoreSub (left.data, right.data));
630 return new BigInteger ((short)-right.sign, CoreSub (right.data, left.data));
633 return new BigInteger (left.sign, CoreAdd (left.data, right.data));
636 public static BigInteger operator* (BigInteger left, BigInteger right)
638 if (left.sign == 0 || right.sign == 0)
639 return new BigInteger (0, ZERO);
641 if (left.data [0] == 1 && left.data.Length == 1) {
644 return new BigInteger ((short)-right.sign, right.data);
647 if (right.data [0] == 1 && right.data.Length == 1) {
650 return new BigInteger ((short)-left.sign, left.data);
653 uint[] a = left.data;
654 uint[] b = right.data;
656 uint[] res = new uint [a.Length + b.Length];
658 for (int i = 0; i < a.Length; ++i) {
663 for (int j = 0; j < b.Length; ++j) {
664 carry = carry + ((ulong)ai) * b [j] + res [k];
665 res[k++] = (uint)carry;
671 res[k++] = (uint)carry;
677 for (m = res.Length - 1; m >= 0 && res [m] == 0; --m) ;
678 if (m < res.Length - 1)
679 res = Resize (res, m + 1);
681 return new BigInteger ((short)(left.sign * right.sign), res);
684 public static BigInteger operator/ (BigInteger dividend, BigInteger divisor)
686 if (divisor.sign == 0)
687 throw new DivideByZeroException ();
689 if (dividend.sign == 0)
693 uint[] remainder_value;
695 DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
698 for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
700 return new BigInteger (0, ZERO);
701 if (i < quotient.Length - 1)
702 quotient = Resize (quotient, i + 1);
704 return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
707 public static BigInteger operator% (BigInteger dividend, BigInteger divisor)
709 if (divisor.sign == 0)
710 throw new DivideByZeroException ();
712 if (dividend.sign == 0)
716 uint[] remainder_value;
718 DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
721 for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
723 return new BigInteger (0, ZERO);
725 if (i < remainder_value.Length - 1)
726 remainder_value = Resize (remainder_value, i + 1);
727 return new BigInteger (dividend.sign, remainder_value);
730 public static BigInteger operator- (BigInteger value)
734 return new BigInteger ((short)-value.sign, value.data);
737 public static BigInteger operator+ (BigInteger value)
742 public static BigInteger operator++ (BigInteger value)
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);
750 return new BigInteger (1, ONE);
754 data = CoreSub (data, 1);
756 data = CoreAdd (data, 1);
758 return new BigInteger (sign, data);
761 public static BigInteger operator-- (BigInteger value)
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);
769 return new BigInteger (-1, ONE);
773 data = CoreAdd (data, 1);
775 data = CoreSub (data, 1);
777 return new BigInteger (sign, data);
780 public static BigInteger operator& (BigInteger left, BigInteger right)
788 uint[] a = left.data;
789 uint[] b = right.data;
793 bool neg_res = (ls == rs) && (ls == -1);
795 uint[] result = new uint [Math.Max (a.Length, b.Length)];
797 ulong ac = 1, bc = 1, borrow = 1;
800 for (i = 0; i < result.Length; ++i) {
807 ac = (uint)(ac >> 32);
816 bc = (uint)(bc >> 32);
822 borrow = word - borrow;
823 word = ~(uint)borrow;
824 borrow = (uint)(borrow >> 32) & 0x1u;
830 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
832 return new BigInteger (0, ZERO);
834 if (i < result.Length - 1)
835 result = Resize (result, i + 1);
837 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
840 public static BigInteger operator| (BigInteger left, BigInteger right)
848 uint[] a = left.data;
849 uint[] b = right.data;
853 bool neg_res = (ls == -1) || (rs == -1);
855 uint[] result = new uint [Math.Max (a.Length, b.Length)];
857 ulong ac = 1, bc = 1, borrow = 1;
860 for (i = 0; i < result.Length; ++i) {
867 ac = (uint)(ac >> 32);
876 bc = (uint)(bc >> 32);
882 borrow = word - borrow;
883 word = ~(uint)borrow;
884 borrow = (uint)(borrow >> 32) & 0x1u;
890 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
892 return new BigInteger (0, ZERO);
894 if (i < result.Length - 1)
895 result = Resize (result, i + 1);
897 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
900 public static BigInteger operator^ (BigInteger left, BigInteger right)
908 uint[] a = left.data;
909 uint[] b = right.data;
913 bool neg_res = (ls == -1) ^ (rs == -1);
915 uint[] result = new uint [Math.Max (a.Length, b.Length)];
917 ulong ac = 1, bc = 1, borrow = 1;
920 for (i = 0; i < result.Length; ++i) {
927 ac = (uint)(ac >> 32);
936 bc = (uint)(bc >> 32);
942 borrow = word - borrow;
943 word = ~(uint)borrow;
944 borrow = (uint)(borrow >> 32) & 0x1u;
950 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
952 return new BigInteger (0, ZERO);
954 if (i < result.Length - 1)
955 result = Resize (result, i + 1);
957 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
960 public static BigInteger operator~ (BigInteger value)
963 return new BigInteger (-1, ONE);
965 uint[] data = value.data;
966 int sign = value.sign;
968 bool neg_res = sign == 1;
970 uint[] result = new uint [data.Length];
972 ulong carry = 1, borrow = 1;
975 for (i = 0; i < result.Length; ++i) {
976 uint word = data [i];
978 carry = ~word + carry;
980 carry = (uint)(carry >> 32);
986 borrow = word - borrow;
987 word = ~(uint)borrow;
988 borrow = (uint)(borrow >> 32) & 0x1u;
994 for (i = result.Length - 1; i >= 0 && result [i] == 0; --i) ;
996 return new BigInteger (0, ZERO);
998 if (i < result.Length - 1)
999 result = Resize (result, i + 1);
1001 return new BigInteger (neg_res ? (short)-1 : (short)1, result);
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)
1008 for (int i = 31; i >= 0; --i) {
1009 uint mask = 1u << i;
1010 if ((word & mask) == mask)
1016 public static BigInteger operator<< (BigInteger value, int shift)
1018 if (shift == 0 || value.sign == 0)
1021 return value >> -shift;
1023 uint[] data = value.data;
1024 int sign = value.sign;
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);
1030 uint[] res = new uint [data.Length + extra_words];
1032 int idx_shift = shift >> 5;
1033 int bit_shift = shift & 0x1F;
1034 int carry_shift = 32 - bit_shift;
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;
1043 return new BigInteger ((short)sign, res);
1046 public static BigInteger operator>> (BigInteger value, int shift)
1048 if (shift == 0 || value.sign == 0)
1051 return value << -shift;
1053 uint[] data = value.data;
1054 int sign = value.sign;
1056 int topMostIdx = BitScanBackward (data [data.Length - 1]);
1057 int idx_shift = shift >> 5;
1058 int bit_shift = shift & 0x1F;
1060 int extra_words = idx_shift;
1061 if (bit_shift > topMostIdx)
1063 int size = data.Length - extra_words;
1067 return new BigInteger (0, ZERO);
1068 return new BigInteger (-1, ONE);
1071 uint[] res = new uint [size];
1072 int carry_shift = 32 - bit_shift;
1074 for (int i = data.Length - 1; i >= idx_shift; --i) {
1075 uint word = data [i];
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;
1083 //Round down instead of toward zero
1085 for (int i = 0; i < idx_shift; i++) {
1086 if (data [i] != 0u) {
1087 var tmp = new BigInteger ((short)sign, res);
1092 if (bit_shift > 0 && (data [idx_shift] << carry_shift) != 0u) {
1093 var tmp = new BigInteger ((short)sign, res);
1098 return new BigInteger ((short)sign, res);
1101 public static bool operator< (BigInteger left, BigInteger right)
1103 return Compare (left, right) < 0;
1106 public static bool operator< (BigInteger left, long right)
1108 return left.CompareTo (right) < 0;
1112 public static bool operator< (long left, BigInteger right)
1114 return right.CompareTo (left) > 0;
1118 [CLSCompliantAttribute (false)]
1119 public static bool operator< (BigInteger left, ulong right)
1121 return left.CompareTo (right) < 0;
1124 [CLSCompliantAttribute (false)]
1125 public static bool operator< (ulong left, BigInteger right)
1127 return right.CompareTo (left) > 0;
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;
1140 public static bool operator<= (long left, BigInteger right)
1142 return right.CompareTo (left) >= 0;
1145 [CLSCompliantAttribute (false)]
1146 public static bool operator<= (BigInteger left, ulong right)
1148 return left.CompareTo (right) <= 0;
1151 [CLSCompliantAttribute (false)]
1152 public static bool operator<= (ulong left, BigInteger right)
1154 return right.CompareTo (left) >= 0;
1157 public static bool operator> (BigInteger left, BigInteger right)
1159 return Compare (left, right) > 0;
1162 public static bool operator> (BigInteger left, long right)
1164 return left.CompareTo (right) > 0;
1167 public static bool operator> (long left, BigInteger right)
1169 return right.CompareTo (left) < 0;
1172 [CLSCompliantAttribute (false)]
1173 public static bool operator> (BigInteger left, ulong right)
1175 return left.CompareTo (right) > 0;
1178 [CLSCompliantAttribute (false)]
1179 public static bool operator> (ulong left, BigInteger right)
1181 return right.CompareTo (left) < 0;
1184 public static bool operator>= (BigInteger left, BigInteger right)
1186 return Compare (left, right) >= 0;
1189 public static bool operator>= (BigInteger left, long right)
1191 return left.CompareTo (right) >= 0;
1194 public static bool operator>= (long left, BigInteger right)
1196 return right.CompareTo (left) <= 0;
1199 [CLSCompliantAttribute (false)]
1200 public static bool operator>= (BigInteger left, ulong right)
1202 return left.CompareTo (right) >= 0;
1205 [CLSCompliantAttribute (false)]
1206 public static bool operator>= (ulong left, BigInteger right)
1208 return right.CompareTo (left) <= 0;
1211 public static bool operator== (BigInteger left, BigInteger right)
1213 return Compare (left, right) == 0;
1216 public static bool operator== (BigInteger left, long right)
1218 return left.CompareTo (right) == 0;
1221 public static bool operator== (long left, BigInteger right)
1223 return right.CompareTo (left) == 0;
1226 [CLSCompliantAttribute (false)]
1227 public static bool operator== (BigInteger left, ulong right)
1229 return left.CompareTo (right) == 0;
1232 [CLSCompliantAttribute (false)]
1233 public static bool operator== (ulong left, BigInteger right)
1235 return right.CompareTo (left) == 0;
1238 public static bool operator!= (BigInteger left, BigInteger right)
1240 return Compare (left, right) != 0;
1243 public static bool operator!= (BigInteger left, long right)
1245 return left.CompareTo (right) != 0;
1248 public static bool operator!= (long left, BigInteger right)
1250 return right.CompareTo (left) != 0;
1253 [CLSCompliantAttribute (false)]
1254 public static bool operator!= (BigInteger left, ulong right)
1256 return left.CompareTo (right) != 0;
1259 [CLSCompliantAttribute (false)]
1260 public static bool operator!= (ulong left, BigInteger right)
1262 return right.CompareTo (left) != 0;
1265 public override bool Equals (object obj)
1267 if (!(obj is BigInteger))
1269 return Equals ((BigInteger)obj);
1272 public bool Equals (BigInteger other)
1274 if (sign != other.sign)
1276 if (data.Length != other.data.Length)
1278 for (int i = 0; i < data.Length; ++i) {
1279 if (data [i] != other.data [i])
1285 public bool Equals (long other)
1287 return CompareTo (other) == 0;
1290 public override string ToString ()
1292 return ToString (10, null);
1295 string ToStringWithPadding (string format, uint radix, IFormatProvider provider)
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;
1305 return "-" + additional + baseStr.Substring (1);
1310 return ToString (radix, provider);
1313 public string ToString (string format)
1315 return ToString (format, null);
1318 public string ToString (IFormatProvider provider)
1320 return ToString (null, provider);
1324 public string ToString (string format, IFormatProvider provider)
1326 if (format == null || format == "")
1327 return ToString (10, provider);
1329 switch (format[0]) {
1336 return ToStringWithPadding (format, 10, provider);
1339 return ToStringWithPadding (format, 16, null);
1341 throw new FormatException (string.Format ("format '{0}' not implemented", format));
1345 static uint[] MakeTwoComplement (uint[] v)
1347 uint[] res = new uint [v.Length];
1350 for (int i = 0; i < v.Length; ++i) {
1352 carry = (ulong)~word + carry;
1354 carry = (uint)(carry >> 32);
1358 uint last = res [res.Length - 1];
1359 int idx = FirstNonFFByte (last);
1361 for (int i = 1; i < idx; ++i)
1362 mask = (mask << 8) | 0xFF;
1364 res [res.Length - 1] = last & mask;
1368 string ToString (uint radix, IFormatProvider provider)
1370 const string characterSet = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
1372 if (characterSet.Length < radix)
1373 throw new ArgumentException ("charSet length less than radix", "characterSet");
1375 throw new ArgumentException ("There is no such thing as radix one notation", "radix");
1379 if (data.Length == 1 && data [0] == 1)
1380 return sign == 1 ? "1" : "-1";
1382 List<char> digits = new List<char> (1 + data.Length * 3 / 10);
1390 dt = MakeTwoComplement (dt);
1391 a = new BigInteger (1, dt);
1396 a = DivRem (a, radix, out rem);
1397 digits.Add (characterSet [(int) rem]);
1400 if (sign == -1 && radix == 10) {
1401 NumberFormatInfo info = null;
1402 if (provider != null)
1403 info = provider.GetFormat (typeof (NumberFormatInfo)) as NumberFormatInfo;
1405 string str = info.NegativeSign;
1406 for (int i = str.Length - 1; i >= 0; --i)
1407 digits.Add (str [i]);
1413 char last = digits [digits.Count - 1];
1414 if (sign == 1 && radix > 10 && (last < '0' || last > '9'))
1419 return new String (digits.ToArray ());
1422 public static BigInteger Parse (string value)
1427 if (!Parse (value, false, out result, out ex))
1433 public static bool TryParse (string value, out BigInteger result)
1436 return Parse (value, true, out result, out ex);
1439 static Exception GetFormatException ()
1441 return new FormatException ("Input string was not in the correct format");
1444 static bool ProcessTrailingWhitespace (bool tryParse, string s, int position, ref Exception exc)
1448 for (int i = position; i < len; i++){
1451 if (c != 0 && !Char.IsWhiteSpace (c)){
1453 exc = GetFormatException ();
1460 static bool Parse (string s, bool tryParse, out BigInteger result, out Exception exc)
1464 bool digits_seen = false;
1471 exc = new ArgumentNullException ("value");
1478 for (i = 0; i < len; i++){
1480 if (!Char.IsWhiteSpace (c))
1486 exc = GetFormatException ();
1490 var info = Thread.CurrentThread.CurrentCulture.NumberFormat;
1492 string negative = info.NegativeSign;
1493 string positive = info.PositiveSign;
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) {
1499 i += negative.Length;
1502 BigInteger val = Zero;
1503 for (; i < len; i++){
1511 if (c >= '0' && c <= '9'){
1512 byte d = (byte) (c - '0');
1517 } else if (!ProcessTrailingWhitespace (tryParse, s, i, ref exc))
1523 exc = GetFormatException ();
1529 else if (sign == -1)
1530 result = new BigInteger (-1, val.data);
1532 result = new BigInteger (1, val.data);
1537 public static BigInteger Min (BigInteger left, BigInteger right)
1540 int rs = right.sign;
1547 int r = CoreCompare (left.data, right.data);
1557 public static BigInteger Max (BigInteger left, BigInteger right)
1560 int rs = right.sign;
1567 int r = CoreCompare (left.data, right.data);
1576 public static BigInteger Abs (BigInteger value)
1578 return new BigInteger ((short)Math.Abs (value.sign), value.data);
1582 public static BigInteger DivRem (BigInteger dividend, BigInteger divisor, out BigInteger remainder)
1584 if (divisor.sign == 0)
1585 throw new DivideByZeroException ();
1587 if (dividend.sign == 0) {
1588 remainder = dividend;
1593 uint[] remainder_value;
1595 DivModUnsigned (dividend.data, divisor.data, out quotient, out remainder_value);
1598 for (i = remainder_value.Length - 1; i >= 0 && remainder_value [i] == 0; --i) ;
1600 remainder = new BigInteger (0, ZERO);
1602 if (i < remainder_value.Length - 1)
1603 remainder_value = Resize (remainder_value, i + 1);
1604 remainder = new BigInteger (dividend.sign, remainder_value);
1607 for (i = quotient.Length - 1; i >= 0 && quotient [i] == 0; --i) ;
1609 return new BigInteger (0, ZERO);
1610 if (i < quotient.Length - 1)
1611 quotient = Resize (quotient, i + 1);
1613 return new BigInteger ((short)(dividend.sign * divisor.sign), quotient);
1616 public static BigInteger Pow (BigInteger value, int exponent)
1619 throw new ArgumentOutOfRangeException("exponent", "exp must be >= 0");
1625 BigInteger result = One;
1626 while (exponent != 0) {
1627 if ((exponent & 1) != 0)
1628 result = result * value;
1632 value = value * value;
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 ();
1644 BigInteger result = One % modulus;
1645 while (exponent.sign != 0) {
1646 if (!exponent.IsEven) {
1647 result = result * value;
1648 result = result % modulus;
1652 value = value * value;
1653 value = value % modulus;
1659 public static BigInteger GreatestCommonDivisor (BigInteger left, BigInteger right)
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);
1670 BigInteger x = new BigInteger (1, left.data);
1671 BigInteger y = new BigInteger (1, right.data);
1675 while (x.data.Length > 1) {
1681 if (x.IsZero) return g;
1683 // TODO: should we have something here if we can convert to long?
1686 // Now we can just do it with single precision. I am using the binary gcd method,
1687 // as it should be faster.
1690 uint yy = x.data [0];
1691 uint xx = (uint)(y % yy);
1695 while (((xx | yy) & 1) == 0) {
1696 xx >>= 1; yy >>= 1; t++;
1699 while ((xx & 1) == 0) xx >>= 1;
1700 while ((yy & 1) == 0) yy >>= 1;
1702 xx = (xx - yy) >> 1;
1704 yy = (yy - xx) >> 1;
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
1713 public static double Log (BigInteger value, Double baseValue)
1715 if (value.sign == -1 || baseValue == 1.0d || baseValue == -1.0d ||
1716 baseValue == Double.NegativeInfinity || double.IsNaN (baseValue))
1719 if (baseValue == 0.0d || baseValue == Double.PositiveInfinity)
1720 return value.IsOne ? 0 : double.NaN;
1722 if (value.sign == 0)
1723 return double.NegativeInfinity;
1725 int length = value.data.Length - 1;
1727 for (int curBit = 31; curBit >= 0; curBit--) {
1728 if ((value.data [length] & (1 << curBit)) != 0) {
1729 bitCount = curBit + length * 32;
1734 long bitlen = bitCount;
1735 Double c = 0, d = 1;
1737 BigInteger testBit = One;
1738 long tempBitlen = bitlen;
1739 while (tempBitlen > Int32.MaxValue) {
1740 testBit = testBit << Int32.MaxValue;
1741 tempBitlen -= Int32.MaxValue;
1743 testBit = testBit << (int)tempBitlen;
1745 for (long curbit = bitlen; curbit >= 0; --curbit) {
1746 if ((value & testBit).sign != 0)
1749 testBit = testBit >> 1;
1751 return (System.Math.Log (c) + System.Math.Log (2) * bitlen) / System.Math.Log (baseValue);
1755 public static double Log (BigInteger value)
1757 return Log (value, Math.E);
1761 public static double Log10 (BigInteger value)
1763 return Log (value, 10);
1766 [CLSCompliantAttribute (false)]
1767 public bool Equals (ulong other)
1769 return CompareTo (other) == 0;
1772 public override int GetHashCode ()
1774 uint hash = (uint)(sign * 0x01010101u);
1776 for (int i = 0; i < data.Length; ++i)
1781 public static BigInteger Add (BigInteger left, BigInteger right)
1783 return left + right;
1786 public static BigInteger Subtract (BigInteger left, BigInteger right)
1788 return left - right;
1791 public static BigInteger Multiply (BigInteger left, BigInteger right)
1793 return left * right;
1796 public static BigInteger Divide (BigInteger dividend, BigInteger divisor)
1798 return dividend / divisor;
1801 public static BigInteger Remainder (BigInteger dividend, BigInteger divisor)
1803 return dividend % divisor;
1806 public static BigInteger Negate (BigInteger value)
1811 public int CompareTo (object obj)
1816 if (!(obj is BigInteger))
1819 return Compare (this, (BigInteger)obj);
1822 public int CompareTo (BigInteger other)
1824 return Compare (this, other);
1827 [CLSCompliantAttribute (false)]
1828 public int CompareTo (ulong other)
1833 return other == 0 ? 0 : -1;
1835 if (data.Length > 2)
1838 uint high = (uint)(other >> 32);
1839 uint low = (uint)other;
1841 return LongCompare (low, high);
1844 int LongCompare (uint low, uint high)
1847 if (data.Length > 1)
1865 public int CompareTo (long other)
1868 int rs = Math.Sign (other);
1871 return ls > rs ? 1 : -1;
1876 if (data.Length > 2)
1881 uint low = (uint)other;
1882 uint high = (uint)((ulong)other >> 32);
1884 int r = LongCompare (low, high);
1891 public static int Compare (BigInteger left, BigInteger right)
1894 int rs = right.sign;
1897 return ls > rs ? 1 : -1;
1899 int r = CoreCompare (left.data, right.data);
1906 static int TopByte (uint x)
1908 if ((x & 0xFFFF0000u) != 0) {
1909 if ((x & 0xFF000000u) != 0)
1913 if ((x & 0xFF00u) != 0)
1918 static int FirstNonFFByte (uint word)
1920 if ((word & 0xFF000000u) != 0xFF000000u)
1922 else if ((word & 0xFF0000u) != 0xFF0000u)
1924 else if ((word & 0xFF00u) != 0xFF00u)
1929 public byte[] ToByteArray ()
1932 return new byte [1];
1934 //number of bytes not counting upper word
1935 int bytes = (data.Length - 1) * 4;
1936 bool needExtraZero = false;
1938 uint topWord = data [data.Length - 1];
1941 //if the topmost bit is set we need an extra
1943 extra = TopByte (topWord);
1944 uint mask = 0x80u << ((extra - 1) * 8);
1945 if ((topWord & mask) != 0) {
1946 needExtraZero = true;
1949 extra = TopByte (topWord);
1952 byte[] res = new byte [bytes + extra + (needExtraZero ? 1 : 0) ];
1955 int end = data.Length - 1;
1956 for (int i = 0; i < end; ++i) {
1957 uint word = data [i];
1959 res [j++] = (byte)word;
1960 res [j++] = (byte)(word >> 8);
1961 res [j++] = (byte)(word >> 16);
1962 res [j++] = (byte)(word >> 24);
1964 while (extra-- > 0) {
1965 res [j++] = (byte)topWord;
1970 int end = data.Length - 1;
1972 uint carry = 1, word;
1974 for (int i = 0; i < end; ++i) {
1976 add = (ulong)~word + carry;
1978 carry = (uint)(add >> 32);
1980 res [j++] = (byte)word;
1981 res [j++] = (byte)(word >> 8);
1982 res [j++] = (byte)(word >> 16);
1983 res [j++] = (byte)(word >> 24);
1986 add = (ulong)~topWord + (carry);
1988 carry = (uint)(add >> 32);
1990 int ex = FirstNonFFByte (word);
1991 bool needExtra = (word & (1 << (ex * 8 - 1))) == 0;
1992 int to = ex + (needExtra ? 1 : 0);
1995 res = Resize (res, bytes + to);
1998 res [j++] = (byte)word;
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);
2016 static byte[] Resize (byte[] v, int len)
2018 byte[] res = new byte [len];
2019 Array.Copy (v, res, Math.Min (v.Length, len));
2023 static uint[] Resize (uint[] v, int len)
2025 uint[] res = new uint [len];
2026 Array.Copy (v, res, Math.Min (v.Length, len));
2030 static uint[] CoreAdd (uint[] a, uint[] b)
2032 if (a.Length < b.Length) {
2041 uint[] res = new uint [bl];
2046 for (; i < sl; i++) {
2047 sum = sum + a [i] + b [i];
2048 res [i] = (uint)sum;
2052 for (; i < bl; i++) {
2054 res [i] = (uint)sum;
2059 res = Resize (res, bl + 1);
2060 res [i] = (uint)sum;
2067 static uint[] CoreSub (uint[] a, uint[] b)
2072 uint[] res = new uint [bl];
2076 for (i = 0; i < sl; ++i) {
2077 borrow = (ulong)a [i] - b [i] - borrow;
2079 res [i] = (uint)borrow;
2080 borrow = (borrow >> 32) & 0x1;
2083 for (; i < bl; i++) {
2084 borrow = (ulong)a [i] - borrow;
2085 res [i] = (uint)borrow;
2086 borrow = (borrow >> 32) & 0x1;
2089 //remove extra zeroes
2090 for (i = bl - 1; i >= 0 && res [i] == 0; --i) ;
2092 res = Resize (res, i + 1);
2098 static uint[] CoreAdd (uint[] a, uint b)
2101 uint[] res = new uint [len];
2105 for (i = 0; i < len; i++) {
2107 res [i] = (uint)sum;
2112 res = Resize (res, len + 1);
2113 res [i] = (uint)sum;
2119 static uint[] CoreSub (uint[] a, uint b)
2122 uint[] res = new uint [len];
2126 for (i = 0; i < len; i++) {
2127 borrow = (ulong)a [i] - borrow;
2128 res [i] = (uint)borrow;
2129 borrow = (borrow >> 32) & 0x1;
2132 //remove extra zeroes
2133 for (i = len - 1; i >= 0 && res [i] == 0; --i) ;
2135 res = Resize (res, i + 1);
2140 static int CoreCompare (uint[] a, uint[] b)
2150 for (int i = al - 1; i >= 0; --i) {
2161 static int GetNormalizeShift(uint value) {
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; }
2173 static void Normalize (uint[] u, int l, uint[] un, int shift)
2178 int rshift = 32 - shift;
2179 for (i = 0; i < l; i++) {
2181 un [i] = (ui << shift) | carry;
2182 carry = ui >> rshift;
2185 for (i = 0; i < l; i++) {
2190 while (i < un.Length) {
2199 static void Unnormalize (uint[] un, out uint[] r, int shift)
2201 int length = un.Length;
2202 r = new uint [length];
2205 int lshift = 32 - shift;
2207 for (int i = length - 1; i >= 0; i--) {
2209 r [i] = (uni >> shift) | carry;
2210 carry = (uni << lshift);
2213 for (int i = 0; i < length; i++) {
2219 const ulong Base = 0x100000000;
2220 static void DivModUnsigned (uint[] u, uint[] v, out uint[] q, out uint[] r)
2226 // Divide by single digit
2233 for (int j = m - 1; j >= 0; j--) {
2237 ulong div = rem / v0;
2242 } else if (m >= n) {
2243 int shift = GetNormalizeShift (v [n - 1]);
2245 uint[] un = new uint [m + 1];
2246 uint[] vn = new uint [n];
2248 Normalize (u, m, un, shift);
2249 Normalize (v, n, vn, shift);
2251 q = new uint [m - n + 1];
2254 // Main division loop
2256 for (int j = m - n; j >= 0; j--) {
2260 rr = Base * un [j + n] + un [j + n - 1];
2261 qq = rr / vn [n - 1];
2262 rr -= qq * vn [n - 1];
2265 // Estimate too big ?
2267 if ((qq >= Base) || (qq * vn [n - 2] > (rr * Base + un [j + n - 2]))) {
2269 rr += (ulong)vn [n - 1];
2277 // Multiply and subtract
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;
2289 t = (long)un [j + n] - b;
2290 un [j + n] = (uint)t;
2292 // Store the calculated value
2296 // Add back vn[0..n] to un[j..j+n]
2301 for (i = 0; i < n; i++) {
2302 c = (ulong)vn [i] + un [j + i] + c;
2303 un [j + i] = (uint)c;
2306 c += (ulong)un [j + n];
2307 un [j + n] = (uint)c;
2311 Unnormalize (un, out r, shift);
2313 q = new uint [] { 0 };