1. 首页
  2. 问题
  3. 快速的bignum平方计算
一尘不染

快速的bignum平方计算

algorithm

为了加快我的bignum划分速度，我需要加快y =x^2bigint的操作，bigint表示为无符号DWORD的动态数组。要清楚：

DWORD x[n+1] = { LSW, ......, MSW };
  • 其中n + 1是已用DWORD的数量
  • 所以数字的价值 x = x[0]+x[1]<<32 + ... x[N]<<32*(n)

问题是： 如何在y = x^2不损失精度的情况下尽快进行计算？ -可以使用 C ++ 和整数运算（带进位的32位）。

我目前的方法是应用乘法，y = x*x并避免多次乘法。

例如：

x = x[0] + x[1]<<32 + ... x[n]<<32*(n)

为了简单起见，让我重写一下：

x = x0+ x1 + x2 + ... + xn

其中index表示数组内部的地址，因此：

y = x*x
y = (x0 + x1 + x2 + ...xn)*(x0 + x1 + x2 + ...xn)
y = x0*(x0 + x1 + x2 + ...xn) + x1*(x0 + x1 + x2 + ...xn) + x2*(x0 + x1 + x2 + ...xn) + ...xn*(x0 + x1 + x2 + ...xn)

y0     = x0*x0
y1     = x1*x0 + x0*x1
y2     = x2*x0 + x1*x1 + x0*x2
y3     = x3*x0 + x2*x1 + x1*x2
...
y(2n-3) = xn(n-2)*x(n  ) + x(n-1)*x(n-1) + x(n  )*x(n-2)
y(2n-2) = xn(n-1)*x(n  ) + x(n  )*x(n-1)
y(2n-1) = xn(n  )*x(n  )

经过仔细观察，很明显几乎所有的对象都xi*xj出现了两次（不是第一个和最后一个），这意味着N*N乘法可以被(N+1)*(N/2)乘法代替。PS，32bit*32bit = 64bit因此每个mul+add操作的结果都将作为处理64+1 bit。

有没有更好的方法来快速计算？我在搜索过程中发现的只是sqrts算法，而不是sqr …

快速平方

!!! 请注意，我代码中的所有数字都首先是MSW，…不是上面的测试中那样（为了简化方程式，首先是LSW，否则将是索引混乱）。

当前的功能性fsqr实现

void arbnum::sqr(const arbnum &x)
{
    // O((N+1)*N/2)
    arbnum c;
    DWORD h, l;
    int N, nx, nc, i, i0, i1, k;
    c._alloc(x.siz + x.siz + 1);
    nx = x.siz - 1;
    nc = c.siz - 1;
    N = nx + nx;
    for (i=0; i<=nc; i++)
        c.dat[i]=0;
    for (i=1; i<N; i++)
        for (i0=0; (i0<=nx) && (i0<=i); i0++)
        {
            i1 = i - i0;
            if (i0 >= i1)
                break;
            if (i1 > nx)
                continue;
            h = x.dat[nx-i0];
            if (!h)
                continue;
            l = x.dat[nx-i1];
            if (!l)
                continue;
            alu.mul(h, l, h, l);
            k = nc - i;
            if (k >= 0)
                alu.add(c.dat[k], c.dat[k], l);
            k--;
            if (k>=0)
                alu.adc(c.dat[k], c.dat[k],h);
            k--;
            for (; (alu.cy) && (k>=0); k--)
                alu.inc(c.dat[k]);
        }
        c.shl(1);
        for (i = 0; i <= N; i += 2)
        {
            i0 = i>>1;
            h = x.dat[nx-i0];
            if (!h)
                continue;
            alu.mul(h, l, h, h);
            k = nc - i;
            if (k >= 0)
                alu.add(c.dat[k], c.dat[k],l);
            k--;
            if (k>=0)
                alu.adc(c.dat[k], c.dat[k], h);
            k--;
            for (; (alu.cy) && (k >= 0); k--)
                alu.inc(c.dat[k]);
        }
        c.bits = c.siz<<5;
        c.exp = x.exp + x.exp + ((c.siz - x.siz - x.siz)<<5) + 1;
        c.sig = sig;
        *this = c;
    }

使用唐津乘法

（感谢卡尔皮斯）

我实现了Karatsuba乘法，但是结果甚至比使用简单O(N^2)乘法都要慢得多，这可能是因为我无法避免避免这种可怕的递归。它的折衷必须是非常大的数字（大于数百个数字）……但是即使如此，仍然存在大量的内存传输。有没有一种方法可以避免递归调用（非递归变量，…几乎所有递归算法都可以通过这种方式完成）。尽管如此，我还是会尝试进行调整，看看会发生什么（避免规范化等，也可能是代码中的一些愚蠢的错误）。无论如何，在解决了Karatsuba的问题后，x*x性能没有太大提高。

优化的唐津乘法

性能测试y = x^2 looped 1000x times, 0.9 < x < 1 ~ 32*98 bits：

x = 0.98765588997654321000000009876... | 98*32 bits
sqr [ 213.989 ms ] ... O((N+1)*N/2) fast sqr
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

x = 0.98765588997654321000... | 195*32 bits
sqr [ 883.01 ms ]
mul1[ 1427.02 ms ]
mul2[ 1089.84 ms ]

x = 0.98765588997654321000... | 389*32 bits
sqr [ 3189.19 ms ]
mul1[ 5553.23 ms ]
mul2[ 3159.07 ms ]

经过Karatsuba的优化后，代码比以前快得多。但是，对于较小的数字，它的速度略小于我的O(N^2)乘法速度的一半。对于更大的数字，以布斯乘法的复杂性给出的比率更快。乘法的阈值大约为32 * 98位，而sqr的大约为32 * 389位，因此，如果输入位的总和超过此阈值，那么Karatsuba乘法将用于加速乘法，而sqr也是如此。

顺便说一句，优化包括：

  • 太大的递归参数将堆垃圾最小化
  • 取而代之的是避免使用任何带进位的bignum语法（+，-）32位ALU。
  • 忽视0*y或x*0或0*0案件
  • 将输入x,y数字的大小重新格式化为2的幂，以避免重新分配
  • 实现模乘以z1 = (x0 + x1)*(y0 + y1)最小化递归

将Schönhage-Strassen乘法修改为sqr实现

我已经测试过使用 FFT 和 NTT 转换来加快sqr计算的速度。结果如下：

  1. 快速傅立叶变换

精度下降，因此需要高精度的复数。实际上，这会大大降低速度，因此不存在加速问题。结果不精确（可能会四舍五入），因此 FFT 无法使用（目前）

  1. NTT

NTT 是有限域 DFT ，因此不会发生精度损失。它需要对无符号整数：modpow, modmul, modadd和进行模运算modsub。

我使用DWORD（32位无符号整数）。由于存在溢出问题， NTT
输入/输出向量大小受到限制！！！对于32位模块化算术，N仅限于此，(2^32)/(max(input[])^2)因此bigint必须将其划分为较小的块（我使用的BYTES最大bigint处理量为

    (2^32)/((2^8)^2) = 2^16 bytes = 2^14 DWORDs = 16384 DWORDs)

该sqr应用仅1xNTT + 1xINTT代替2xNTT + 1xINTT乘法，但 NTT
使用太慢和阈值数大小对我实施的实际使用（过大mul和也sqr）。

甚至有可能超过溢出限制，因此应使用64位模块化算术，这会进一步降低速度。因此， NTT 对我而言也无法使用。

一些测量：

a = 0.98765588997654321000 | 389*32 bits
looped 1x times
sqr1[ 3.177 ms ] fast sqr
sqr2[ 720.419 ms ] NTT sqr
mul1[ 5.588 ms ] simpe mul
mul2[ 3.172 ms ] karatsuba mul
mul3[ 1053.382 ms ] NTT mul

我的实现：

void arbnum::sqr_NTT(const arbnum &x)
{
    // O(N*log(N)*(log(log(N)))) - 1x NTT
    // Schönhage-Strassen sqr
    // To prevent NTT overflow: n <= 48K * 8 bit -> result siz <= 12K * 32 bit -> x.siz + y.siz <= 12K!!!
    int i, j, k, n;
    int s = x.sig*x.sig, exp0 = x.exp + x.exp - ((x.siz+x.siz)<<5) + 2;
    i = x.siz;
    for (n = 1; n < i; n<<=1)
        ;
    if (n + n > 0x3000) {
        _error(_arbnum_error_TooBigNumber);
        zero();
        return;
    }
    n <<= 3;
    DWORD *xx, *yy, q, qq;
    xx = new DWORD[n+n];
    #ifdef _mmap_h
    if (xx)
        mmap_new(xx, (n+n) << 2);
    #endif
    if (xx==NULL) {
        _error(_arbnum_error_NotEnoughMemory);
        zero();
        return;
    }
    yy = xx + n;

    // Zero padding (and split DWORDs to BYTEs)
    for (i--, k=0; i >= 0; i--)
    {
        q = x.dat[i];
        xx[k] = q&0xFF; k++; q>>=8;
        xx[k] = q&0xFF; k++; q>>=8;
        xx[k] = q&0xFF; k++; q>>=8;
        xx[k] = q&0xFF; k++;
    }
    for (;k<n;k++)
        xx[k] = 0;

    //NTT
    fourier_NTT ntt;

    ntt.NTT(yy,xx,n);    // init NTT for n

    // Convolution
    for (i=0; i<n; i++)
        yy[i] = modmul(yy[i], yy[i], ntt.p);

    //INTT
    ntt.INTT(xx, yy);

    //suma
    q=0;
    for (i = 0, j = 0; i<n; i++) {
        qq = xx[i];
        q += qq&0xFF;
        yy[n-i-1] = q&0xFF;
        q>>=8;
        qq>>=8;
        q+=qq;
    }

    // Merge WORDs to DWORDs and copy them to result
    _alloc(n>>2);
    for (i = 0, j = 0; i<siz; i++)
    {
        q  =(yy[j]<<24)&0xFF000000; j++;
        q |=(yy[j]<<16)&0x00FF0000; j++;
        q |=(yy[j]<< 8)&0x0000FF00; j++;
        q |=(yy[j]    )&0x000000FF; j++;
        dat[i] = q;
    }

    #ifdef _mmap_h
    if (xx)
        mmap_del(xx);
    #endif
    delete xx;
    bits = siz<<5;
    sig = s;
    exp = exp0 + (siz<<5) - 1;
        // _normalize();
    }

结论

对于较小的数字，这是我的快速sqr方法的最佳选择，在阈值 Karatsuba乘法之后更好。但是我仍然认为应该忽略一些琐碎的事情。还有其他想法吗？

NTT优化

经过大量优化（主要是 NTT）后：堆栈溢出问题[模块化算法和NTT（有限域DFT）优化。

一些值已更改：

a = 0.98765588997654321000 | 1553*32bits
looped 10x times
mul2[ 28.585 ms ] Karatsuba mul
mul3[ 26.311 ms ] NTT mul

因此现在，在大约1500 * 32位的阈值之后， NTT 乘法最终比 Karatsuba 更快。

发现一些测量结果和错误

a = 0.99991970486 | 1553*32 bits
looped: 10x
sqr1[  58.656 ms ] fast sqr
sqr2[  13.447 ms ] NTT sqr
mul1[ 102.563 ms ] simpe mul
mul2[  28.916 ms ] Karatsuba mul Error
mul3[  19.470 ms ] NTT mul

我发现我的 Karatsuba （上/下）流到bignum 的每个段的 LSBDWORD。研究完成后，我将更新代码…

而且，经过进一步 NTT 优化阈值改变，所以对于 NTT SQR 它310*32 bits = 9920 bits的 操作数
，以及用于 NTT MUL 它1396*32 bits = 44672 bits的 结果 （操作数的位的总和）。

@greybeard修复了Karatsuba代码

//---------------------------------------------------------------------------
void arbnum::_mul_karatsuba(DWORD *z, DWORD *x, DWORD *y, int n)
{
    // Recursion for Karatsuba
    // z[2n] = x[n]*y[n];
    // n=2^m
    int i;
    for (i=0; i<n; i++)
        if (x[i]) {
            i=-1;
            break;
        } // x==0 ?

    if (i < 0)
        for (i = 0; i<n; i++)
            if (y[i]) {
                i = -1;
                break;
            } // y==0 ?

    if (i >= 0) {
        for (i = 0; i < n + n; i++)
            z[i]=0;
            return;
        } // 0.? = 0

    if (n == 1) {
        alu.mul(z[0], z[1], x[0], y[0]);
        return;
    }

    if (n< 1)
        return;
    int n2 = n>>1;
    _mul_karatsuba(z+n, x+n2, y+n2, n2);                         // z0 = x0.y0
    _mul_karatsuba(z  , x   , y   , n2);                         // z2 = x1.y1
    DWORD *q = new DWORD[n<<1], *q0, *q1, *qq;
    BYTE cx,cy;
    if (q == NULL) {
        _error(_arbnum_error_NotEnoughMemory);
        return;
    }
    #define _add { alu.add(qq[i], q0[i], q1[i]); for (i--; i>=0; i--) alu.adc(qq[i], q0[i], q1[i]); } // qq = q0 + q1 ...[i..0]
    #define _sub { alu.sub(qq[i], q0[i], q1[i]); for (i--; i>=0; i--) alu.sbc(qq[i], q0[i], q1[i]); } // qq = q0 - q1 ...[i..0]
    qq = q;
    q0 = x + n2;
    q1 = x;
    i = n2 - 1;
    _add;
    cx = alu.cy; // =x0+x1

    qq = q + n2;
    q0 = y + n2;
    q1 = y;
    i = n2 - 1;
    _add;
    cy = alu.cy; // =y0+y1

    _mul_karatsuba(q + n, q + n2, q, n2);                       // =(x0+x1)(y0+y1) mod ((2^N)-1)

    if (cx) {
        qq = q + n;
        q0 = qq;
        q1 = q + n2;
        i = n2 - 1;
        _add;
        cx = alu.cy;
    }// += cx*(y0 + y1) << n2

    if (cy) {
        qq = q + n;
        q0 = qq;
        q1 = q;
        i = n2 -1;
        _add;
        cy = alu.cy;
    }// +=cy*(x0+x1)<<n2

    qq = q + n;  q0 = qq; q1 = z + n; i = n - 1; _sub;  // -=z0
    qq = q + n;  q0 = qq; q1 = z;     i = n - 1; _sub;  // -=z2
    qq = z + n2; q0 = qq; q1 = q + n; i = n - 1; _add;  // z1=(x0+x1)(y0+y1)-z0-z2

    DWORD ccc=0;

    if (alu.cy)
        ccc++;    // Handle carry from last operation
    if (cx || cy)
        ccc++;    // Handle carry from before last operation
    if (ccc)
    {
        i = n2 - 1;
        alu.add(z[i], z[i], ccc);
        for (i--; i>=0; i--)
            if (alu.cy)
                alu.inc(z[i]);
            else
                break;
    }

    delete[] q;
    #undef _add
    #undef _sub
    }

//---------------------------------------------------------------------------
void arbnum::mul_karatsuba(const arbnum &x, const arbnum &y)
{
    // O(3*(N)^log2(3)) ~ O(3*(N^1.585))
    // Karatsuba multiplication
    //
    int s = x.sig*y.sig;
    arbnum a, b;
    a = x;
    b = y;
    a.sig = +1;
    b.sig = +1;
    int i, n;
    for (n = 1; (n < a.siz) || (n < b.siz); n <<= 1)
        ;
    a._realloc(n);
    b._realloc(n);
    _alloc(n + n);
    for (i=0; i < siz; i++)
        dat[i]=0;
    _mul_karatsuba(dat, a.dat, b.dat, n);
    bits = siz << 5;
    sig = s;
    exp = a.exp + b.exp + ((siz-a.siz-b.siz)<<5) + 1;
    //    _normalize();
    }
//---------------------------------------------------------------------------

我的arbnum电话号码代表：

// dat is MSDW first ... LSDW last
DWORD *dat; int siz,exp,sig,bits;
  • dat[siz]是螳螂。LSDW表示最低有效的DWORD。
  • exp 是MSB的指数 dat[0]
  • 尾数中存在第一个非零位！！！
    // |-----|---------------------------|---------------|------|

    // | sig | MSB mantisa LSB | exponent | bits |
    // |-----|---------------------------|---------------|------|
    // | +1 | 0.(0 … 0) | 2^0 | 0 | +zero
    // | -1 | 0.(0 … 0) | 2^0 | 0 | -zero
    // |-----|---------------------------|---------------|------|
    // | +1 | 1.(dat[0] … dat[siz-1]) | 2^exp | n | +number
    // | -1 | 1.(dat[0] … dat[siz-1]) | 2^exp | n | -number
    // |-----|---------------------------|---------------|------|
    // | +1 | 1.0 | 2^+0x7FFFFFFE | 1 | +infinity
    // | -1 | 1.0 | 2^+0x7FFFFFFE | 1 | -infinity
    // |-----|---------------------------|---------------|------|


阅读 306

收藏
2020-07-28

共1个答案

一尘不染

如果我正确地理解了您的算法，似乎位数O(n^2)在哪里n。

您看过Karatsuba算法吗？它使用分而治之的方法加快了乘法速度。可能值得一看。

2020-07-28