algorithm primes enumeration factors

algorithm - ¿Enumerar factores de un número directamente en orden ascendente sin ordenar?



primes enumeration (4)

[Estoy publicando esta respuesta solo como una formalidad para completar. Ya he elegido la respuesta de otra persona como la respuesta aceptada.]

Descripción general del algoritmo. Al buscar la forma más rápida de generar una lista de factores en la memoria (valores sin signo de 64 bits en mi caso), me decidí por un algoritmo híbrido que implementa una clasificación bidimensional, que aprovecha el conocimiento interno del ordenar claves (es decir, son solo enteros y, por lo tanto, se pueden calcular). El método específico es algo más cercano a un "ProxMapSort" pero con dos niveles de teclas (mayor, menor) en lugar de uno solo. La clave principal es simplemente el logaritmo en base 2 del valor. La clave menor es el número mínimo de dígitos más significativos del valor necesario para producir una extensión razonable en la segunda capa de cubos. Los factores se producen primero en una matriz de trabajo temporal de factores no clasificados. A continuación, se analiza su distribución y se asigna y completa una matriz de índices de depósito. Finalmente, los factores se almacenan directamente en su lugar en la matriz ordenada final, utilizando la ordenación por inserción. La gran mayoría de los depósitos tienen solo 1, 2 o 3 valores. Se dan ejemplos en el código fuente, que se adjunta al final de esta respuesta.

Complejidad espacial. La utilización de la memoria es aproximadamente 4 veces mayor que la de una solución basada en Quicksort, pero en realidad esto es bastante insignificante, ya que la memoria máxima utilizada en el peor de los casos (para la entrada de 64 bits) es de 5,5 MB, de los cuales 4,0 MB se mantienen por solo un pequeño puñado de milisegundos.

Complejidad de tiempo de ejecución. El rendimiento es mucho mejor que una solución codificada a mano basada en Quicksort: para números con un número no trivial de factores, es irregular aproximadamente 2.5 veces más rápido. En mis 3.4 GHz. Intel i7, produce 184.320 factores de 18.401.055.938.125.660.800 en orden ordenado en 0.0052 segundos, o alrededor de 96 ciclos de reloj por factor, o alrededor de 35 millones de factores por segundo.

Gráficos La memoria y el rendimiento en tiempo de ejecución se perfilaron para los 47,616 representantes canónicos de las clases de equivalencia de firmas primas de números hasta 2⁶⁴ – 1. Estos son los llamados "números altamente factorizables" en el espacio de búsqueda de 64 bits.

El tiempo de ejecución total es ~ 2.5 veces mejor que una solución basada en Quicksort para recuentos de factores no triviales, que se muestra a continuación en este diagrama de registro-registro:

El número de factores ordenados producidos por segundo es esencialmente la inversión de lo anterior. El rendimiento por factor disminuye después del punto óptimo de aproximadamente 2000 factores, pero no mucho. El rendimiento se ve afectado por los tamaños de caché L1, L2 y L3, así como por el recuento de factores primos únicos del número que se factoriza, que aumenta aproximadamente con el logaritmo del valor de entrada.

El uso máximo de memoria es una línea recta en este gráfico log-log, ya que es proporcional al logaritmo en base 2 del número de factores. Tenga en cuenta que el uso máximo de memoria es solo por un período de tiempo muy breve; Las matrices de trabajo de corta duración se descartan en milisegundos. Después de descartar los arreglos temporales, lo que queda es la lista final de factores, que es el mismo uso mínimo que se ve en la solución basada en Quicksort.

Código fuente. A continuación se adjunta un programa de prueba de concepto en el lenguaje de programación C (específicamente, C11). Se ha probado en x86-64 con Clang / LLVM, aunque también debería funcionar bien con GCC.

/*============================================================================== DESCRIPTION This is a small proof-of-concept program to test the idea of "sorting" factors using a form of bucket sort. The method is essentially a 2D version of ProxMapSort that has tuned for vast, nonlinear distributions using two keys (major, minor) rather than one. The major key is simply the floor of the base-2 logarithm of the value, and the minor key is derived from the most significant bits of the value. INPUT Input is given on the command line, either as a single argument giving the number to be factored or an even number of arguments giving the 2-tuples that comprise the prime-power factorization of the desired number. For example, the number 75600 = 2^4 x 3^3 x 5^2 x 7 can be given by the following list of arguments: 2 4 3 3 5 2 7 1 Note: If a single number is given, it will require factoring to produce its prime-power factorization. Since this is just a small test program, a very crude factoring method is used that is extremely fast for small prime factors but extremely slow for large prime factors. This is actually fine, because the largest factor lists occur with small prime factors anyway, and it is the production of large factor lists at which this program aims to be proficient. It is simply not interesting to be fast at producing the factor list of a number like 17293823921105882610 = 2 x 3 x 5 x 576460797370196087, because it has only 32 factors. Numbers with tens or hundreds of thousands of factors are much more interesting. OUTPUT Results are written to standard output. A list of factors in ascending order is produced, followed by runtime (in microseconds) required to generate the list (not including time to print it). STATISTICS Bucket size statistics for the 47616 canonical representatives of the prime signature equivalence classes of 64-bit numbers: ============================================================== Bucket size Total count of factored Total count of b numbers needing size b buckets of size b -------------------------------------------------------------- 1 47616 (100.0%) 514306458 (76.2%) 2 47427 (99.6%) 142959971 (21.2%) 3 43956 (92.3%) 16679329 (2.5%) 4 27998 (58.8%) 995458 (0.1%) 5 6536 (13.7%) 33427 (<0.1%) 6 400 (0.8%) 729 (<0.1%) 7 12 (<0.1%) 18 (<0.1%) -------------------------------------------------------------- ~ 47616 (100.0%) 674974643 (100.0%) -------------------------------------------------------------- Thus, no 64-bit number (of the input set) ever requires more than 7 buckets, and the larger the bucket size the less frequent it is. This is highly desirable. Note that although most numbers need at least 1 bucket of size 5, the vast majority of buckets (99.9%) are of size 1, 2, or 3, meaning that insertions are extremely efficient. Therefore, the use of insertion sort for the buckets is clearly the right choice and is arguably optimal for performance. AUTHOR Todd Lehman 2015/05/08 */ #include <inttypes.h> #include <limits.h> #include <stdbool.h> #include <stdlib.h> #include <stdio.h> #include <stdarg.h> #include <string.h> #include <time.h> #include <math.h> #include <assert.h> typedef unsigned int uint; typedef uint8_t uint8; typedef uint16_t uint16; typedef uint32_t uint32; typedef uint64_t uint64; #define ARRAY_CAPACITY(x) (sizeof(x) / sizeof((x)[0])) //----------------------------------------------------------------------------- // This structure is sufficient to represent the prime-power factorization of // all 64-bit values. The field names ω and Ω are dervied from the standard // number theory functions ω(n) and Ω(n), which count the number of unique and // non-unique prime factors of n, respectively. The field name d is derived // from the standard number theory function d(n), which counts the number of // divisors of n, including 1 and n. // // The maximum possible value here of ω is 15, which occurs for example at // n = 7378677391061896920 = 2^3 x 3^2 x 5 x 7 x 11 x 13 x 17 x 19 x 23 x 29 // 31 x 37 x 41 x 43 x 47, which has 15 unique prime factors. // // The maximum possible value of Ω here is 63, which occurs for example at // n = 2^63 and n = 2^62 x 3, both of which have 63 non-unique prime factors. // // The maximum possible value of d here is 184320, which occurs at // n = 18401055938125660800 = 2^7 x 3^4 x 5^2 x 7^2 x 11 x 13 x 17 x 19 x 23 x // 29 x 31 x 37 x 41. // // Maximum possible exponents when exponents are sorted in decreasing order: // // Index Maximum Bits Example of n // ----- ------- ---- -------------------------------------------- // 0 63 6 (2)^63 // 1 24 5 (2*3)^24 // 2 13 4 (2*3*5)^13 // 3 8 4 (2*3*5*7)^8 // 4 5 3 (2*3*5*7*11)^5 // 5 4 3 (2*3*5*7*11*13)^4 // 6 3 2 (2*3*5*7*11*13*17)^3 // 7 2 2 (2*3*5*7*11*13*17*19)^2 // 8 2 2 (2*3*5*7*11*13*17*19*23)^2 // 9 1 1 (2*3*5*7*11*13*17*19*23*29)^1 // 10 1 1 (2*3*5*7*11*13*17*19*23*29*31)^1 // 11 1 1 (2*3*5*7*11*13*17*19*23*29*31*37)^1 // 12 1 1 (2*3*5*7*11*13*17*19*23*29*31*37*41)^1 // 13 1 1 (2*3*5*7*11*13*17*19*23*29*31*37*41*43)^1 // 14 1 1 (2*3*5*7*11*13*17*19*23*29*31*37*41*43*47)^1 // ----- ------- ---- -------------------------------------------- // 15 63 37 // #pragma pack(push, 8) typedef struct { uint8 e[16]; // Exponents. uint64 p[16]; // Primes in increasing order. uint8 ω; // Count of prime factors without multiplicity. uint8 Ω; // Count of prime factors with multiplicity. uint32 d; // Count of factors of n, including 1 and n. uint64 n; // Value of n on which all other fields of this struct depend. } PrimePowerFactorization; // 176 bytes with 8-byte packing #pragma pack(pop) #define MAX_ω 15 #define MAX_Ω 63 //----------------------------------------------------------------------------- // Fatal error: print error message and abort. void fatal_error(const char *format, ...) { va_list args; va_start(args, format); vfprintf(stderr, format, args); exit(1); } //----------------------------------------------------------------------------- // Compute 64-bit 2-adic integer inverse. uint64 uint64_inv(const uint64 x) { assert(x != 0); uint64 y = 1; for (uint i = 0; i < 6; i++) // 6 = log2(log2(2**64)) = log2(64) y = y * (2 - (x * y)); return y; } //------------------------------------------------------------------------------ // Compute 2 to arbitrary power. This is just a portable and abstract way to // write a left-shift operation. Note that the use of the UINT64_C macro here // is actually required, because the result of 1U<<x is not guaranteed to be a // 64-bit result; on a 32-bit compiler, 1U<<32 is 0 or is undefined. static inline uint64 uint64_pow2(x) { return UINT64_C(1) << x; } //------------------------------------------------------------------------------ // Deduce native word size (int, long, or long long) for 64-bit integers. // This is needed for abstracting certain compiler-specific intrinsic functions. #if UINT_MAX == 0xFFFFFFFFFFFFFFFFU #define UINT64_IS_U #elif ULONG_MAX == 0xFFFFFFFFFFFFFFFFUL #define UINT64_IS_UL #elif ULLONG_MAX == 0xFFFFFFFFFFFFFFFFULL #define UINT64_IS_ULL #else //error "Unable to deduce native word size of 64-bit integers." #endif //------------------------------------------------------------------------------ // Define abstracted intrinsic function for counting leading zeros. Note that // the value is well-defined for nonzero input but is compiler-specific for // input of zero. #if defined(UINT64_IS_U) && __has_builtin(__builtin_clz) #define UINT64_CLZ(x) __builtin_clz(x) #elif defined(UINT64_IS_UL) && __has_builtin(__builtin_clzl) #define UINT64_CLZ(x) __builtin_clzl(x) #elif defined(UINT64_IS_ULL) && __has_builtin(__builtin_clzll) #define UINT64_CLZ(x) __builtin_clzll(x) #else #undef UINT64_CLZ #endif //------------------------------------------------------------------------------ // Compute floor of base-2 logarithm y = log_2(x), where x > 0. Uses fast // intrinsic function if available; otherwise resorts to hand-rolled method. static inline uint uint64_log2(uint64 x) { assert(x > 0); #if defined(UINT64_CLZ) return 63 - UINT64_CLZ(x); #else #define S(k) if ((x >> k) != 0) { y += k; x >>= k; } uint y = 0; S(32); S(16); S(8); S(4); S(2); S(1); return y; #undef S #endif } //------------------------------------------------------------------------------ // Compute major key, given a nonzero number. The major key is simply the // floor of the base-2 logarithm of the number. static inline uint major_key(const uint64 n) { assert(n > 0); uint k1 = uint64_log2(n); return k1; } //------------------------------------------------------------------------------ // Compute minor key, given a nonzero number, its major key, k1, and the // bit-size b of major bucket k1. The minor key, k2, is is computed by first // removing the most significant 1-bit from the number, because it adds no // information, and then extracting the desired number of most significant bits // from the remainder. For example, given the number n=1463 and a major bucket // size of b=6 bits, the keys are computed as follows: // // Step 0: Given number n = 0b10110110111 = 1463 // // Step 1: Compute major key: k1 = floor(log_2(n)) = 10 // // Step 2: Remove high-order 1-bit: n'' = 0b0110110111 = 439 // // Step 3: Compute minor key: k2 = n'' >> (k1 - b) // = 0b0110110111 >> (10 - 6) // = 0b0110110111 >> 4 // = 0b011011 // = 27 static inline uint minor_key(const uint64 n, const uint k1, const uint b) { assert(n > 0); assert(k1 >= 0); assert(b > 0); const uint k2 = (uint)((n ^ uint64_pow2(k1)) >> (k1 - b)); return k2; } //------------------------------------------------------------------------------ // Raw unsorted factor. #pragma push(pack, 4) typedef struct { uint64 n; // Value of factor. uint32 k1; // Major key. uint32 k2; // Minor key. } UnsortedFactor; #pragma pop(pack) //------------------------------------------------------------------------------ // Compute sorted list of factors, given a prime-power factorization. static uint64 memory_usage; uint64 *compute_factors(const PrimePowerFactorization ppf) { memory_usage = 0; if (ppf.n == 0) return NULL; uint64 *sorted_factors = calloc(ppf.d, sizeof(*sorted_factors)); if (!sorted_factors) fatal_error("Failed to allocate array of %"PRIu32" factors.", ppf.d); memory_usage += ppf.d * sizeof(*sorted_factors); UnsortedFactor *unsorted_factors = malloc(ppf.d * sizeof(*unsorted_factors)); if (!unsorted_factors) fatal_error("Failed to allocate array of %"PRIu32" factors.", ppf.d); memory_usage += ppf.d * sizeof(*unsorted_factors); // These arrays are indexed by the major key of a number. uint32 major_counts[64]; // Counts of factors in major buckets. uint32 major_spans[64]; // Counts rounded up to power of 2. uint32 major_bits[64]; // Base-2 logarithm of bucket size. uint32 major_indexes[64]; // Indexes into minor array. memset(major_counts, 0, sizeof(major_counts)); memset(major_spans, 0, sizeof(major_spans)); memset(major_bits, 0, sizeof(major_bits)); memset(major_indexes, 0, sizeof(major_indexes)); // --- Step 1: Produce unsorted list of factors from prime-power // factorization. At the same time, count groups of factors by their // major keys. { // This array is for counting in the multi-radix number system dictated by // the exponents of the prime-power factorization. An invariant is that // e[i] <= ppf.e[i] for all i (0 < i <ppf.ω). uint8 e[MAX_ω]; for (uint i = 0; i < ppf.ω; i++) e[i] = 0; // Initialize inverse-prime-powers. This array allows for division by // p[i]**e[i] extremely quickly in the main loop below. Note that 2-adic // inverses are not defined for even numbers (of which 2 is the only prime), // so powers of 2 must be handled specially. uint64 pe_inv[MAX_ω]; for (uint i = 0; i < ppf.ω; i++) { uint64 pe = 1; for (uint j = 1; j <= ppf.e[i]; j++) pe *= ppf.p[i]; pe_inv[i] = uint64_inv(pe); } uint64 n = 1; // Current factor accumulator. for (uint k = 0; k < ppf.d; k++) // k indexes into unsorted_factors[]. { //printf("unsorted_factors[%u] = %"PRIu64" j = %u/n", k, n, j); assert(ppf.n % n == 0); unsorted_factors[k].n = n; uint k1 = major_key(n); assert(k1 < ARRAY_CAPACITY(major_counts)); unsorted_factors[k].k1 = k1; major_counts[k1] += 1; // Increment the remainder of the multi-radix number e[]. for (uint i = 0; i < ppf.ω; i++) { if (e[i] == ppf.e[i]) // Carrying is occurring. { if (ppf.p[i] == 2) n >>= ppf.e[i]; // Divide n by 2 ** ppf.e[i]. else n *= pe_inv[i]; // Divide n by ppf.p[i] ** ppf.e[i]. e[i] = 0; } else // Carrying is not occurring. { n *= ppf.p[i]; e[i] += 1; break; } } } assert(n == 1); // n always cycles back to 1, not to ppf.n. assert(unsorted_factors[ppf.d-1].n == ppf.n); } // --- Step 2: Define the major bits array, the major spans array, the major // index array, and count the total spans. uint32 total_spans = 0; { uint32 k = 0; for (uint k1 = 0; k1 < ARRAY_CAPACITY(major_counts); k1++) { uint32 count = major_counts[k1]; uint32 bits = (count <= 1)? count : uint64_log2(count - 1) + 1; major_bits[k1] = bits; major_spans[k1] = (count > 0)? (UINT32_C(1) << bits) : 0; major_indexes[k1] = k; k += major_spans[k1]; } total_spans = k; } // --- Step 3: Allocate and populate the minor counts array. Note that it // must be initialized to zero. uint32 *minor_counts = calloc(total_spans, sizeof(*minor_counts)); if (!minor_counts) fatal_error("Failed to allocate array of %"PRIu32" counts.", total_spans); memory_usage += total_spans * sizeof(*minor_counts); for (uint k = 0; k < ppf.d; k++) { const uint64 n = unsorted_factors[k].n; const uint k1 = unsorted_factors[k].k1; const uint k2 = minor_key(n, k1, major_bits[k1]); assert(k2 < major_spans[k1]); unsorted_factors[k].k2 = k2; minor_counts[major_indexes[k1] + k2] += 1; } // --- Step 4: Define the minor indexes array. // // NOTE: Instead of allocating a separate array, the earlier-allocated array // of minor indexes is simply repurposed here using an alias. uint32 *minor_indexes = minor_counts; // Alias the array for repurposing. { uint32 k = 0; for (uint i = 0; i < total_spans; i++) { uint32 count = minor_counts[i]; // This array is the same array... minor_indexes[i] = k; // ...as this array. k += count; } } // --- Step 5: Populate the sorted factors array. Note that the array must // be initialized to zero earlier because values of zero are used // as sentinels in the bucket lists. for (uint32 i = 0; i < ppf.d; i++) { uint64 n = unsorted_factors[i].n; const uint k1 = unsorted_factors[i].k1; const uint k2 = unsorted_factors[i].k2; // Insert factor into bucket using insertion sort (which happens to be // extremely fast because we know the bucket sizes are always very small). uint32 k; for (k = minor_indexes[major_indexes[k1] + k2]; sorted_factors[k] != 0; k++) { assert(k < ppf.d); if (sorted_factors[k] > n) { uint64 t = sorted_factors[k]; sorted_factors[k] = n; n = t; } } sorted_factors[k] = n; } // --- Step 6: Validate array of sorted factors. { for (uint32 k = 1; k < ppf.d; k++) { if (sorted_factors[k] == 0) fatal_error("Produced a factor of 0 at index %"PRIu32".", k); if (ppf.n % sorted_factors[k] != 0) fatal_error("Produced non-factor %"PRIu64" at index %"PRIu32".", sorted_factors[k], k); if (sorted_factors[k-1] == sorted_factors[k]) fatal_error("Duplicate factor %"PRIu64" at index %"PRIu32".", sorted_factors[k], k); if (sorted_factors[k-1] > sorted_factors[k]) fatal_error("Out-of-order factors %"PRIu64" and %"PRIu64" " "at indexes %"PRIu32" and %"PRIu32".", sorted_factors[k-1], sorted_factors[k], k-1, k); } } free(minor_counts); free(unsorted_factors); return sorted_factors; } //------------------------------------------------------------------------------ // Compute prime-power factorization of a 64-bit value. Note that this function // is designed to be fast *only* for numbers with very simple factorizations, // e.g., those that produce large factor lists. Do not attempt to factor // large semiprimes with this function. (The author does know how to factor // large numbers efficiently; however, efficient factorization is beyond the // scope of this small test program.) PrimePowerFactorization compute_ppf(const uint64 n) { PrimePowerFactorization ppf; if (n == 0) { ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 }; } else if (n == 1) { ppf = (PrimePowerFactorization){ .p = { 1 }, .e = { 1 }, .ω = 1, .Ω = 1, .d = 1, .n = 1 }; } else { ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = n }; uint64 m = n; for (uint64 p = 2; p * p <= m; p += 1 + (p > 2)) { if (m % p == 0) { assert(ppf.ω <= MAX_ω); ppf.p[ppf.ω] = p; ppf.e[ppf.ω] = 0; while (m % p == 0) { m /= p; ppf.e[ppf.ω] += 1; } ppf.d *= (1 + ppf.e[ppf.ω]); ppf.Ω += ppf.e[ppf.ω]; ppf.ω += 1; } } if (m > 1) { assert(ppf.ω <= MAX_ω); ppf.p[ppf.ω] = m; ppf.e[ppf.ω] = 1; ppf.d *= 2; ppf.Ω += 1; ppf.ω += 1; } } return ppf; } //------------------------------------------------------------------------------ // Parse prime-power factorization from a list of ASCII-encoded base-10 strings. // The values are assumed to be 2-tuples (p,e) of prime p and exponent e. // Primes must not exceed 2^64 - 1. Exponents must not exceed 2^8 - 1. The // constructed value must not exceed 2^64 - 1. PrimePowerFactorization parse_ppf(const uint pairs, const char *const values[]) { assert(pairs <= MAX_ω); PrimePowerFactorization ppf; if (pairs == 0) { ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 }; } else { ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = 1 }; for (uint i = 0; i < pairs; i++) { ppf.p[i] = (uint64)strtoumax(values[(i*2)+0], NULL, 10); ppf.e[i] = (uint8)strtoumax(values[(i*2)+1], NULL, 10); // Validate prime value. if (ppf.p[i] < 2) // (Ideally this would actually do a primality test.) fatal_error("Factor %"PRIu64" is invalid.", ppf.p[i]); // Accumulate count of unique prime factors. if (ppf.ω > UINT8_MAX - 1) fatal_error("Small-omega overflow at factor %"PRIu64"^%"PRIu8".", ppf.p[i], ppf.e[i]); ppf.ω += 1; // Accumulate count of total prime factors. if (ppf.Ω > UINT8_MAX - ppf.e[i]) fatal_error("Big-omega wverflow at factor %"PRIu64"^%"PRIu8".", ppf.p[i], ppf.e[i]); ppf.Ω += ppf.e[i]; // Accumulate total divisor count. if (ppf.d > UINT32_MAX / (1 + ppf.e[i])) fatal_error("Divisor count overflow at factor %"PRIu64"^%"PRIu8".", ppf.p[i], ppf.e[i]); ppf.d *= (1 + ppf.e[i]); // Accumulate value. for (uint8 k = 1; k <= ppf.e[i]; k++) { if (ppf.n > UINT64_MAX / ppf.p[i]) fatal_error("Value overflow at factor %"PRIu64".", ppf.p[i]); ppf.n *= ppf.p[i]; } } } return ppf; } //------------------------------------------------------------------------------ // Main control. Parse command line and produce list of factors. int main(const int argc, const char *const argv[]) { PrimePowerFactorization ppf; uint values = (uint)argc - 1; // argc is always guaranteed to be at least 1. if (values == 1) { ppf = compute_ppf((uint64)strtoumax(argv[1], NULL, 10)); } else { if (values % 2 != 0) fatal_error("Odd number of arguments (%u) given.", values); uint pairs = values / 2; ppf = parse_ppf(pairs, &argv[1]); } // Run for (as close as possible to) a fixed amount of time, tallying the // elapsed CPU time. uint64 iterations = 0; double cpu_time = 0.0; const double cpu_time_limit = 0.05; while (cpu_time < cpu_time_limit) { clock_t clock_start = clock(); uint64 *factors = compute_factors(ppf); clock_t clock_end = clock(); cpu_time += (double)(clock_end - clock_start) / (double)CLOCKS_PER_SEC; if (++iterations == 1) { for (uint32 i = 0; i < ppf.d; i++) printf("%"PRIu64"/n", factors[i]); } if (factors) free(factors); } // Print the average amount of CPU time required for each iteration. uint mem_scale = (memory_usage >= 1e9)? 9: (memory_usage >= 1e6)? 6: (memory_usage >= 1e3)? 3: 0; char *mem_units = (mem_scale == 9)? "GB": (mem_scale == 6)? "MB": (mem_scale == 3)? "KB": "B"; printf("%"PRIu64" %"PRIu32" factors %.6f ms %.3f ns/factor %.3f %s/n", ppf.n, ppf.d, cpu_time/iterations * 1e3, cpu_time/iterations * 1e9 / (double)(ppf.d? ppf.d : 1), (double)memory_usage / pow(10, mem_scale), mem_units); return 0; }

¿Existe un algoritmo eficiente para enumerar los factores de un número n , en orden ascendente, sin ordenar? Por "eficiente" quiero decir:

  1. El algoritmo evita una búsqueda de divisores por fuerza bruta comenzando con la factorización de potencia primaria de n .

  2. La complejidad de tiempo de ejecución del algoritmo es O ( d log₂ d ) o mejor, donde d es el recuento de divisores de n .

  3. La complejidad espacial del algoritmo es O ( d ).

  4. El algoritmo evita una operación de clasificación. Es decir, los factores se producen en orden en lugar de producirse fuera de orden y luego se ordenan. Aunque enumerar usando un enfoque recursivo simple y luego ordenar es O ( d log₂ d ), hay un costo muy feo para todos los accesos de memoria involucrados en la clasificación.

Un ejemplo trivial es n = 360 = 2³ × 3² × 5, que tiene d = 24 factores: {1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 18, 20, 24, 30, 36, 40, 45, 60, 72, 90, 120, 180, 360}.

Un ejemplo más serio es n = 278282512406132373381723386382308832000 = 2⁸ × 3⁴ × 5³ × 7² × 11² × 13² × 17 × 19 × 23 × 29 × 31 × 37 × 41 × 43 × 47 × 53 × 59 × 61 × 67 × 71 × 73 × 79, que tiene d = 318504960 factores (¡obviamente demasiados para enumerarlos aquí!). Este número, por cierto, tiene el mayor número de factores de cualquier número hasta 2 ^ 128.

Podría jurar que vi una descripción de dicho algoritmo hace unas semanas, con código de muestra, pero ahora parece que no puedo encontrarlo en ningún lado. Usó algún truco de magia para mantener una lista de índices progenitores en la lista de salida para cada factor primo. (Actualización: Estaba confundiendo la generación de factores con los números de Hamming, que funcionan de manera similar).

Actualizar

Terminé usando una solución que es O ( d ) en tiempo de ejecución, tiene una carga de memoria extremadamente baja creando la salida de O ( d ) en el lugar, y es significativamente más rápido que cualquier otro método que conozco. He publicado esta solución como respuesta, con código fuente C. Es una versión simplificada y altamente optimizada de un hermoso algoritmo presentado aquí en Haskell por otro colaborador, Will Ness. Seleccioné la respuesta de Will como la respuesta aceptada, ya que proporcionó una solución muy elegante que cumplía con todos los requisitos establecidos originalmente.

En resumen: extraiga repetidamente el siguiente factor más pequeño de un montón, luego empuje hacia atrás cada múltiplo de ese factor que todavía es un factor de n. Use un truco para evitar que surjan duplicados, de modo que el tamaño del montón nunca exceda d. La complejidad temporal es O (kd log d), donde k es el número de factores primos distintos.

La propiedad clave que utilizamos es que si x e y son factores de n con y = x * p para algún factor p> = 2, es decir, si los factores primos de x son un subconjunto propio de los factores primos de y - entonces x <y. Esto significa que siempre es seguro generar x antes de que cualquier y se inserte en el montón. El montón se usa efectivamente solo para comparar factores donde ninguno de los dos es un subconjunto del otro.

Un primer intento: los duplicados ralentizan las cosas

Será útil describir primero un algoritmo que produce la respuesta correcta, pero también produce muchos duplicados:

  1. Establecer prev = NULL.
  2. Inserte 1 en un montón H.
  3. Extraiga la parte superior del montón t de H. Si el montón está vacío, deténgase.
  4. Si t == prev, entonces goto 3. [EDITAR: Corregido]
  5. Salida t.
  6. Establecer prev = t.
  7. Para cada factor primo distinto p de n:
    • Si n% (t * p) == 0 (es decir, si t * p sigue siendo un factor de n), presione t * p sobre H.
  8. Ir a 3.

El único problema con el algoritmo anterior es que puede generar el mismo factor muchas veces. Por ejemplo, si n = 30, entonces el factor 15 se generará como hijo del factor 5 (multiplicando por el factor primo 3), y también como hijo del factor 3 (multiplicando por 5). Una forma de evitar esto es notar que cualquier duplicado debe leerse en un bloque contiguo cuando alcanza la parte superior del montón, por lo que simplemente puede verificar si la parte superior del montón es igual al valor recién extraído y mantener extrayéndolo y descartándolo si es así. Pero es posible un mejor enfoque:

Matar duplicados en la fuente

¿De cuántas maneras se puede generar un factor x? Primero considere el caso en el que x no contiene factores primos con multiplicidad> 1. En ese caso, si contiene m factores primos distintos, entonces hay m-1 factores "principales" que lo generarán como un "elemento secundario" en el caso anterior. algoritmo: cada uno de estos padres consiste en un subconjunto de m-1 de los factores primos m, siendo el factor primo restante el que se agrega al niño. (Si x tiene un factor primo con multiplicidad> 1, entonces hay de hecho m padres.) Si tuviéramos una forma de decidir exactamente cuál de estos padres es el "elegido" que realmente genera x como niño, y esta regla dio como resultado una prueba que se podría aplicar a cada padre en el momento en que se abre el padre , luego podríamos evitar crear duplicados en primer lugar.

Podemos usar la siguiente regla: para cualquier x dada, elija el padre potencial y que no tiene el factor m más grande de x. Esto hace una regla simple: un padre y produce un hijo x si y solo si x = y * p para alguna p mayor o igual a cualquier factor primo que ya esté en y. Esto es fácil de comprobar: simplemente recorra los factores primos en orden decreciente, generando hijos para cada uno, hasta que lleguemos a un factor primo que ya divide y. En el ejemplo anterior, el padre 3 producirá 15, pero el padre 5 no (porque 3 <5), por lo que 15 de hecho se produce solo una vez. Para n = 30, el árbol completo se ve así:

1 /|/ 2 3 5 /| / 6 10 15 | 30

Observe que cada factor se genera exactamente una vez.

El nuevo algoritmo libre de duplicados es el siguiente:

  1. Inserte 1 en un montón H.
  2. Extraiga la parte superior del montón t de H. Si el montón está vacío, deténgase.
  3. Salida t.
  4. Para cada factor primo distinto p de n en orden decreciente :
    • Si n% (t * p) == 0 (es decir, si t * p sigue siendo un factor de n), presione t * p sobre H.
    • Si t% p == 0 (es decir, si t ya contiene p como factor), deténgase.
  5. Ir a 2.

Esta respuesta proporciona una implementación en C que creo que es la más rápida y eficiente en memoria.

Descripción general del algoritmo. Este algoritmo se basa en el enfoque de fusión ascendente introducido por Will Ness en otra respuesta , pero se simplifica aún más para que las listas que se fusionan en realidad nunca existan en ningún lugar de la memoria. El elemento principal de cada lista se prepara y mantiene en una pequeña matriz, mientras que todos los demás elementos de las listas se construyen sobre la marcha según sea necesario. Este uso de "listas fantasmas" —figuras de la imaginación del código en ejecución— reduce en gran medida la huella de la memoria, así como el volumen de accesos a la memoria, tanto de lectura como de escritura, y también mejora la localidad espacial, que a su vez aumenta significativamente la velocidad del algoritmo Los factores en cada nivel se escriben directamente en su lugar de descanso final en la matriz de salida, en orden.

La idea básica es producir los factores usando la inducción matemática en la factorización de potencia primaria. Por ejemplo, para producir los factores de 360, los factores de 72 se calculan y luego se multiplican por las potencias relevantes de 5, en este caso {1,5}; para producir los factores de 72, los factores de 8 se calculan y luego se multiplican por las potencias relevantes de 3, en este caso {1,3,9}; Para producir los factores de 8, el caso base 1 se multiplica por las potencias relevantes de 2, en este caso {1,2,4,8}. Por lo tanto, en cada paso inductivo, se calcula un producto cartesiano entre conjuntos de factores existentes y conjuntos de potencias primarias, y los resultados se reducen a enteros mediante multiplicación.

A continuación se muestra una ilustración para el número 360. De izquierda a derecha hay celdas de memoria; una fila representa un paso de tiempo. El tiempo se representa como un flujo vertical.

Complejidad espacial. Las estructuras de datos temporales para producir los factores son extremadamente pequeñas: solo se utiliza el espacio O (log₂ ( n )), donde n es el número cuyos factores se están generando. Por ejemplo, en la implementación de 128 bits de este algoritmo en C, un número como 333,939,014,887,358,848,058,068,063,658,770,598,400 (cuyo logaritmo de base 2 es ≈127.97) asigna 5.1 GB para almacenar la lista de sus 318,504,960 factores, pero usa solo 360 bytes de sobrecarga adicional para producir esa lista. Como máximo, se necesitan aproximadamente 5 KB de sobrecarga para cualquier número de 128 bits.

Complejidad de tiempo de ejecución. El tiempo de ejecución depende completamente de los exponentes de la factorización de potencia principal (por ejemplo, la firma principal ). Para obtener mejores resultados, los exponentes más grandes deben procesarse primero y los exponentes más pequeños al final, de modo que el tiempo de ejecución esté dominado en las etapas finales por fusiones de baja complejidad, que en la práctica a menudo resultan ser fusiones binarias. El tiempo de ejecución asintótico es O ( c ( n ) d ( n )), donde d ( n ) es el recuento divisor de n y donde c ( n ) depende constantemente de la firma principal de n . Es decir, cada clase de equivalencia asociada con una firma principal tiene una constante diferente. Las firmas primas dominadas por pequeños exponentes tienen constantes más pequeñas; Las firmas primas dominadas por grandes exponentes tienen constantes más grandes. Por lo tanto, la complejidad del tiempo de ejecución está agrupada por la firma principal.

Gráficos El rendimiento en tiempo de ejecución se perfiló en 3.4 GHz. Intel Core i7 para un muestreo de 66,591 valores de n que tienen factores d ( n ) para d ( n ) únicos de hasta 160 millones. El mayor valor de n perfilado para este gráfico fue 14,550,525,518,294,259,162,294,162,737,840,640,000, produciendo 159,744,000 factores en 2.35 segundos.

El número de factores ordenados producidos por segundo es esencialmente la inversión de lo anterior. La agrupación por firma principal es evidente en los datos. El rendimiento también se ve afectado por los tamaños de caché L1, L2 y L3, así como por las limitaciones físicas de RAM.

Código fuente. A continuación se adjunta un programa de prueba de concepto en el lenguaje de programación C (específicamente, C11). Se ha probado en x86-64 con Clang / LLVM, aunque también debería funcionar bien con GCC.

/*============================================================================== DESCRIPTION This is a small proof-of-concept program to test the idea of generating the factors of a number in ascending order using an ultra-efficient sortless method. INPUT Input is given on the command line, either as a single argument giving the number to be factored or an even number of arguments giving the 2-tuples that comprise the prime-power factorization of the desired number. For example, the number 75600 = 2^4 x 3^3 x 5^2 x 7 can be given by the following list of arguments: 2 4 3 3 5 2 7 1 Note: If a single number is given, it will require factoring to produce its prime-power factorization. Since this is just a small test program, a very crude factoring method is used that is extremely fast for small prime factors but extremely slow for large prime factors. This is actually fine, because the largest factor lists occur with small prime factors anyway, and it is the production of large factor lists at which this program aims to be proficient. It is simply not interesting to be fast at producing the factor list of a number like 17293823921105882610 = 2 x 3 x 5 x 576460797370196087, because it has only 32 factors. Numbers with tens or hundreds of thousands of factors are much more interesting. OUTPUT Results are written to standard output. A list of factors in ascending order is produced, followed by runtime required to generate the list (not including time to print it). AUTHOR Todd Lehman 2015/05/10 */ //----------------------------------------------------------------------------- #include <inttypes.h> #include <limits.h> #include <stdbool.h> #include <stdlib.h> #include <stdio.h> #include <stdarg.h> #include <string.h> #include <ctype.h> #include <time.h> #include <math.h> #include <assert.h> //----------------------------------------------------------------------------- typedef unsigned int uint; typedef uint8_t uint8; typedef uint16_t uint16; typedef uint32_t uint32; typedef uint64_t uint64; typedef __uint128_t uint128; #define UINT128_MAX (uint128)(-1) #define UINT128_MAX_STRLEN 39 //----------------------------------------------------------------------------- #define ARRAY_CAPACITY(x) (sizeof(x) / sizeof((x)[0])) //----------------------------------------------------------------------------- // This structure encode a single prime-power pair for the prime-power // factorization of numbers, for example 3 to the 4th power. #pragma pack(push, 8) typedef struct { uint128 p; // Prime. uint8 e; // Power (exponent). } PrimePower; // 24 bytes using 8-byte packing #pragma pack(pop) //----------------------------------------------------------------------------- // Prime-power factorization structure. // // This structure is sufficient to represent the prime-power factorization of // all 128-bit values. The field names ω and Ω are dervied from the standard // number theory functions ω(n) and Ω(n), which count the number of unique and // non-unique prime factors of n, respectively. The field name d is derived // from the standard number theory function d(n), which counts the number of // divisors of n, including 1 and n. // // The maximum possible value here of ω is 26, which occurs at // n = 232862364358497360900063316880507363070 = 2 x 3 x 5 x 7 x 11 x 13 x 17 x // 19 x 23 x 29 x 31 x 37 x 41 x 43 x 47 x 53 x 59 x 61 x 67 x 71 x 73 x 79 x // 83 x 89 x 97 x 101, which has 26 unique prime factors. // // The maximum possible value of Ω here is 127, which occurs at n = 2^127 and // n = 2^126 x 3, both of which have 127 non-unique prime factors. // // The maximum possible value of d here is 318504960, which occurs at // n = 333939014887358848058068063658770598400 = 2^9 x 3^5 x 5^2 x 7^2 x 11^2 x // 13^2 x 17 x 19 x 23 x 29 x 31 x 37 x 41 x 43 x 47 x 53 x 59 x 61 x 67 x 71 x // 73 x 79. // #pragma pack(push, 8) typedef struct { PrimePower f[32]; // Primes and their exponents. uint8 ω; // Count of prime factors without multiplicity. uint8 Ω; // Count of prime factors with multiplicity. uint32 d; // Count of factors of n, including 1 and n. uint128 n; // Value of n on which all other fields depend. } PrimePowerFactorization; // 656 bytes using 8-byte packing #pragma pack(pop) #define MAX_ω 26 #define MAX_Ω 127 //----------------------------------------------------------------------------- // Fatal error: print error message and abort. void fatal_error(const char *format, ...) { va_list args; va_start(args, format); vfprintf(stderr, format, args); exit(1); } //------------------------------------------------------------------------------ uint128 uint128_from_string(const char *const str) { assert(str != NULL); uint128 n = 0; for (int i = 0; isdigit(str[i]); i++) n = (n * 10) + (uint)(str[i] - ''0''); return n; } //------------------------------------------------------------------------------ void uint128_to_string(const uint128 n, char *const strbuf, const uint strbuflen) { assert(strbuf != NULL); assert(strbuflen >= UINT128_MAX_STRLEN + 1); // Extract digits into string buffer in reverse order. uint128 a = n; char *s = strbuf; do { *(s++) = ''0'' + (uint)(a % 10); a /= 10; } while (a != 0); *s = ''/0''; // Reverse the order of the digits. uint l = strlen(strbuf); for (uint i = 0; i < l/2; i++) { char t = strbuf[i]; strbuf[i] = strbuf[l-1-i]; strbuf[l-1-i] = t; } // Verify result. assert(uint128_from_string(strbuf) == n); } //------------------------------------------------------------------------------ char *uint128_to_static_string(const uint128 n, const uint i) { static char str[2][UINT128_MAX_STRLEN + 1]; assert(i < ARRAY_CAPACITY(str)); uint128_to_string(n, str[i], ARRAY_CAPACITY(str[i])); return str[i]; } //------------------------------------------------------------------------------ // Compute sorted list of factors, given a prime-power factorization. uint128 *compute_factors(const PrimePowerFactorization ppf) { const uint128 n = ppf.n; const uint d = (uint)ppf.d; const uint ω = (uint)ppf.ω; if (n == 0) return NULL; uint128 *factors = malloc((d + 1) * sizeof(*factors)); if (!factors) fatal_error("Failed to allocate array of %u factors.", d); uint128 *const factors_end = &factors[d]; // --- Seed the factors[] array. factors_end[0] = 0; // Dummy value to simplify looping in bottleneck code. factors_end[-1] = 1; // Seed value. if (n == 1) return factors; // --- Iterate over all prime factors. uint range = 1; for (uint i = 0; i < ω; i++) { const uint128 p = ppf.f[i].p; const uint e = ppf.f[i].e; // --- Initialize phantom input lists and output list. assert(e < 128); assert(range < d); uint128 *restrict in[128]; uint128 pe[128], f[128]; for (uint j = 0; j <= e; j++) { in[j] = &factors[d - range]; pe[j] = (j == 0)? 1 : pe[j-1] * p; f[j] = pe[j]; } uint active_list_count = 1 + e; range *= 1 + e; uint128 *restrict out = &factors[d - range]; // --- Merge phantom input lists to output list, until all input lists are // extinguished. while (active_list_count > 0) { if (active_list_count == 1) { assert(out == in[0]); while (out != factors_end) *(out++) *= pe[0]; in[0] = out; } else if (active_list_count == 2) { // This section of the code is the bottleneck of the entire factor- // producing algorithm. Other portions need to be fast, but this // *really* needs to be fast; therefore, it has been highly optimized. // In fact, it is by far most frequently the case here that pe[0] is 1, // so further optimization is warranted in this case. uint128 f0 = f[0], f1 = f[1]; uint128 *in0 = in[0], *in1 = in[1]; const uint128 pe0 = pe[0], pe1 = pe[1]; if (pe[0] == 1) { while (true) { if (f0 < f1) { *(out++) = f0; f0 = *(++in0); if (in0 == factors_end) break; } else { *(out++) = f1; f1 = *(++in1) * pe1; } } } else { while (true) { if (f0 < f1) { *(out++) = f0; f0 = *(++in0) * pe0; if (in0 == factors_end) break; } else { *(out++) = f1; f1 = *(++in1) * pe1; } } } f[0] = f0; f[1] = f1; in[0] = in0; in[1] = in1; } else if (active_list_count == 3) { uint128 f0 = f[0], f1 = f[1], f2 = f[2]; uint128 *in0 = in[0], *in1 = in[1], *in2 = in[2]; const uint128 pe0 = pe[0], pe1 = pe[1], pe2 = pe[2]; while (true) { if (f0 < f1) { if (f0 < f2) { *(out++) = f0; f0 = *(++in0) * pe0; if (in0 == factors_end) break; } else { *(out++) = f2; f2 = *(++in2) * pe2; } } else { if (f1 < f2) { *(out++) = f1; f1 = *(++in1) * pe1; } else { *(out++) = f2; f2 = *(++in2) * pe2; } } } f[0] = f0; f[1] = f1, f[2] = f2; in[0] = in0; in[1] = in1, in[2] = in2; } else if (active_list_count >= 3) { while (true) { // Chose the smallest multiplier. uint k_min = 0; uint128 f_min = f[0]; for (uint k = 0; k < active_list_count; k++) if (f[k] < f_min) { f_min = f[k]; k_min = k; } // Write the output factor, advance the input pointer, and // produce a new factor in the array f[] of list heads. *(out++) = f_min; f[k_min] = *(++in[k_min]) * pe[k_min]; if (in[k_min] == factors_end) { assert(k_min == 0); break; } } } // --- Remove the newly emptied phantom input list. Note that this is // guaranteed *always* to be the first remaining non-empty list. assert(in[0] == factors_end); for (uint j = 1; j < active_list_count; j++) { in[j-1] = in[j]; pe[j-1] = pe[j]; f[j-1] = f[j]; } active_list_count -= 1; } assert(out == factors_end); } // --- Validate array of sorted factors. #ifndef NDEBUG { for (uint k = 0; k < d; k++) { if (factors[k] == 0) fatal_error("Produced a factor of 0 at index %u.", k); if (n % factors[k] != 0) fatal_error("Produced non-factor %s at index %u.", uint128_to_static_string(factors[k], 0), k); if ((k > 0) && (factors[k-1] == factors[k])) fatal_error("Duplicate factor %s at index %u.", uint128_to_static_string(factors[k], 0), k); if ((k > 0) && (factors[k-1] > factors[k])) fatal_error("Out-of-order factors %s and %s at indexes %u and %u.", uint128_to_static_string(factors[k-1], 0), uint128_to_static_string(factors[k], 1), k-1, k); } } #endif return factors; } //------------------------------------------------------------------------------ // Print prime-power factorization of a number. void print_ppf(const PrimePowerFactorization ppf) { printf("%s = ", uint128_to_static_string(ppf.n, 0)); if (ppf.n == 0) { printf("0"); } else { for (uint i = 0; i < ppf.ω; i++) { if (i > 0) printf(" x "); printf("%s", uint128_to_static_string(ppf.f[i].p, 0)); if (ppf.f[i].e > 1) printf("^%"PRIu8"", ppf.f[i].e); } } printf("/n"); } //------------------------------------------------------------------------------ int compare_powers_ascending(const void *const pf1, const void *const pf2) { const PrimePower f1 = *((const PrimePower *)pf1); const PrimePower f2 = *((const PrimePower *)pf2); return (f1.e < f2.e)? -1: (f1.e > f2.e)? +1: 0; // Not an error; duplicate exponents are common. } //------------------------------------------------------------------------------ int compare_powers_descending(const void *const pf1, const void *const pf2) { const PrimePower f1 = *((const PrimePower *)pf1); const PrimePower f2 = *((const PrimePower *)pf2); return (f1.e < f2.e)? +1: (f1.e > f2.e)? -1: 0; // Not an error; duplicate exponents are common. } //------------------------------------------------------------------------------ int compare_primes_ascending(const void *const pf1, const void *const pf2) { const PrimePower f1 = *((const PrimePower *)pf1); const PrimePower f2 = *((const PrimePower *)pf2); return (f1.p < f2.p)? -1: (f1.p > f2.p)? +1: 0; // Error; duplicate primes must never occur. } //------------------------------------------------------------------------------ int compare_primes_descending(const void *const pf1, const void *const pf2) { const PrimePower f1 = *((const PrimePower *)pf1); const PrimePower f2 = *((const PrimePower *)pf2); return (f1.p < f2.p)? +1: (f1.p > f2.p)? -1: 0; // Error; duplicate primes must never occur. } //------------------------------------------------------------------------------ // Sort prime-power factorization. void sort_ppf(PrimePowerFactorization *const ppf, const bool primes_major, // Best false const bool primes_ascending, // Best false const bool powers_ascending) // Best false { int (*compare_primes)(const void *, const void *) = primes_ascending? compare_primes_ascending : compare_primes_descending; int (*compare_powers)(const void *, const void *) = powers_ascending? compare_powers_ascending : compare_powers_descending; if (primes_major) { mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_powers); mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_primes); } else { mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_primes); mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_powers); } } //------------------------------------------------------------------------------ // Compute prime-power factorization of a 128-bit value. Note that this // function is designed to be fast *only* for numbers with very simple // factorizations, e.g., those that produce large factor lists. Do not attempt // to factor large semiprimes with this function. (The author does know how to // factor large numbers efficiently; however, efficient factorization is beyond // the scope of this small test program.) PrimePowerFactorization compute_ppf(const uint128 n) { PrimePowerFactorization ppf; if (n == 0) { ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 }; } else if (n == 1) { ppf = (PrimePowerFactorization){ .f[0] = { .p = 1, .e = 1 }, .ω = 1, .Ω = 1, .d = 1, .n = 1 }; } else { ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = n }; uint128 m = n; for (uint128 p = 2; p * p <= m; p += 1 + (p > 2)) { if (m % p == 0) { assert(ppf.ω <= MAX_ω); ppf.f[ppf.ω].p = p; ppf.f[ppf.ω].e = 0; while (m % p == 0) { m /= p; ppf.f[ppf.ω].e += 1; } ppf.d *= (1 + ppf.f[ppf.ω].e); ppf.Ω += ppf.f[ppf.ω].e; ppf.ω += 1; } } if (m > 1) { assert(ppf.ω <= MAX_ω); ppf.f[ppf.ω].p = m; ppf.f[ppf.ω].e = 1; ppf.d *= 2; ppf.Ω += 1; ppf.ω += 1; } } return ppf; } //------------------------------------------------------------------------------ // Parse prime-power factorization from a list of ASCII-encoded base-10 strings. // The values are assumed to be 2-tuples (p,e) of prime p and exponent e. // Primes must not exceed 2^128 - 1 and must not be repeated. Exponents must // not exceed 2^8 - 1, but can of course be repeated. The constructed value // must not exceed 2^128 - 1. PrimePowerFactorization parse_ppf(const uint pairs, const char *const values[]) { assert(pairs <= MAX_ω); PrimePowerFactorization ppf; if (pairs == 0) { ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 }; } else { ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = 1 }; for (uint i = 0; i < pairs; i++) { ppf.f[i].p = uint128_from_string(values[(i*2)+0]); ppf.f[i].e = (uint8)strtoumax(values[(i*2)+1], NULL, 10); // Validate prime value. if (ppf.f[i].p < 2) // (Ideally this would actually do a primality test.) fatal_error("Factor %s is invalid.", uint128_to_static_string(ppf.f[i].p, 0)); // Accumulate count of unique prime factors. if (ppf.ω > UINT8_MAX - 1) fatal_error("Small-omega overflow at factor %s^%"PRIu8".", uint128_to_static_string(ppf.f[i].p, 0), ppf.f[i].e); ppf.ω += 1; // Accumulate count of total prime factors. if (ppf.Ω > UINT8_MAX - ppf.f[i].e) fatal_error("Big-omega wverflow at factor %s^%"PRIu8".", uint128_to_static_string(ppf.f[i].p, 0), ppf.f[i].e); ppf.Ω += ppf.f[i].e; // Accumulate total divisor count. if (ppf.d > UINT32_MAX / (1 + ppf.f[i].e)) fatal_error("Divisor count overflow at factor %s^%"PRIu8".", uint128_to_static_string(ppf.f[i].p, 0), ppf.f[i].e); ppf.d *= (1 + ppf.f[i].e); // Accumulate value. for (uint8 k = 1; k <= ppf.f[i].e; k++) { if (ppf.n > UINT128_MAX / ppf.f[i].p) fatal_error("Value overflow at factor %s.", uint128_to_static_string(ppf.f[i].p, 0)); ppf.n *= ppf.f[i].p; } } } return ppf; } //------------------------------------------------------------------------------ // Main control. Parse command line and produce list of factors. int main(const int argc, const char *const argv[]) { bool primes_major = false; bool primes_ascending = false; bool powers_ascending = false; PrimePowerFactorization ppf; // --- Parse prime-power sort specifier (if present). uint value_base = 1; uint value_count = (uint)argc - 1; if ((argc > 1) && (argv[1][0] == ''-'')) { static const struct { char *str; bool primes_major, primes_ascending, powers_ascending; } sort_options[] = { // Sorting criteria: // ---------------------------------------- { "ep", 0,0,0 }, // Exponents descending, primes descending { "Ep", 0,0,1 }, // Exponents ascending, primes descending { "eP", 0,1,0 }, // Exponents descending, primes ascending { "EP", 0,1,1 }, // Exponents ascending, primes ascending { "p", 1,0,0 }, // Primes descending (exponents irrelevant) { "P", 1,1,0 }, // Primes ascending (exponents irrelevant) }; bool valid = false; for (uint i = 0; i < ARRAY_CAPACITY(sort_options); i++) { if (strcmp(&argv[1][1], sort_options[i].str) == 0) { primes_major = sort_options[i].primes_major; primes_ascending = sort_options[i].primes_ascending; powers_ascending = sort_options[i].powers_ascending; valid = true; break; } } if (!valid) fatal_error("Bad sort specifier: /"%s/"", argv[1]); value_base += 1; value_count -= 1; } // --- Prime factorization from either a number or a raw prime factorization. if (value_count == 1) { uint128 n = uint128_from_string(argv[value_base]); ppf = compute_ppf(n); } else { if (value_count % 2 != 0) fatal_error("Odd number of arguments (%u) given.", value_count); uint pairs = value_count / 2; ppf = parse_ppf(pairs, &argv[value_base]); } // --- Sort prime factorization by either the default or the user-overridden // configuration. sort_ppf(&ppf, primes_major, primes_ascending, powers_ascending); print_ppf(ppf); // --- Run for (as close as possible to) a fixed amount of time, tallying the // elapsed CPU time. uint128 iterations = 0; double cpu_time = 0.0; const double cpu_time_limit = 0.10; uint128 memory_usage = 0; while (cpu_time < cpu_time_limit) { clock_t clock_start = clock(); uint128 *factors = compute_factors(ppf); clock_t clock_end = clock(); cpu_time += (double)(clock_end - clock_start) / (double)CLOCKS_PER_SEC; memory_usage = sizeof(*factors) * ppf.d; if (++iterations == 0) //1) { for (uint32 i = 0; i < ppf.d; i++) printf("%s/n", uint128_to_static_string(factors[i], 0)); } if (factors) free(factors); } // --- Print the average amount of CPU time required for each iteration. uint memory_scale = (memory_usage >= 1e9)? 9: (memory_usage >= 1e6)? 6: (memory_usage >= 1e3)? 3: 0; char *memory_units = (memory_scale == 9)? "GB": (memory_scale == 6)? "MB": (memory_scale == 3)? "KB": "B"; printf("%s %"PRIu32" factors %.6f ms %.3f ns/factor %.3f %s/n", uint128_to_static_string(ppf.n, 0), ppf.d, cpu_time/iterations * 1e3, cpu_time/iterations * 1e9 / (double)(ppf.d? ppf.d : 1), (double)memory_usage / pow(10, memory_scale), memory_units); return 0; }

Podemos fusionar flujos de múltiplos, producidos para que no haya duplicados en primer lugar.

Comenzando con la lista [1] , para cada factor primo único p , multiplicamos la lista por p iterativamente k veces (donde k es la multiplicidad de p ) para obtener k nuevas listas y fusionarlas con esa lista entrante.

En Haskell

ordfactors n = -- <----------------------<---<---<-----/ foldr g [1] -- g (19,1) (g (7,1) (g (3,2) (g (2,3) [1]))) . reverse -- [(19,1),(7,1),(3,2),(2,3)] ^ . primePowers -- [(2,3),(3,2),(7,1),(19,1)] | $ n -- 9576 --->--->--->----/ where g (p,k) xs = mergeAsTree . take (k+1) -- take 3 [a,b,c,d..] = [a,b,c] . iterate (map (p*)) -- iterate f x = [x, f x, f(f x),...] $ xs {- g (2,3) [1] = let a0 = [1] a1 = map (2*) a0 -- [2] a2 = map (2*) a1 -- [4] a3 = map (2*) a2 -- [8] in mergeAsTree [ a0, a1, a2, a3 ] -- xs2 = [1,2,4,8] g (3,2) xs2 = let b0 = xs2 -- [1,2,4,8] b1 = map (3*) b0 -- [3,6,12,24] b2 = map (3*) b1 -- [9,18,36,72] in mergeAsTree [ b0, b1, b2 ] -- xs3 = [1,2,3,4,6,8,9,12,...] g (7,1) xs3 = mergeAsTree [ xs3, map (7*) xs3 ] -- xs4 = [1,2,3,4,6,7,8,9,...] g (19,1) xs4 = mergeAsTree [ xs4, map (19*) xs4 ] -} mergeAsTree xs = foldi merge [] xs -- [a,b] --> merge a b -- [a,b,c,d] --> merge (merge a b) (merge c d) foldi f z [] = z foldi f z (x:xs) = f x (foldi f z (pairs f xs)) pairs f (x:y:t) = f x y : pairs f t pairs _ t = t

foldi organiza las merge binarias (que presuponen que no hay duplicados en las secuencias que se fusionan, para una operación simplificada) en forma de árbol para mayor velocidad; mientras que foldr funciona de manera lineal.

primePowers n | n > 1 = -- 9576 => [(2,3),(3,2),(7,1),(19,1)] map (/xs -> (head xs,length xs)) -- ^ . group -- [[2,2,2],[3,3],[7],[19]] | . factorize -- [2,2,2,3,3,7,19] | $ n -- 9576 --->--->--->---/

group es una función estándar que agrupa a los iguales adyacentes en una lista, y factorize es un algoritmo conocido de división de prueba que produce factores primos de un número en orden no decreciente:

factorize n | n > 1 = go n (2:[3,5..]) -- 9576 = 2*2*2*3*3*7*19 where -- produce prime factors only go n (d:t) | d*d > n = [n] | r == 0 = d : go q (d:t) | otherwise = go n t where (q,r) = quotRem n d factorize _ = []

(.) es un operador de composición funcional, encadenando la salida de su función de argumento a su derecha como entrada en la de su izquierda. Se puede leer en voz alta (y $ ) como "de".