Ответ 1
Если бы он преподавал еще больше, а Треугольник Pascal не прошел, потому что это привело бы к еще большему количеству операций. К счастью, операция mod может быть перенесена под PI, поэтому вам не нужно использовать большую int, но вместо 64-битной арифметики (или 32-битный modmul).
PI(ai+aj) mod p == PI((ai+aj)mod p) mod p ... 1<=i<j<=n
поэтому наивное решение С++ (где p<2^16) ваша задача требует вместо этого 64-битных переменных (к которым у меня нет доступа к простым термам).
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)
therms для 3 почти то же самое, поэтому мы можем это использовать. подсчитайте, сколько из тех же a[i] там подсчитывает частичный PI для одного из них. поместите это значение по счетчику и для каждого экземпляра умножите также на a[i]^instance здесь С++ пример:
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') с помощью:
- sort
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
В скобках указано время в [мс], 0: означает наивный подход, 1: означает отсортированную и частичную RLE-разложение PI. Последнее значение является результатом для 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] full RLE (a[i]) -> Разделение RLE (a[i]+a[j])
Я объединил все подходы выше и создал гораздо более быструю версию. Алгоритм таков:
-
RLE encode
a[i]просто создайте гистограмму
a[i]по сортировке ведра вO(n), а затем упакуйте в кодированиеvalue[n'],count[n'], так что в массиве нет нулей. Это довольно быстро. -
преобразовать RLE (
a[i]) в RLE (a[i]+a[j])просто создайте счетчик каждого терма
a[i]+a[j]в окончательном PI, подобном разложению RLE (a[i]+a[j]), но вO(N'.N')без какой-либо временной операции. Да, это квадратично, но наn'<=nи с очень маленьким постоянным временем. Но эта часть является узким местом... -
вычислить modpi из RLE (
a[i]+a[j])Это просто
modmul/modpowвO(n')наибольшее постоянное время, но низкая сложность, так что все еще очень быстро.
Код С++ для этого:
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. Здесь С++:
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.
Здесь comaprison к подходу 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 сек:).
Здесь несколько раз, чтобы сравнить с дополнительными быстрыми методами свертки, я попытался использовать Karatsuba и FastSQR из Быстрое вычисление квадратов bignum, чтобы избежать использования 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 делает это бесполезным для этой задачи.