Вычисление phi(k) для 1
https://stackoverflow.com/questions/1024640
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Учитывая большое N, мне нужно выполнить итерацию по всем phi (k) таким образом, чтобы 1 < k < N :
- временная сложность должна быть
O(N logN) - сложность памяти должна быть ниже
O(N)(поскольку значения N будут составлять около 1012)
Возможно ли это?Если да, то каким образом?
Решение
Это может быть сделано со сложностью памяти O (Sqrt (N)) и сложностью процессора O (N * Log(Log(N))) с оптимизированным оконным решетом Eratosthenes, как реализовано в примере кода ниже.
Поскольку язык не был указан, и поскольку я не знаю Python, я реализовал его в VB.net , однако я могу преобразовать его в C #, если вам это нужно.
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
Обратите внимание, что при O (N * Log(Log(N))) эта процедура в среднем разлагает каждое число на O(Log(Log(N))), что намного, намного быстрее, чем самые быстрые алгоритмы факторизации с одним N, приведенные в некоторых ответах здесь.Фактически, при N = 10 ^ 12 это 2400 в разы быстрее!
Я протестировал эту процедуру на своем ноутбуке Intel Core 2 с частотой 2 ГГц, и она вычисляет более 3 000 000 значений Phi () в секунду.При такой скорости вам потребуется около 4 дней, чтобы вычислить 10 ^ 12 значений.Я также протестировал его на корректность до 100 000 000 без каких-либо ошибок.Он основан на 64-битных целых числах, поэтому все, что не превышает 2 ^ 63 (10 ^ 19), должно быть точным (хотя и слишком медленным для кого-либо).
У меня также есть Visual Studio WinForm (также VB.net) для ее запуска / тестирования, которую я могу предоставить, если вы этого хотите.
Дайте мне знать, если у вас возникнут какие-либо вопросы.
Как и просили в комментариях, я добавил ниже версию кода на C #.Однако, поскольку в настоящее время я занят некоторыми другими проектами, у меня нет времени конвертировать его самостоятельно, поэтому я воспользовался одним из онлайн-сайтов по преобразованию VB в C # (http://www.carlosag.net/tools/codetranslator/).Так что имейте в виду, что это было преобразовано автоматически, и у меня еще не было времени протестировать или проверить это самостоятельно.
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))));
}
}
Другие советы
Никто не нашел более быстрого способа вычислить phi (k) (он же, Постоянная функция Эйлера), чем путем первого нахождения простых множителей k.Лучшие математики мира потратили много циклов процессора на решение этой задачи с 1977 года, поскольку поиск более быстрого способа решения этой проблемы привел бы к ослаблению Алгоритм с открытым ключом RSA.(Как открытый, так и закрытый ключ в RSA вычисляются на основе phi(n), где n - произведение двух больших простых чисел.)
Вычисление phi (k) должно быть выполнено с использованием простой факторизации k, что является единственным разумным способом сделать это.Если вам нужно уточнить этот вопрос, википедия содержит формулу.
Если вам теперь нужно вычислить все простые делители каждого числа от 1 до большого N, вы умрете от старости, прежде чем увидите какой-либо результат, поэтому я бы поступил наоборот, т. Е.постройте все числа ниже N, используя их возможные простые множители, т. е.все простые числа меньше или равны N.
Таким образом, ваша проблема будет похожа на вычисление всех делителей числа, только вы не знаете, какое максимальное количество раз вы можете заранее найти определенное простое число при факторизации.Настройка итератора, первоначально написанного Тимом Питерсом в списке python (что-то Я писал в блоге о...) чтобы включить функцию totient, возможная реализация на python, которая выдает k пар phi (k), могла бы быть следующей:
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
Если вам нужна помощь в вычислении всех простых множителей ниже N, я также написал об этом в блоге...Имейте в виду, однако, что вычисление всех простых чисел ниже 1012 это само по себе довольно примечательный подвиг...
Это из проекта Эйлер 245?Я помню этот вопрос, и я отказался от него.
Самый быстрый способ, который я нашел для вычисления totient, состоял в том, чтобы умножить простые множители (p-1) вместе, учитывая, что k не имеет повторяющихся множителей (чего никогда не было, если я правильно помню проблему).
Таким образом, для вычисления коэффициентов, вероятно, было бы лучше использовать gmpy.next_prime - следующий прайм или пьекм (факторизация эллиптической кривой).
Вы также могли бы просеять основные факторы, как предлагает Джейми.Для чисел до 1012, максимальный простой коэффициент ниже 1 миллиона, что должно быть разумным.
Если вы запомните факторизацию, это может еще больше ускорить вашу функцию phi.
Для такого рода задач я использую итератор, который возвращает для каждого целого числа m < N список простых чисел < кврт(N), которые делят m.Для реализации такого итератора я использую массив A длины R где R > sqrt(N).В каждой точке массива A содержит список простых чисел, которые делят целые числа m..m+R-1.То есть. A[m % R] содержит простые числа , делящиеся m.Каждое простое p находится ровно в одном списке, т. е.в A[m % R] для наименьшего целого числа в диапазоне m .. m+R-1 , которое делится на p.При генерации следующего элемента итератора просто список в A[m % R] возвращается.Затем список простых чисел удаляется из A[m % R] и каждое простое число p добавляется к A[(m+p) % R].
Со списком простых чисел < кврт(N) разделяющий m легко найти факторизацию m, поскольку существует не более одного простого числа , большего , чем sqrt(N).
Этот метод имеет сложность O(N журнал(log(N))) в предположении, что все операции, включая операции со списком, занимают O(1).Потребность в памяти равна O(sqrt(N)).
К сожалению, здесь есть некоторые постоянные накладные расходы, поэтому я искал более элегантный способ генерировать значения phi (n), но так и не добился успеха.
Вот эффективный генератор python.Предостережение заключается в том, что это не дает результатов по порядку.Она основана на https://stackoverflow.com/a/10110008/412529 .
Сложность памяти равна O (log (N)), поскольку в ней одновременно должен храниться только список простых множителей для одного числа.
Сложность процессора едва ли сверхлинейна, что-то вроде 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)
Я думаю, вы можете пойти другим путем.Вместо того чтобы разлагать целое число k на множители, чтобы получить phi (k), вы можете попытаться сгенерировать все целые числа от 1 до N из простых чисел и быстро получить phi (k).Например, если Pn является максимальным простым числом, которое меньше N, вы можете сгенерировать все целые числа меньше N с помощью
P1 i 1 * P2 i 2 * ...* Pn i n где каждый яj выполнить от 0 до [log (N) / log (Pj)] ([] - функция floor).
Таким образом, вы можете быстро получить phi (k) без дорогостоящей простой факторизации и по-прежнему выполнять итерации по всем k от 1 до N (не по порядку, но я думаю, что порядок вас не волнует).
Просейте ингредиенты, чтобы 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)))))))))
Это разлагает N = PQ, где P & Q - простые числа.
Работает довольно хорошо в Elixir или Erlang.
Вы можете попробовать различные генераторы для вашей псевдослучайной последовательности. x*x + 1 широко используется.
Эта линия: defp f0(x, n), do: rem((x * x) + 1, n)
Другие возможные моменты улучшения:лучшее или альтернативное нод, рем и абс функции
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
