VERSION 1.0 CLASS
BEGIN
  MultiUse = -1  'True
END
Attribute VB_Name = "素数"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = False
Attribute VB_Exposed = False
Const SNmax = 100000                                            '作成する素数の数
Private Slist() As Long                                         '素数リスト
Dim Bobj As ベルヌーイ                                          'ベルヌーイクラスオブジェクト

'初期設定で素数リストを作成
Private Sub Class_Initialize()
    ReDim Slist(1 To SNmax, 1 To 1)
    Slist(1, 1) = 2: Slist(2, 1) = 3: cnt& = 2
    For i& = 5 To &H7FFFFFFF Step 2
        For k& = 2 To cnt&: sk& = Slist(k&, 1)
            If (i& Mod sk&) = 0 Then Exit For
            If i& < sk& ^ 2 Then
                cnt& = cnt& + 1
                Slist(cnt&, 1) = i&
                Exit For
            End If
        Next k&: If SNmax <= cnt& Then Exit For
    Next i&
End Sub

'n番目の素数
Public Function Nsosuu(n As Long) As Long
    If n = 0 Then Nsosuu = SNmax: Exit Function                 'n=0:素数の数
    If (n < 0) Or (n > SNmax) Then Nsosuu = -1: Exit Function   'nが異常
    Nsosuu = Slist(n, 1)                                        'n番目の素数を返却
End Function

'素因数分解
Public Function Soinsuu(n As LongLong) As String
    If n < 2 Then Exit Function
    m^ = n: d$ = ""
    For k& = 1 To SNmax: s^ = Slist(k&, 1): Ss% = 0
        Do Until m^ Mod s^: m^ = m^ \ s^: Ss% = Ss% + 1: Loop
        Select Case Ss%
        Case 0
        Case 1:     d$ = d$ & " * " & s^
        Case Else:  d$ = d$ & " * " & s^ & "^" & Ss%
        End Select: If m^ < s^ ^ 2 Then Exit For
    Next k&: If m^ > 1 Then d$ = d$ & " * " & m^
    Soinsuu = Mid(d$, 4)
End Function

'Moebius関数 [μ(n)]        '-1:素因数が奇数,0:同一素因が有る,1:素因数が偶数
Public Function Myu(n As Integer) As Integer
Static ms(&H7FFF) As Byte
    If ms(n) Then Myu = ms(n) - 2:    Exit Function
    If n < 2 Then Myu = 1: ms(n) = 3: Exit Function
    s$ = Me.Soinsuu((n)): a% = 0: m% = 0
    If InStr(s$, "^") Then Myu = 0: ms(n) = 2: Exit Function
    Do: a% = InStr(a% + 1, s$, "*"): If a% Then m% = m% + 1 Else Exit Do
    Loop: If m% Mod 2 Then Myu = 1: ms(n) = 3 Else Myu = -1: ms(n) = 1
End Function

'x以下の素数個数π(x) x=素数の場合:-0.5
Public Function Spi(x As Double, Optional f As Boolean = False) As Double
    If f Then
            s# = 0
            For m% = 1 To 20
                p# = Me.Myu(m%) * Me.Lpi(x ^ (1 / m%)) / m%
                If p# = 0 Then Exit For
                s# = s# + p#
            Next m%: Spi = s#
    Else:   If x < 2 Then Spi = 0: Exit Function
            r& = WorksheetFunction.XMatch(x, Slist, -1, 2)
            If x = Slist(r&, 1) Then Spi = r& - 0.5 _
                                Else Spi = r&
    End If
End Function

'Π(x)
Public Function Lpi(x As Double) As Double
    s# = 0
    For n% = 1 To 20
        p# = Me.Spi(x ^ (1 / n%)) / n%
        If p# = 0 Then Exit For
        s# = s# + p#
    Next n%
    Lpi = s#
End Function

'対数積分(対数の逆数を0〜xまで積分)  l=0:Li(x)
Public Function li(x As Double, Optional l As Double = 1.04516378011749) As Double
Const Gm = 0.577215664901533 - 1.04516378011749         'γ-li(2)
    lx# = Log(x): s# = Gm + l + Log(Abs(lx#)): m# = 1
    For k% = 1 To 100
        m# = m# * lx# / k%
        s# = s# + m# / k%
        If Abs(m#) < 1E-16 Then Exit For
    Next k%: li = s#
End Function

'Riemann関数：x以下の素数の数(近似値)
Public Function Riemann(x As Double) As Double
Static sn(170) As Double
    If x < 0.000000000000001 Then Exit Function
    lx# = Log(x): lxn# = 1: s# = 1
    For n% = 1 To 140: lxn# = lxn# * lx#
        If sn(n%) < 1 Then sn(n%) = n% * Me.Zeta(n% + 1) * WorksheetFunction.Fact(n%)
        r# = lxn# / sn(n%)
        s# = s# + r#: If Abs(r# / s#) < 0.000000000000001 Then Exit For
    Next n%: If Abs(r# / s#) < 0.0000000001 Then Riemann = s#
End Function

'x以下の素数の数π(x) z:ζ(x)の非自明零点(0.5±zi),n:使用する零点数
Public Function SpiZero(x As Double, Optional z As Range = Nothing, _
                                     Optional ByVal n As Integer = -1) As Double
    s# = Me.Riemann(x)
    For i% = 1 To 50
        r# = Me.Riemann(x ^ (-2 * i%))
        s# = s# - r#: If Abs(r#) < 0.001 Then Exit For
    Next i%
    If z Is Nothing Then SpiZero = s#: Exit Function
    If n < 0 Then n = z.Rows.Count
    For i% = 1 To n
        a# = z(i%, 1): alx# = a# * Log(x): ar# = Atn(a# / 0.5)
        mn% = Log(x) / Log(2) + 1: t# = 0
        For j% = 1 To mn%
            t# = t# + Myu(j%) * x ^ (1 / 2 / j%) * Cos(alx# / j% - ar#)
        Next j%
        s# = s# - 2 * t# / Sqr(0.25 + a# ^ 2) / Log(x)
    Next i%
    SpiZero = s#
End Function

'リーマンのゼータ関数ζ(x)
Public Function Zeta(x As Double, Optional z0 As Double = 1) As Double
Const n0 = 20, n2 = 1 / n0 / n0
    If Bobj Is Nothing Then Set Bobj = New ベルヌーイ
    cs# = z0
    For n% = 2 To n0:   cs# = cs# + n% ^ (-x): Next n%
    cm# = n0 ^ (1 - x): cs# = cs# + Bobj.Bn(0) * cm# / (x - 1)
                        cs# = cs# + Bobj.Bn(1) * cm# / n0
    cm# = cm# * x * n2: cs# = cs# + Bobj.Bn(2) / 2 * cm#
    For n% = 4 To 80 Step 2
        cm# = cm# * (x + n% - 3) * (x + n% - 2) * n2
        e# = Bobj.Bn(n%) / WorksheetFunction.Fact(n%) * cm#
        cs# = cs# + e#: If Abs(e#) < 1E-16 Then Exit For
    Next n%
    Zeta = cs#
End Function

'非自明の零点虚数部T以下の零点数n   零点はd=0.5
Public Function ZeroN(t As Double, Optional d As Double = 0.5) As Double
Const PAI2 = 3.14159265358979 * 2
    ZeroN = (t * (Log(t / PAI2) - 1) + 1 / 24 / t) / PAI2 + 7 / 8 + d
End Function