c++ - BigInteger 数字实现和性能

Question

我用 C++ 编写了一个 BigInteger 类，它应该能够对任何大小的所有数字进行运算。目前，我正在尝试通过比较现有算法来实现一种非常快速的乘法方法，并测试它们最适合的数字数量，但我遇到了非常意想不到的结果。我尝试对 500 位数字进行 20 次乘法运算，并为它们计时。这是结果:

karatsuba:
  14.178 seconds

long multiplication:
  0.879 seconds

维基百科告诉我

It follows that, for sufficiently large n, Karatsuba's algorithm will perform fewer shifts and single-digit additions than longhand multiplication, even though its basic step uses more additions and shifts than the straightforward formula. For small values of n, however, the extra shift and add operations may make it run slower than the longhand method. The point of positive return depends on the computer platform and context. As a rule of thumb, Karatsuba is usually faster when the multiplicands are longer than 320–640 bits.

因为我的数字至少有 1500 位长，所以这是非常出乎意料的，因为维基百科说 karatsuba 应该运行得更快。我相信我的问题可能出在我的加法算法上，但我不知道如何让它更快，因为它已经在 O(n) 中运行了。我将在下面发布我的代码，以便您可以检查我的实现。我会把不相关的部分排除在外。
我也在想，也许我使用的结构不是最好的。我用小端表示每个数据段。因此，例如，如果我将数字 123456789101112 存储到长度为 3 的数据段中，它将如下所示:

{112,101,789,456,123}

所以这就是为什么我现在要问实现 BigInteger 类的最佳结构和最佳方法是什么？为什么 karatsuba 算法比长乘法算法慢？

这是我的代码:(我很抱歉太长了)

using namespace std;

bool __longmult=true;
bool __karatsuba=false;

struct BigInt {
public:
    vector <int> digits;

    BigInt(const char * number) {
        //constructor is not relevant   
    }
    BigInt() {}

    void BigInt::operator = (BigInt a) {
        digits=a.digits;
    }

    friend BigInt operator + (BigInt,BigInt);
    friend BigInt operator * (BigInt,BigInt);

    friend ostream& operator << (ostream&,BigInt);
};

BigInt operator + (BigInt a,BigInt b) {
    if (b.digits.size()>a.digits.size()) {
        a.digits.swap(b.digits); //make sure a has more or equal amount of digits than b
    }
    int carry=0;

    for (unsigned int i=0;i<a.digits.size();i++) {
        int sum;
        if (i<b.digits.size()) {
            sum=b.digits[i]+a.digits[i]+carry;
        } else if (carry==1) {
            sum=a.digits[i]+carry;
        } else {
            break; // if carry is 0 and no more digits in b are left then we are done already
        }

        if (sum>=1000000000) {
            a.digits[i]=sum-1000000000;
            carry=1;
        } else {
            a.digits[i]=sum;
            carry=0;
        }
    }

    if (carry) {
        a.digits.push_back(1);
    }

    return a;
}

BigInt operator * (BigInt a,BigInt b) {
    if (__longmult) {
        BigInt res;
        for (unsigned int i=0;i<b.digits.size();i++) {
            BigInt temp;
            temp.digits.insert(temp.digits.end(),i,0); //shift to left for i 'digits'

            int carry=0;
            for (unsigned int j=0;j<a.digits.size();j++) {
                long long prod=b.digits[i];
                prod*=a.digits[j];
                prod+=carry;
                int t=prod%1000000000;
                temp.digits.push_back(t);
                carry=(prod-t)/1000000000;
            }
            if (carry>0) {
                temp.digits.push_back(carry);
            }
            res+=temp;
        }
        return res;
    } else if (__karatsuba) {
        BigInt res;
        BigInt a1,a0,b1,b0;
        assert(a.digits.size()>0 && b.digits.size()>0);
        while (a.digits.size()!=b.digits.size()) { //add zeroes for equal size
            if (a.digits.size()>b.digits.size()) {
                b.digits.push_back(0);
            } else {
                a.digits.push_back(0);
            }
        }

        if (a.digits.size()==1) {
            long long prod=a.digits[0];
            prod*=b.digits[0];

            res=prod;//conversion from long long to BigInt runs in constant time
            return res;

        } else {
            for (unsigned int i=0;i<a.digits.size();i++) {
                if (i<(a.digits.size()+(a.digits.size()&1))/2) { //split the number in 2 equal parts
                    a0.digits.push_back(a.digits[i]);
                    b0.digits.push_back(b.digits[i]);
                } else {
                    a1.digits.push_back(a.digits[i]);
                    b1.digits.push_back(b.digits[i]);
                }
            }
        }

        BigInt z2=a1*b1;
        BigInt z0=a0*b0;
        BigInt z1 = (a1 + a0)*(b1 + b0) - z2 - z0;

        if (z2==0 && z1==0) {
            res=z0;
        } else if (z2==0) {
            z1.digits.insert(z1.digits.begin(),a0.digits.size(),0);
            res=z1+z0;
        } else {
            z1.digits.insert(z1.digits.begin(),a0.digits.size(),0);
            z2.digits.insert(z2.digits.begin(),2*a0.digits.size(),0);
            res=z2+z1+z0;
        }

        return res;
    }
}

int main() {
    clock_t start, end;

    BigInt a("984561231354629875468546546534125215534125215634987498548489456125215421563498749854848945612385663498749854848945612521542156349874985484894561238561698774565123165221393856169877456512316552156349874985484894561238561698774565123165221392213935215634987498548489456123856169877456512316522139521563498749854848945612385616987745651231652213949651465123151354686324848945612385616987745651231652213949651465123151354684132319321005482265341252156349874985484894561252154215634987498548489456123856264596162131");
    BigInt b("453412521563498749853412521563498749854848945612521542156349874985484894561238565484894561252154215634987498548489456123856848945612385616935462987546854521563498749854848945653412521563498749854848945612521542156349874985484894561238561238754579785616987745651231652213965465341235215634987495215634987498548489456123856169877456512316522139854848774565123165223546298754685465465341235215634987498548354629875468546546534123521563498749854844139496514651231513546298754685465465341235215634987498548435468");

    __longmult=false;
    __karatsuba=true;

    start=clock();
    for (int i=0;i<20;i++) {
        a*b;
    }
    end=clock();
    printf("\nTook %f seconds\n", (double)(end-start)/CLOCKS_PER_SEC);

    __longmult=true;
    __karatsuba=false;

    start=clock();
    for (int i=0;i<20;i++) {
        a*b;
    }
    end=clock();
    printf("\nTook %f seconds\n", (double)(end-start)/CLOCKS_PER_SEC);

    return 0;
}

Answer 1

1. 你使用 std::vector

   对于您的数字，请确保其中没有不必要的重新分配。所以在操作前分配空间以避免它。我也不使用它，所以我不知道数组范围检查速度变慢。

   检查你是否不移动它!!!这是 O(N) ... 即插入到第一个位置...

2. 优化您的实现

   here你可以找到我的实现优化和未优化的比较

   x=0.98765588997654321000000009876... | 98*32 bits...
   mul1[ 363.472 ms ]... O(N^2) classic multiplication
   mul2[ 349.384 ms ]... O(3*(N^log2(3))) optimized karatsuba multiplication
   mul3[ 9345.127 ms]... O(3*(N^log2(3))) unoptimized karatsuba multiplication 
   

   我的 Karatsuba 实现阈值大约是 3100 位 ... ~ 944 位!!!代码越优化，情人阈值就越高。

   
   尝试从函数操作数中删除不必要的数据

   //BigInt operator + (BigInt a,BigInt b)
   BigInt operator + (const BigInt &a,const BigInt &b)
   

   这样你就不会在每次 + 调用时在堆上创建另一个 a,b 拷贝，而且速度更快:

   mul(BigInt &ab,const BigInt &a,const BigInt &b) // ab = a*b
   

3. Schönhage-Strassen 乘法

   这是基于FFT 或NTT 的。我的阈值很大...... ~ 49700bits ...... ~ 15000digits 所以如果你不打算使用这么大的数字那么忘记它。实现也在上面的链接中。

   
   here是我的NTT实现(尽可能优化)

4. 总结

   无论您使用小端还是大端都没有关系，但您应该以不使用插入操作的方式对您的操作进行编码。

   
   您对数字使用十进制基数很慢，因为您需要使用除法和模运算。如果您选择base as power of 2，那么只需位运算就足够了并且它还会从代码中删除许多速度最慢的if语句。如果您需要基数作为 10 的幂，那么在某些情况下您可以使用最大的，这可以将 div、mod 减少到几个减法

   2^32 = 4 294 967 296 ... int = +/- 2147483648
   base = 1 000 000 000
   
   //x%=base
   while (x>=base) x-=base;
   

   在某些平台上，最大循环数是 2^32/base 或 2^31/base 这比取模更快，而且基数越大，您需要的操作就越少，但要注意溢出!!!