给定一个大的 N，我需要遍历所有 phi(k) 使得 1 < k < N ：
- 时间复杂度必须是
O(N logN) - memory-complexity 必须是 sub
O(N)（因为 N 的值将在 10 12左右）
是否可以？如果是这样，怎么做？
algorithm - 计算 phi(k) 为 1
给定一个大的 N，我需要遍历所有 phi(k) 使得 1 < k < N ：
- 时间复杂度必须是
O(N logN)
- memory-complexity 必须是 sub
O(N)（因为 N 的值将在 10 12左右）
是否可以？如果是这样，怎么做？
PythonPower
问问题
17926 次
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的所有素数本身就是一项了不起的壮举......
4
9 回答
21
这可以通过内存复杂度 O(Sqrt(N)) 和 CPU 复杂度 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 = 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/代码翻译器/）。所以请注意，这是自动转换的，我还没有时间自己测试或检查它。
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))));
}
}
于 2009-07-16T00:30:49.203 回答
4
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的所有素数本身就是一项了不起的壮举......
于 2009-06-21T21:34:57.373 回答
1
这是来自 Project Euler 245 的吗？我记得那个问题，我已经放弃了。
我发现计算 totient 的最快方法是将素数 (p-1) 相乘，因为 k 没有重复的因数（如果我没记错问题，情况就不会如此）。
因此，对于计算因子，最好使用gmpy.next_prime或pyecm（椭圆曲线分解）。
您也可以按照 Jaime 的建议筛选主要因素。对于高达 10 12的数字，最大素因数低于 100 万应该是合理的。
如果你记住分解，它可以进一步加快你的 phi 函数。
于 2009-06-21T21:56:46.240 回答
1
对于这类问题，我使用了一个迭代器，它为每个整数m < N返回除 m的素数 < sqrt( N ) 列表。为了实现这样一个迭代器，我使用了一个长度为R的数组A ，其中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]。
用一个素数列表 < sqrt( N ) 除m，很容易找到m的因式分解，因为最多有一个素数大于 sqrt( N )。
该方法的复杂度为 O( N log(log( N )))，假设包括列表操作在内的所有操作都需要 O(1)。内存要求为 O(sqrt( N ))。
不幸的是，这里有一些恒定的开销，因此我一直在寻找一种更优雅的方式来生成值 phi(n)，但因此我没有成功。
于 2009-06-23T06:45:43.250 回答
1
这是一个高效的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)
于 2012-12-24T04:34:10.410 回答
0
我想你可以反过来。您可以尝试从素数生成从 1 到 N 的所有整数并快速获得 phi(k)，而不是分解整数 k 以获得 phi(k)。例如，如果 P n是小于 N 的最大素数，则可以通过以下方式生成所有小于 N 的整数
P 1 i 1 * P 2 i 2 * ... * P n i n其中每个 i j从 0 运行到 [log (N) / log (P j )] （[] 是底函数）。
这样，您无需昂贵的素数分解即可快速获得 phi(k)，并且仍然遍历 1 和 N 之间的所有 k（不是按顺序，但我认为您不关心顺序）。
于 2009-07-12T22:56:42.150 回答
0
将 totients 筛分到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)))))))))
于 2012-05-06T14:23:34.510 回答
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
于 2016-06-19T02:32:14.930 回答
给定一个大的 N，我需要遍历所有 phi(k) 使得 1 < k < N ：
- 时间复杂度必须是
O(N logN) - memory-complexity 必须是 sub
O(N)（因为 N 的值将在 10 12左右）
是否可以？如果是这样，怎么做？
PythonPower
问问题
17926 次
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的所有素数本身就是一项了不起的壮举......
9 回答
21
这可以通过内存复杂度 O(Sqrt(N)) 和 CPU 复杂度 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 = 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/代码翻译器/）。所以请注意，这是自动转换的，我还没有时间自己测试或检查它。
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))));
}
}
于 2009-07-16T00:30:49.203 回答
4
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的所有素数本身就是一项了不起的壮举......
于 2009-06-21T21:34:57.373 回答
1
这是来自 Project Euler 245 的吗？我记得那个问题，我已经放弃了。
我发现计算 totient 的最快方法是将素数 (p-1) 相乘，因为 k 没有重复的因数（如果我没记错问题，情况就不会如此）。
因此，对于计算因子，最好使用gmpy.next_prime或pyecm（椭圆曲线分解）。
您也可以按照 Jaime 的建议筛选主要因素。对于高达 10 12的数字，最大素因数低于 100 万应该是合理的。
如果你记住分解，它可以进一步加快你的 phi 函数。
于 2009-06-21T21:56:46.240 回答
1
对于这类问题，我使用了一个迭代器，它为每个整数m < N返回除 m的素数 < sqrt( N ) 列表。为了实现这样一个迭代器，我使用了一个长度为R的数组A ，其中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]。
用一个素数列表 < sqrt( N ) 除m，很容易找到m的因式分解，因为最多有一个素数大于 sqrt( N )。
该方法的复杂度为 O( N log(log( N )))，假设包括列表操作在内的所有操作都需要 O(1)。内存要求为 O(sqrt( N ))。
不幸的是，这里有一些恒定的开销，因此我一直在寻找一种更优雅的方式来生成值 phi(n)，但因此我没有成功。
于 2009-06-23T06:45:43.250 回答
1
这是一个高效的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)
于 2012-12-24T04:34:10.410 回答
0
我想你可以反过来。您可以尝试从素数生成从 1 到 N 的所有整数并快速获得 phi(k)，而不是分解整数 k 以获得 phi(k)。例如，如果 P n是小于 N 的最大素数，则可以通过以下方式生成所有小于 N 的整数
P 1 i 1 * P 2 i 2 * ... * P n i n其中每个 i j从 0 运行到 [log (N) / log (P j )] （[] 是底函数）。
这样，您无需昂贵的素数分解即可快速获得 phi(k)，并且仍然遍历 1 和 N 之间的所有 k（不是按顺序，但我认为您不关心顺序）。
于 2009-07-12T22:56:42.150 回答
0
将 totients 筛分到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)))))))))
于 2012-05-06T14:23:34.510 回答
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
于 2016-06-19T02:32:14.930 回答