c++ - numeros - factorizacion prima de 64
Fuerza bruta, factorización prima de un solo hilo (7)
Para su consideración está la siguiente función que se puede usar para (relativamente rápidamente) factorizar un entero sin signo de 64 bits en sus factores primos. Tenga en cuenta que la factorización no es probabalística (es decir, es exacta). El algoritmo ya es lo suficientemente rápido como para encontrar que un número es primo o tiene pocos factores muy grandes en cuestión de segundos, en el hardware moderno.
La pregunta: ¿Se pueden realizar mejoras al algoritmo presentado, al tiempo que se mantiene de un solo hilo, de modo que pueda factorizar (de forma arbitraria) enteros de 64 bits sin signo muy grandes más rápidamente, preferiblemente sin utilizar un enfoque probabalístico (por ejemplo, Miller-Rabin) para determinar la primalidad?
// system specific typedef for ulong should go here (or use boost::uint64_t)
typedef unsigned __int64 ulong;
typedef std::vector<ulong> ULongVector;
// Caller needs to pass in an empty factors vector
void GetFactors(ULongVector &factors, ulong num)
{
// Num has to be at least 2 to contain "prime" factors
if (num<2)
return;
ulong workingNum=num;
ulong nextOffset=2; // Will be used to skip multiples of 3, later
// Factor out factors of 2
while (workingNum%2==0)
{
factors.push_back(2);
workingNum/=2;
}
// Factor out factors of 3
while (workingNum%3==0)
{
factors.push_back(3);
workingNum/=3;
}
// If all of the factors were 2s and 3s, done...
if (workingNum==1)
return;
// sqrtNum is the (inclusive) upper bound of our search for factors
ulong sqrtNum=(ulong) sqrt(double(workingNum+0.5));
// Factor out potential factors that are greate than or equal to 5
// The variable n represents the next potential factor to be tested
for (ulong n=5;n<=sqrtNum;)
{
// Is n a factor of the current working number?
if (workingNum%n==0)
{
// n is a factor, so add it to the list of factors
factors.push_back(n);
// Divide current working number by n, to get remaining number to factor
workingNum/=n;
// Check if we''ve found all factors
if (workingNum==1)
return;
// Recalculate the new upper bound for remaining factors
sqrtNum=(ulong) sqrt(double(workingNum+0.5));
// Recheck if n is a factor of the new working number,
// in case workingNum contains multiple factors of n
continue;
}
// n is not or is no longer a factor, try the next odd number
// that is not a multiple of 3
n+=nextOffset;
// Adjust nextOffset to be an offset from n to the next non-multiple of 3
nextOffset=(nextOffset==2UL ? 4UL : 2UL);
}
// Current workingNum is prime, add it as a factor
factors.push_back(workingNum);
}
Gracias
Edición: he añadido aún más comentarios. La razón por la que se pasa un vector por referencia, es para permitir que el vector se reutilice entre llamadas y evitar las asignaciones dinámicas. La razón por la que el vector no se vacía en la función, es para permitir el impar de agregar los factores "num" actuales a los factores que ya están en el vector.
La función en sí no es bonita y puede ser refactora, pero la pregunta es sobre cómo hacer que el algoritmo sea más rápido. Entonces, por favor, no hay sugerencias sobre cómo hacer que la función sea más bonita, legible, o C ++ ish. Eso es juego de niños. Mejorar este algoritmo, para que pueda encontrar los factores (probados) más rápido, es la parte difícil.
Actualización: Potatoswatter tiene algunas soluciones excelentes hasta ahora, asegúrese de revisar su solución MMX cerca de la parte inferior, también.
Compare tal enfoque con un tamiz (pre-generado). El módulo es caro, por lo que ambos enfoques esencialmente hacen dos cosas: generan factores potenciales y realizan operaciones de módulo. Cualquiera de los dos programas debería generar razonablemente un nuevo factor candidato en menos ciclos que los módulos, por lo que cualquiera de los dos programas está enlazado a módulo.
El enfoque dado filtra una proporción constante de todos los enteros, es decir, los múltiplos de 2 y 3, o el 75%. Uno de cada cuatro números (como se indica) se utiliza como argumento para el operador de módulo. Lo llamaré un filtro de salto.
Por otro lado, un tamiz usa solo primos como argumentos para el operador de módulo, y la diferencia promedio entre primos sucesivos se rige por el teorema del número primo para que sea 1 / ln (N). Por ejemplo, e ^ 20 tiene poco menos de 500 millones, por lo que los números que superan los 500 millones tienen menos de un 5% de probabilidad de ser mejores. Si se consideran todos los números hasta 2 ^ 32, el 5% es una buena regla general.
Por lo tanto, un tamiz pasará 5 veces menos tiempo en operaciones div como su filtro de omisión. El siguiente factor a considerar es la velocidad a la que el tamiz produce números primos, es decir, los lee de la memoria o del disco. Si obtener una prima es más rápido que 4 div s, entonces el tamiz es más rápido. De acuerdo con mis tablas, el rendimiento de mi Core2 es a lo sumo uno por cada 12 ciclos. Estos serán problemas difíciles de división, por lo que vamos a presupuestar de manera conservadora 50 ciclos por primo. Para un procesador de 2.5 GHz, eso es 20 nanosegundos.
En 20 ns, un disco duro de 50 MB / s puede leer aproximadamente un byte. La solución simple es utilizar 4 bytes por primo, por lo que la unidad será más lenta. Pero, podemos ser más inteligentes. Si queremos codificar todos los números primos en orden, simplemente podemos codificar sus diferencias. Una vez más, la diferencia esperada es 1 / ln (N). Además, todos son iguales, lo que ahorra un bit extra. Y nunca son cero, lo que hace que la extensión a una codificación multibyte sea gratuita. Por lo tanto, al usar un byte por primo, las diferencias hasta 512 se pueden almacenar en un byte, lo que nos lleva hasta 303371455241 de acuerdo con ese artículo de Wikipedia .
Por lo tanto, dependiendo del disco duro, una lista almacenada de números primos debe ser aproximadamente igual en velocidad al verificar la primalidad. Si puede almacenarse en la RAM (es de 203 MB, por lo que las ejecuciones subsiguientes probablemente afectarán al caché del disco), entonces el problema desaparece por completo, ya que la velocidad del FSB generalmente difiere de la velocidad del procesador en un factor menor que el ancho del FSB en bytes - Es decir, el FSB puede transferir más de una prima por ciclo. Entonces, el factor de mejora es la reducción en las operaciones de división, es decir, cinco veces. Esto se confirma por los resultados experimentales a continuación.
Por supuesto, entonces hay multihilo. Se pueden asignar rangos de primos o de candidatos filtrados por omisión a diferentes subprocesos, haciendo que ambos enfoques sean vergonzosamente paralelos. No hay optimizaciones que no impliquen aumentar el número de circuitos divisores paralelos, a menos que de alguna manera elimines el módulo.
Aquí hay un programa así. Está templado para que pueda agregar bignums.
/*
* multibyte_sieve.cpp
* Generate a table of primes, and use it to factorize numbers.
*
* Created by David Krauss on 10/12/10.
*
*/
#include <cmath>
#include <bitset>
#include <limits>
#include <memory>
#include <fstream>
#include <sstream>
#include <iostream>
#include <iterator>
#include <stdint.h>
using namespace std;
char const primes_filename[] = "primes";
enum { encoding_base = (1<< numeric_limits< unsigned char >::digits) - 2 };
template< typename It >
unsigned decode_gap( It &stream ) {
unsigned gap = static_cast< unsigned char >( * stream ++ );
if ( gap ) return 2 * gap; // only this path is tested
gap = ( decode_gap( stream )/2-1 ) * encoding_base; // deep recursion
return gap + decode_gap( stream ); // shallow recursion
}
template< typename It >
void encode_gap( It &stream, uint32_t gap ) {
unsigned len = 0, bytes[4];
gap /= 2;
do {
bytes[ len ++ ] = gap % encoding_base;
gap /= encoding_base;
} while ( gap );
while ( -- len ) { // loop not tested
* stream ++ = 0;
* stream ++ = bytes[ len + 1 ];
}
* stream ++ = bytes[ 0 ];
}
template< size_t lim >
void generate_primes() {
auto_ptr< bitset< lim / 2 > > sieve_p( new bitset< lim / 2 > );
bitset< lim / 2 > &sieve = * sieve_p;
ofstream out_f( primes_filename, ios::out | ios::binary );
ostreambuf_iterator< char > out( out_f );
size_t count = 0;
size_t last = sqrtl( lim ) / 2 + 1, prev = 0, x = 1;
for ( ; x != last; ++ x ) {
if ( sieve[ x ] ) continue;
size_t n = x * 2 + 1; // translate index to number
for ( size_t m = x + n; m < lim/2; m += n ) sieve[ m ] = true;
encode_gap( out, ( x - prev ) * 2 );
prev = x;
}
for ( ; x != lim / 2; ++ x ) {
if ( sieve[ x ] ) continue;
encode_gap( out, ( x - prev ) * 2 );
prev = x;
}
cout << prev * 2 + 1 << endl;
}
template< typename I >
void factorize( I n ) {
ifstream in_f( primes_filename, ios::in | ios::binary );
if ( ! in_f ) {
cerr << "Could not open primes file./n"
"Please generate it with ''g'' command./n";
return;
}
while ( n % 2 == 0 ) {
n /= 2;
cout << "2 ";
}
unsigned long factor = 1;
for ( istreambuf_iterator< char > in( in_f ), in_end; in != in_end; ) {
factor += decode_gap( in );
while ( n % factor == 0 ) {
n /= factor;
cout << factor << " ";
}
if ( n == 1 ) goto finish;
}
cout << n;
finish:
cout << endl;
}
int main( int argc, char *argv[] ) {
if ( argc != 2 ) goto print_help;
unsigned long n;
if ( argv[1][0] == ''g'' ) {
generate_primes< (1ul<< 32) >();
} else if ( ( istringstream( argv[1] ) >> n ).rdstate() == ios::eofbit )
factorize( n );
} else goto print_help;
return 0;
print_help:
cerr << "Usage:/n/t" << argv[0] << " <number> -- factorize number./n"
"/t" << argv[0] << " g -- generate primes file in current directory./n";
}
Rendimiento en un MacBook Pro de 2.2 GHz:
dkrauss$ time ./multibyte_sieve g
4294967291
real 2m8.845s
user 1m15.177s
sys 0m2.446s
dkrauss$ time ./multibyte_sieve 18446743721522234449
4294967231 4294967279
real 0m5.405s
user 0m4.773s
sys 0m0.458s
dkrauss$ time ./mike 18446743721522234449
4294967231 4294967279
real 0m25.147s
user 0m24.170s
sys 0m0.096s
Este código es bastante más lento, y estoy bastante seguro de que entiendo por qué. No es increíblemente más lento, pero definitivamente más lento en el rango de 10-20%. La división no se debe hacer una vez por cada bucle, pero la única forma de hacerlo es llamar a sqrt o algo similar.
// system specific typedef for ulong should go here (or use boost::uint64_t)
typedef std::vector<ulong> ULongVector;
void GetFactors(ULongVector &factors, ulong num)
{
if (num<2)
return;
ulong workingNum=num;
ulong nextOffset=2;
while (workingNum%2==0)
{
factors.push_back(2);
workingNum/=2;
}
while (workingNum%3==0)
{
factors.push_back(3);
workingNum/=3;
}
ulong n = 5;
while ((workingNum != 1) && ((workingNum / n) >= n)) {
// Is workingNum divisible by n?
if (workingNum%n==0)
{
// n is a factor!
// so is the number multiplied by n to get workingNum
// Insert n into the list of factors
factors.push_back(n);
// Divide working number by n
workingNum/=n;
} else {
n+=nextOffset;
nextOffset=(nextOffset==2UL ? 4UL : 2UL);
}
}
if (workingNum != 1) {
// workingNum is prime, add it to the list of factors
factors.push_back(workingNum);
}
}
Incorporando algunas ideas de Omnifarious, además de otras mejoras:
// system specific typedef for ulong should go here (or use boost::uint64_t)
typedef unsigned __int64 ulong;
typedef std::vector<ulong> ULongVector;
// Caller needs to pass in an empty factors vector
void GetFactors(ULongVector &factors, ulong num)
{
if (num<2)
return;
ulong workingNum=num;
// Factor out factors of 2
while (workingNum%2==0)
{
factors.push_back(2);
workingNum/=2;
}
// Factor out factors of 3
while (workingNum%3==0)
{
factors.push_back(3);
workingNum/=3;
}
if (workingNum==1)
return;
// Factor out factors >=5
ulong nextOffset=2;
char nextShift = 1;
ulong n = 5;
ulong nn = 25;
do {
// Is workingNum divisible by n?
if (workingNum%n==0)
{
// n is a factor!
// so is the number multiplied by n to get workingNum
// Insert n into the list of factors
factors.push_back(n);
// Divide working number by n
workingNum/=n;
// Test for done...
if (workingNum==1)
return;
// Try n again
}
else {
nn += (n << (nextShift+1)) + (1<<(nextShift*2)); // (n+b)^2 = n^2 + 2*n*b + b*2
n += nextOffset;
nextOffset ^= 6;
nextShift ^= 3;
// invariant: nn == n*n
if (n & 0x100000000LL) break; // careful of integer wraparound in n^2
}
} while (nn <= workingNum);
// workingNum is prime, add it to the list of factors
factors.push_back(workingNum);
}
La generalización natural es calcular las omisiones de página utilizando números primos más conocidos que solo 2 y 3. Me gusta 2, 3, 5, 7, 11, para un período de patrón de 2310 (eh, buen número). Y tal vez más, pero tiene rendimientos decrecientes: un gráfico de tiempos de ejecución podría establecer exactamente dónde la precomputación comienza a tener un impacto negativo, pero, por supuesto, depende del número de números a factorizar ...
Hah, les dejo los detalles de codificación a ustedes, amigos. :-)
Salud y salud,
- Alf
Mi otra respuesta es bastante larga y bastante diferente de esta, así que aquí hay algo más.
En lugar de simplemente filtrar los múltiplos de los dos primeros primos, o codificar todos los primos relevantes en un byte cada uno, este programa filtra los múltiplos de todos los primos que caben en ocho bits, específicamente del 2 al 211. Así que en lugar de pasar el 33% de Números, esto pasa alrededor del 10% al operador de la división.
Los números primos se mantienen en tres registros SSE, y sus módulos con el contador de ejecución se mantienen en otros tres. Si el módulo de cualquier primo con el contador es igual a cero, el contador no puede ser primo. Además, si cualquier módulo es igual a uno, entonces (contador + 2) no puede ser primo, etc., hasta (contador + 30). Los números pares se ignoran y las desviaciones como +3, +6 y +5 se omiten. El procesamiento de vectores permite actualizar dieciséis variables de tamaño byte a la vez.
Después de lanzar un fregadero de cocina lleno de micro optimizaciones (pero nada más específico de la plataforma que una directiva en línea), obtuve un aumento de rendimiento de 1.78x (24 s frente a 13.4 s en mi computadora portátil). Si se utiliza una biblioteca bignum (incluso una muy rápida), la ventaja es mayor. Vea a continuación una versión más legible, previa a la optimización.
/*
* factorize_sse.cpp
* Filter out multiples of the first 47 primes while factorizing a number.
*
* Created by David Krauss on 10/14/10.
*
*/
#include <cmath>
#include <sstream>
#include <iostream>
#include <xmmintrin.h>
using namespace std;
inline void remove_factor( unsigned long &n, unsigned long factor ) __attribute__((always_inline));
void remove_factor( unsigned long &n, unsigned long factor ) {
while ( n % factor == 0 ) {
n /= factor;
cout << factor << " ";
}
}
int main( int argc, char *argv[] ) {
unsigned long n;
if ( argc != 2
|| ( istringstream( argv[1] ) >> n >> ws ).rdstate() != ios::eofbit ) {
cerr << "Usage: " << argv[0] << " <number>/n";
return 1;
}
int primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211 };
for ( int *p = primes; p < primes + sizeof primes/sizeof *primes; ++ p ) {
remove_factor( n, * p );
}
//int histo[8] = {}, total = 0;
enum { bias = 15 - 128 };
__m128i const prime1 = _mm_set_epi8( 21, 21, 21, 22, 22, 26, 26, 17, 19, 23, 29, 31, 37, 41, 43, 47 ),
prime2 = _mm_set_epi8( 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127 ),
prime3 = _mm_set_epi8( 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211 ),
vbias = _mm_set1_epi8( bias ),
v3 = _mm_set1_epi8( 3+bias ), v5 = _mm_set1_epi8( 5+bias ), v6 = _mm_set1_epi8( 6+bias ), v8 = _mm_set1_epi8( 8+bias ),
v9 = _mm_set1_epi8( 9+bias ), v11 = _mm_set1_epi8( 11+bias ), v14 = _mm_set1_epi8( 14+bias ), v15 = _mm_set1_epi8( 15+bias );
__m128i mod1 = _mm_add_epi8( _mm_set_epi8( 3, 10, 17, 5, 16, 6, 19, 8, 9, 11, 14, 15, 18, 20, 21, 23 ), vbias ),
mod2 = _mm_add_epi8( _mm_set_epi8( 26, 29, 30, 33, 35, 36, 39, 41, 44, 48, 50, 51, 53, 54, 56, 63 ), vbias ),
mod3 = _mm_add_epi8( _mm_set_epi8( 65, 68, 69, 74, 75, 78, 81, 83, 86, 89, 90, 95, 96, 98, 99, 105 ), vbias );
for ( unsigned long factor = 1, limit = sqrtl( n ); factor <= limit + 30; factor += 30 ) {
if ( n == 1 ) goto done;
// up to 2^32, distribution of number candidates produced (0 up to 7) is
// 0.010841 0.0785208 0.222928 0.31905 0.246109 0.101023 0.0200728 0.00145546
unsigned candidates[8], *cand_pen = candidates;
* cand_pen = 6;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v3 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v3 ), _mm_cmpeq_epi8( mod3, v3 ) ) ) );
* cand_pen = 10;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v5 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v5 ), _mm_cmpeq_epi8( mod3, v5 ) ) ) );
* cand_pen = 12;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v6 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v6 ), _mm_cmpeq_epi8( mod3, v6 ) ) ) );
* cand_pen = 16;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v8 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v8 ), _mm_cmpeq_epi8( mod3, v8 ) ) ) );
* cand_pen = 18;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v9 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v9 ), _mm_cmpeq_epi8( mod3, v9 ) ) ) );
* cand_pen = 22;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v11 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v11 ), _mm_cmpeq_epi8( mod3, v11 ) ) ) );
* cand_pen = 28;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v14 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v14 ), _mm_cmpeq_epi8( mod3, v14 ) ) ) );
* cand_pen = 30;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v15 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v15 ), _mm_cmpeq_epi8( mod3, v15 ) ) ) );
/*++ total;
++ histo[ cand_pen - candidates ];
cout << "( ";
while ( cand_pen != candidates ) cout << factor + * -- cand_pen << " ";
cout << ")" << endl; */
mod1 = _mm_sub_epi8( mod1, _mm_set1_epi8( 15 ) ); // update residuals
__m128i mask1 = _mm_cmplt_epi8( mod1, _mm_set1_epi8( 1+bias ) );
mask1 = _mm_and_si128( mask1, prime1 ); // residual goes to zero or negative?
mod1 = _mm_add_epi8( mask1, mod1 ); // combine reset into zero or negative
mod2 = _mm_sub_epi8( mod2, _mm_set1_epi8( 15 ) );
__m128i mask2 = _mm_cmplt_epi8( mod2, _mm_set1_epi8( 1+bias ) );
mask2 = _mm_and_si128( mask2, prime2 );
mod2 = _mm_add_epi8( mask2, mod2 );
mod3 = _mm_sub_epi8( mod3, _mm_set1_epi8( 15 ) );
__m128i mask3 = _mm_cmplt_epi8( mod3, _mm_set1_epi8( 1+bias ) );
mask3 = _mm_and_si128( mask3, prime3 );
mod3 = _mm_add_epi8( mask3, mod3 );
if ( cand_pen - candidates == 0 ) continue;
remove_factor( n, factor + candidates[ 0 ] );
if ( cand_pen - candidates == 1 ) continue;
remove_factor( n, factor + candidates[ 1 ] );
if ( cand_pen - candidates == 2 ) continue;
remove_factor( n, factor + candidates[ 2 ] );
if ( cand_pen - candidates == 3 ) continue;
remove_factor( n, factor + candidates[ 3 ] );
if ( cand_pen - candidates == 4 ) continue;
remove_factor( n, factor + candidates[ 4 ] );
if ( cand_pen - candidates == 5 ) continue;
remove_factor( n, factor + candidates[ 5 ] );
if ( cand_pen - candidates == 6 ) continue;
remove_factor( n, factor + candidates[ 6 ] );
}
cout << n;
done:
/*cout << endl;
for ( int hx = 0; hx < 8; ++ hx ) cout << (float) histo[hx] / total << " ";*/
cout << endl;
}
.
dkrauss$ /usr/local/bin/g++ main.cpp -o factorize_sse -O3 --profile-use
dkrauss$ time ./factorize_sse 18446743721522234449
4294967231 4294967279
real 0m13.437s
user 0m13.393s
sys 0m0.011s
A continuación se muestra el primer borrador de los anteriores. Optimizaciones incluidas
- Haga que el contador cíclico se fusione de forma incondicional (evite una rama).
- Obtenga ILP desenrollando el bucle por un factor de 15, aumentando el paso a 30.
- Inspirado en su optimización.
- 30 parece ser un punto dulce, ya que elimina los múltiplos de 2, 3 y 5 de forma gratuita.
- Los intervalos de tiempo entre 5 y 15 pueden tener varios múltiplos en una zancada, así que ponga varias copias en una fase variada en el vector.
-
remove_factorfactor deremove_factor. - Cambie las llamadas condicionales e impredecibles
remove_factora escrituras de matrices no ramificadas. -
remove_factorcompletamente el bucle final con la llamadaremove_factory asegúrese de que la función esté siempre en línea.- Elimine la iteración final desenrollada ya que siempre hay un múltiplo de 7 entre los candidatos.
- Agrega otro vector que contenga todos los números primos restantes que sean lo suficientemente pequeños.
- Haga más espacio agregando un sesgo a los contadores y agregue otro vector. Ahora solo hay seis primos más que podrían filtrarse sin tener que saltar hasta 16 bits, y también me he quedado sin registros: el bucle necesita 3 vectores primos, 3 vectores de módulo, 8 constantes para buscar y una constante para cada uno. Incrementar por y hacer la verificación de rango. Eso hace 16.
- La ganancia es mínima (pero positiva) en esta aplicación, pero el propósito original de esta técnica era filtrar números primos para el tamiz en la otra respuesta. Así que estad atentos…
Versión legible:
/*
* factorize_sse.cpp
* Filter out multiples of the first 17 primes while factorizing a number.
*
* Created by David Krauss on 10/14/10.
*
*/
#include <cmath>
#include <sstream>
#include <iostream>
#include <xmmintrin.h>
using namespace std;
int main( int argc, char *argv[] ) {
unsigned long n;
if ( argc != 2
|| ( istringstream( argv[1] ) >> n >> ws ).rdstate() != ios::eofbit ) {
cerr << "Usage: " << argv[0] << " <number>/n";
return 1;
}
int primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59 };
for ( int *p = primes; p < primes + sizeof primes/sizeof *primes; ++ p ) {
while ( n % * p == 0 ) {
n /= * p;
cout << * p << " ";
}
}
if ( n != 1 ) {
__m128i mod = _mm_set_epi8( 1, 2, 3, 5, 6, 8, 9, 11, 14, 15, 18, 20, 21, 23, 26, 29 );
__m128i const prime = _mm_set_epi8( 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59 ),
one = _mm_set1_epi8( 1 );
for ( unsigned long factor = 1, limit = sqrtl( n ); factor < limit; ) {
factor += 2;
__m128i mask = _mm_cmpeq_epi8( mod, one ); // residual going to zero?
mod = _mm_sub_epi8( mod, one ); // update other residuals
if ( _mm_movemask_epi8( mask ) ) {
mask = _mm_and_si128( mask, prime ); // reset cycle if going to zero
mod = _mm_or_si128( mask, mod ); // combine reset into zeroed position
} else while ( n % factor == 0 ) {
n /= factor;
cout << factor << " ";
if ( n == 1 ) goto done;
}
}
cout << n;
}
done:
cout << endl;
}
No estoy seguro de cuán efectivos serían, pero en lugar de
while (workingNum%2==0)
Podrías hacerlo
while (workingNum & 1 == 0)
No estoy seguro de si gcc o msvc (o cualquier compilador que esté usando) es lo suficientemente inteligente como para cambiar la expresión workingNum% 2, pero es probable que esté haciendo la división y mirando el módulo en edx ...
Mi siguiente sugerencia podría ser completamente innecesaria dependiendo de su compilador, pero podría intentar poner workingNum / = 3 antes de la llamada al método. G ++ podría ser lo suficientemente inteligente como para ver la división innecesaria y simplemente usar el cociente en eax (también podría hacer esto dentro de su bucle más grande). O, un enfoque más completo (pero doloroso) sería ensamblar en línea el siguiente código.
while (workingNum%3==0)
{
factors.push_back(3);
workingNum/=3;
}
El compilador probablemente está traduciendo la operación modular en una división y luego mirando el módulo en edx. El problema es que estás realizando la división nuevamente (y dudo que el compilador esté viendo que simplemente realizaste la división implícitamente en la condición del bucle). Entonces, podrías ensamblar esto en línea. Esto presenta dos cuestiones:
1) El método de llamada para push_back (3). Esto podría alterar los registros, haciendo esto completamente innecesario.
2) Obtener un registro para workingNum, pero esto puede determinarse mediante una verificación modular inicial (para forzarlo en% eax), o en el momento actual, estará / debería estar en eax.
Puede escribir el bucle como (suponiendo que workingNum está en eax, y esta es la sintaxis AT&T de 32 bits, solo porque no sé el ensamblado de 64 bits o la sintaxis de Intel)
asm( "
movl $3, %ebx
WorkNumMod3Loop:
movl %eax, %ecx # just to be safe, backup workingNUm
movl $0, %edx # zero out edx
idivl $3 # divide by 3. quotient in eax, remainder in edx
cmpl $0, %edx # compare if it''s 0
jne AfterNumMod3Loop # if 0 is the remainder, jump out
# no need to perform division because new workingNum is already in eax
#factors.push_back(3) call
je WorkNumMod3Loop
AfterNumMod3Loop:
movl %ecx, %eax"
);
Debes mirar la salida de ensamblaje para esos bucles.Es posible que su compilador ya esté haciendo estas optimizaciones, pero lo dudo. No me sorprendería si poner el número de trabajo / = n antes de la llamada al método mejora un poco el rendimiento en algunos casos.
El método de factorización de Fermat es simple y rápido para encontrar pares de factores primos grandes siempre que lo detenga antes de que vaya demasiado lejos y se vuelva lento. Sin embargo, en mis pruebas en números aleatorios, tales casos han sido demasiado raros para ver alguna mejora.
... sin utilizar un enfoque probabalístico (por ejemplo, Miller-Rabin) para determinar la primalidad
Con una distribución uniforme, el 75% de sus entradas necesitarán mil millones de iteraciones de bucle, por lo que vale la pena gastar un millón de operaciones en técnicas menos deterministas, incluso si obtiene una respuesta no concluyente y tiene que volver a la división de prueba.
He encontrado que la variación Brent del Método Rho de Pollard es muy buena, aunque más complicada de codificar y comprender. El mejor ejemplo que he visto es de este foro de discusión . El método se basa en la suerte, pero a menudo ayuda a valer la pena.
La prueba de primalidad de Miller-Rabin es en realidad determinista hasta aproximadamente 10 ^ 15, lo que puede ahorrarle el problema de una búsqueda infructuosa.
Probé unas pocas docenas de variaciones y me decidí por lo siguiente para factorizar los valores int64:
- División de prueba sobre pequeños factores. (Yo uso los primeros 8000 primos precomputados.)
- 10 intentos con Rho de Pollard, cada uno con 16 iteraciones
- División de prueba a sqrt (n).
Tenga en cuenta que el Rho de Pollard encuentra factores que no son necesariamente primos, por lo que la recursión se puede usar para factorizarlos.