Расчет phi (k) для 1 <k <N

Question

Учитывая большое N, мне нужно перебрать все фи (k) так, чтобы 1

• сложность времени должна быть O(N logN)
• сложность памяти должна быть меньше O(N) (поскольку значения N будут около 10 ¹² )

Возможно ли это? Если да, то как?

Answer 1

Это можно сделать с помощью сложности памяти 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 раз быстрее!

Я протестировал эту процедуру на моем 2-ГГц ноутбуке Intel Core 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))));
    }
}

Answer 2

Никто не нашел более быстрого способа вычисления phi (k) (иначе, общая функция Эйлера ), чем сначала найти простые множители k. Лучшие математики мира с 1977 года бросили в процесс много циклов ЦП, поскольку поиск более быстрого способа решения этой проблемы может привести к слабости алгоритма RSA . (И открытый, и закрытый ключи в RSA рассчитываются на основе phi (n), где n - произведение двух больших простых чисел.)

Answer 3

Вычисление 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, я также написал об этом в блоге ... Имейте в виду, хотя вычисление всех простых чисел ниже 10 ¹² сам по себе довольно замечательный подвиг ...

Answer 4

Вот эффективный генератор Python. Предостережение в том, что оно не дает результатов по порядку. Он основан на 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)

Answer 5

Для такого рода проблем я использую итератор, который возвращает для каждого целого числа m список простых чисел N ), которые делят m. Для реализации такого итератора я использую массив A длиной R , где 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].

Со списком простых чисел N ), разделяющих m , легко найти факторизацию m , поскольку существует не более одного простого числа больше чем sqrt ( N ).

Этот метод имеет сложность O ( N log (log ( N ))) в предположении, что все операции, включая операции со списком, принимают O (1). Требуемый объем памяти составляет O (sqrt ( N )).

К сожалению, здесь есть некоторые постоянные издержки, поэтому я искал более элегантный способ генерирования значений phi (n), но так как я не добился успеха.

Answer 6

Это из Project Euler 245? Я помню этот вопрос, и я отказался от него.

Самый быстрый способ, который я нашел для вычисления коэффициента, состоял в том, чтобы умножить простые множители (p-1) вместе, учитывая, что k не имеет повторяющихся факторов (что никогда не было так, если я правильно помню проблему).

Так что для расчета факторов, вероятно, было бы лучше использовать gmpy.next_prime или pyecm (факторизация эллиптической кривой).

Вы также можете определить основные факторы, как предполагает Хайме. Для чисел до 10 ¹² максимальный простой коэффициент меньше 1 миллиона, что должно быть разумным.

Если вы запомните факторизации, это может ускорить вашу функцию phi.

Answer 7

Это множитель N = PQ, где P & Q простые.

Работает довольно хорошо, в эликсире или эрланге.

Вы можете попробовать разные генераторы для вашей псевдослучайной последовательности. 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

Answer 8

Сито 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)))))))))

Answer 9

Я думаю, вы можете пойти другим путем. Вместо того, чтобы разложить целое число k, чтобы получить phi (k), вы можете попытаться сгенерировать все целые числа от 1 до N из простых чисел и быстро получить phi (k). Например, если P _n - максимальное простое число, которое меньше N, вы можете сгенерировать все целые числа, меньшие N, с помощью

P ₁ ^{i ₁} * P ₂ ^{i ₂} *. .. * P _n ^{i _n} , где каждый i _j работает от 0 до [log (N) / log (P * 1025) * j )] ([] - функция пола).

Таким образом, вы можете быстро получить phi (k) без дорогой простой факторизации и по-прежнему перебирать все k между 1 и N (не по порядку, но я думаю, что вы не заботитесь о порядке).