给定一个大的N，我需要迭代所有phi(k)，使得 1 < k < N:
- 时间复杂度必须为
O(N logN) - 内存复杂度必须是子
O(N)（因为N的值将约为10^12）
是否可能？如果可能，如何实现？
计算1<k<N时的phi(k)
14
- PythonPower
7
6那是一个欧拉计划的题目吗？ - Marc Gravell
1请参考“有多少个小于N的数与N互质？”：https://dev59.com/jXNA5IYBdhLWcg3wUL1Y，关于如何快速计算phi(k)的同一问题。 - ShreevatsaR
1这不是一个欧拉计划题目，但它是由其中一个问题引发的，我现在已经解决了。 - PythonPower
2ShreevatsaR：那个问题是关于孤立地计算 Phi(k)。这个问题是关于按顺序计算所有 Phi(k)，直到某个 N。这是一个相关但仍然显着不同的问题。 - RBarryYoung
1顺序不重要；"quickly" 是 O(n log n) 的成员。 - PythonPower
显示剩余2条评论
9个回答
21
这可以通过优化的窗口埃拉托色尼筛法实现，其内存复杂度为O(Sqrt(N))，CPU复杂度为O(N * Log(Log(N)))。以下是代码示例。
由于未指定编程语言且我不熟悉Python，因此我使用了VB.net实现，但如果需要，我可以将其转换为C#。
请注意，在 O(N * Log(Log(N))) 的时间复杂度下，这个程序可以在平均 O(Log(Log(N))) 的时间内对每个数字进行因式分解，比某些答案中提到的最快的单个 N 因式分解算法要快得多。事实上，在 N = 10^12 时，它快了2400倍！
我已经在我的 2Ghz 英特尔 Core 2 笔记本电脑上测试过这个程序，并且每秒计算超过 3,000,000 个 Phi() 值。以这个速度计算 10^12 个值需要大约 4 天时间。我还在没有任何错误的情况下对其进行了正确性测试，最高测试值为 100,000,000。它基于 64 位整数，因此精确度可以达到 2^63 (10^19)，尽管速度对任何人来说都太慢了。
如果您需要，我还有一个 Visual Studio WinForm（也是 VB.net）用于运行/测试该程序。
如果您有任何问题，请告诉我。
由于未指定编程语言且我不熟悉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倍！
我已经在我的 2Ghz 英特尔 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))));
}
}
- RBarryYoung
4
非常感谢您提供如此清晰易读的代码。同时也感谢您提供了一些有趣的基准测试，尽管我预计与优化后的C相比，VB.net可能会慢得多；我期望10^12个数字的运行时间可以缩短到一天或更短。 - PythonPower
是的，我同意，使用高度优化的中级语言如C，将1天的速度提高4倍似乎是合理的。真正的挑战将是如何处理每秒10百万个Totient值而不会显著减慢速度。 - RBarryYoung
你好，我正在进行一个研究项目，希望能够实现你上面所描述的算法。但是，我对VB.net不熟悉。是否可以给我发送伪代码，或者用C语言实现它？请告诉我。谢谢。 - nicole
@nicole 我已经更新了我的答案，加入了自动生成的转换代码。我现在没有时间测试它，但是事件声明看起来好像没有正确转换。 - RBarryYoung
5
- Bill the Lizard
4
你说的没错，但只适用于(k)的单个值。我不知道有什么理由相信O(Compute(phi(k)))|0<k<=N == N*O(Compute(phi(k)))。实际上，乍一看，像优化的窗口埃拉托色尼筛法这样简单的东西就证明了另外一种情况。 - RBarryYoung
没有人会争论 O(Compute(phi(k))) | 0 < k <= N == N * O(Compute(phi(k))) 这一点。问题仍然是要找到极大数 k 的质因数。 - Bill the Lizard
2在这个领域，10^12（或2^40）被认为是适中的，不算大或非常大。问题并不是找到k的质因数，而是找到所有k≤N的质因数，这可以比用于分解RSA公钥等事物的例程快几千倍来完成。 - RBarryYoung
好的。我本意是说“非常大量的k”，这与之前完全不同。 - Bill the Lizard
4
计算phi(k)必须使用k的质因数分解，这是唯一合理的方法。如果您需要恢复记忆，请参考wikipedia上的公式。
现在，如果您需要计算1到大N之间每个数字的所有质因子，您将在看到任何结果之前死去，所以我建议采用另一种方式，即使用它们可能的质因数构建所有小于N的数字，即所有小于或等于N的质数。
因此，您的问题将类似于计算数字的所有除数，只是您事先不知道在因式分解中可能找到某个质数的最大次数。通过调整最初由Tim Peters编写的Python列表迭代器（我已经在博客中介绍过...）以包括欧拉函数，可以得到以下可能的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的所有质数本身就是一项相当了不起的成就...
- Jaime
2
1
这里有一个高效的Python生成器。但是需要注意的是，它不会按顺序产生结果。它基于https://dev59.com/BWkw5IYBdhLWcg3wJHKF#10110008。
内存复杂度为O(log(N))，因为它只需要一次为单个数字存储质因数列表。
CPU复杂度仅略微超线性，大约是O(N log log N)。
内存复杂度为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)
- jnnnnn
1
这是来自欧拉计划245题吗？我记得那道题，而且我已经放弃了。
我发现计算totient的最快方法是将质因数(p-1)相乘，假设k没有重复的因子（如果我正确地记得问题的话，这从未是个问题）。
因此，对于计算因子，最好使用gmpy.next_prime或pyecm（椭圆曲线分解）。
您也可以像Jaime建议的那样筛选质因数。对于不超过1012的数字，最大质因数低于100万，这应该是合理的。
如果您记忆化因式分解，它可以进一步加速您的phi函数。
- Unknown
1
对于这类问题，我使用一个迭代器，为每个整数sqrt(N)。在每个点上，数组A包含整数m..m+R-1的质因数列表。即A[m%R]包含除m外的质因数。每个质数p只在一个列表中，即在范围m..m+R-1中可被p整除的最小整数所在的A[m%R]中。当生成迭代器的下一个元素时，只需返回A[m%R]中的列表。然后从A[m%R]中删除质数列表，并将每个质数p附加到A[(m+p)%R]中。
如果有一个小于sqrt(N)的质数列表可以整除m，那么很容易找到m的因式分解，因为最多只有一个大于sqrt(N)的质数。
在假设所有操作（包括列表操作）都需要O(1)的情况下，该方法的复杂度为O(N log(log(N)))。内存需求为O(sqrt(N))。
不幸的是，这里存在一些常量开销，因此我正在寻找更优雅的方法来生成phi(n)的值，但目前还没有成功。
- Accipitridae
3
谢谢您的帖子。我正在努力计算数组A中的值。如果N=20且R=5，那么A会是什么呢？干杯！ - PythonPower
非常酷的方法。但是你不需要将质数列表排序（或至少分组）以便从中计算Phi(m)吗？你如何保持它们排序，或者每个m都必须进行排序？ - RBarryYoung
很抱歉我没有早些回答。我不需要对质数列表进行排序。所有必要的是，当迭代器返回小于sqrt（N）的质因子列表<m的A[m％R]包含这些质因子。为了得到phi，仍然需要分解m，但是鉴于我们拥有所有质因子，这很容易。我最近实现了这种方法：计算n <10 ^ 12的所有phi（n）需要17小时在Core 2上完成。我可能比筛法做更多的操作，但我的方法似乎更加缓存友好，因为我使用的内存少于5MB。 - Accipitridae
0
这个因式分解N = PQ，其中P和Q是质数。
在Elixir或Erlang中运行得非常好。
您可以尝试不同的生成器来生成伪随机序列。通常使用x * x + 1。
这一行：defp f0(x, n), do: rem((x * x) + 1, n)
其他可能改进的点：更好或替代的gcd，rem和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
- Charles Okwuagwu
0
筛选小于等于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)))))))))
- user448810
0
我认为你可以采用另一种方法。不必将整数k分解为phi(k)，而是尝试从质数生成所有1到N的整数，并快速获得phi(k)。例如，如果Pn是小于N的最大质数，则可以通过以下方式生成小于N的所有整数：
P1 i 1 * P2 i 2 * ... * Pn i n，其中每个ij从0到[log (N) / log (Pj)]([]是向下取整函数)。
这样，您可以快速获得phi(k)，而无需进行昂贵的质因数分解，并且仍然可以迭代1到N之间的所有k（不按顺序，但我认为您不关心顺序）。
- gineer
1
关于内存，根据素数定理，n ~ N/ln(N)。因此，一个良好的实现应该具有O(n/log n)的内存使用率。在这里查看素数定理：http://en.wikipedia.org/wiki/Prime_number_theorem - gineer
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接