计算岛(K)1< K< N [英] Calculating phi(k) for 1<k<N

查看:99
本文介绍了计算岛(K)1< K< N的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

给出一个大的N,我需要遍历所有岛(K),使得1< K< ñ快。由于N的值将是大约10 12 ,重要的是该存储器的复杂性是子为O(n)。

Given a large N, I need to iterate through all phi(k) such that 1 < k < N quickly. Since the values of N will be around 1012, it is important that the memory complexity is sub O(n).

这可能吗?如果是这样,怎么样?

Is it possible? If so, how?

推荐答案

这可与内存复杂度为O(SQRT(N))和CPU的复杂度为O完成的(N *日志(日志(N)))与优化的窗口埃拉托色尼的筛,如code例如实现如下。

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.

由于没有语言被指定,当我不知道Python的,我已经在VB.net实现它,但是我可以,如果你需要转换为C#。

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

请注意,在O(N *日志(日志(N))),该程序被分解每一个数字在澳平均(日志(日志(N))),这是很多,比最快的单个n分解得更快算法通过选址这里的一些答复。事实上,在N = 10 ^ 12它的 2400 快倍!

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!

我已经测试这个程序在我的2GHz的英特尔酷睿2笔记本电脑和它计算每秒超过300万岛()值。在这个速度下,它会带你4天左右来计算10 ^ 12个值。我还测试了它的正确性高达亿没有任何错误。它是基于64位整数,所以任何达到2 ^ 63(10 ^ 19)应该是准确的(尽管任何人都太慢了)。

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

我也有一个Visual Studio的WinForm(也VB.net)运行/测试它,我可以提供,如果你想要它。

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.

这篇关于计算岛(K)1&LT; K&LT; N的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
相关文章
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆