verguenza - tag de preguntas incomodas
¡Optimízame!(C, rendimiento)-seguimiento de la pregunta (16)
Gracias a algunos usuarios de stackOverflow muy útiles en Bit twiddling: ¿qué bit está configurado? , He construido mi función (publicada al final de la pregunta).
Cualquier sugerencia, incluso pequeñas sugerencias, sería apreciada. Con suerte, mejorará mi código, pero al menos debería enseñarme algo. :)
Visión de conjunto
Esta función se llamará al menos 10 13 veces, y posiblemente tan seguido como 10 15 . Es decir, este código se ejecutará durante meses con toda probabilidad, por lo que cualquier sugerencia de rendimiento sería útil.
Esta función representa el 72-77% del tiempo del programa, en función del perfil y alrededor de una docena de ejecuciones en diferentes configuraciones (optimizando ciertos parámetros no relevantes aquí).
Por el momento, la función se ejecuta en un promedio de 50 relojes. No estoy seguro de cuánto se puede mejorar esto, pero estaría encantado de verlo correr en 30.
Observación clave
Si en algún punto del cálculo puede decir que el valor que será devuelto será pequeño (valor exacto negociable, por ejemplo, debajo de un millón) , puede cancelar anticipadamente . Solo me interesan los grandes valores.
Así es como espero ahorrar la mayor cantidad de tiempo, en lugar de otras micro optimizaciones (¡aunque también son bienvenidas!).
Información de rendimiento
- smallprimes es una matriz de bits (64 bits); en promedio se establecerán unos 8 bits, pero podrían ser tan pocos como 0 o tanto como 12.
- q generalmente será distinto de cero. (Observe que la función sale temprano si q y las primas son cero).
- rys a menudo serán 0. Si q es cero, rys también lo serán; si r es cero, s también lo será.
- Como dice el comentario al final, nu suele ser 1 al final, así que tengo un caso especial eficiente para él.
- Los cálculos debajo del caso especial pueden parecer desbordamiento de riesgo, pero a través del modelado apropiado, he demostrado que, por mi opinión, esto no ocurrirá, así que no se preocupe por ese caso.
- Las funciones no definidas aquí (ugcd, minuu, star, etc.) ya se han optimizado; ninguno toma mucho tiempo para correr. pr es una pequeña matriz (todo en L1). Además, todas las funciones llamadas aquí son funciones puras .
- Pero si realmente te importa ... ugcd es el gcd , minuu es el mínimo, vals es el número de 0 binarios finales, __builtin_ffs es la ubicación del binario 1 más a la izquierda, la estrella es (n-1) >> vals (n- 1), pr es una matriz de primos de 2 a 313.
- Los cálculos se están realizando actualmente en un Phenom II 920 x4, aunque las optimizaciones para i7 o Woodcrest siguen siendo de interés (si obtengo tiempo de cálculo en otros nodos).
- Me complacería responder cualquier pregunta que tenga sobre la función o sus componentes.
Lo que realmente hace
Agregado en respuesta a una solicitud. No necesita leer esta parte.
La entrada es un número impar n con 1 <n <4282250400097. Las otras entradas proporcionan la factorización del número en este sentido particular:
smallprimes & 1 se establece si el número es divisible entre 3, smallprimes & 2 se establece si el número es divisible entre 5, lowprimes & 4 se establece si el número es divisible entre 7, lowprimes & 8 se establece si el número es divisible entre 11, etc. hasta el máximo un bit significativo que representa 313. Un número divisible por el cuadrado de un primo no se representa de manera diferente a un número divisible solo por ese número. (De hecho, se pueden descartar múltiplos de cuadrados; en la etapa de preprocesamiento en otra función, los múltiplos de cuadrados de primos <= lim tienen valores pequeños yq se ponen a 0 para que se eliminen, donde el valor óptimo de lim se determina por experimentación. )
q, r y s representan factores más grandes del número. Cualquier factor restante (que puede ser mayor que la raíz cuadrada del número, o si s es distinto de cero puede ser incluso menor) se puede encontrar dividiendo los factores de n.
Una vez que todos los factores se recuperan de esta manera, se cuenta el número de bases, 1 <= b <n, a las que n es un fuerte pseudoprima , usando una fórmula matemática que se explica mejor con el código.
Mejoras hasta el momento
- Empujó la prueba de salida temprana hacia arriba. Esto claramente ahorra trabajo, así que hice el cambio.
- Las funciones apropiadas ya están en línea, por lo que
__attribute__ ((inline))no hace nada. Curiosamente, marcar lasbasesfunciones principales y algunos de los ayudantes con__attribute ((hot))dañó el rendimiento en casi un 2% y no puedo entender por qué (pero es reproducible con más de 20 pruebas). Entonces no hice ese cambio. Del mismo modo,__attribute__ ((const)), en el mejor de los casos, no ayudó. Estaba más que ligeramente sorprendido por esto.
Código
ulong bases(ulong smallprimes, ulong n, ulong q, ulong r, ulong s)
{
if (!smallprimes & !q)
return 0;
ulong f = __builtin_popcountll(smallprimes) + (q > 1) + (r > 1) + (s > 1);
ulong nu = 0xFFFF; // "Infinity" for the purpose of minimum
ulong nn = star(n);
ulong prod = 1;
while (smallprimes) {
ulong bit = smallprimes & (-smallprimes);
ulong p = pr[__builtin_ffsll(bit)];
nu = minuu(nu, vals(p - 1));
prod *= ugcd(nn, star(p));
n /= p;
while (n % p == 0)
n /= p;
smallprimes ^= bit;
}
if (q) {
nu = minuu(nu, vals(q - 1));
prod *= ugcd(nn, star(q));
n /= q;
while (n % q == 0)
n /= q;
} else {
goto BASES_END;
}
if (r) {
nu = minuu(nu, vals(r - 1));
prod *= ugcd(nn, star(r));
n /= r;
while (n % r == 0)
n /= r;
} else {
goto BASES_END;
}
if (s) {
nu = minuu(nu, vals(s - 1));
prod *= ugcd(nn, star(s));
n /= s;
while (n % s == 0)
n /= s;
}
BASES_END:
if (n > 1) {
nu = minuu(nu, vals(n - 1));
prod *= ugcd(nn, star(n));
f++;
}
// This happens ~88% of the time in my tests, so special-case it.
if (nu == 1)
return prod << 1;
ulong tmp = f * nu;
long fac = 1 << tmp;
fac = (fac - 1) / ((1 << f) - 1) + 1;
return fac * prod;
}
Asegúrate de que tus funciones estén en línea. Si están fuera de línea, la sobrecarga puede sumarse, especialmente en el primer ciclo while. La mejor manera de estar seguro es examinar el conjunto.
¿Has probado
star( pr[__builtin_ffsll(bit)] )yvals( pr[__builtin_ffsll(bit)] - 1)? Eso cambiaría un trabajo simple para una búsqueda en matriz, pero podría valer la pena si las tablas son lo suficientemente pequeñas.No calcule
fhasta que realmente lo necesite (cerca del final, después de su salida anticipada). Puedes reemplazar el código en BASES_END por algo como
BASES_END:
ulong addToF = 0;
if (n > 1) {
nu = minuu(nu, vals(n - 1));
prod *= ugcd(nn, star(n));
addToF = 1;
}
// ... early out if nu == 1...
// ... compute f ...
f += addToF;
Espero que ayude.
¿ smallprimes pasar una matriz de números primos en lugar de dividirlos en smallprimes , q , smallprimes ? Como no sé lo que hace el código externo, probablemente esté equivocado, pero existe la posibilidad de que también tenga una función para convertir algunos números primos en un mapa de bits de smallprimes , y dentro de esta función, convierta el mapa de bits a una matriz de primos, con eficacia. Además, parece que haces un procesamiento idéntico para elementos de smallprimes , q , r y s . Debería ahorrarle una pequeña cantidad de procesamiento por llamada.
Además, parece saber que los primos pasados dividen n . ¿Tienes suficiente conocimiento sobre el poder de cada primo que divide n? Puede ahorrar mucho tiempo si puede eliminar la operación del módulo pasando esa información a esta función. En otras palabras, si n es pow(p_0,e_0)*pow(p_1,e_1)*...*pow(p_k,e_k)*n_leftover , y si sabes más acerca de estos e_i sy n_leftover , pasarlos significa muchas cosas que no tienes que hacer en esta función.
Puede haber una manera de descubrir n_leftover (la parte no factorizada de n ) con menos número de operaciones de módulo, pero es solo una corazonada, por lo que puede necesitar experimentar un poco con él. La idea es usar gcd para eliminar los factores conocidos de n repetidamente hasta que se deshaga de todos los factores primos conocidos. Déjame dar un código casi-c:
factors=p_0*p_1*...*p_k*q*r*s;
n_leftover=n/factors;
do {
factors=gcd(n_leftover, factors);
n_leftover = n_leftover/factors;
} while (factors != 1);
No estoy del todo seguro de que esto sea mejor que el código que tiene, y mucho menos las sugerencias mod / div combinadas que puede encontrar en otras respuestas, pero creo que vale la pena intentarlo. Siento que será una victoria, especialmente para los números con un gran número de pequeños factores primos.
¿El código en su publicación principal es la versión optimizada? If yes, there is still too many divide operations which greatly eat CPU cycles.
This code is overexecute innecessarily a bit
if (!smallprimes & !q)
return 0;
change to logical and &&
if (!smallprimes && !q)
return 0;
will make it short circuited faster without eveluating q
And the following code
ulong bit = smallprimes & (-smallprimes);
ulong p = pr[__builtin_ffsll(bit)];
which is used to find the last set bit of smallprimes. Why don''t you use the simpler way
ulong p = pr[__builtin_ctz(smallprimes)];
Another culprit for decreased performance maybe too many program branching. You may consider changing to some other less-branch or branch-less equivalents
¿Has intentado utilizar la optimización guiada por perfil?
Compile y vincule el programa con la opción -fprofile-generate , luego ejecute el programa sobre un conjunto de datos representativos (digamos, un día de cómputo).
Luego vuelva a compilar y vincularlo con la -fprofile-use lugar.
¿Intentó una versión de búsqueda de tabla del primer ciclo while? Puede dividir las smallprimes en valores de 4 16 bits, buscar su contribución y fusionarlas. Pero tal vez necesites los efectos secundarios.
1) Haría que el compilador escupe el ensamblaje que genera e intente y deduzca si lo que hace es lo mejor que puede hacer ... y si detecta problemas, cambie el código para que el conjunto se vea mejor. De esta forma, también puede asegurarse de que las funciones que espera que estén en línea (como estrella y val) estén realmente en línea. (Es posible que tenga que agregar pragma, o incluso convertirlos en macros)
2) Es genial que pruebes esto en una máquina multinúcleo, pero este ciclo es singlethreaded. Supongo que hay una función paraguas que divide la carga en unos pocos hilos para que se usen más núcleos.
3) Es difícil sugerir aceleraciones si lo que la función real intenta calcular no está claro. Por lo general, las aceleraciones más impresionantes no se logran con el twiddling de bits, pero con un cambio en el algoritmo. Así que un poco de comentarios pueden ayudar; ^)
4) Si realmente desea una velocidad de 10 * o más, consulte CUDA o openCL que le permite ejecutar programas C en su hardware de gráficos. ¡Brilla con funciones como estas!
5) Estás haciendo montones de módulo y divide justo después uno del otro. En C esto es 2 comandos separados (primero ''/'' y luego ''%''). Sin embargo, en el ensamblaje, este es 1 comando: ''DIV'' o ''IDIV'' que devuelve el resto y el cociente de una vez:
B.4.75 IDIV: Signed Integer Divide
IDIV r/m8 ; F6 /7 [8086]
IDIV r/m16 ; o16 F7 /7 [8086]
IDIV r/m32 ; o32 F7 /7 [386]
IDIV performs signed integer division. The explicit operand provided is the divisor; the dividend and destination operands are implicit, in the following way:
For IDIV r/m8, AX is divided by the given operand; the quotient is stored in AL and the remainder in AH.
For IDIV r/m16, DX:AX is divided by the given operand; the quotient is stored in AX and the remainder in DX.
For IDIV r/m32, EDX:EAX is divided by the given operand; the quotient is stored in EAX and the remainder in EDX.
Por lo tanto, requerirá un ensamblaje en línea, pero supongo que habrá una aceleración significativa ya que hay algunos lugares en su código que pueden beneficiarse de esto.
Estás pasando la factorización completa de n, por lo que estás factorizando enteros consecutivos y luego usando los resultados de esa factorización aquí. Me parece que podría beneficiarse al hacer algo de esto al momento de encontrar los factores.
Por cierto, tengo un código muy rápido para encontrar los factores que estás usando sin hacer ninguna división. Es un poco como un tamiz, pero produce factores de números consecutivos muy rápidamente. Puede encontrarlo y publicarlo si cree que puede ser útil.
edit tuvo que volver a crear el código aquí:
#include #define SIZE (1024*1024) //must be 2^n #define MASK (SIZE-1) typedef struct { int p; int next; } p_type; p_type primes[SIZE]; int sieve[SIZE]; void init_sieve() { int i,n; int count = 1; primes[1].p = 3; sieve[1] = 1; for (n=5;SIZE>n;n+=2) { int flag = 0; for (i=1;count>=i;i++) { if ((n%primes[i].p) == 0) { flag = 1; break; } } if (flag==0) { count++; primes[count].p = n; sieve[n>>1] = count; } } } int main() { int ptr,n; init_sieve(); printf("init_done/n"); // factor odd numbers starting with 3 for (n=1;1000000000>n;n++) { ptr = sieve[n&MASK]; if (ptr == 0) //prime { // printf("%d is prime",n*2+1); } else //composite { // printf ("%d has divisors:",n*2+1); while(ptr!=0) { // printf ("%d ",primes[ptr].p); sieve[n&MASK]=primes[ptr].next; //move the prime to the next number it divides primes[ptr].next = sieve[(n+primes[ptr].p)&MASK]; sieve[(n+primes[ptr].p)&MASK] = ptr; ptr = sieve[n&MASK]; } } // printf("/n"); } return 0; }
La función init crea una base de factores e inicializa el tamiz. Esto lleva unos 13 segundos en mi computadora portátil. Luego, se toman en cuenta todos los números de hasta mil millones o se determina que son primos en otros 25 segundos. Los números inferiores a SIZE nunca se informan como primos porque tienen 1 factor en la base de factores, pero eso podría modificarse.
La idea es mantener una lista vinculada para cada entrada en el tamiz. Los números se tienen en cuenta simplemente sacando sus factores de la lista vinculada. A medida que se extraen, se insertan en la lista para el siguiente número que será divisible por ese primo. Esto también es muy caché. El tamaño del tamiz debe ser mayor que el primo más grande en la base del factor. Como es, este tamiz podría alcanzar hasta 2 ** 40 en aproximadamente 7 horas, lo que parece ser su objetivo (excepto que n necesita tener 64 bits).
Su algoritmo podría fusionarse en esto para hacer uso de los factores tal como se identifican en lugar de los bits de empaque y primos grandes en variables para pasar a su función. O su función podría cambiarse para tomar la lista enlazada (podría crear un enlace ficticio para pasar los números primos fuera de la base del factor).
Espero eso ayude.
Por cierto, esta es la primera vez que publico este algoritmo públicamente.
Intente reemplazar este patrón (para r yq también):
n /= p;
while (n % p == 0)
n /= p;
Con este:
ulong m;
...
m = n / p;
do {
n = m;
m = n / p;
} while ( m * p == n);
En mis pruebas limitadas, obtuve una pequeña aceleración (10%) al eliminar el módulo.
Además, si p, q o r son constantes, el compilador reemplazará las divisiones por multiplicaciones. Si hay pocas opciones para p, qor, o si algunas son más frecuentes, puede obtener algo especializando la función para esos valores.
Pareces estar perdiendo mucho tiempo haciendo divisiones por los factores. Es mucho más rápido reemplazar una división por una multiplicación por el recíproco de divisor (división: ~ 15-80 ( ! ) Ciclos, dependiendo del divisor, multiplicación: ~ 4 ciclos), IF , por supuesto, puede calcular previamente los recíprocos.
Si bien esto parece poco probable que sea posible con q , r , s - debido al rango de esos valores, es muy fácil de hacer con p , que siempre proviene de la matriz pequeña, estática pr [] . Precalcular los recíprocos de esos primos y almacenarlos en otra matriz. Luego, en lugar de dividir por p , multiplica por el recíproco tomado de la segunda matriz. (O crea una única serie de estructuras).
Ahora, obtener el resultado de división exacto mediante este método requiere algunos trucos para compensar los errores de redondeo. Encontrará los detalles sangrientos de esta técnica en este documento , en la página 138.
EDITAR:
Después de consultar Hacker''s Delight (un excelente libro, BTW) sobre el tema, parece que puedes hacerlo aún más rápido explotando el hecho de que todas las divisiones en tu código son exactas (es decir, el resto es cero).
Parece que para cada divisor d que es impar y base B = 2 word_size , existe un único inverso multiplicativo d⃰ que cumple las condiciones: d⃰ < B y d·d⃰ ≡ 1 (mod B) . Para cada x que es un múltiplo exacto de d , esto implica x/d ≡ x·d⃰ (mod B) . Lo que significa que puedes simplemente reemplazar una división con una multiplicación, sin correcciones, controles, problemas de redondeo, lo que sea. (Las pruebas de estos teoremas se pueden encontrar en el libro.) ¡ Tenga en cuenta que este inverso multiplicativo no necesita ser igual al recíproco como lo define el método anterior!
¿Cómo verificar si una x dada es un múltiplo exacto de d , es decir, x mod d = 0 ? ¡Fácil! x mod d = 0 iff x·d⃰ mod B ≤ ⌊(B-1)/d⌋ . Tenga en cuenta que este límite superior se puede precalcular.
Entonces, en el código:
unsigned x, d;
unsigned inv_d = mulinv(d); //precompute this!
unsigned limit = (unsigned)-1 / d; //precompute this!
unsigned q = x*inv_d;
if(q <= limit)
{
//x % d == 0
//q == x/d
} else {
//x % d != 0
//q is garbage
}
Suponiendo que la matriz pr[] convierte en una matriz de struct prime :
struct prime {
ulong p;
ulong inv_p; //equal to mulinv(p)
ulong limit; //equal to (ulong)-1 / p
}
el while(smallprimes) en su código se convierte en:
while (smallprimes) {
ulong bit = smallprimes & (-smallprimes);
int bit_ix = __builtin_ffsll(bit);
ulong p = pr[bit_ix].p;
ulong inv_p = pr[bit_ix].inv_p;
ulong limit = pr[bit_ix].limit;
nu = minuu(nu, vals(p - 1));
prod *= ugcd(nn, star(p));
n *= inv_p;
for(;;) {
ulong q = n * inv_p;
if (q > limit)
break;
n = q;
}
smallprimes ^= bit;
}
Y para la función mulinv() :
ulong mulinv(ulong d) //d needs to be odd
{
ulong x = d;
for(;;)
{
ulong tmp = d * x;
if(tmp == 1)
return x;
x *= 2 - tmp;
}
}
Tenga en cuenta que puede reemplazar ulong con cualquier otro tipo sin firmar, solo use el mismo tipo de forma consistente.
Las pruebas, por qué s y cómo están todos disponibles en el libro. Una lectura recomendada con entusiasmo :-).
Pequeña optimización pero:
ulong f;
ulong nn;
ulong nu = 0xFFFF; // "Infinity" for the purpose of minimum
ulong prod = 1;
if (!smallprimes & !q)
return 0;
// no need to do this operations before because of the previous return
f = __builtin_popcountll(smallprimes) + (q > 1) + (r > 1) + (s > 1);
nn = star(n);
Por cierto: debe editar su publicación para agregar star() y otras funciones que usa la definición
Primero un poco quisquilloso ;-) deberías tener más cuidado con los tipos que estás usando. En algunos lugares, parece suponer que ulong tiene 64 bits de ancho, use uint64_t allí. Y también para todos los demás tipos, reconsidere cuidadosamente lo que espera de ellos y use el tipo apropiado.
La optimización que pude ver es la división de enteros. Tu código hace eso mucho, esta es probablemente la cosa más cara que estás haciendo. La división de enteros pequeños ( uint32_t ) puede ser mucho más eficiente que la de los grandes. En particular para uint32_t hay una instrucción de ensamblador que hace división y módulo de una vez, llamado divl .
Si usa los tipos apropiados, su compilador puede hacer todo eso por usted. Pero será mejor que verifique el ensamblador (opción -S a gcc) como alguien ya dijo. De lo contrario, es fácil incluir algunos pequeños fragmentos de ensamblador aquí y allá. Encontré algo así en algún código mío:
register uint32_t a asm("eax") = 0;
register uint32_t ret asm("edx") = 0;
asm("divl %4"
: "=a" (a), "=d" (ret)
: "0" (a), "1" (ret), "rm" (divisor));
Como puedes ver, usa registros especiales eax y edx y cosas así ...
Si las divisiones que está utilizando son con números que no se conocen en tiempo de compilación, pero se usan con frecuencia en tiempo de ejecución (dividiendo por el mismo número muchas veces), entonces sugeriría usar la biblioteca libdivide , que básicamente implementa en tiempo de ejecución el optimizaciones que los compiladores hacen para compilar constantes de tiempo (usando máscaras de cambios, etc.). Esto puede proporcionar un gran beneficio. También evitando usar x% y == 0 para algo como z = x / y, z * y == x como ergosys sugirió ergosys anteriormente también debería tener una mejora medible.
Si su compilador admite atributos de función GCC, puede marcar sus funciones puras con este atributo:
ulong star(ulong n) __attribute__ ((const));
Este atributo indica al compilador que el resultado de la función depende solo de su (s) argumento (s). Esta información puede ser utilizada por el optimizador.
¿Hay algún motivo por el que haya vals() código abierto vals() lugar de usar __builtin_ctz() ?
Si va a salir de inmediato (!smallprimes&!q) ¿por qué no hacer esa prueba incluso antes de llamar a la función, y guardar la llamada de la llamada sobrecarga?
Además, parece que efectivamente tiene 3 funciones diferentes que son lineales, excepto para el ciclo lowprimes. bases1(s,n,q) , bases2(s,n,q,r) y bases3(s,n,q,r,s) .
Puede ser una victoria crearlos realmente como 3 funciones separadas sin las ramas y los gotos, y llamar al apropiado:
if (!(smallprimes|q)) { r = 0;}
else if (s) { r = bases3(s,n,q,r,s);}
else if (r) { r = bases2(s,n,q,r); }
else { r = bases1(s,n,q);
Esto sería más efectivo si el procesamiento previo ya le ha dado al código de llamada un cierto "conocimiento" de qué función ejecutar y usted no tiene que probarlo.
Todavía no está claro qué es lo que está buscando. Muy frecuentemente los problemas teóricos de los números permiten grandes aceleraciones al derivar propiedades matemáticas que las soluciones deben satisfacer.
Si de hecho está buscando los enteros que maximizan el número de no testigos para la prueba de MR (es decir, oeis.org/classic/A141768 que menciona), entonces podría ser posible usar que el número de no testigos no puede ser mayor. que phi (n) / 4 y que los enteros para los que tienen muchos no testigos son el producto de dos números primos de la forma
(k + 1) * (2k + 1)
o son números de Carmichael con 3 factores primos. Pensaría que, por encima de un límite, todos los enteros en la secuencia tienen esta forma y que es posible verificar esto demostrando un límite superior para los testigos de todos los demás enteros. Por ejemplo, los enteros con 4 o más factores siempre tienen como máximo phi (n) / 8 no testigos. Resultados similares pueden derivarse de su fórmula para la cantidad de bases para otros enteros.
En cuanto a las micro-optimizaciones: cada vez que sepa que un número entero es divisible por algún cociente, entonces es posible reemplazar la división por una multiplicación con el inverso del módulo del cociente 2 ^ 64. Y las pruebas n% q == 0 pueden ser reemplazadas por una prueba
n * inverse_q <max_q,
donde inverse_q = q ^ (- 1) mod 2 ^ 64 y max_q = 2 ^ 64 / q. Obviamente inverse_q y max_q necesitan ser precalculados, para ser eficientes, pero ya que estás usando un tamiz, supongo que esto no debería ser un obstáculo.
solo un pensamiento, pero tal vez usar las opciones de optimización de tu compilador sería útil, si no lo has hecho aún. Otro pensamiento sería que si el dinero no es un problema, podría usar el compilador Intel C / C ++, suponiendo que usa un procesador Intel. También supongo que otros fabricantes de procesadores (AMD, etc.) tendrían compiladores similares