brickner
Reputation: 6585
Calculating the square of BigInteger
I'm using .NET 4's System.Numerics.BigInteger structure.
I need to calculate the square (x2) of very large numbers - millions of decimal digits.
If x is a BigInteger, What is the time complexity of:
x*x;
or
BigInteger.Pow(x,2);
?
How can multiply such big numbers in the fastest way using .NET 4 BigInteger? Is there an implementation for Schönhage–Strassen algorithm?
Upvotes: 5
Views: 3530
Answers (4)
saeid
Reputation: 21
First of All,
System.Numerics.BigIntegerdoes NOT use the [Karatsuba algorithm] with O(n 0.5 ) and it uses standard schoolbook multiplication O(n 2 ).By this code you can Multiple two 30,000 bit (approximately 9000 decimal digit) in Just 1.4 millisecond.
public void benchMark() { Xint U, V,Temp; int n = 30000; while (n > 29000) { U = RND(n << 1); //_______________________ sw.Restart(); Temp = U * U; sw.Stop(); label7.Text = Convert.ToString("Micro " + sw.Elapsed.TotalMilliseconds + " ms"); //_______________________ } n>>=1; } public static Xint MTP(Xint U, Xint V) { return MTP(U, V, Xint.Max(U.Sign * U, V.Sign * V).ToByteArray().Length << 3); } public static Xint MTP(Xint U, Xint V, int n) { if (n <= 3000) return U * V; if (n <= 6000) return TC2(U, V, n); if (n <= 10000) return TC3(U, V, n); if (n <= 40000) return TC4(U, V, n); return TC2P(U, V, n); } private static Xint MTPr(Xint U, Xint V, int n) { if (n <= 3000) return U * V; if (n <= 6000) return TC2(U, V, n); if (n <= 10000) return TC3(U, V, n); return TC4(U, V, n); } private static Xint TC2(Xint U1, Xint V1, int n) { n >>= 1; Xint Mask = (Xint.One << n) - 1; Xint U0 = U1 & Mask; U1 >>= n; Xint V0 = V1 & Mask; V1 >>= n; Xint P0 = MTPr(U0, V0, n); Xint P2 = MTPr(U1, V1, n); return ((P2 << n) + (MTPr(U0 + U1, V0 + V1, n) - (P0 + P2)) << n) + P0; } private static Xint TC3(Xint U2, Xint V2, int n) { n = (int)((long)(n) * 0x55555556 >> 32); // n /= 3; Xint Mask = (Xint.One << n) - 1; Xint U0 = U2 & Mask; U2 >>= n; Xint U1 = U2 & Mask; U2 >>= n; Xint V0 = V2 & Mask; V2 >>= n; Xint V1 = V2 & Mask; V2 >>= n; Xint W0 = MTPr(U0, V0, n); Xint W4 = MTPr(U2, V2, n); Xint P3 = MTPr((((U2 << 1) + U1) << 1) + U0, (((V2 << 1) + V1 << 1)) + V0, n); U2 += U0; V2 += V0; Xint P2 = MTPr(U2 - U1, V2 - V1, n); Xint P1 = MTPr(U2 + U1, V2 + V1, n); Xint W2 = (P1 + P2 >> 1) - (W0 + W4); Xint W3 = W0 - P1; W3 = ((W3 + P3 - P2 >> 1) + W3) / 3 - (W4 << 1); Xint W1 = P1 - (W4 + W3 + W2 + W0); return ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0; } private static Xint TC4(Xint U3, Xint V3, int n) { n >>= 2; Xint Mask = (Xint.One << n) - 1; Xint U0 = U3 & Mask; U3 >>= n; Xint U1 = U3 & Mask; U3 >>= n; Xint U2 = U3 & Mask; U3 >>= n; Xint V0 = V3 & Mask; V3 >>= n; Xint V1 = V3 & Mask; V3 >>= n; Xint V2 = V3 & Mask; V3 >>= n; Xint W0 = MTPr(U0, V0, n); // 0 U0 += U2; U1 += U3; V0 += V2; V1 += V3; Xint P1 = MTPr(U0 + U1, V0 + V1, n); // 1 Xint P2 = MTPr(U0 - U1, V0 - V1, n); // -1 U0 += 3 * U2; U1 += 3 * U3; V0 += 3 * V2; V1 += 3 * V3; Xint P3 = MTPr(U0 + (U1 << 1), V0 + (V1 << 1), n); // 2 Xint P4 = MTPr(U0 - (U1 << 1), V0 - (V1 << 1), n); // -2 Xint P5 = MTPr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2), V0 + 12 * V2 + ((V1 + 12 * V3) << 2), n); // 4 Xint W6 = MTPr(U3, V3, n); // inf Xint W1 = P1 + P2; Xint W4 = (((((P3 + P4) >> 1) - (W1 << 1)) / 3 + W0) >> 2) - 5 * W6; Xint W2 = (W1 >> 1) - (W6 + W4 + W0); P1 = P1 - P2; P4 = P4 - P3; Xint W5 = ((P1 >> 1) + (5 * P4 + P5 - W0 >> 4) - ((((W6 << 4) + W4) << 4) + W2)) / 45; W1 = ((P4 >> 2) + (P1 << 1)) / 3 + (W5 << 2); Xint W3 = (P1 >> 1) - (W1 + W5); return ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0; } private static Xint TC2P(Xint A, Xint B, int n) { n >>= 1; Xint Mask = (Xint.One << n) - 1; Xint[] U = new Xint[3]; U[0] = A & Mask; A >>= n; U[2] = A; U[1] = U[0] + A; Xint[] V = new Xint[3]; V[0] = B & Mask; B >>= n; V[2] = B; V[1] = V[0] + B; Xint[] P = new Xint[3]; Parallel.For(0, 3, (int i) => P[i] = MTPr(U[i], V[i], n)); return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0]; } private static long current_n; private static Xint Product(int n) { if (n > 2) return MTP(Product(n - (n >>= 1)), Product(n)); if (n > 1) return (current_n += 2) * (current_n += 2); return current_n += 2; } public static Xint[] DQR(Xint A, Xint B) { int n = bL(B); int m = bL(A) - n; if (m <= 6000) return SmallDivRem(A, B); int signA = A.Sign; A *= signA; int signB = B.Sign; B *= signB; Xint[] QR = new Xint[2]; if (m <= n) QR = D21(A, B, n); else { Xint Mask = (Xint.One << n) - 1; m /= n; Xint[] U = new Xint[m]; int i = 0; for (; i < m; i++) { U[i] = A & Mask; A >>= n; } QR = D21(A, B, n); A = QR[0]; for (i--; i >= 0; i--) { QR = D21(QR[1] << n | U[i], B, n); A = A << n | QR[0]; } QR[0] = A; } QR[0] *= signA * signB; QR[1] *= signA; return QR; } private static Xint[] SmallDivRem(Xint A, Xint B) { Xint[] QR = new Xint[2]; QR[0] = Xint.DivRem(A, B, out QR[1]); return QR; } private static Xint[] D21(Xint A, Xint B, int n) { if (n <= 6000) return SmallDivRem(A, B); int s = n & 1; A <<= s; B <<= s; n++; n >>= 1; Xint Mask = (Xint.One << n) - 1; Xint B1 = B >> n; Xint B2 = B & Mask; Xint[] QR1 = D32(A >> (n << 1), A >> n & Mask, B, B1, B2, n); Xint[] QR2 = D32(QR1[1], A & Mask, B, B1, B2, n); QR2[0] |= QR1[0] << n; QR2[1] >>= s; return QR2; } private static Xint[] D32(Xint A12, Xint A3, Xint B, Xint B1, Xint B2, int n) { Xint[] QR = new Xint[2]; if ((A12 >> n) != B1) QR = D21(A12, B1, n); else { QR[0] = (Xint.One << n) - 1; QR[1] = A12 + B1 - (B1 << n); } QR[1] = (QR[1] << n | A3) - MTP(QR[0], B2, n); while (QR[1] < 0) { QR[0] -= 1; QR[1] += B; } return QR; } public static Xint SQ(Xint U) { return SQ(U, U.Sign * U.ToByteArray().Length << 3); } public static Xint SQ(Xint U, int n) { if (n <= 700) return U * U; if (n <= 3000) return Xint.Pow(U, 2); if (n <= 6000) return SQ2(U, n); if (n <= 10000) return SQ3(U, n); if (n <= 40000) return SQ4(U, n); return SQ2P(U, n); } private static Xint SQr(Xint U, int n) { if (n <= 3000) return Xint.Pow(U, 2); if (n <= 6000) return SQ2(U, n); if (n <= 10000) return SQ3(U, n); return SQ4(U, n); } private static Xint SQ2(Xint U1, int n) { n >>= 1; Xint U0 = U1 & ((Xint.One << n) - 1); U1 >>= n; Xint P0 = SQr(U0, n); Xint P2 = SQr(U1, n); return ((P2 << n) + (SQr(U0 + U1, n) - (P0 + P2)) << n) + P0; } private static Xint SQ3(Xint U2, int n) { n = (int)((long)(n) * 0x55555556 >> 32); Xint Mask = (Xint.One << n) - 1; Xint U0 = U2 & Mask; U2 >>= n; Xint U1 = U2 & Mask; U2 >>= n; Xint W0 = SQr(U0, n); Xint W4 = SQr(U2, n); Xint P3 = SQr((((U2 << 1) + U1) << 1) + U0, n); U2 += U0; Xint P2 = SQr(U2 - U1, n); Xint P1 = SQr(U2 + U1, n); Xint W2 = (P1 + P2 >> 1) - (W0 + W4); Xint W3 = W0 - P1; W3 = ((W3 + P3 - P2 >> 1) + W3) / 3 - (W4 << 1); Xint W1 = P1 - (W4 + W3 + W2 + W0); return ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0; } private static Xint SQ4(Xint U3, int n) { n >>= 2; Xint Mask = (Xint.One << n) - 1; Xint U0 = U3 & Mask; U3 >>= n; Xint U1 = U3 & Mask; U3 >>= n; Xint U2 = U3 & Mask; U3 >>= n; Xint W0 = SQr(U0, n); // 0 U0 += U2; U1 += U3; Xint P1 = SQr(U0 + U1, n); // 1 Xint P2 = SQr(U0 - U1, n); // -1 U0 += 3 * U2; U1 += 3 * U3; Xint P3 = SQr(U0 + (U1 << 1), n); // 2 Xint P4 = SQr(U0 - (U1 << 1), n); // -2 Xint P5 = SQr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2), n); // 4 Xint W6 = SQr(U3, n); // inf Xint W1 = P1 + P2; Xint W4 = (((((P3 + P4) >> 1) - (W1 << 1)) / 3 + W0) >> 2) - 5 * W6; Xint W2 = (W1 >> 1) - (W6 + W4 + W0); P1 = P1 - P2; P4 = P4 - P3; Xint W5 = ((P1 >> 1) + (5 * P4 + P5 - W0 >> 4) - ((((W6 << 4) + W4) << 4) + W2)) / 45; W1 = ((P4 >> 2) + (P1 << 1)) / 3 + (W5 << 2); Xint W3 = (P1 >> 1) - (W1 + W5); return ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0; } private static Xint SQ2P(Xint A, int n) { n >>= 1; Xint[] U = new Xint[3]; U[0] = A & ((Xint.One << n) - 1); A >>= n; U[2] = A; U[1] = U[0] + A; Xint[] P = new Xint[3]; Parallel.For(0, 3, (int i) => P[i] = SQr(U[i], n)); return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0]; } public static Xint POW(Xint X, int y) { if (y > 1) return ((y & 1) == 1) ? SQ(POW(X, y >> 1)) * X : SQ(POW(X, y >> 1)); if (y == 1) return X; return 1; } public static Xint[] SR(Xint A) { return SR(A, bL(A)); } public static Xint[] SR(Xint A, int n) { if (n < 53) return SR52(A); int m = n >> 2; Xint Mask = (Xint.One << m) - 1; Xint A0 = A & Mask; A >>= m; Xint A1 = A & Mask; A >>= m; Xint[] R = SR(A, n - (m << 1)); Xint[] D = DQR(R[1] << m | A1, R[0] << 1); R[0] = (R[0] << m) + D[0]; R[1] = (D[1] << m | A0) - SQ(D[0], m); if (R[1] < 0) { R[0] -= 1; R[1] += (R[0] << 1) | 1; } return R; } private static Xint[] SR52(Xint A) { double a = (double)A; long q = (long)Math.Sqrt(a); long r = (long)(a) - q * q; Xint[] QR = { q, r }; return QR; } public static Xint SRO(Xint A) { return SRO(A, bL(A)); } public static Xint SRO(Xint A, int n) { if (n < 53) return (int)Math.Sqrt((double)A); Xint[] R = SROr(A, n, 1); return R[0]; } private static Xint[] SROr(Xint A, int n, int rc) // rc=recursion counter { if (n < 53) return SR52(A); int m = n >> 2; Xint Mask = (Xint.One << m) - 1; Xint A0 = A & Mask; A >>= m; Xint A1 = A & Mask; A >>= m; Xint[] R = SROr(A, n - (m << 1), rc + 1); Xint[] D = DQR((R[1] << m) | A1, R[0] << 1); R[0] = (R[0] << m) + D[0]; rc--; if (rc != 0) { R[1] = (D[1] << m | A0) - SQ(D[0], m); if (R[1] < 0) { R[0] -= 1; R[1] += (R[0] << 1) | 1; } return R; } n = (bL(D[0]) << 1) - bL(D[1] << m | A0); if (n < 0) return R; if (n > 1) { R[0] -= 1; return R; } int shift = (bL(D[0]) - 31) << 1; long d0 = (int)(D[0] >> (shift >> 1)); long d = (long)((D[1] >> (shift - m)) | (A0 >> shift)) - d0 * d0; if (d < 0) { R[0] -= 1; return R; } if (d > d0 << 1) return R; R[0] -= (1 - (((D[1] << m) | A0) - SQ(D[0], m)).Sign) >> 1; return R; } public static int bL(Xint U) { byte[] bytes = (U.Sign * U).ToByteArray(); int i = bytes.Length - 1; return (i << 3) + bitLengthMostSignificantByte(bytes[i]); } private static int bitLengthMostSignificantByte(byte b) { return b < 08 ? b < 02 ? b < 01 ? 0 : 1 : b < 04 ? 2 : 3 : b < 32 ? b < 16 ? 4 : 5 : b < 64 ? 6 : 7; } public static int fL2(int i) { return i < 1 << 15 ? i < 1 << 07 ? i < 1 << 03 ? i < 1 << 01 ? i < 1 << 00 ? -1 : 00 : i < 1 << 02 ? 01 : 02 : i < 1 << 05 ? i < 1 << 04 ? 03 : 04 : i < 1 << 06 ? 05 : 06 : i < 1 << 11 ? i < 1 << 09 ? i < 1 << 08 ? 07 : 08 : i < 1 << 10 ? 09 : 10 : i < 1 << 13 ? i < 1 << 12 ? 11 : 12 : i < 1 << 14 ? 13 : 14 : i < 1 << 23 ? i < 1 << 19 ? i < 1 << 17 ? i < 1 << 16 ? 15 : 16 : i < 1 << 18 ? 17 : 18 : i < 1 << 21 ? i < 1 << 20 ? 19 : 20 : i < 1 << 22 ? 21 : 22 : i < 1 << 27 ? i < 1 << 25 ? i < 1 << 24 ? 23 : 24 : i < 1 << 26 ? 25 : 26 : i < 1 << 29 ? i < 1 << 28 ? 27 : 28 : i < 1 << 30 ? 29 : 30; } private static int seed; public static Xint RND(int n) { if (n < 2) return n; if (seed == int.MaxValue) seed = 0; else seed++; Random rand = new Random(seed); byte[] bytes = new byte[(n + 15) >> 3]; rand.NextBytes(bytes); int i = bytes.Length - 1; bytes[i] = 0; n = (i << 3) - n; i--; bytes[i] >>= n; bytes[i] |= (byte)(128 >> n); return new Xint(bytes); } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++ } }
Upvotes: 2
Alexandre C.
Reputation: 56976
A quite simple method to implement is based on FFT. Since multiplying two numbers amounts to perform a convolution of their coefficients followed by a pass to eliminate carries, you should be able to do the convolution in O(n log n) operations via FFT methods (n = number of digits).
Numerical recipes has a chapter on this. This is definitely faster than divide and conquer methods, like Karatsuba, for such big numbers.
Upvotes: 2
Cloudanger
Reputation: 9524
You can use a C# wrapper for GNU MP Bignum Library, which is probably as fast as you can get. For pure C# implementation you could try IntX.
Fastest multiplication algorithm is actually Fürer's algorithm, but I haven't found any implementations for it.
Upvotes: 1
Joey
Reputation: 354714
That depends on how large your numbers are. As the Wikipedia page tells you:
In practice the Schönhage–Strassen algorithm starts to outperform older methods such as Karatsuba and Toom–Cook multiplication for numbers beyond 2215 to 2217 (10,000 to 40,000 decimal digits).
System.Numerics.BigInteger uses the Karatsuba algorithm or standard schoolbook multiplication, depending on the size of the numbers. Karatsuba has a time complexity of O(n log2 3). But if your numbers are smaller than above quoted figures, then you likely won't see much speedup from implementing Schönhage–Strassen.
As for Pow() this itself squares the number at least once during its calculations (and it does that by simply doing num * num – so I think this won't be more efficient, either.
Upvotes: 7
Related Questions
- Calculate square root of a BigInteger (System.Numerics.BigInteger)
- BigInteger Power calculation with double base C#
- BigInteger source code
- Complexity of checking whether a BigInteger is a perfect square
- How to solve BigInteger incorrect result?
- BigInteger.Pow(BigInteger, BigInteger)?
- Multiplying BigInteger by Double
- Calculating very large whole numbers
- how to get the BigInteger to the pow Double in C#?
- Working with BigIntegers in C#