为 1<k<N 计算 phi(k)

Question

给定一个很大的 N，我需要遍历所有 phi(k) 使得 1

• 时间复杂度必须为O(N logN)
• memory-complexity 必须低于 O(N)（因为 N 的值约为 10¹²）

有可能吗？如果有，怎么做？

Answer 1

这可以通过内存复杂度 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 分解要快得多这里的一些回复所定位的算法。事实上，在 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），如果您需要，我可以提供。

如果您有任何问题，请告诉我。

---

根据 cmets 的要求，我在下面添加了 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))));
    }
}

Answer 2

没有人找到比首先找到 k 的质因数更快的方法来计算 phi(k)（又名Euler's totient function）。自 1977 年以来，世界上最优秀的数学家已经为这个问题投入了许多 CPU 周期，因为找到一种更快的方法来解决这个问题会在RSA public-key algorithm 中造成一个弱点。 （RSA中的公钥和私钥都是根据phi(n)计算的，其中n是两个大素数的乘积。）

Answer 3

必须使用 k 的素数分解来计算 phi(k)，这是唯一合理的方法。如果您需要对此进行复习，wikipedia carries the formula。

如果您现在必须计算 1 到大 N 之间的每个数字的所有素数除数，那么您将在看到任何结果之前死于老年，所以我会反过来，即构建所有低于 N 的数字, 使用它们可能的素数因子，即所有小于或等于 N 的素数。

因此，您的问题将类似于计算一个数字的所有除数，只是您事先不知道在分解中可以找到某个素数的最大次数是多少。调整最初由 Tim Peters 在 python 列表（I've blogged about...）上编写的迭代器以包含 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 的所有素数方面需要帮助，我也有 blogged about it... 请记住，虽然计算低于 10¹² 的所有素数本身就非常了不起壮举...

Answer 4

这是来自 Project Euler 245 的吗？我记得那个问题，我已经放弃了。

我发现计算 totient 的最快方法是将素数 (p-1) 相乘，因为 k 没有重复的因数（如果我没记错的话，情况就不会如此）。

所以对于计算因子，最好使用gmpy.next_prime 或pyecm（椭圆曲线分解）。

您也可以按照 Jaime 的建议筛选主要因素。对于高达 10¹² 的数字，最大素因数低于 100 万应该是合理的。

如果你记住因式分解，它可以进一步加快你的 phi 函数。

Answer 5

对于这类问题，我使用了一个迭代器，它为每个整数 m 返回除 m 的素数 N) 列表。为了实现这样的迭代器，我使用了一个长度为 R 的数组 A 其中 R > sqrt(N) .在每个点，数组 A 包含除以整数 m .. m+R-1 的素数列表。 IE。 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(log(N)))，假设包括列表操作在内的所有操作都需要O(1)。内存需求为 O(sqrt(N))。

不幸的是，这里有一些恒定的开销，因此我一直在寻找一种更优雅的方法来生成值 phi(n)，但是我没有成功。

Answer 6

这是一个高效的 p​​ython 生成器。需要注意的是，它不会按顺序产生结果。它基于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)

Answer 7

我认为你可以反过来。您可以尝试从素数生成从 1 到 N 的所有整数并快速获得 phi(k)，而不是分解整数 k 以获得 phi(k)。例如，如果 P_n 是小于 N 的最大素数，则可以通过

生成所有小于 N 的整数

P₁^{i ₁} * P₂^{i ₂} * ... * P_n^{i _n} 其中每个 i_j 从 0 运行到 [ log (N) / log (P_j)] （[] 是底函数）。

这样，您可以快速获得 phi(k)，而无需进行昂贵的素因数分解，并且仍然遍历 1 和 N 之间的所有 k（不是按顺序，但我认为您不关心顺序）。

Answer 8

将全部筛分到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)))))))))

Answer 9

这分解了 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