计算成对和模 10^9+7 乘积的替代方法比 O(N^2) 快

Question

给定一个大小为N的整数数组A，我想计算

这是过去一次校际编程竞赛中的一个问题。我们必须编写一个程序来解决最多 5 个该问题的实例，其中 N ≤ 200,000 并且每个 a_{i ≤ 200,000，在 20 秒的运行时间限制内。显然，O(N²) 解决方案会超过时间限制。根据editorial，预期的解决方案涉及使用快速傅里叶变换的多项式乘法。我正在寻找比没有 FFT（也不是 NTT）的天真的 O(N²) 算法更快地解决这个问题的替代算法。这个问题有没有简单优雅的解决方案？}

已知事实：

mod 可以在产品内部“分发”，因为 (x*y) % m = ((x%m) * (y%m)) % m

更新： 这是比赛期间的输入/输出测试用例文件：如果它在 20 秒内通过，它将被接受。 输入：https://www.dropbox.com/s/nw81hts9rniter5/algol.in?dl=0 输出：https://www.dropbox.com/s/kpa7wit35xr4xm4/algol.out?dl=0

Answer 1

已经给了它更多的教导，帕斯卡三角是不行的，因为它会导致更多的操作。幸运的是 mod 操作可以在 PI 下移动，因此您不需要使用 big int 而是使用 64 位算法（或 32 位 modmul）。

PI(ai+aj) mod p == PI((ai+aj)mod p) mod p ... 1<=i<j<=n

太天真了 C++ 解决方案是（p<2^16）你的任务需要 64 位变量而不是（我无法在简单的 therms 中访问）。

DWORD modpi(DWORD *a,int n,DWORD p)
    {
    int i,j;
    DWORD x=1;
    for (i=0;i<n;i++)
     for (j=i+1;j<n;j++)
        {
        x*=a[i]+a[j];
        x%=p;
        }
    return x;
    }

现在p 比max(a[i]) 大很多，所以你可以改变：

x%=p;

与：

while (x>=p) x-=p;

但在如今的 CPU 上，这甚至更慢。

这种方法还是太慢了（~280 ms for n=10000）。如果我们对这些值重新排序（排序），那么事情就会突然变得更好。数组中的每个值多次导致简化，因为它的部分和几乎相同。例如：

a[] = { 2,3,3,4 }
x = (2+3).(2+3).(2+4)
  . (3+3).(3+4)
  . (3+4)

3 的温度几乎相同，因此我们可以使用它。计算有多少相同的a[i]，然后计算其中单个的部分 PI。通过计数为它提供动力，并为每个实例乘以 a[i]^instance 这里 C++ 示例：

DWORD modpi1(DWORD *a,int n,DWORD p)
    {
    int i,j,k;
    DWORD x,y,z;
    sort_asc(a,n);
    for (x=1,i=0;i<n;i++)
        {
        // count the values in k
        for (k=1;(i+1<n)&&(a[i]==a[i+1]);i++,k++);
        // compute partial modPI y
        for (y=1,j=i+1;j<n;j++)
            {
            y*=a[i]+a[j];
            y%=p;
            }
        // add the partial modPI y^k;
        for (j=0;j<k;j++)
            {
            x*=y;
            x%=p;
            }
        // add powers for instances of a[i]
        for (k--;k;k--)
         for (j=0;j<k;j++)
            {
            x*=a[i]+a[i];
            x%=p;
            }
        }
    return x;
    }

这为数组中每个多次出现的值提供了一些加速。但是由于您的数组与其中可能的数字一样大，因此不要期望太多。对于均匀随机数据，max(a[i])~=n 和 快速排序 的速度略低于 50%。但是如果你使用像 MSalters 这样的素数分解，答案表明你可能会得到真正的加速，因为重复率应该比~1 高得多，但这需要大量的工作来处理方程。

此代码为O(N.N')，其中N' 是a[] 中不同值的计数。您还可以通过以下方式将其进一步增强为O(N'.N')：

1. 按桶排序O(n)或快速排序O(n.log(n))对a[i]进行排序

2. 执行 RLE（运行长度编码）O(n)

3. 帐户也计入部分总和O(n'.n') 其中n'<=n

   素数分解应该只是将n'<= n 更改为n' <<< n。

这里使用 快速排序 完整的 32 位 modmul 进行一些测量（使用 32 位 x86 asm，这会大大降低我的编译器的速度）。随机数据凡max(a[i])~=n：

n=1000;
[   4.789 ms]0: 234954047
[   3.044 ms]1: 234954047
n=10000;
[ 510.544 ms]0: 629694784
[ 330.876 ms]1: 629694784
n=20000;
[2126.041 ms]0: 80700577
[1350.966 ms]1: 80700577

括号中是以 [ms] 为单位的时间，0: 表示幼稚方法，1: 表示 PI 的排序和部分 RLE 分解。最后一个值是p=1000000009的结果

如果这还不够，那么使用 DFT/NTT 分开我认为没有其他可能的加速。

[Edit1] a[i] 的完整 RLE 分解

//---------------------------------------------------------------------------
const DWORD p=1000000009;
const int n=10000;
const int m=n;
DWORD a[n];
//---------------------------------------------------------------------------
DWORD modmul(DWORD a,DWORD b,DWORD p)
    {
    DWORD _a,_b;
    _a=a;
    _b=b;
    asm {
        mov    eax,_a
        mov    ebx,_b
        mul    ebx   // H(edx),L(eax) = eax * ebx
        mov    ebx,p
        div    ebx   // eax = H(edx),L(eax) / ebx
        mov    _a,edx// edx = H(edx),L(eax) % ebx
        }
    return _a;
    }
//---------------------------------------------------------------------------
DWORD modpow(DWORD a,DWORD b,DWORD p)
    {   // b is not mod(p) !
    int i;
    DWORD d=1;
    for (i=0;i<32;i++)
        {
        d=modmul(d,d,p);
        if (DWORD(b&0x80000000)) d=modmul(d,a,p);
        b<<=1;
        }
    return d;
    }
//---------------------------------------------------------------------------
DWORD modpi(DWORD *a,int n,DWORD p)
    {
    int i,j,k;
    DWORD x,y;
    DWORD *value=new DWORD[n+1];// RLE value
    int   *count=new int[n+1];  // RLE count
    // O(n) bucket sort a[] -> count[] because max(a[i])<=n
    for (i=0;i<=n;i++) count[i]=0;
    for (i=0;i< n;i++) count[a[i]]++;
    // O(n) RLE packing value[n],count[n]
    for (i=0,j=0;i<=n;i++)
     if (count[i])
        {
        value[j]=    i;
        count[j]=count[i];
        j++;
        } n=j;
    // compute the whole PI to x
    for (x=1,i=0;i<n;i++)
        {
        // compute partial modPI value[i]+value[j] to y
        for (y=1,j=i+1;j<n;j++)
         for (k=0;k<count[j];k++)
          y=modmul(y,value[i]+value[j],p);
        // add the partial modPI y^count[j];
        x=modmul(x,modpow(y,count[i],p),p);
        // add powers for instances of value[i]
        for (j=0,k=1;k<count[i];k++) j+=k;
        x=modmul(x,modpow(value[i]+value[i],j,p),p);
        }
    delete[] value;
    delete[] count;
    return x;
    }
//---------------------------------------------------------------------------

这甚至更快，因为它在O(n) 和 RLE 在O(n) 中排序，所以这导致O(N'.N')。如果有的话，您可以利用更高级的modmul,modpow 例程。但是对于值的均匀分布，这仍然没有接近可用速度。

[edit2] a[i]+a[j] 的完整 RLE 分解

DWORD modpi(DWORD *a,int n,DWORD p) // full RLE(a[i]+a[j]) O(n'.n') n' <= 2n
    {
    int i,j;
    DWORD x,y;
    DWORD nn=(n+1)*2;
    int   *count=new int[nn+1]; // RLE count
    // O(n^2) bucket sort a[] -> count[] because max(a[i]+a[j])<=nn
    for (i=0;i<=nn;i++) count[i]=0;
    for (i=0;i<n;i++)
     for (j=i+1;j<n;j++)
      count[a[i]+a[j]]++;
    // O(n') compute the whole PI to x
    for (x=1,y=0;y<=nn;y++)
     if (count[y])
      x=modmul(x,modpow(y,count[y],p),p);
    delete[] count;
    return x;
    }
//---------------------------------------------------------------------------

在接近理想的时间时，这甚至更快，但仍然相差几个数量级。

n=20000
[3129.710 ms]0: 675975480 // O(n^2) naive
[2094.998 ms]1: 675975480 // O(n'.n) partial RLE decomposition of a[i] , n'<= n
[2006.689 ms]2: 675975480 // O(n'.n') full RLE decomposition of a[i] , n'<= n
[ 729.983 ms]3: 675975480 // T(c0.n^2+c1.n') full RLE decomposition of a[i]+a[j] , n'<= 2n , c0 <<< c1

[Edit3]完整的RLE(a[i])->RLE(a[i]+a[j])分解

我结合了上述所有方法并创建了更快的版本。算法是这样的：

1. RLE 编码a[i]

   只需在O(n) 中通过桶排序创建a[i] 的直方图，然后打包到编码value[n'],count[n']，因此数组中不存在零。这非常快。

2. 将 RLE(a[i]) 转换为 RLE(a[i]+a[j])

   只需在最终 PI 中创建每个 a[i]+a[j] therm 的计数，类似于 RLE(a[i]+a[j]) 分解，但在 O(n'.n') 中不需要任何时间要求的操作。是的，这是二次的，但在 n'<=n 上并且 非常小 恒定时间。但这部分是瓶颈...

3. 从 RLE(a[i]+a[j]) 计算 modpi

   这很简单 modmul/modpow 在 O(n') 最大的常数时间但复杂度低所以仍然非常快。

C++ 代码：

DWORD modpi(DWORD *a,int n,DWORD p) // T(c0.n+c1.n'.n'+c2.n'') full RLE(a[i]->a[i]+a[j]) n' <= n , n'' <= 2n , c0 <<< c1 << c2
    {
    int i,j,k;
    DWORD x,y;
    DWORD nn=(n+1)*2;
    DWORD *rle_iv =new DWORD[ n+1]; // RLE a[i] value
    int   *rle_in =new int[ n+1];   // RLE a[i] count
    int   *rle_ij=new int[nn+1];    // RLE (a[i]+a[j]) count
    // O(n) bucket sort a[] -> rle_i[] because max(a[i])<=n
    for (i=0;i<=n;i++) rle_in[i]=0;
    for (i=0;i<n;i++)  rle_in[a[i]]++;
    for (x=0,i=0;x<=n;x++)
     if (rle_in[x])
        {
        rle_iv[i]=       x;
        rle_in[i]=rle_in[x];
        i++;
        } n=i;
    // O(n'.n') convert rle_iv[]/in[] to rle_ij[]
    for (i=0;i<=nn;i++) rle_ij[i]=0;
    for (i=0;i<n;i++)
        {
        rle_ij[rle_iv[i]+rle_iv[i]]+=(rle_in[i]*(rle_in[i]-1))>>1; // 1+2+3+...+(rle_iv[i]-1)
        for (j=i+1;j<n;j++)
         rle_ij[rle_iv[i]+rle_iv[j]]+=rle_in[i]*rle_in[j];
        }
    // O(n') compute the whole PI to x
    for (x=1,y=0;y<=nn;y++)
     if (rle_ij[y])
      x=modmul(x,modpow(y,rle_ij[y],p),p);
    delete[] rle_iv;
    delete[] rle_in;
    delete[] rle_ij;
    return x;
    }

和对比测量：

n=10000
[ 751.606 ms] 814157062 O(n^2) naive
[ 515.944 ms] 814157062 O(n'.n) partial RLE(a[i]) n' <= n
[ 498.840 ms] 814157062 O(n'.n') full RLE(a[i]) n' <= n
[ 179.896 ms] 814157062 T(c0.n^2+c1.n') full RLE(a[i]+a[j]) n' <= 2n , c0 <<< c1
[  66.695 ms] 814157062 T(c0.n+c1.n'.n'+c2.n'') full RLE(a[i]->a[i]+a[j]) n' <= n , n'' <= 2n , c0 <<< c1 << c2
n=20000
[ 785.177 ms] 476588184 T(c0.n^2+c1.n') full RLE(a[i]+a[j]) n' <= 2n , c0 <<< c1
[ 255.503 ms] 476588184 T(c0.n+c1.n'.n'+c2.n'') full RLE(a[i]->a[i]+a[j]) n' <= n , n'' <= 2n , c0 <<< c1 << c2
n=100000
[6158.516 ms] 780587335 T(c0.n+c1.n'.n'+c2.n'') full RLE(a[i]->a[i]+a[j]) n' <= n , n'' <= 2n , c0 <<< c1 << c2

最后一次是这种方法。加倍 n 将运行时间乘以 cca 4 次。所以对于n=200000，我的设置运行时间约为 24 秒。

[Edit4] 我的NTT 比较方法

我知道你想避免 FFT，但我仍然认为这有利于比较。 32 位 NTT 就可以了。因为它仅应用于直方图，该直方图仅由几位宽且大部分等于1 的指数组成，即使在n=200000 上也能防止溢出。这里C++来源：

DWORD modpi(DWORD *a,int n,int m,DWORD p) // O(n.log(n) RLE(a[i])+NTT convolution
    {
    int i,z;
    DWORD x,y;
    for (i=1;i<=m;i<<=1); m=i<<1;   // m power of 2 > 2*(n+1)
    #ifdef _static_arrays
    m=2*M;
    DWORD rle[2*M];                 // RLE a[i]
    DWORD con[2*M];                 // convolution c[i]
    DWORD tmp[2*M];                 // temp
    #else
    DWORD *rle =new DWORD[m];       // RLE a[i]
    DWORD *con =new DWORD[m];       // convolution c[i]
    DWORD *tmp =new DWORD[m];       // temp
    #endif
    fourier_NTT ntt;
    // O(n) bucket sort a[] -> rle[] because max(a[i])<=n
    for (i=0;i<m;i++) rle[i]=0.0;
    for (i=0;i<n;i++) rle[a[i]]++;

    // O(m.log(m)) NTT convolution
    for (i=0;i<m;i++) con[i]=rle[i];
    ntt.NTT(tmp,con,m);
    for (i=0;i<m;i++) tmp[i]=ntt.modmul(tmp[i],tmp[i]);
    ntt.iNTT(con,tmp,m);
    // O(n') compute the whole PI to x
    for (x=1,i=0;i<m;i++)
        {
        z=con[i];
        if (int(i&1)==0) z-=int(rle[(i+1)>>1]);
        z>>=1; y=i;
        if ((y)&&(z)) x=modmul(x,modpow(y,z,p),p);
        }
    #ifdef _static_arrays
    #else
    delete[] rle;
    delete[] con;
    delete[] tmp;
    #endif
    return x;
    }

你可以忽略_static_arrays（处理它，因为它没有定义）它只是为了更简单的调试。当心卷积ntt.modmul 不适用于p 的任务，而是使用NTTs modulo ！！！如果您想绝对确定这适用于更高的n 或使用 64 位 NTT 的不同数据分布。

这里比较 Edit3 方法：

n=200000
[24527.645 ms] 863132560 O(m^2) RLE(a[i]) -> RLE(a[i]+a[j]) m <= n
[  754.409 ms] 863132560 O(m.log(m)) RLE(a[i])+NTT

如您所见，我距离估计的 ~24 秒并不太远 :)。

这里有时与我尝试使用来自Fast bignum square computation 的 Karatsuba 和 FastSQR 的其他快速卷积方法进行比较，以避免使用 FFT/NTT：

n=10000
[ 749.033 ms] 149252794 O(n^2)        naive
[1077.618 ms] 149252794 O(n'^2)       RLE(a[i])+fast_sqr32
[ 568.510 ms] 149252794 O(n'^1.585)   RLE(a[i])+Karatsuba32
[  65.805 ms] 149252794 O(n'^2)       RLE(a[i]) -> RLE(a[i]+a[j])
[  53.833 ms] 149252794 O(n'.log(n')) RLE(a[i])+FFT
[  34.129 ms] 149252794 O(n'.log(n')) RLE(a[i])+NTT
n=20000
[3084.546 ms] 365847531 O(n^2)        naive
[4311.491 ms] 365847531 O(n'^2)       RLE(a[i])+fast_sqr32
[1672.769 ms] 365847531 O(n'^1.585)   RLE(a[i])+Karatsuba32
[ 238.725 ms] 365847531 O(n'^2)       RLE(a[i]) -> RLE(a[i]+a[j])
[ 115.047 ms] 365847531 O(n'.log(n')) RLE(a[i])+FFT
[  71.587 ms] 365847531 O(n'.log(n')) RLE(a[i])+NTT
n=40000
[12592.250 ms] 347013745 O(n^2)        naive
[17135.248 ms] 347013745 O(n'^2)       RLE(a[i])+fast_sqr32
[5172.836 ms] 347013745 O(n'^1.585)   RLE(a[i])+Karatsuba32
[ 951.256 ms] 347013745 O(n'^2)       RLE(a[i]) -> RLE(a[i]+a[j])
[ 242.918 ms] 347013745 O(n'.log(n')) RLE(a[i])+FFT
[ 152.553 ms] 347013745 O(n'.log(n')) RLE(a[i])+NTT

遗憾的是，Karatsuba 的开销太大，因此阈值高于 n=200000，使其对这项任务毫无用处。

Answer 2

由于ai <= 200.000 和N<=200.000，总共可能有 40.000.000.000 个术语，但您知道ai + aj <= 400.000。最多可以有 400.000 个唯一术语。这已经好 5 个数量级了。

但是，这些项中的大多数都不是质数；在 400.000 之下只有约 40.000 个素数。您最终可能会得到每个单独术语的更高的多重性，但这没什么大不了的。即使对于大 X，计算 (prime^N) 模 1000000007 也足够快。

您可以合理地预先计算所有数字

这种方法实现了加速，因为我们延迟了乘法，而是计算通过查找找到的小素因数。当我们需要进行乘法运算时，我们已经有了一系列指数，并且可以使用重复平方来有效地减少它们。

当“众所周知的事实”是质数分解很难时，我们使用质数分解作为加速可能是违反直觉的。但这是可能的，因为每一项都很小，而且我们反复需要相同的因式分解。

[编辑] 从 cmets 看来，要弄清楚 ai+aj 的多重性似乎很困难，因为您只能计算 i<j 所在的项。但这不是问题。计算所有项 ai+aj 的重数，并除以 2，因为 aj+i==ai+aj。这仅对i==j 的对角线是错误的。这可以通过在除以 2 之前添加所有项 ai+ai 的多重性来解决。

例如：a={1 2 3}，要考虑的术语是{1+1, 1+2, 1+3, 2+2, 2+3, 3+3} [三角形]。 4 的重数是 2（通过 1+3 和 2+2）。相反，请考虑{1+1, 1+2, 1+3, 2+1, 2+2, 2+3, 3+1, 3+2, 3+3} [正方形] + {1+1, 2+2, 3+3} [对角线]。 4 的重数现在是 4（1+3,2+2,3+1 和 2+2），除以 2 得到正确的结果。

由于a[] 的顺序对于方形变体不再重要，您可以对其使用计数排序。例如。给定{4,5,6,5}，我们得到4:1, 5:2, 6:1。因此 10 的重数是4+6:1, 5+5:2, 6+4:1