Cálculo de phi (k) para 1 < k < N
https://stackoverflow.com/questions/1024640
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPregunta
Dado un N grande, necesito iterar a través de todo phi (k) de modo que 1 < k < N:
- la complejidad del tiempo debe ser
O (N logN) - la complejidad de la memoria debe ser sub
O (N)(ya que los valores de N serán alrededor de 10 12 )
¿Es posible? Si es así, ¿cómo?
Solución
Esto se puede hacer con la complejidad de la memoria O (Sqrt (N)) y la complejidad de la CPU O (N * Log (Log (N))) con un tamiz de Eratóstenes con ventana optimizado, como se implementa en el siguiente ejemplo de código.
Como no se especificó ningún idioma y como no conozco Python, lo he implementado en VB.net, sin embargo, puedo convertirlo a C # si lo necesita.
Imports System.Math
Public Class TotientSerialCalculator
'Implements an extremely efficient Serial Totient(phi) calculator '
' This implements an optimized windowed Sieve of Eratosthenes. The'
' window size is set at Sqrt(N) both to optimize collecting and '
' applying all of the Primes below Sqrt(N), and to minimize '
' window-turning overhead. '
' '
' CPU complexity is O( N * Log(Log(N)) ), which is virtually linear.'
' '
' MEM Complexity is O( Sqrt(N) ). '
' '
' This is probalby the ideal combination, as any attempt to further '
'reduce memory will almost certainly result in disproportionate increases'
'in CPU complexity, and vice-versa. '
Structure NumberFactors
Dim UnFactored As Long 'the part of the number that still needs to be factored'
Dim Phi As Long 'the totient value progressively calculated'
' (equals total numbers less than N that are CoPrime to N)'
'MEM = 8 bytes each'
End Structure
Private ReportInterval As Long
Private PrevLast As Long 'the last value in the previous window'
Private FirstValue As Long 'the first value in this windows range'
Private WindowSize As Long
Private LastValue As Long 'the last value in this windows range'
Private NextFirst As Long 'the first value in the next window'
'Array that stores all of the NumberFactors in the current window.'
' this is the primary memory consumption for the class and it'
' is 16 * Sqrt(N) Bytes, which is O(Sqrt(N)).'
Public Numbers() As NumberFactors
' For N=10^12 (1 trilion), this will be 16MB, which should be bearable anywhere.'
'(note that the Primes() array is a secondary memory consumer'
' at O(pi(Sqrt(N)), which will be within 10x of O(Sqrt(N)))'
Public Event EmitTotientPair(ByVal k As Long, ByVal Phi As Long)
'===== The Routine To Call: ========================'
Public Sub EmitTotientPairsToN(ByVal N As Long)
'Routine to Emit Totient pairs {k, Phi(k)} for k = 1 to N'
' 2009-07-14, RBarryYoung, Created.'
Dim i As Long
Dim k As Long 'the current number being factored'
Dim p As Long 'the current prime factor'
'Establish the Window frame:'
' note: WindowSize is the critical value that controls both memory'
' usage and CPU consumption and must be SQRT(N) for it to work optimally.'
WindowSize = Ceiling(Sqrt(CDbl(N)))
ReDim Numbers(0 To WindowSize - 1)
'Initialize the first window:'
MapWindow(1)
Dim IsFirstWindow As Boolean = True
'adjust this to control how often results are show'
ReportInterval = N / 100
'Allocate the primes array to hold the primes list:'
' Only primes <= SQRT(N) are needed for factoring'
' PiMax(X) is a Max estimate of the number of primes <= X'
Dim Primes() As Long, PrimeIndex As Long, NextPrime As Long
'init the primes list and its pointers'
ReDim Primes(0 To PiMax(WindowSize) - 1)
Primes(0) = 2 '"prime" the primes list with the first prime'
NextPrime = 1
'Map (and Remap) the window with Sqrt(N) numbers, Sqrt(N) times to'
' sequentially map all of the numbers <= N.'
Do
'Sieve the primes across the current window'
PrimeIndex = 0
'note: cant use enumerator for the loop below because NextPrime'
' changes during the first window as new primes <= SQRT(N) are accumulated'
Do While PrimeIndex < NextPrime
'get the next prime in the list'
p = Primes(PrimeIndex)
'find the first multiple of (p) in the current window range'
k = PrevLast + p - (PrevLast Mod p)
Do
With Numbers(k - FirstValue)
.UnFactored = .UnFactored \ p 'always works the first time'
.Phi = .Phi * (p - 1) 'Phi = PRODUCT( (Pi-1)*Pi^(Ei-1) )'
'The loop test that follows is probably the central CPU overhead'
' I believe that it is O(N*Log(Log(N)), which is virtually O(N)'
' ( for instance at N = 10^12, Log(Log(N)) = 3.3 )'
Do While (.UnFactored Mod p) = 0
.UnFactored = .UnFactored \ p
.Phi = .Phi * p
Loop
End With
'skip ahead to the next multiple of p: '
'(this is what makes it so fast, never have to try prime factors that dont apply)'
k += p
'repeat until we step out of the current window:'
Loop While k < NextFirst
'if this is the first window, then scan ahead for primes'
If IsFirstWindow Then
For i = Primes(NextPrime - 1) + 1 To p ^ 2 - 1 'the range of possible new primes'
'Dont go beyond the first window'
If i >= WindowSize Then Exit For
If Numbers(i - FirstValue).UnFactored = i Then
'this is a prime less than SQRT(N), so add it to the list.'
Primes(NextPrime) = i
NextPrime += 1
End If
Next
End If
PrimeIndex += 1 'move to the next prime'
Loop
'Now Finish & Emit each one'
For k = FirstValue To LastValue
With Numbers(k - FirstValue)
'Primes larger than Sqrt(N) will not be finished: '
If .UnFactored > 1 Then
'Not done factoring, must be an large prime factor remaining: '
.Phi = .Phi * (.UnFactored - 1)
.UnFactored = 1
End If
'Emit the value pair: (k, Phi(k)) '
EmitPhi(k, .Phi)
End With
Next
're-Map to the next window '
IsFirstWindow = False
MapWindow(NextFirst)
Loop While FirstValue <= N
End Sub
Sub EmitPhi(ByVal k As Long, ByVal Phi As Long)
'just a placeholder for now, that raises an event to the display form'
' periodically for reporting purposes. Change this to do the actual'
' emitting.'
If (k Mod ReportInterval) = 0 Then
RaiseEvent EmitTotientPair(k, Phi)
End If
End Sub
Public Sub MapWindow(ByVal FirstVal As Long)
'Efficiently reset the window so that we do not have to re-allocate it.'
'init all of the boundary values'
FirstValue = FirstVal
PrevLast = FirstValue - 1
NextFirst = FirstValue + WindowSize
LastValue = NextFirst - 1
'Initialize the Numbers prime factor arrays'
Dim i As Long
For i = 0 To WindowSize - 1
With Numbers(i)
.UnFactored = i + FirstValue 'initially equal to the number itself'
.Phi = 1 'starts at mulplicative identity(1)'
End With
Next
End Sub
Function PiMax(ByVal x As Long) As Long
'estimate of pi(n) == {primes <= (n)} that is never less'
' than the actual number of primes. (from P. Dusart, 1999)'
Return (x / Log(x)) * (1.0 + 1.2762 / Log(x))
End Function
End Class
Tenga en cuenta que en O (N * Log (Log (N))), esta rutina está factorizando cada número en O (Log (Log (N))) en promedio, que es mucho, mucho más rápido que la factorización N más rápida algoritmos ubicados por algunas de las respuestas aquí. De hecho, en N = 10 ^ 12 es 2400 veces más rápido
He probado esta rutina en mi laptop Intel Core 2 de 2Ghz y calcula más de 3,000,000 de valores Phi () por segundo. A esta velocidad, le tomaría aproximadamente 4 días calcular 10 ^ 12 valores. También lo probé para corregir hasta 100,000,000 sin ningún error. Se basa en enteros de 64 bits, por lo que cualquier cosa hasta 2 ^ 63 (10 ^ 19) debe ser precisa (aunque demasiado lenta para cualquier persona).
También tengo un Visual Studio WinForm (también VB.net) para ejecutarlo / probarlo, que puedo proporcionar si lo desea.
Avíseme si tiene alguna pregunta.
Como se solicitó en los comentarios, he agregado a continuación una versión C # del código. Sin embargo, debido a que actualmente estoy en medio de otros proyectos, no tengo tiempo para convertirlo yo mismo, así que he usado uno de los sitios de conversión de VB a C # en línea ( http://www.carlosag.net/tools/codetranslator/ ). Tenga en cuenta que esto se convirtió automáticamente y todavía no he tenido tiempo de probarlo o comprobarlo.
using System.Math;
public class TotientSerialCalculator {
// Implements an extremely efficient Serial Totient(phi) calculator '
// This implements an optimized windowed Sieve of Eratosthenes. The'
// window size is set at Sqrt(N) both to optimize collecting and '
// applying all of the Primes below Sqrt(N), and to minimize '
// window-turning overhead. '
// '
// CPU complexity is O( N * Log(Log(N)) ), which is virtually linear.'
// '
// MEM Complexity is O( Sqrt(N) ). '
// '
// This is probalby the ideal combination, as any attempt to further '
// reduce memory will almost certainly result in disproportionate increases'
// in CPU complexity, and vice-versa. '
struct NumberFactors {
private long UnFactored; // the part of the number that still needs to be factored'
private long Phi;
}
private long ReportInterval;
private long PrevLast; // the last value in the previous window'
private long FirstValue; // the first value in this windows range'
private long WindowSize;
private long LastValue; // the last value in this windows range'
private long NextFirst; // the first value in the next window'
// Array that stores all of the NumberFactors in the current window.'
// this is the primary memory consumption for the class and it'
// is 16 * Sqrt(N) Bytes, which is O(Sqrt(N)).'
public NumberFactors[] Numbers;
// For N=10^12 (1 trilion), this will be 16MB, which should be bearable anywhere.'
// (note that the Primes() array is a secondary memory consumer'
// at O(pi(Sqrt(N)), which will be within 10x of O(Sqrt(N)))'
//NOTE: this part looks like it did not convert correctly
public event EventHandler EmitTotientPair;
private long k;
private long Phi;
// ===== The Routine To Call: ========================'
public void EmitTotientPairsToN(long N) {
// Routine to Emit Totient pairs {k, Phi(k)} for k = 1 to N'
// 2009-07-14, RBarryYoung, Created.'
long i;
long k;
// the current number being factored'
long p;
// the current prime factor'
// Establish the Window frame:'
// note: WindowSize is the critical value that controls both memory'
// usage and CPU consumption and must be SQRT(N) for it to work optimally.'
WindowSize = Ceiling(Sqrt(double.Parse(N)));
object Numbers;
this.MapWindow(1);
bool IsFirstWindow = true;
ReportInterval = (N / 100);
// Allocate the primes array to hold the primes list:'
// Only primes <= SQRT(N) are needed for factoring'
// PiMax(X) is a Max estimate of the number of primes <= X'
long[] Primes;
long PrimeIndex;
long NextPrime;
// init the primes list and its pointers'
object Primes;
-1;
Primes[0] = 2;
// "prime" the primes list with the first prime'
NextPrime = 1;
// Map (and Remap) the window with Sqrt(N) numbers, Sqrt(N) times to'
// sequentially map all of the numbers <= N.'
for (
; (FirstValue <= N);
) {
PrimeIndex = 0;
// note: cant use enumerator for the loop below because NextPrime'
// changes during the first window as new primes <= SQRT(N) are accumulated'
while ((PrimeIndex < NextPrime)) {
// get the next prime in the list'
p = Primes[PrimeIndex];
// find the first multiple of (p) in the current window range'
k = (PrevLast
+ (p
- (PrevLast % p)));
for (
; (k < NextFirst);
) {
// With...
UnFactored;
p;
// always works the first time'
(Phi
* (p - 1));
while (// TODO: Warning!!!! NULL EXPRESSION DETECTED...
) {
(UnFactored % p);
UnFactored;
(Phi * p);
}
// skip ahead to the next multiple of p: '
// (this is what makes it so fast, never have to try prime factors that dont apply)'
k = (k + p);
// repeat until we step out of the current window:'
}
// if this is the first window, then scan ahead for primes'
if (IsFirstWindow) {
for (i = (Primes[(NextPrime - 1)] + 1); (i
<= (p | (2 - 1))); i++) {
// the range of possible new primes'
// TODO: Warning!!! The operator should be an XOR ^ instead of an OR, but not available in CodeDOM
// Dont go beyond the first window'
if ((i >= WindowSize)) {
break;
}
if ((Numbers[(i - FirstValue)].UnFactored == i)) {
// this is a prime less than SQRT(N), so add it to the list.'
Primes[NextPrime] = i;
NextPrime++;
}
}
}
PrimeIndex++;
// move to the next prime'
}
// Now Finish & Emit each one'
for (k = FirstValue; (k <= LastValue); k++) {
// With...
// Primes larger than Sqrt(N) will not be finished: '
if ((Numbers[(k - FirstValue)].UnFactored > 1)) {
// Not done factoring, must be an large prime factor remaining: '
(Numbers[(k - FirstValue)].Phi * (Numbers[(k - FirstValue)].UnFactored - 1).UnFactored) = 1;
Numbers[(k - FirstValue)].Phi = 1;
}
// Emit the value pair: (k, Phi(k)) '
this.EmitPhi(k, Numbers[(k - FirstValue)].Phi);
}
// re-Map to the next window '
IsFirstWindow = false;
this.MapWindow(NextFirst);
}
}
void EmitPhi(long k, long Phi) {
// just a placeholder for now, that raises an event to the display form'
// periodically for reporting purposes. Change this to do the actual'
// emitting.'
if (((k % ReportInterval)
== 0)) {
EmitTotientPair(k, Phi);
}
}
public void MapWindow(long FirstVal) {
// Efficiently reset the window so that we do not have to re-allocate it.'
// init all of the boundary values'
FirstValue = FirstVal;
PrevLast = (FirstValue - 1);
NextFirst = (FirstValue + WindowSize);
LastValue = (NextFirst - 1);
// Initialize the Numbers prime factor arrays'
long i;
for (i = 0; (i
<= (WindowSize - 1)); i++) {
// With...
// initially equal to the number itself'
Phi = 1;
// starts at mulplicative identity(1)'
}
}
long PiMax(long x) {
// estimate of pi(n) == {primes <= (n)} that is never less'
// than the actual number of primes. (from P. Dusart, 1999)'
return ((x / Log(x)) * (1 + (1.2762 / Log(x))));
}
}
Otros consejos
Nadie ha encontrado una forma más rápida de calcular phi (k) (también conocido como Función totient de Euler ) que al encontrar primero los factores primos de k. Los mejores matemáticos del mundo han lanzado muchos ciclos de CPU al problema desde 1977, ya que encontrar una forma más rápida de resolver este problema crearía una debilidad en el Algoritmo de clave pública RSA . (Tanto la clave pública como la privada en RSA se calculan en función de phi (n), donde n es el producto de dos números primos grandes).
El cálculo de phi (k) debe hacerse utilizando la factorización prima de k, que es la única forma sensata de hacerlo. Si necesita una actualización sobre eso, wikipedia lleva la fórmula .
Si ahora tiene que calcular todos los divisores primos de cada número entre 1 y una N grande, morirá de vejez antes de ver algún resultado, por lo que iría al revés, es decir, construiría todos los números debajo de N , utilizando sus posibles factores primos, es decir, todos los primos menores o iguales a N.
Por lo tanto, su problema será similar al cálculo de todos los divisores de un número, solo que no sabe cuál es el número máximo de veces que puede encontrar un cierto primo en la factorización de antemano. Ajustar un iterador originalmente escrito por Tim Peters en la lista de Python (algo He publicado un blog sobre ...) para incluir la función totient, una posible implementación en python que produce pares k, phi (k) podría ser la siguiente:
def composites(factors, N) :
"""
Generates all number-totient pairs below N, unordered, from the prime factors.
"""
ps = sorted(set(factors))
omega = len(ps)
def rec_gen(n = 0) :
if n == omega :
yield (1,1)
else :
pows = [(1,1)]
val = ps[n]
while val <= N :
pows += [(val, val - pows[-1][0])]
val *= ps[n]
for q, phi_q in rec_gen(n + 1) :
for p, phi_p in pows :
if p * q > N :
break
else :
yield p * q, phi_p * phi_q
for p in rec_gen() :
yield p
Si necesita ayuda para calcular todos los factores primos por debajo de N, también tengo escribió en un blog al respecto ... Tenga en cuenta, aunque calcular todos los números primos por debajo de 10 12 es en sí mismo una hazaña notable. .
¿Es esto del Proyecto Euler 245? Recuerdo esa pregunta y me di por vencida.
La forma más rápida que encontré para calcular totient fue multiplicar los factores primos (p-1) juntos, dado que k no tiene factores repetidos (lo cual nunca fue el caso si recuerdo el problema correctamente).
Entonces, para calcular factores, probablemente sería mejor usar gmpy.next_prime o pyecm (factorización de curva elíptica).
También podría tamizar los factores primos como sugiere Jaime. Para números de hasta 10 12 , el factor primo máximo es inferior a 1 millón, lo que debería ser razonable.
Si memoriza factorizaciones, podría acelerar aún más su función phi.
Para este tipo de problemas, estoy usando un iterador que devuelve para cada entero m < N la lista de primos < sqrt ( N ) que divide m. Para implementar dicho iterador, estoy usando una matriz A de longitud R donde R > sqrt ( N ). En cada punto, la matriz A contiene una lista de números primos que dividen los enteros m .. m + R -1. Es decir. A [ m % R ] contiene números primos que dividen m . Cada primo p está exactamente en una lista, es decir, en A [ m % R ] para el entero más pequeño en el rango m .. m + R -1 que es divisible por p . Al generar el siguiente elemento del iterador, simplemente se devuelve la lista en A [ m % R ]. Luego, la lista de primos se elimina de A [ m % R ] y cada primo p se agrega a A [( m + p )% R].
Con una lista de primos < sqrt ( N ) dividiendo m es fácil encontrar la factorización de m , ya que hay como máximo un primo mayor que sqrt ( N ).
Este método tiene complejidad O ( N log (log ( N ))) bajo el supuesto de que todas las operaciones, incluidas las operaciones de lista, toman O (1). El requisito de memoria es O (sqrt ( N )).
Desafortunadamente, hay una sobrecarga constante aquí, por lo tanto, estaba buscando una forma más elegante de generar los valores phi (n), pero no he tenido éxito.
Aquí hay un generador de Python eficiente. La advertencia es que no produce los resultados en orden. Se basa en https://stackoverflow.com/a/10110008/412529 .
La complejidad de la memoria es O (log (N)) ya que solo tiene que almacenar una lista de factores primos para un solo número a la vez.
La complejidad de la CPU es apenas superlineal, algo así como O (N log log N).
def totientsbelow(N):
allprimes = primesbelow(N+1)
def rec(n, partialtot=1, min_p = 0):
for p in allprimes:
if p > n:
break
# avoid double solutions such as (6, [2,3]), and (6, [3,2])
if p < min_p: continue
yield (p, p-1, [p])
for t, tot2, r in rec(n//p, partialtot, min_p = p): # uses integer division
yield (t*p, tot2 * p if p == r[0] else tot2 * (p-1), [p] + r)
for n, t, factors in rec(N):
yield (n, t)
Creo que puedes ir al revés. En lugar de factorizar un entero k para obtener phi (k), puede intentar generar todos los enteros de 1 a N a partir de números primos y obtener phi (k) rápidamente. Por ejemplo, si P n es el primo máximo que es menor que N, puede generar todos los enteros menores que N por
P 1 i 1 * P 2 i 2 * ... * P n i n donde cada i j va de 0 a [ log (N) / log (P j )] ([] es la función de piso).
De esa manera, puede obtener phi (k) rápidamente sin una costosa factorización prima y aún recorrer todo k entre 1 y N (no en orden, pero creo que no le importa el orden).
Tamice los clientes para n :
(define (totients n)
(let ((tots (make-vector (+ n 1))))
(do ((i 0 (+ i 1))) ((< n i))
(vector-set! tots i i))
(do ((i 2 (+ i 1))) ((< n i) tots)
(when (= i (vector-ref tots i))
(vector-set! tots i (- i 1))
(do ((j (+ i i) (+ i j))) ((< n j))
(vector-set! tots j
(* (vector-ref tots j) (- 1 (/ i)))))))))
Esto factoriza N = PQ, donde P & amp; Q son primos.
Funciona bastante bien, en Elixir o Erlang.
Puede probar diferentes generadores para su secuencia pseudoaleatoria. x * x + 1 se usa comúnmente.
Esta línea: defp f0 (x, n), do: rem ((x * x) + 1, n)
Otros posibles puntos de mejora: mejores o alternativas gcd , rem y abs funciones
defmodule Factorizer do
def factorize(n) do
t = System.system_time
x = pollard(n, 2_000_000, 2_000_000)
y = div(n, x)
p = min(x, y)
q = max(x, y)
t = System.system_time - t
IO.puts "
Factorized #{n}: into [#{p} , #{q}] in #{t} μs
"
{p, q}
end
defp gcd(a,0), do: a
defp gcd(a,b), do: gcd(b,rem(a,b))
defp pollard(n, a, b) do
a = f0(a, n)
b = f0(f0(b, n), n)
p = gcd(abs(b - a), n)
case p > 1 do
true -> p
false -> pollard(n, a, b)
end
end
defp f0(x, n), do: rem((x * x) + 1, n)
end
