当前位置： 文江博客话题详情

计算 BigInteger 的平方

发布于 2024-09-06 10:36:50 字数 539 浏览 4 评论 0原文

我正在使用 .NET 4 的 System.Numerics.BigInteger 结构。

我需要计算非常大的数字的平方 (x2) - 数百万十进制数字。

如果 x 是 BigInteger，则：

x*x;

或 的

BigInteger.Pow(x,2);

时间复杂度是多少？

如何使用 .NET 4 BigInteger 以最快的方式乘以如此大的数字？是否有 Schönhage–Strassen 算法 的实现？

原文

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?

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（4）

葬花如无物 2024-09-13 10:36:51

这取决于你的人数有多大。正如维基百科页面告诉您的那样：

实际上，对于 2215 到 22< 的数字，Schönhage–Strassen 算法的性能开始优于 Karatsuba 和 Toom–Cook 乘法等旧方法。 support>17（10,000 到 40,000 位十进制数字）。

System.Numerics.BigInteger 使用 Karatsuba 算法 或标准教科书乘法，取决于数字的大小。 Karasuba 的时间复杂度为 O(n log2 3)。但如果您的数字小于上面引用的数字，那么您可能不会看到实施 Schönhage-Strassen 带来的加速。

至于 Pow() ，它本身在计算过程中至少对数字进行一次平方（并且它通过简单地执行 num * num 来实现这一点 - 所以我认为这不会是也更有效率。

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.

回复 原文
违心° 2024-09-13 10:36:51

一种非常简单的实现方法是基于 FFT。由于将两个数字相乘相当于执行其系数的卷积，然后通过一次消除进位，因此您应该能够通过 FFT 方法（n = 位数）以 O(n log n) 运算进行卷积。

数值食谱有一章介绍这一点。对于如此大的数字来说，这绝对比分而治之的方法更快，比如唐叶。

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.

回复 原文
乱了心跳 2024-09-13 10:36:51

  • 首先，System.Numerics.BigInteger 不使用 [Karatsuba
    算法]，O(n 0.5 )，它使用标准教科书乘法 O(n 2 )。

  • 通过此代码，您可以在短短 1.4 毫秒内乘以两个 30,000 位（大约 9000 个十进制数字）。

     public void benchMark()
     {
         Xint U、V、温度；
         整数 n = 30000；
         而（n>29000）
         {
             U = RND(n << 1);
             //_______________________
             sw.Restart();
             温度=U*U；
             sw.Stop();
             label7.Text = Convert.ToString("Micro " + sw.Elapsed.TotalMilliseconds + " ms");
             //_______________________

         }
      n>>=1；
     }

     公共静态Xint MTP（Xint U，Xint V）
     {
         return MTP(U, V, Xint.Max(U.Sign * U, V.Sign * V).ToByteArray().Length << 3);
     }
     公共静态Xint MTP（Xint U，Xint V，int n）
     {
         如果 (n <= 3000) 返回 U * V；
         if (n <= 6000) 返回 TC2(U, V, n);
         if (n <= 10000) 返回 TC3(U, V, n);
         if (n <= 40000) 返回 TC4(U, V, n);
         返回TC2P（U，V，n）；
     }
     私有静态 Xint MTPr(Xint U, Xint V, int n)
     {
         如果 (n <= 3000) 返回 U * V；
         if (n <= 6000) 返回 TC2(U, V, n);
         if (n <= 10000) 返回 TC3(U, V, n);
         返回TC4（U，V，n）；
     }
     私有静态 Xint TC2(Xint U1, Xint V1, int n)
     {
         n>>=1；
         Xint Mask = (Xint.One << n) - 1;
         Xint U0 = U1 &面具; U1>>=n；
         Xint V0 = V1 &面具; 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;
     }
     私有静态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 &面具; U2>>=n；
         Xint U1 = U2 &面具; U2>>=n；
         Xint V0 = V2 &面具; V2>>=n；
         Xint V1 = V2 &面具; 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);
         返回 ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
     }
     私有静态Xint TC4（Xint U3，Xint V3，int n）
     {
         n>>=2；
         Xint Mask = (Xint.One << n) - 1;
         Xint U0 = U3 &面具; U3>>=n；
         Xint U1 = U3 &面具; U3>>=n；
         Xint U2 = U3 &面具; U3>>=n；
         Xint V0 = V3 &面具; V3>>=n；
         Xint V1 = V3 &面具; V3>>=n；
         Xint V2 = V3 &面具; 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); // 信息

         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)；
         返回 ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << ＜n)+W0；
     }
     私有静态Xint TC2P（Xint A，Xint B，int n）
     {
         n>>=1；
         Xint Mask = (Xint.One << n) - 1;
         Xint[] U = 新Xint[3];
         U[0]=A&面具; A>>=n； U[2] = A； U[1] = U[0] + A；
         Xint[] V = 新Xint[3];
         V[0]=B&面具; B>>=n； V[2] = B； V[1] = V[0] + B；
         Xint[] P = 新Xint[3];
         Parallel.For(0, 3, (int i) => P[i] = MTPr(U[i], V[i], n));
         返回 ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0];
     }

     私有静态长current_n；
    私有静态Xint产品（int n）
     {
         如果 (n > 2) 返回 MTP(Product(n - (n >>= 1)), Product(n));
         if (n > 1) return (current_n += 2) * (current_n += 2);
         返回当前_n += 2；
     }

     公共静态Xint [] DQR（Xint A，Xint B）
     {
         int n = bL(B);
         int m = bL(A) - n;
         if (m <= 6000) 返回 SmallDivRem(A, B);
         int signA = A.Sign; A *= 符号A;
         int SignB = B.Sign; B *= 符号B；
         Xint[] QR = 新Xint[2];
         如果 (m <= n) QR = D21(A, B, n)；
         别的
         {
             Xint Mask = (Xint.One << n) - 1;
             米/=n；
             Xint[] U = 新Xint[m];
             整数 i = 0;
             对于 (; i < m; i++)
             {
                 U[i] = A &面具;
                 A>>=n；
             }
             QR = D21(A, B, n);
             A = QR[0];
             对于 (i--; i >= 0; i--)
             {
                 QR = D21(QR[1] << n | U[i], B, n);
                 A = A << n | QR[0]；
             }
             QR[0] = A；
         }
         QR[0] *= 符号A * 符号B;
         QR[1] *= 符号A；
         返回二维码；
     }
     私有静态 Xint[] SmallDivRem(Xint A, Xint B)
     {
         Xint[] QR = 新Xint[2];
         QR[0] = Xint.DivRem(A, B, 输出 QR[1]);
         返回二维码；
     }
     私有静态 Xint[] D21(Xint A, Xint B, int n)
     {
         if (n <= 6000) 返回 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&面具;
         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；
         返回QR2；
     }
     私有静态Xint [] D32（Xint A12，Xint A3，Xint B，Xint B1，Xint B2，int n）
     {
         Xint[] QR = 新Xint[2];
         if ((A12 >> n) != B1) QR = D21(A12, B1, n);
         别的
         {
             QR[0] = (Xint.One << n) - 1;
             QR[1]=A12+B1-(B1＜＜n)；
         }
         QR[1]=(QR[1]＜＜n|A3)-MTP(QR[0]，B2，n)；
         而（QR[1] < 0）
         {
             QR[0] -= 1；
             QR[1] += B；
         }
         返回二维码；
     }

     公共静态Xint SQ(Xint U)
     {
         return SQ(U, U.Sign * U.ToByteArray().Length << 3);
     }
     公共静态Xint SQ（Xint U，int n）
     {
         如果 (n <= 700) 返回 U * U；
         if (n <= 3000) 返回 Xint.Pow(U, 2);
         if (n <= 6000) 返回 SQ2(U, n);
         if (n <= 10000) 返回 SQ3(U, n);
         if (n <= 40000) 返回 SQ4(U, n);
         返回 SQ2P(U, n);
     }
     私有静态 Xint SQr(Xint U, int n)
     {
         if (n <= 3000) 返回 Xint.Pow(U, 2);
         if (n <= 6000) 返回 SQ2(U, n);
         if (n <= 10000) 返回 SQ3(U, n);
         返回 SQ4(U, n);
     }
     私有静态 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);
         返回 ((P2 << n) + (SQr(U0 + U1, n) - (P0 + P2)) << n) + P0;
     }
     私有静态 Xint SQ3(Xint U2, int n)
     {
         n = (int)((long)(n) * 0x55555556 >> 32);
         Xint Mask = (Xint.One << n) - 1;
         Xint U0 = U2 &面具; U2>>=n；
         Xint U1 = U2 &面具; 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);
         返回 ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
     }
     私有静态 Xint SQ4(Xint U3, int n)
     {
         n>>=2；
         Xint Mask = (Xint.One << n) - 1;
         Xint U0 = U3 &面具; U3>>=n；
         Xint U1 = U3 &面具; U3>>=n；
         Xint U2 = U3 &面具; 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); // 信息
         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)；
         返回 ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << ＜n)+W0；
     }
     私有静态Xint SQ2P（Xint A，int n）
     {
         n>>=1；
         Xint[] U = 新Xint[3];
         U[0]=A& ((Xint.One << n) - 1); A>>=n； U[2] = A； U[1] = U[0] + A；
         Xint[] P = 新Xint[3];
         Parallel.For(0, 3, (int i) => P[i] = SQr(U[i], n));
         返回 ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0];
     }

     公共静态Xint POW（Xint X，int y）
     {
         if (y > 1) return ((y & 1) == 1) ? SQ(POW(X, y >> 1)) * X : SQ(POW(X, y >> 1));
         如果（y==1）返回X；
         返回1；
     }

     公共静态Xint[] SR(Xint A)
     {
         返回 SR(A, bL(A));
     }
     公共静态Xint [] SR（Xint A，int n）
     {
         如果 (n < 53) 返回 SR52(A)；
         int m = n >>> 2；
         Xint Mask = (Xint.One << m) - 1;
         Xint A0=A&面具; A>>=m；
         Xint A1=A&面具; 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)；
         如果（R[1] < 0）
         {
             R[0] -= 1；
             R[1] += (R[0] << 1) | 1;
         }
         返回R；
     }
     私有静态Xint[] SR52(Xint A)
     {
         双a = (双)A;
         长 q = (长)Math.Sqrt(a);
         长 r = (长)(a) - q * q;
         Xint[] QR = { q, r };
         返回二维码；
     }

     公共静态Xint SRO（Xint A）
     {
         返回 SRO(A, bL(A));
     }
     公共静态Xint SRO（Xint A，int n）
     {
         if (n < 53) return (int)Math.Sqrt((double)A);
         Xint[] R = SROr(A, n, 1);
         返回R[0]；
     }
     private static Xint[] SROr(Xint A, int n, int rc) // rc=递归计数器
     {
         如果 (n < 53) 返回 SR52(A)；
         int m = n >>> 2；
         Xint Mask = (Xint.One << m) - 1;
         Xint A0 = A&面具; A>>=m；
         Xint A1=A&面具; 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--;
         如果（rc！= 0）
         {
             R[1]=(D[1]＜＜m|A0)-SQ(D[0]，m)；
             如果（R[1] < 0）
             {
                 R[0] -= 1；
                 R[1] += (R[0] << 1) | 1;
             }
             返回R；
         }
         n = (bL(D[0]) << 1) - bL(D[1] << m | A0);
         如果 (n < 0) 返回 R；
         如果（n>1）
         {
             R[0] -= 1；
             返回R；
         }
         int 移位 = (bL(D[0]) - 31) << 1;
         long d0 = (int)(D[0] >> (shift >> 1));
         长 d = (长)((D[1] >> (移位 - m)) | (A0 >> 移位)) - d0 * d0;
         如果 (d < 0)
         {
             R[0] -= 1；
             返回R；
         }
         如果(d＞d0＜＜1)返回R；
         R[0]-＝(1-(((D[1]＜＜m)|A0)-SQ(D[0]，m)).Sign)>> 1;
         返回R；
     }

     公共静态 int bL(Xint U)
     {
         byte[] 字节 = (U.Sign * U).ToByteArray();
         int i = bytes.Length - 1;
         return (i << 3) + bitLengthMostSignificantByte(bytes[i]);
     }
     私有静态 int bitLengthMostSignificantByte(字节 b)
     {
         返回 b < 08 ？ b＜ 02 ? b＜ 01 ？ 0：1：
                                  b＜ 04 ? 2：3：
                         b＜ 32？ b＜ 16？ 4：5：
                                  b＜ 64？ 6：7；
     }

     公共静态 int fL2(int i)
     {
         返回
         我< 1<<; 15？我< 1<<; 07 ?我< 1<<; 03 ?我< 1<<; 01 ？我< 1<<; 00 ? -1：00：
                                                                 我< 1<<; 02 ? 01：02：
                                                   我< 1<<; 05 ?我< 1<<; 04 ? 03：04：
                                                                 我< 1<<; 06 ？ 05 : 06 :
                                     我< 1<<; 11？我< 1<<; 09 ?我< 1<<; 08 ？ 07:08:
                                                                 我< 1<<; 10？ 09:10:
                                                   我< 1<<; 13？我< 1<<; 12？ 11:12:
                                                                 我< 1<<； 14？ 13：14：
                       我< 1<<; 23？我< 1<<; 19 ？我< 1<<; 17 ？我< 1<<; 16？ 15：16：
                                                                 我< 1<<; 18？ 17：18：
                                                   我< 1<<; 21？我< 1<<; 20？ 19:20:
                                                                 我< 1<<; 22？ 21:22:
                                     我< 1<<; 27？我< 1<<; 25？我< 1<<; 24？ 23：24：
                                                                 我< 1<<; 26？ 25：26：
                                                   我< 1<<; 29？我< 1<<; 28？ 27：28：
                                                                 我< 1<<; 30？ 29:30；
     }

     私有静态 int 种子；
     公共静态Xint RND（int n）
     {
         如果 (n < 2) 返回 n；
         if (seed == int.MaxValue)seed = 0;否则种子++；
         随机兰特=新随机（种子）；
         byte[] 字节 = 新字节[(n + 15)>> 3];
         rand.NextBytes(字节);
         int i = bytes.Length - 1;
         字节[i] = 0;
         n = (i << 3) - n；
         我 - ;
         字节[i]>>=n；
         bytes[i] |= (字节)(128 >> n);
         返回新的Xint（字节）；
     }
 //++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++


  }
}

  • First of All,System.Numerics.BigInteger does 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);
     }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++


  }
}
回复 原文
贪了杯 2024-09-13 10:36:51

您可以使用 GNU MP Bignum 库的 C# 包装器，它的速度可能与您一样快可以得到。对于纯 C# 实现，您可以尝试 IntX。

最快的乘法算法实际上是Fürer的算法，但我还没有找到它的任何实现。

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.

回复 原文
~没有更多了~

关于作者

酒浓于脸红

暂无简介

0 文章
0 评论
24 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

lioqio

文章 0 评论 0

Single

文章 0 评论 0

禾厶谷欠

文章 0 评论 0

alipaysp_2zg8elfGgC

文章 0 评论 0

qq_N6d4X7

文章 0 评论 0

放低过去

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文