Calcolo di phi (k) per 1 < k < N
https://stackoverflow.com/questions/1024640
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianDomanda
Dato un N grande, devo scorrere attraverso tutto il phi (k) in modo tale che 1 < k < N:
- la complessità temporale deve essere
O (N logN) - la complessità della memoria deve essere sub
O (N)(poiché i valori di N saranno circa 10 12 )
È possibile? In tal caso, come?
Soluzione
Questo può essere fatto con la complessità della memoria O (Sqrt (N)) e la complessità della CPU O (N * Log (Log (N))) con un setaccio a finestre ottimizzato di Eratostene, come implementato nell'esempio di codice seguente.
Poiché non è stato specificato alcun linguaggio e poiché non conosco Python, l'ho implementato in VB.net, tuttavia posso convertirlo in C # se ne hai bisogno.
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
Nota che su O (N * Log (Log (N))), questa routine considera in media ogni numero su O (Log (Log (N))) che è molto, molto più veloce della fattorizzazione N singola più veloce algoritmi situati in alcune delle risposte qui. In effetti, a N = 10 ^ 12 è 2400 volte più veloce!
Ho testato questa routine sul mio laptop Intel Core 2 da 2 Ghz e calcola oltre 3.000.000 di valori Phi () al secondo. A questa velocità, occorrerebbero circa 4 giorni per calcolare 10 ^ 12 valori. Ho anche testato la correttezza fino a 100.000.000 senza errori. Si basa su numeri interi a 64 bit, quindi qualsiasi cosa fino a 2 ^ 63 (10 ^ 19) dovrebbe essere accurata (anche se troppo lenta per chiunque).
Ho anche un Visual Studio WinForm (anche VB.net) per eseguirlo / testarlo, che posso fornire se lo desideri.
Fammi sapere se hai domande.
Come richiesto nei commenti, ho aggiunto di seguito una versione C # del codice. Tuttavia, poiché attualmente sono nel mezzo di alcuni altri progetti, non ho tempo per convertirli da solo, quindi ho usato uno dei siti di conversione online da VB a C # ( http://www.carlosag.net/tools/codetranslator/ ). Quindi, tieni presente che questo è stato convertito automaticamente e non ho ancora avuto il tempo di testarlo o controllarlo da solo.
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))));
}
}
Altri suggerimenti
Nessuno ha trovato un modo più veloce per calcolare phi (k) (aka, Funzione totient di Eulero ) che trovando prima i fattori primi di k. I migliori matematici del mondo hanno lanciato molti cicli CPU al problema dal 1977, poiché trovare un modo più rapido per risolvere questo problema creerebbe un punto debole nella rel <
Il calcolo di phi (k) deve essere fatto usando la scomposizione in fattori primi di k, che è l'unico modo sensato di farlo. Se hai bisogno di un aggiornamento su questo, wikipedia porta la formula .
Se ora devi calcolare tutti i divisori primi di ogni numero tra 1 e una N grande, morirai di vecchiaia prima di vedere qualsiasi risultato, quindi farei il contrario, cioè costruirò tutti i numeri sotto N , usando i loro possibili fattori primi, ovvero tutti i numeri primi inferiori o uguali a N.
Il tuo problema sarà quindi simile al calcolo di tutti i divisori di un numero, solo tu non sai qual è il numero massimo di volte in cui potresti trovare un certo numero primo nella fattorizzazione in anticipo. Ottimizzare un iteratore originariamente scritto da Tim Peters nell'elenco Python (qualcosa Ho scritto un blog su ...) per includere la funzione totient, una possibile implementazione in python che produce coppie k, phi (k) potrebbe essere la seguente:
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
Se hai bisogno di aiuto per calcolare tutti i fattori primi inferiori a N, ho anche ha scritto un blog al riguardo ... Tieni presente, sebbene il calcolo di tutti i numeri primi inferiori a 10 12 sia di per sé un'impresa davvero notevole. .
Questo proviene dal Project Euler 245? Ricordo quella domanda e mi sono arreso.
Il modo più veloce che ho trovato per il calcolo del totiente è stato quello di moltiplicare i fattori primi (p-1) insieme, dato che k non ha fattori ripetuti (che non è mai stato il caso se ricordo correttamente il problema).
Quindi, per calcolare i fattori, sarebbe probabilmente meglio usare gmpy.next_prime o pyecm (fattorizzazione con curva ellittica).
Potresti anche setacciare i fattori primi come suggerisce Jaime. Per numeri fino a 10 12 , il fattore primo massimo è inferiore a 1 milione, il che dovrebbe essere ragionevole.
Se memorizzi le fattorizzazioni, potresti velocizzare ulteriormente la tua funzione phi.
Per questo tipo di problemi sto usando un iteratore che ritorna per ogni numero intero m < N l'elenco dei numeri primi < sqrt ( N ) che divide m. Per implementare un simile iteratore sto usando un array A di lunghezza R dove R > sqrt ( N ). Ad ogni punto l'array A contiene un elenco di numeri primi che dividono numeri interi m .. m + R -1. Cioè A [ m % R ] contiene numeri primi che dividono m . Ogni p primo è esattamente in un elenco, ovvero in A [ m % R ] per il numero intero più piccolo in l'intervallo m .. m + R -1 che è divisibile per p . Quando si genera l'elemento successivo dell'iteratore, viene restituito semplicemente l'elenco in A [ m % R ]. Quindi l'elenco di numeri primi viene rimosso da A [ m % R ] e ogni primo p viene aggiunto a A [( m + p )% R].
Con un elenco di numeri primi < sqrt ( N ) che divide m è facile trovare la fattorizzazione di m , dato che esiste al massimo un numero primo maggiore di sqrt ( N ).
Questo metodo ha complessità O ( N log (log ( N ))) presupponendo che tutte le operazioni, incluse quelle dell'elenco, assumano O (1). Il requisito di memoria è O (sqrt ( N )).
Sfortunatamente c'è un sovraccarico costante qui, quindi stavo cercando un modo più elegante per generare i valori phi (n), ma quindi non ho avuto successo.
Ecco un efficiente generatore di Python. L'avvertenza è che non produce i risultati in ordine. Si basa su https://stackoverflow.com/a/10110008/412529 .
La complessità della memoria è O (log (N)) in quanto deve solo memorizzare un elenco di fattori primi per un singolo numero alla volta.
La complessità della CPU è appena superlineare, qualcosa come 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)
Penso che tu possa fare il contrario. Invece di fattorizzare un numero intero k per ottenere phi (k), puoi provare a generare tutti i numeri interi da 1 a N dai numeri primi e ottenere rapidamente phi (k). Ad esempio, se P n è il massimo primo che è inferiore a N, puoi generare tutti gli interi inferiori a N di
P 1 i 1 * P 2 i 2 * ... * P n i n dove ogni i j viene eseguito da 0 a [ log (N) / log (P j )] ([] è la funzione floor).
In questo modo, puoi ottenere phi (k) rapidamente senza costose fattorizzazioni prime e ancora scorrere tutte le k tra 1 e N (non in ordine ma penso che non ti interessi all'ordine).
Setaccia i totient su 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)))))))))
Questo fattorizza N = PQ, dove P & amp; Q sono primi.
Funziona abbastanza bene, in elisir o Erlang.
Puoi provare diversi generatori per la tua sequenza pseudo-casuale. x * x + 1 è comunemente usato.
Questa riga: defp f0 (x, n), do: rem ((x * x) + 1, n)
Altri possibili punti di miglioramento: funzioni gcd , rem e abs migliori o alternative
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
