Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
437 views
in Technique[技术] by (71.8m points)

algorithm - Calculating phi(k) for 1<k<N

Given a large N, I need to iterate through all phi(k) such that 1 < k < N :

  • time-complexity must be O(N logN)
  • memory-complexity must be sub O(N) (since the values of N will be around 1012)

Is it possible? If so, how?

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

This can be done with Memory complexity O(Sqrt(N)) and CPU complexity O(N * Log(Log(N))) with an optimized windowed Sieve of Eratosthenes, as implemented in the code example below.

As no language was specified and as I do not know Python, I have implemented it in VB.net, however I can convert it to C# if you need that.

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

Note that at O(N * Log(Log(N))), this routine is factoring each number at O(Log(Log(N))) on average which is much, much faster than the fastest single N factorization algorithms sited by some of the replies here. In fact, at N = 10^12 it is 2400 times faster!

I have tested this routine on my 2Ghz Intel Core 2 laptop and it computes over 3,000,000 Phi() values per second. At this speed, it would take you about 4 days to compute 10^12 values. I have also tested it for correctness up to 100,000,000 without any errors. It is based in 64-bit integers, so anything up to 2^63 (10^19) should be accurate (though too slow for anyone).

I also have a Visual Studio WinForm (also VB.net) for running/testing it, that I can provide if you want it.

Let me know if you have any questions.


As requested in the comments, I have added below a C# version of the code. However, because I am currently in the middle of some other projects, I do not have time to convert it myself, so I have used one of the online VB to C# conversion sites (http://www.carlosag.net/tools/codetranslator/). So be aware that this was auto-converted and I have not had time to test or check it myself yet.

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))

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

1.4m articles

1.4m replys

5 comments

57.0k users

...