حساب phi(k) لـ 1<k<N
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))) باستخدام منخل إراتوستينس ذو النوافذ المُحسّنة، كما تم تنفيذه في مثال التعليمات البرمجية أدناه.
نظرًا لعدم تحديد لغة ولأنني لا أعرف لغة 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) (المعروف أيضًا باسم، دالة أويلر totient) من خلال إيجاد العوامل الأولية لـ k أولاً.لقد طرح أفضل علماء الرياضيات في العالم العديد من دورات وحدة المعالجة المركزية لحل هذه المشكلة منذ عام 1977، نظرًا لأن العثور على طريقة أسرع لحل هذه المشكلة من شأنه أن يخلق ضعفًا في حل هذه المشكلة. خوارزمية المفتاح العام RSA.(يتم حساب كل من المفتاح العام والخاص في RSA على أساس phi(n)، حيث n هو منتج عددين أوليين كبيرين.)
حساب phi(k) يجب أن يتم باستخدام التحليل الأولي لـ k، وهي الطريقة المعقولة الوحيدة للقيام بذلك.إذا كنت بحاجة إلى تجديد معلوماتك عن ذلك، ويكيبيديا تحمل الصيغة.
إذا كان عليك الآن حساب جميع المقسومات الأولية لكل رقم يقع بين 1 وN كبير، فسوف تموت بسبب الشيخوخة قبل أن ترى أي نتيجة، لذلك سأذهب في الاتجاه المعاكس، أي.بناء جميع الأعداد تحت N، وذلك باستخدام العوامل الأولية المحتملة، أي.جميع الأعداد الأولية أقل من أو تساوي N.
وبالتالي فإن مشكلتك ستكون مشابهة لحساب جميع مقسومات الرقم، إلا أنك لا تعرف ما هو الحد الأقصى لعدد المرات التي قد تجد فيها عددًا أوليًا معينًا في التحليل مسبقًا.التغيير والتبديل في مكرر كتبه في الأصل تيم بيترز في قائمة بايثون (شيء ما لقد قمت بالتدوين عن...) لتضمين دالة 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؟أتذكر هذا السؤال، وقد تخليت عنه.
أسرع طريقة وجدتها لحساب المجموع هي ضرب العوامل الأولية (p-1) معًا، نظرًا لأن k لا تحتوي على عوامل متكررة (وهو ما لم يحدث أبدًا إذا كنت أتذكر المشكلة بشكل صحيح).
لذا، بالنسبة لحساب العوامل، ربما يكون من الأفضل استخدامها gmpy.next_prime أو pyecm (تحليل المنحنى الإهليلجي).
يمكنك أيضًا غربلة العوامل الأولية كما يقترح خايمي.للأعداد حتى 1012, ، الحد الأقصى للعامل الأولي أقل من مليون والذي يجب أن يكون معقولاً.
إذا حفظت العوامل، فقد يؤدي ذلك إلى تسريع وظيفة phi بشكل أكبر.
بالنسبة لهذا النوع من المشكلات، أستخدم مكررًا يقوم بإرجاع كل عدد صحيح م <ن قائمة الأعداد الأولية <sqrt(ن) التي تقسم م.لتنفيذ مثل هذا المكرر أستخدم مصفوفة أ من الطول ر أين ر > سرت (ن).في كل نقطة المصفوفة أ يحتوي على قائمة الأعداد الأولية التي تقسم الأعداد الصحيحة م ..م+ر-1.أي. أ[م % ر] يحتوي على تقسيم الأعداد الأولية م.كل رئيس الوزراء ص موجود في قائمة واحدة بالضبط، أي.في أ[م % ر] لأصغر عدد صحيح في النطاق م .. م+ر-1 يقبل القسمة على ص.عند إنشاء العنصر التالي من المكرر، ما عليك سوى إدراج القائمة فيه أ[م % ر] يتم إرجاع.ثم تتم إزالة قائمة الأعداد الأولية من أ[م % ر] ويتم إلحاق كل رئيس الوزراء ل أ[(م+ص)% ر].
مع قائمة الأعداد الأولية <sqrt(ن) تقسيم م فمن السهل العثور على التخصيم م, ، نظرًا لوجود عدد أولي واحد أكبر من sqrt(ن).
هذه الطريقة لها تعقيد O(ن سجل (سجل (ن))) على افتراض أن جميع العمليات بما في ذلك عمليات القائمة تأخذ O(1).متطلبات الذاكرة هي O(sqrt(ن)).
يوجد لسوء الحظ بعض الحمل الثابت هنا، ومن ثم كنت أبحث عن طريقة أكثر أناقة لتوليد القيم phi(n)، لكنني لم أنجح.
إليك مولد بايثون فعال.التحذير هو أنه لا يعطي النتائج بالترتيب.تعتمد على 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) بسرعة.على سبيل المثال، إذا كان Pن هو الحد الأقصى الأولي الذي هو أقل من N، يمكنك إنشاء جميع الأعداد الصحيحة أقل من N
ص1 أنا 1 * ص2 أنا 2 * ...* صن أنا ن حيث كل طي يعمل من 0 إلى [سجل (N) / سجل (Pي)] ([] هي دالة الأرضية).
بهذه الطريقة، يمكنك الحصول على phi(k) بسرعة بدون تحليل أولي باهظ الثمن ولا تزال تتكرر خلال كل k بين 1 وN (ليس بالترتيب ولكن أعتقد أنك لا تهتم بالترتيب).
غربال totients ن:
(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)
نقاط التحسين المحتملة الأخرى:أفضل أو بديل gcd, rem و عضلات المعدة المهام
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
