以O(N ^ 2)速度计算成对总和mod 10 ^ 9 + 7的乘积的替代方法
给定大小为N的整数数组A,我想计算
这是过去大学间编程竞赛中的一个问题.我们必须写这将解决了这个问题的5实例程序,与ñ ≤200,000,每一个我 ≤20万辆,20秒运行时间限制内.显然,O(N 2)解决方案将超过时限.根据社论,预期的解决方案涉及使用快速傅里叶变换的多项式乘法.我正在寻找可以比没有FFT(也不是NTT)的朴素O(N 2)算法更快地解决这个问题的替代算法.这个问题有没有简单而优雅的解决方案？
已知事实:
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
如果给它更多的教导，帕斯卡三角形是不行的，因为它会导致更多的操作。幸运的是，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对于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')通过以下方式进一步增强此功能：
a[i]按桶排序O(n)或快速排序排序O(n.log(n))执行 RLE（行程编码）
O(n)帐户也计入部分总和，
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] 完整的 RLE 分解a[i]
//---------------------------------------------------------------------------
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] 完整的 RLE 分解a[i]+a[j]
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]) 分解
我结合了上述所有方法并创建了更快的版本。算法是这样的：
RLE编码
a[i]只需创建一个
a[i]按桶排序的直方图O(n)，然后打包到编码中value[n'],count[n']，这样数组中就不存在零。这速度相当快。将 RLE(
a[i]) 转换为 RLE(a[i]+a[j])只需在最终 PI 中创建每个
a[i]+a[j]热值的计数，类似于 RLE(a[i]+a[j]) 分解，但不需要O(n'.n')任何时间要求较高的操作。是的，这是二次的，但是在非常小的常数时间n'<=n上。但这部分是瓶颈......a[i]+a[j]从 RLE( )计算 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运行时间乘以 cca4倍。因此，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，而是用于 NTT 模运算！如果您想绝对确定这可以在更高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 秒并不太远:)。
这里有时与其他快速卷积方法进行比较，我尝试使用快速 bignum 平方计算中的 Karasuba 和 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使其无法完成此任务。
| 归档时间： |
|
| 查看次数： |
715 次 |
| 最近记录： |