代码之家  ›  专栏  ›  技术社区  ›  brickner

计算大整数的平方

biginteger complexity-theory .net-4.0 math c#
  •  5
  • brickner  · 技术社区  · 14 年前

    我用的是.NET 4 System.Numerics.BigInteger structure .

    我需要计算平方(x 二 )数量非常大- 数以百万计的十进制数字 .

    如果 x 是一个 BigInteger ,什么是时间复杂性: 

    x*x;

    或 

    BigInteger.Pow(x,2);

    ?

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

    5 回复  |  直到 10 年前
        1
  •  7
  •   Joey Gumbo    14 年前

    这取决于你的数字有多大。正如维基百科网页上所说:

    在实践中,sch_¶nhage_ strassen算法开始优于旧方法,如karatsuba和toom_ cook multiplex for numbers beyond 2。 二 十五 到2 二 十七 (10000到40000位小数)。

    System.Numerics.BigInteger 使用 Karatsuba algorithm 或者标准的课本乘法,取决于数字的大小。Karatsuba的时间复杂度为O( n γ测井 二 α3 )但如果你的数字小于以上引用的数字,那么你可能不会看到实施sch_nhage_strassen的速度加快。

    至于 Pow() 在计算的过程中,这个数字本身至少平方一次(它通过简单地做 num * num 所以我认为这也不会更有效率。

        2
  •  3
  •   Tim Buktu    13 年前

    这里有一个sch_¶nhage-strassen实现:

    https://github.com/tbuktu/ntru/blob/master/src/main/java/net/sf/ntru/arith/Sch%C3%B6nhageStrassen.java

        3
  •  2
  •   Alexandre C.    14 年前

    一个非常简单的实现方法是基于FFT的。由于将两个数字相乘,对其系数进行卷积,然后通过一次消除进位,因此您应该能够通过FFT方法(n=位数)在O(n log n)操作中进行卷积。

    数字食谱有一章是关于这个的。对于如此庞大的数字来说,这绝对比像空竹这样的分而治之的方法要快。

        4
  •  2
  •   saeid    10 年前
    • 首先, System.Numerics.BigInteger 不使用[空竹] 算法]带O(N) 零点五 )它使用标准的教科书乘法 二 )

    • 通过此代码,您可以在1.4毫秒内多出两个30000位(约9000位十进制数字)。 

          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);
    }
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      } }

        5
  •  1
  •   Cloudanger    14 年前

    你可以使用 C# wrapper 对于GNU MP Bignum库,这可能是最快的。对于纯C实现,您可以尝试 IntX .

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

    推荐文章
    zadspecial  ·  为什么大人物们不加入呢?
    6 年前
    user3596335  ·  如何计算大整数以N为底的对数?
    6 年前
    Ras  ·  如何将解密后的BigInteger转换回字符串?
    6 年前
    pacoverflow  ·  为什么Perl模运算符使用bignum而不是Math::BigInt处理大整数?
    6 年前
    babakahn  ·  如何重载这个使用数组的bigint类的+运算符?C++
    6 年前
    Aaron  ·  Python v 3.5.2和v 2.7.12中的结果不同,但v 2.7.12是否正确?
    7 年前
    datv  ·  如何将两个BigInteger变量作为无符号变量进行比较?
    7 年前
    krltos  ·  在使用Hibernate查询填充的列表中循环时发生ClassCastException
    7 年前
    Darth7urtle  ·  如何创建一个新的BigInteger变量并生成一个返回BigInteger的方法?[副本]
    7 年前
    ArkhamWarfare  ·  Java Fibonocci序列BigInteger
    7 年前
    关于   移动版
    代码之家 - 一站式码农服务社区
    沪ICP备11025650号