Calcul de phi (k) pour 1
https://stackoverflow.com/questions/1024640
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
Etant donné un grand N, je dois parcourir tous les phi (k) tels que 1 < k < N:
- la complexité temporelle doit être
O (N logN) - la complexité de la mémoire doit être inférieure à
O (N)(puisque les valeurs de N seront d'environ 10 12 )
Est-ce possible? Si oui, comment?
La solution
Cela peut être fait avec la complexité de la mémoire O (Sqrt (N)) et la complexité de la CPU O (N * Log (Log (N))) avec un tamis d’Ératosthène optimisé, comme implémenté dans l’exemple de code ci-dessous.
Comme aucun langage n'a été spécifié et que je ne connais pas Python, je l'ai implémenté dans VB.net, mais je peux le convertir en C # si vous en avez besoin.
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
Notez qu’en O (N * Log (Log (N))), cette routine factorise chaque nombre en O (Log (Log (N))) en moyenne, ce qui est beaucoup, beaucoup plus rapide que la factorisation à N simple la plus rapide. algorithmes situés par certaines des réponses ici. En fait, à N = 10 ^ 12, 2400 fois est plus rapide!
J'ai testé cette routine sur mon ordinateur portable 2Ghz Intel Core 2 et elle calcule plus de 3 000 000 de valeurs Phi () par seconde. À cette vitesse, il vous faudrait environ 4 jours pour calculer 10 ^ 12 valeurs. Je l'ai également testé pour l'exactitude jusqu'à 100.000.000 sans aucune erreur. Il est basé sur des entiers 64 bits, donc tout ce qui va jusqu’à 2 ^ 63 (10 ^ 19) devrait être précis (bien que trop lent pour quiconque).
J'ai également un fichier Visual Studio WinForm (également VB.net) pour l'exécuter / le tester, que je peux vous fournir si vous le souhaitez.
Faites-moi savoir si vous avez des questions.
Comme demandé dans les commentaires, j'ai ajouté ci-dessous une version C # du code. Cependant, comme je suis actuellement au milieu de certains autres projets, je n’ai pas le temps de le convertir moi-même; j’ai donc utilisé l’un des sites de conversion de VB en C # en ligne ( http://www.carlosag.net/tools/codetranslator/ ). Sachez donc que cela a été converti automatiquement et que je n’ai pas encore eu le temps de le tester ou de le vérifier moi-même.
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))));
}
}
Autres conseils
Personne n'a trouvé de moyen plus rapide de calculer phi (k) ( Fonction totale d'Euler ) que d’abord en trouvant les facteurs premiers de k. Les meilleurs mathématiciens du monde ont jeté de nombreux cycles de processeurs sur le problème depuis 1977, car trouver un moyen plus rapide de résoudre ce problème créerait une faiblesse dans le Algorithme de clé publique RSA . (La clé publique et la clé privée dans RSA sont calculées sur la base de phi (n), où n est le produit de deux grands nombres premiers.)
Le calcul de phi (k) doit être effectué à l'aide de la factorisation première de k, qui est la seule façon sensée de le faire. Si vous avez besoin d'une mise à jour à ce sujet, wikipedia contient la formule .
Si vous devez maintenant calculer tous les diviseurs premiers de chaque nombre compris entre 1 et un grand N, vous mourrez de vieillesse avant de voir un résultat, alors je ferais l'inverse, c'est-à-dire que tous les nombres sont inférieurs à N , en utilisant leurs facteurs premiers possibles, c’est-à-dire tous les nombres premiers inférieurs ou égaux à N.
Votre problème ressemblera donc au calcul de tous les diviseurs d’un nombre, si ce n’est que vous ne savez pas quel est le nombre maximal de fois où vous pouvez trouver un certain nombre premier dans la factorisation. Modification d'un itérateur écrit à l'origine par Tim Peters sur la liste python (quelque chose de que j'ai blogué à propos de ...) pour inclure la fonction totient, une implémentation possible en python donnant k, phi (k) paires pourrait être la suivante:
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 vous avez besoin d'aide pour calculer tous les facteurs premiers inférieurs à N, j'ai aussi a blogué à ce sujet ... N'oubliez pas que calculer tous les nombres premiers inférieurs à 10 12 est en soi un exploit remarquable. .
S'agit-il du projet Euler 245? Je me souviens de cette question et j’y ai renoncé.
Le moyen le plus rapide que j'ai trouvé pour calculer le total était de multiplier les facteurs premiers (p-1), étant donné que k n'a pas de facteurs répétés (ce qui n'a jamais été le cas si je me souviens bien du problème).
Donc, pour calculer les facteurs, il serait probablement préférable d'utiliser gmpy.next_prime ou pyecm (factorisation de la courbe elliptique).
Vous pouvez également trier les facteurs premiers comme le suggère Jaime. Pour des nombres allant jusqu'à 10 12 , le facteur premier maximum est inférieur à 1 million, ce qui devrait être raisonnable.
Si vous mémorisez des factorisations, cela pourrait accélérer encore plus votre fonction phi.
Pour ces types de problèmes, j’utilise un itérateur qui renvoie pour chaque entier m < N la liste des nombres premiers < sqrt ( N ) qui divise m. Pour implémenter un tel itérateur, j'utilise un tableau A de longueur R où R > sqrt ( N ). À chaque point, le tableau A contient la liste des nombres premiers qui divisent les entiers m .. m + R -1. C'est à dire. A [ m % R ] contient des nombres premiers divisant m . Chaque prime p se trouve dans une seule liste, c'est-à-dire dans A [ m % R ] pour le plus petit entier de la plage m .. m + R -1 divisible par p . Lors de la génération de l'élément suivant de l'itérateur, la liste dans A [ m % R ] est renvoyée. Ensuite, la liste des nombres premiers est supprimée de A [ m % R ] et chaque nombre premier p est ajouté à A . [( m + p )% R].
Avec une liste de nombres premiers < sqrt ( N ) divisant m , il est facile de trouver la factorisation de m , car il existe au plus un nombre premier plus grand que sqrt ( N ).
Cette méthode a une complexité O ( N log (log ( N ))) en supposant que toutes les opérations, y compris les opérations de liste, prennent O (1). La mémoire requise est O (sqrt ( N )).
Malheureusement, il y a des frais généraux constants ici, c'est pourquoi je cherchais un moyen plus élégant de générer les valeurs phi (n), mais je n'ai pas réussi.
Voici un générateur de python efficace. La mise en garde est que cela ne donne pas les résultats dans l'ordre. Il est basé sur https://stackoverflow.com/a/10110008/412529 .
La complexité de la mémoire est O (log (N)) car il ne doit stocker qu'une liste de facteurs premiers pour un nombre à la fois.
La complexité du processeur est à peine superlinéaire, quelque chose comme 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)
Je pense que vous pouvez faire l’inverse. Au lieu de factoriser un entier k pour obtenir phi (k), vous pouvez essayer de générer tous les entiers de 1 à N à partir de nombres premiers et obtenir rapidement phi (k). Par exemple, si P n est l’amorce maximale inférieure à N, vous pouvez générer tous les entiers inférieurs à N en
.P 1 i 1 * P 2 i 2 * ... * P n i n où chaque i j va de 0 à [ log (N) / log (P j )] ([] est la fonction de plancher).
De cette façon, vous pouvez obtenir rapidement phi (k) sans une factorisation première coûteuse et continuer à parcourir tous les k compris entre 1 et N (pas dans l’ordre mais je pense que l’ordre vous importe peu).
Tamiser les totients sur 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)))))))))
Ceci factorise N = PQ, où P & amp; Q sont premiers.
Fonctionne très bien, à Elixir ou à Erlang.
Vous pouvez essayer différents générateurs pour votre séquence pseudo-aléatoire. x * x + 1 est couramment utilisé.
Cette ligne: defp f0 (x, n), faites: rem ((x * x) + 1, n)
Autres points d'amélioration possibles: fonctions gcd , rem et abs abs
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
