Вычисление phi (k) для 1 <k <N

Question

Учитывая большой N, мне нужно выполнить итерацию по всем phi (k) таким образом, чтобы 1 < k < N быстро. Так как значения N будут около 10 , важно, чтобы сложность памяти была sub O (n).

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

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 это раз быстрее!

Я протестировал эту процедуру на своем ноутбуке Intel Core 2 на 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, что является единственным разумным способом ее выполнения. Если вам нужен обновитель, wikipedia carries the formula.

Если вам нужно вычислить все простые делители каждого числа между 1 и большим N, вы умрете от старости, прежде чем увидеть какой-либо результат, поэтому я бы пошел наоборот, то есть построил все числа ниже N, используя их возможные простые множители, т. Е. Все простые числа, меньшие или равные N.

Ваша проблема, таким образом, будет аналогична вычислению всех делителей числа, только вы не знаете, какое максимальное количество раз вы можете заранее найти некоторое преимущество в факторизации. Подстройка итератор первоначально написанный Тимом Петерс в списке питона (то I've blogged about ...), чтобы включить функцию totient, возможную реализацию в питона, который дает к, Фи (к) пар может быть следующим:

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 3

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

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

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

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

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

Answer 4

Никто не нашел более быстрый способ вычислить phi (k) (aka, Euler's totient function), чем сначала найти простые множители k. Лучшие математики мира бросили много циклов процессора в проблеме с 1977 года, так как поиск более быстрого способа решения этой проблемы создаст слабость в RSA public-key algorithm. (Как публичный, так и закрытый ключ в RSA рассчитываются на основе phi (n), где n - произведение двух больших простых чисел.)

Answer 5

Для таких проблем я использую итератор, который возвращает для каждого целого числа m < N Список простых чисел < sqrt (N), которые разделяют m. Для реализации такого итератора я использую массив длины R где R > SQRT (N).В каждой точке массив A содержит список простых чисел, которые делят целые числа m .. m + R -1. То есть A [m% R] содержит простые числа, делящие м. Каждое простое р находится точно в одном списке, то есть в [м% R] для наименьшее целое число в диапазоне м .. м + R -1, что делится на p. При генерации следующего элемента итератора просто список в A [m% R] возвращается. Затем список простых чисел удаляются из [м% R] и каждый простой р добавляется к [(м + р)% R].

В списке простых чисел < SQRT (N), разделяющей м, легко найти разложение м, так как есть более одного премьера больше, чем SQRT (N).

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

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

Answer 6

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

P ^я * Р ^я * ... * Р _п^{я _п} где каждый я _J запуска от 0 до [журнал (N)/журнал (Р _J)] ([] - функция пола).

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

Answer 7

просеять totients к н:

(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 8

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

Это дофакторизовывается N = PQ, где Р & Q являются простыми.

Хорошо работает в Эликсире или Эрланге.

Вы можете попробовать различные генераторы для вашей псевдослучайной последовательности. x*x + 1 обычно используется.

Эта линия: defp f0(x, n), do: rem((x * x) + 1, n)

Другие возможные точки улучшения: лучше или альтернативный НОД, бэр и абс функции

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