algorithm - vestir - Forma alternativa de calcular el producto de las sumas por pares mod 10 ^ 9+7 más rápido que O(N ^ 2)
se requiere un algoritmo para determinar de n cantidades cuantas son cero (2)
Como ai <= 200.000 y N<=200.000 , puede haber 40.000.000.000 términos en total, pero usted sabe que ai + aj <= 400.000 . Puede haber como máximo 400.000 términos únicos . Eso ya es 5 órdenes de magnitud mejor.
Sin embargo, la mayoría de estos términos no son primos; Hay solo ~ 40.000 primos por debajo de 400.000. Puede terminar con una multiplicidad algo mayor de cada término individual, pero eso no es un gran problema. El cálculo (prime ^ N) módulo 1000000007 es lo suficientemente rápido incluso para grandes X.
Puede pre-calcular razonablemente la factorización de todos los números <= 400.000 y obtener los números primos <= 400.000 como un efecto secundario libre.
Este método logra una aceleración porque demoramos la multiplicación y, en cambio, contamos los factores primos pequeños encontrados a través de una búsqueda. Para cuando tengamos que hacer las multiplicaciones, tendremos una serie de exponentes y podremos usar la cuadratura repetida para reducirlos de manera eficiente.
Tal vez sea contrario a la intuición que usemos la factorización prima como una aceleración, cuando el "hecho bien conocido" es que la factorización prima es difícil. Pero esto es posible porque cada término es pequeño y repetidamente necesitamos la misma factorización.
[editar] De los comentarios, parece que averiguar la multiplicidad de ai+aj es difícil, ya que solo puedes contar los términos donde i<j . Pero eso no es un problema. Cuente la multiplicidad de todos los términos ai + aj, y divida por dos ya que aj + i == ai + aj. Esto solo es incorrecto para la diagonal donde i==j . Esto se arregla agregando la multiplicidad de todos los términos ai + ai antes de dividir por 2.
Ej: a={1 2 3} , los términos a considerar son {1+1, 1+2, 1+3, 2+2, 2+3, 3+3} [triángulo]. La multiplicidad de 4 es 2 (a través de 1 + 3 y 2 + 2). En su lugar, considere {1+1, 1+2, 1+3, 2+1, 2+2, 2+3, 3+1, 3+2, 3+3} [cuadrado] + {1+1, 2+2, 3+3} [diagonal]. La multiplicidad de 4 es ahora 4 (1 + 3,2 + 2,3 + 1 y 2 + 2), se divide por 2 para obtener el resultado correcto.
Y como el orden de a[] ya no importa para la variante cuadrada, puede usar una ordenación de conteo en ella. Por ejemplo, dado {4,5,6,5} , obtenemos 4:1, 5:2, 6:1 . Así, la multiplicidad de 10 es 4+6:1, 5+5:2, 6+4:1
Dada una matriz A de enteros de tamaño N , quiero calcular
Esto fue un problema en una pasada competencia de programación interuniversitaria. Tuvimos que escribir un programa que resolviera hasta 5 instancias de este problema, con N ≤ 200,000 y cada una i ≤ 200,000, dentro de un límite de tiempo de ejecución de 20 segundos. Claramente, una solución de O ( N 2 ) superaría el límite de tiempo. Según el editorial , la solución pretendida consistía en la multiplicación polinomial utilizando una transformada rápida de Fourier. Estoy buscando algoritmos alternativos que resolverían esto más rápido que el algoritmo ingenuo O ( N 2 ) sin FFT (ni NTT). ¿Hay alguna solución simple y elegante para este problema?
Hechos conocidos:
mod podría ''distribuirse'' dentro del producto porque (x * y)% m = ((x% m) * (y% m))% m
Actualización: aquí está el archivo de prueba / entrada durante el concurso: si pasa esto en 20 segundos, se aceptará. Entrada: https://www.dropbox.com/s/nw81hts9rniter5/algol.in?dl=0 Salida: https://www.dropbox.com/s/kpa7wit35xr4xm4/algol.out?dl=0
Le había dado algo más de enseñado y el triángulo de Pascal es un no ir porque conduciría a más operaciones. Afortunadamente, la operación de modificación se puede mover debajo de la PI, por lo que no necesita usar aritmética de 64 bits sino int (o modmul de 32 bits).
PI(ai+aj) mod p == PI((ai+aj)mod p) mod p ... 1<=i<j<=n
por lo que la solución ingenua de C ++ es (donde p<2^16 ) su tarea requiere variables de 64 bits (a las que no tengo acceso en términos simples).
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;
}
Ahora la p es mucho más grande que el max(a[i]) por lo que podría cambiar:
x%=p;
con:
while (x>=p) x-=p;
Pero en la actualidad la CPU es incluso más lenta.
Aún así, este enfoque es demasiado lento ( ~280 ms para n=10000 ). Si reordenamos los valores (los ordenamos), de repente, las cosas mejoraron mucho. Cada valor que se encuentra en la matriz varias veces conduce a la simplificación, ya que su suma parcial es casi la misma. Por ejemplo:
a[] = { 2,3,3,4 }
x = (2+3).(2+3).(2+4)
. (3+3).(3+4)
. (3+4)
Las termias para 3 son casi las mismas, así que podemos usar eso. cuente cuántos de los mismos a[i] hay, luego cuente el PI parcial para uno de ellos. alimente esto por el recuento y para cada instancia, multiplíquese también por a[i]^instance aquí Ejemplo de 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;
}
Esto le da un poco de aceleración para cada aparición múltiple de valor en la matriz. Pero como su matriz es tan grande como los números posibles en ella, entonces no espere demasiado de esto. Para datos uniformemente aleatorios donde max(a[i])~=n y ordenación rápida es la velocidad, un poco menos del 50%. Pero si usa la descomposición principal como la respuesta de MSalters, sugiere que podría obtener una aceleración real porque entonces la tasa de repetición debería ser mucho mayor que ~1 pero eso requeriría mucho trabajo para manejar las ecuaciones.
Este código es O(N.N'') donde N'' es el conteo de valores distintos en a[] . También puede mejorar esto más a O(N''.N'') mediante:
- ordene
a[i]por cubeta clasifiqueO(n)oO(n.log(n))rápidamenteO(n.log(n)) - hacer RLE (codificación de longitud de ejecución)
O(n) la cuenta también cuenta a la suma parcial
O(n''.n'')donden''<=nLa descomposición principal debería cambiar
n''<= nan'' <<< n.
Aquí algunas mediciones usando modmul de 32 bits de ordenamiento rápido (usando 32bit x86 asm que ralentiza las cosas considerablemente en mi compilador). Datos aleatorios donde 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
Entre paréntesis está el tiempo en [ms], 0: significa enfoque ingenuo, 1: significa ordenación y descomposición parcial de RLE de PI. El último valor es el resultado para p=1000000009
Si eso todavía no es suficiente, a parte de usar DFT / NTT, no veo otra aceleración posible.
[Edit1] descomposición RLE completa de 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;
}
//---------------------------------------------------------------------------
Esto es incluso un poco más rápido, ya que clasifica en O(n) y RLE en O(n) por lo que resulta en O(N''.N'') . Puedes aprovechar las rutinas modmul,modpow más avanzadas si tienes alguna. Pero para una distribución uniforme de los valores, esto todavía no está cerca de las velocidades utilizables.
[edit2] descomposición RLE completa de 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;
}
//---------------------------------------------------------------------------
Y esto es incluso más rápido acercándose a los tiempos deseables, pero aún quedan pocas magnitudes.
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] ) descomposición
Combine todos los enfoques anteriores y creo una versión mucho más rápida. El algoritmo es así:
RLE codifica
a[i]simplemente cree un histograma de
a[i]clasificando en cubeta enO(n)y luego empácelo alvalue[n''],count[n'']codificaciónvalue[n''],count[n'']para que no haya ningún cero en la matriz. Esto es bastante rápido.convertir RLE (
a[i]) a RLE (a[i]+a[j])simplemente cree un recuento de cada
a[i]+a[j]term en la PI final similar a la descomposición de RLE (a[i]+a[j]) pero enO(n''.n'')sin ningún tiempo queO(n''.n'')operación. Sí, esto es cuadrático pero enn''<=ny con un tiempo constante muy pequeño . Pero esta parte es el cuello de botella ...calcular el modpi desde RLE (
a[i]+a[j])Esto es simple
modmul/modpowenO(n'')mayor tiempo constante pero de baja complejidad, por lo que aún es muy rápido.
El código C ++ para esto:
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;
}
Y medidas de comparación:
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
Los últimos tiempos son para este método. Duplicar n multiplica el tiempo de ejecución por cca 4 veces. así que para n=200000 el tiempo de ejecución es de alrededor de 24 segundos en mi configuración.
[Edit4] mi enfoque NTT para la comparación
Sé que quieres evitar la FFT pero aún así creo que esto es bueno para la comparación. El 32bit NTT está bien para esto. Porque se aplica solo en el histograma, que consiste solo en exponentes que tienen solo unos pocos bits de ancho y en su mayoría son iguales a 1 que evita los desbordamientos incluso en n=200000 . Aquí la fuente de 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;
}
Puede ignorar los _static_arrays (manejarlo ya que no está definido) es solo para una depuración más simple. ¡Cuidado con la convolución ntt.modmul no funciona con las tareas p sino con el módulo NTTs! Si desea estar absolutamente seguro de que esto funcionaría en n altas o diferentes distribuciones de datos, use 64bit NTT.
Aquí se relaciona con el enfoque 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
Como puede ver, no estaba muy lejos de los ~ 24 segundos estimados :).
Aquí algunas veces para comparar con los métodos adicionales de convolución rápida que probé con el uso de Karatsuba y FastSQR del cálculo de bignum cuadrado rápido para evitar el uso de 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
Lamentablemente, la sobrecarga de Karatsuba es demasiado grande, por lo que el umbral está por encima de n=200000 por lo que es inútil para esta tarea.