计算 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)) 和 CPU 复杂度 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 = 10^12 时，它是 2400 快几倍！
我在我的 2Ghz Intel Core 2 笔记本电脑上测试了这个例程，它每秒计算超过 3,000,000 个 Phi() 值。按照这个速度，计算 10^12 个值大约需要 4 天。我还测试了它的正确性高达 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 年以来，世界上最好的数学家已经投入了许多 CPU 周期来解决这个问题，因为找到更快的方法来解决这个问题会造成计算能力的弱点。 RSA公钥算法. 。（RSA中的公钥和私钥都是基于phi(n)计算的，其中n是两个大素数的乘积。）
phi(k) 的计算必须使用 k 的素因数分解来完成，这是唯一合理的方法。如果您需要复习一下， 维基百科有这个公式.
如果你现在必须计算 1 到大 N 之间每个数字的所有素因数，那么在看到任何结果之前你就会老死，所以我会采取相反的方法，即使用其可能的质因数构建 N 以下的所有数字，即所有小于或等于 N 的素数。
因此，你的问题将类似于计算一个数字的所有除数，只是你不知道在因式分解中可以找到某个素数的最大次数是多少。调整最初由 Tim Peters 在 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 以下的所有素因数，我也有 在博客上谈论了它...请记住，虽然计算 10 以下的所有素数12 本身就是一项了不起的壮举......
这是来自欧拉 245 计划吗？我记得这个问题，我已经放弃了。
我发现计算 totient 的最快方法是将质因数 (p-1) 相乘，前提是 k 没有重复的因数（如果我没记错问题，则绝不会出现这种情况）。
因此，为了计算因子，最好使用 gmpy.next_prime 或者 派克姆 （椭圆曲线分解）。
您还可以按照 Jaime 的建议筛选主要因素。对于 10 以内的数字12, ，最大素因数低于 100 万，应该是合理的。
如果你记住因式分解，它可以进一步加速你的 phi 函数。
对于这类问题，我使用一个迭代器来返回每个整数 米<N 素数列表 < sqrt(氮）除以 m。为了实现这样的迭代器，我使用了一个数组 A 长度 右 在哪里 右 > 开方（氮）。数组中的每一点 A 包含可除整数的素数列表 米..米+R-1。IE。 A[米 % 右] 包含素数除法 米. 。每个素数 p 恰好在一个列表中，即在 A[米 % 右] 为范围内的最小整数 米 .. 米+右-1 可以整除 p. 。当生成迭代器的下一个元素时，只需列表中 A[米 % 右] 返回。然后从素数列表中删除 A[米 % 右] 并且每个质数 p 都附加到 A[(米+p)%R]。
具有素数列表 < sqrt(氮) 除 米 很容易求出因式分解 米, ，因为最多有一个大于 sqrt(氮).
该方法的复杂度为O(氮 日志（日志（氮))) 假设所有操作（包括列表操作）都需要 O(1)。内存需求为 O(sqrt(氮)).
不幸的是，这里存在一些恒定的开销，因此我一直在寻找一种更优雅的方法来生成值 phi(n)，但因此我没有成功。
这是一个高效的 python 生成器。需要注意的是，它不会按顺序产生结果。它是基于 https://stackoverflow.com/a/10110008/412529 .
内存复杂度为 O(log(N))，因为它一次只需存储单个数字的质因数列表。
CPU 复杂度只是勉强超线性，类似于 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)
我认为你可以反过来。您可以尝试从素数生成从 1 到 N 的所有整数并快速获得 phi(k)，而不是分解整数 k 来获得 phi(k)。例如，如果 Pn 是小于 N 的最大素数，您可以通过以下方式生成所有小于 N 的整数
磷1 我 1 * P2 我 2 * ...* Pn 我 n 其中每个我j 从0运行到[log(N)/log(Pj)]（[]是下取整函数）。
这样，您可以快速获得 phi(k) ，而无需昂贵的素因数分解，并且仍然迭代 1 和 N 之间的所有 k （不按顺序，但我认为您不关心顺序）。
将totents过筛至 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)
其他可能的改进点：更好或替代方案 最大CD, 雷姆 和 腹肌 功能
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
不隶属于 StackOverflow
