Attribute VB_Name = "多倍長整数"
'多倍長整数ルーチン集（1億進数[10進数8桁]整数）

Const Ketasuu = 1000                            '10進数桁数
Const Ncnt = Ketasuu \ 8                        '配列数
Const BASE = 100000000                          '1億
Const BASE1 = BASE - 1                          '最大数=1億-1

'多倍長整数形式の定義
Public Type Lng
    Pm As Byte                                  '0:正数、1:負数
    Dt(0 To Ncnt) As Long                       '1億進数(10進数8桁)の配列,Dt(0)はデータ数
End Type

Dim MdSv As Lng                                 '余り保存
Dim Sobj As 素数                                '素数クラスのオブジェクト

'------- ExcelのWorkSheetでユーザー定義型が使用できない為、文字列に変換 -----------
' オプションのfは数字の区切記号","を付加するかを指定

'多倍長整数_加算 x+y (WorkSheet用)
Public Function LngAddStr(x As String, y As String, Optional f As Boolean = True) As String
    LngAddStr = Lng2Str(LngAdd(Str2Lng(x), Str2Lng(y)), f)
End Function

'多倍長整数_減算 x-y (WorkSheet用)
Public Function LngSubStr(x As String, y As String, Optional f As Boolean = True) As String
    LngSubStr = Lng2Str(LngSub(Str2Lng(x), Str2Lng(y)), f)
End Function

'多倍長整数_乗算 x*y (WorkSheet用)
Public Function LngMulStr(x As String, y As String, Optional f As Boolean = True) As String
    LngMulStr = Lng2Str(LngMul(Str2Lng(x), Str2Lng(y)), f)
End Function

'多倍長整数_除算 x/y (WorkSheet用) 出力:(商,余り)の配列
Public Function LngDivStr(x As String, y As String, Optional f As Boolean = True) As Variant
Dim a As Lng
    a = LngDiv(Str2Lng(x), Str2Lng(y))
    LngDivStr = Array(Lng2Str(a, f), Lng2Str(MdSv, f))
End Function

'多倍長整数_合同式 x mod m (WorkSheet用)
Public Function LngModStr(x As String, m As String, Optional f As Boolean = True) As String
    LngModStr = Lng2Str(LngMod(Str2Lng(x), Str2Lng(m)), f)
End Function

'多倍長整数_最大公約数 GCD(x,y) (WorkSheet用)   ユークリッドの互除法
Public Function LngGcdStr(x As String, y As String, Optional f As Boolean = True) As String
    LngGcdStr = Lng2Str(LngGcd(Str2Lng(x), Str2Lng(y)), f)
End Function

'多倍長整数_最小公倍数 LCM(x,y) (WorkSheet用)   x*y/GCD(x,y)
Public Function LngLcmStr(x As String, y As String, Optional f As Boolean = True) As String
    LngLcmStr = Lng2Str(LngLcm(Str2Lng(x), Str2Lng(y)), f)
End Function

'オイラー関数φ(n) (WorkSheet用)  nと互いに素の数
Public Function LngFaiStr(n As String) As Long
    LngFaiStr = LngFai(Str2Lng(n))
End Function

'オイラー関数φ(n) (検証用)  nと互いに素の数
Public Function LongFai(n As Long) As Long
Dim nl As Lng, c1 As Lng
    If n > 10000 Then Exit Function
    If n = 1 Then LongFai = 1: Exit Function
    c& = 0: nl = LngLongSet(n): c1 = LngLongSet(1)
    For i& = 1 To n - 1     '互いに素(GCD(n,i&)=1)をカウント
        If LngZeroChk(LngSub(LngGcd(nl, LngLongSet(i&)), c1)) Then c& = c& + 1
    Next i&
    LongFai = c&
End Function

'多倍長整数_累乗 x^n (WorkSheet用)
Public Function LngPowStr(x As String, n As Long, Optional f As Boolean = True) As String
    LngPowStr = Lng2Str(LngPow(Str2Lng(x), n), f)
End Function

'多倍長整数_合同累乗 x^p mod m (WorkSheet用)
Public Function LngModPowStr(x As String, p As String, m As String, Optional f As Boolean = True) As String
    LngModPowStr = Lng2Str(LngModPow(Str2Lng(x), Str2Lng(p), Str2Lng(m)), f)
End Function

'多倍長整数_階乗 n! (WorkSheet用)
Public Function LngFactStr(n As Long, Optional f As Boolean = True) As String
    LngFactStr = Lng2Str(LngFact(n), f)
End Function

'多倍長整数_順列 nPk=n!/(n-k)! (WorkSheet用)
Public Function LngPermutStr(n As Long, k As Long) As String
    LngPermutStr = Lng2Str(LngPermut(n, k))
End Function

'多倍長整数_組合 nCk=nPk/k!=n!/(n-k)!/k! (WorkSheet用)
Public Function LngCombinStr(n As Long, k As Long) As String
    LngCombinStr = Lng2Str(LngCombin(n, k))
End Function

'多倍長整数素因数分解 (WorkSheet用) & 素数判定(先頭が＝:素数／→:合成数)
Public Function LngSoinsuuStr(n As String, Optional f As Boolean = True, _
                                           Optional m As Long = 7) As String
Dim nl As Lng
    t! = Timer                                                                  '現在時刻
    nl = Str2Lng(n): s$ = LngSoinsuu(nl, f)                                     '素因数分解
    If (InStr(s$, " * ") = 0) And (InStr(s$, "^") = 0) Then
'           If LngSosuuTest(nl, m) Then s$ = "＝ " & s$ _
'                                  Else s$ = "→ ？(素数リスト不足)"            '素数検査
            s$ = "＝ " & s$
    Else:   s$ = "→ " & s$                                                     '＝:素数 / →:合成数
    End If
    LngSoinsuuStr = s$
    t! = Timer - t!: If t! > 1 Then Debug.Print Format(t!, "0.0秒"), Len(n), s$ '計算時間
End Function

'素因数分解の検証 (WorkSheet用)
Public Function LngSoinsuuChk(n As String, Optional f As Boolean = True) As String
Dim s As Lng:
    s = LngLongSet(1)
    For i% = 3 To 1023
        a1% = InStr(i%, n, " * ")                                               '区切り記号
        If a1% = 0 Then s1$ = Mid(n, i%) _
                   Else s1$ = Mid(n, i%, a1% - i%)                              '素因数取り出し
        a2% = InStr(s1$, "^")                                                   '同一素数が有るか？
        If a2% = 0 Then s = LngMul(s, Str2Lng(s1$)) _
                   Else s = LngMul(s, LngPow(Str2Lng(Left(s1$, a2% - 1)), _
                                             Val(Mid(s1$, a2% + 1))))           '数値変換
        If a1% = 0 Then Exit For                                                '終了判定
        i% = a1% + 2                                                            '次文字位置
    Next i%
    LngSoinsuuChk = Lng2Str(s, f)                                               '文字変換
End Function

'フェルマーテストによる素数検査   m:検査回数
Public Function LngSosuuTestStr(n As String, Optional m As Long = 7) As Boolean
    t! = Timer
    LngSosuuTestStr = LngSosuuTest(Str2Lng(n), m)
    t! = Timer - t!: If t! > 0.1 Then Debug.Print Format(t!, "0.00秒"), n       '計算時間
End Function

'Lucasテスト Mpが素数か？
Public Function LucasTest(p As Long) As Boolean
Dim mp As Lng, xi As Lng, c2 As Lng
    xi = LngLongSet(4):   c2 = LngLongSet(2)
    mp = LngSub(LngPow(c2, p), LngLongSet(1))
    For i& = 2 To p - 1
        xi = LngMod(LngSub(LngMul(xi, xi), c2), mp)
    Next i&
    LucasTest = LngZeroChk(xi)
End Function

'---------- 素数クラス -----------------

'n番目の素数
Public Function Nsosuu(n As Long) As Long
    If Sobj Is Nothing Then Set Sobj = New 素数                                 '素数リスト生成(素数クラス)
    Nsosuu = Sobj.Nsosuu(n)
End Function

'Moebius関数 [μ(n)]        '-1:素因数が奇数,0:同一素因が有る,1:素因数が偶数
Public Function Myu(n As Integer) As Integer
    If Sobj Is Nothing Then Set Sobj = New 素数                                 '素数リスト生成(素数クラス)
    Myu = Sobj.Myu(n)
End Function

'x以下の素数の数 π(x)   素数は-0.5
Public Function Spi(x As Double) As Double
    If Sobj Is Nothing Then Set Sobj = New 素数                                 '素数リスト生成(素数クラス)
    Spi = Sobj.Spi(x)
End Function

'Π(x)   Σ(1/n)π(x^(1/n))
Public Function Lpi(x As Double) As Double
    If Sobj Is Nothing Then Set Sobj = New 素数                                 '素数リスト生成(素数クラス)
    Lpi = Sobj.Lpi(x)
End Function

'対数積分 Li(x)   ∫dx/Ln(x)
Public Function li(x As Double) As Double
    If Sobj Is Nothing Then Set Sobj = New 素数                                 '素数リスト生成(素数クラス)
    li = Sobj.li(x, 0)
End Function

'リーマン関数 R(x)
Public Function Riemann(x As Double) As Double
    If Sobj Is Nothing Then Set Sobj = New 素数                                 '素数リスト生成(素数クラス)
    Riemann = Sobj.Riemann(x)
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
    If Sobj Is Nothing Then Set Sobj = New 素数                                 '素数リスト生成(素数クラス)
    SpiZero = Sobj.SpiZero(x, z, n)
End Function

'リーマンのゼータ関数 ζ(x)
Public Function Zeta(x As Double, Optional z0 As Double = 1) As Double
    If Sobj Is Nothing Then Set Sobj = New 素数                                 '素数リスト生成(素数クラス)
    Zeta = Sobj.Zeta(x, z0)
End Function

'非自明の零点虚数部t以下の零点数n   零点は+0.5
Public Function ZeroN(t As Double, Optional d As Double = 0.5) As Double
    If Sobj Is Nothing Then Set Sobj = New 素数                                 '素数リスト生成(素数クラス)
    ZeroN = Sobj.ZeroN(t, d)
End Function

'--------------- VBA専用の内部関数 -------------------------

'文字列→多倍長整数_変換
Private Function Str2Lng(x As String) As Lng
Dim lx As Lng
    s$ = Trim(x): If Left(s$, 1) = "+" Then s$ = Mid(s$, 2)                     '不要文字削除
    For i% = 1 To Ketasuu
        a% = InStr(s$, ","): If a% = 0 Then Exit For                            '","が有るか？
        s$ = Left(s$, a% - 1) & Mid(s$, a% + 1)                                 '区切記号","を削除
    Next i%
    For i% = 1 To Ncnt                                                          '10進8桁を1億進数に変換
        s1$ = Right(s$, 8)
        If Left(s1$, 1) = "-" Then                                              'マイナス処理
            If Len(s1$) > 1 Then lx.Dt(i%) = Val(Mid(s1$, 2))                   '"-"を削除
            lx.Pm = 1: lx.Dt(0) = i%: Exit For                                  '符号とデータ数設定
        End If
        lx.Dt(i%) = Val(s1$)                                                    'i%桁目設定
        If Len(s$) <= 8 Then lx.Dt(0) = i%: Exit For                            '残りが8文字以下で終了
        s$ = Left(s$, Len(s$) - 8)                                              '最後の8文字削除
    Next i%
    Str2Lng = lx
End Function

'多倍長整数→文字列_変換     f=True:3桁毎に","を挿入
Private Function Lng2Str(x As Lng, Optional f As Boolean = True) As String
    s$ = ""
    For i% = 1 To x.Dt(0) - 1                                                   '1億進数を10進8桁に変換
        s$ = Right("00000000" & x.Dt(i%), 8) & s$
    Next i%
    s$ = x.Dt(i%) & s$
    If f Then
        For i% = Len(s$) - 3 To 1 Step -3                                       '区切記号","を付加
            s$ = Left(s$, i%) & "," & Mid(s$, i% + 1)
        Next i%
    End If
    If x.Pm Then s$ = "-" & s$                                                  '符号を付加
    Lng2Str = s$
End Function

'多倍長整数_加算 x+y
Private Function LngAdd(x As Lng, y As Lng) As Lng
Dim a As Lng
    k& = 0
    If x.Pm = y.Pm Then
            For i% = 1 To Ncnt
                s& = x.Dt(i%) + y.Dt(i%) + k&                                   '加算
                If s& > BASE1 Then k& = 1: s& = s& - BASE _
                              Else k& = 0                                       '繰り上げ判定
                a.Dt(i%) = s&
                If (x.Dt(0) < i%) And (y.Dt(0) < i%) Then e% = i%: Exit For
            Next i%
    Else:   For i% = 1 To Ncnt
                s& = x.Dt(i%) - y.Dt(i%) - k&                                   '減算
                If s& < 0 Then k& = 1: s& = s& + BASE _
                          Else k& = 0                                           '繰り下げ判定
                a.Dt(i%) = s&
                If (x.Dt(0) < i%) And (y.Dt(0) < i%) Then e% = i%: Exit For
            Next i%: a.Pm = k&                                                  '符号設定
            If k& Then                                                          '結果が負か？
                For i% = 1 To e%                                                '補数処理
                    a.Dt(i%) = BASE1 - a.Dt(i%) + k&
                    If a.Dt(i%) > BASE1 Then a.Dt(i%) = a.Dt(i%) - BASE: k& = 1 _
                                        Else k& = 0                             '繰り上げ判定
                Next i%
            End If
    End If
    For i% = e% To 1 Step -1
        If a.Dt(i%) > 0 Then Exit For                                           'データ数a.Dt(0)設定
    Next i%: a.Dt(0) = i%
    If x.Pm = 1 Then a.Pm = 1 - a.Pm                                            '符号設定
    LngAdd = a
End Function

'多倍長整数_減算 x-y
Private Function LngSub(x As Lng, y As Lng) As Lng
Dim yy As Lng
    yy = y: yy.Pm = 1 - y.Pm                                                    'yの符号を反転
    LngSub = LngAdd(x, yy)
End Function

'多倍長整数_乗算 x*y
Private Function LngMul(x As Lng, y As Lng) As Lng
Dim a As Lng
    For i% = 1 To y.Dt(0): For j% = 1 To x.Dt(0)
        k% = i% + j%                                                            'k%:計算桁位置
        b^ = CLngLng(a.Dt(k% - 1)) + CLngLng(x.Dt(j%)) * CLngLng(y.Dt(i%))      'Long * Long -> LongLong
        a.Dt(k% - 1) = CLng(b^ Mod BASE)                                        '下位桁設定
        a.Dt(k) = CLng(CLngLng(a.Dt(k)) + b^ \ BASE)                            '上位桁設定
    Next j%, i%
    a.Pm = IIf(x.Pm <> y.Pm, 1, 0)                                              '符号設定
    For i% = k% + 1 To 1 Step -1                                                'データ数設定
        If a.Dt(i%) > 0 Then a.Dt(0) = i%: Exit For
    Next
    LngMul = a
End Function

'多倍長整数_除算 x/y         ※筆算と同じ方法
Private Function LngDiv(x As Lng, y As Lng) As Lng
Dim a As Lng, x1 As Lng, xs As Lng
    a.Pm = IIf(x.Pm <> y.Pm, 1, 0)                                              '返却値の符号設定
    For k% = x.Dt(0) - 1 To 0 Step -1
        x1 = LngSiftR(x, k%)                                                    '被除数xを右にk%桁シフト
        If LngSub(x1, y).Pm = 0 Then Exit For                                   '被除数xが除数yより大きくなるまで繰返し
    Next k%: k% = k% + 1
    a.Dt(0) = k%                                                                '返却値のデータ数設定
    If y.Dt(0) > 30 Then n% = y.Dt(0) - 20 Else n% = 0                          'Doble化オーバーフロー対策
    y0# = LngDbl(y, n%)                                                         '除数yの浮動小数点化
    For i% = k% To 1 Step -1
        x0# = LngDbl(x1, n%): xs = x1                                           '被除数xの浮動小数点化
        a.Dt(i%) = Int(x0# / y0#)                                               '返却値のi%桁目のデータ設定(近似値)
        x1 = LngSub(x1, LngMul(LngLongSet(a.Dt(i%)), y))                        '余り
        If x1.Pm = 1 Then _
            a.Dt(i%) = a.Dt(i%) - 1: _
            x1 = LngSub(xs, LngMul(LngLongSet(a.Dt(i%)), y))                    '余りがマイナスで再計算
        If LngSub(x1, y).Pm = 0 Then _
            a.Dt(i%) = a.Dt(i%) + 1: _
            x1 = LngSub(xs, LngMul(LngLongSet(a.Dt(i%)), y))                    '余りが除数より大きく再計算
        If i% = 1 Then Exit For                                                 '終了判定(x.Dt(0)破壊防止)
        x1 = LngSiftL(x1, x.Dt(i% - 1))                                         '被除数xを左にシフト
    Next i%
    MdSv = x1                                                                   '余りを保存
    LngDiv = a                                                                  '商を返却
End Function

'n桁右シフト
Private Function LngSiftR(x As Lng, n As Integer) As Lng
Dim a As Lng
    If n < 1 Then LngSiftR = x: Exit Function
    For i% = n + 1 To x.Dt(0)
        a.Dt(i% - n) = x.Dt(i%)
    Next i%
    a.Dt(0) = x.Dt(0) - n
    LngSiftR = a
End Function

'左シフト & 下位桁にd挿入
Private Function LngSiftL(x As Lng, d As Long) As Lng
Dim a As Lng
    If d > BASE1 Then Exit Function
    If d < 0 Then LngSiftL = x: Exit Function
    a.Dt(1) = d
    For i% = 1 To x.Dt(0)
        a.Dt(i% + 1) = x.Dt(i%)
    Next i%
    a.Dt(0) = x.Dt(0) + 1
    LngSiftL = a
End Function

'Long->Lng
Private Function LngLongSet(x As Long) As Lng
Dim a As Lng
    If x > BASE1 Then Exit Function
    If x < 0 Then a.Pm = 1
    a.Dt(1) = x
    a.Dt(0) = 1
    LngLongSet = a
End Function

'Lng->Double        ※近似値
Private Function LngDbl(x As Lng, Optional n As Integer = 0) As Double
Dim a As Lng
    a = LngSiftR(x, n)
    xf% = a.Dt(0):  x0# = a.Dt(xf%)
    If xf% > 1 Then x0# = x0# * BASE + a.Dt(xf% - 1)
    If xf% > 2 Then x0# = x0# * BASE + a.Dt(xf% - 2)
    If xf% > 3 Then x0# = x0# * BASE ^ (xf% - 3)
    If x.Pm = 1 Then x0# = -x0#
    LngDbl = x0#
End Function

'多倍長整数の零チェック
Private Function LngZeroChk(x As Lng) As Boolean
    LngZeroChk = (x.Dt(0) = 1 And x.Dt(1) = 0) Or (x.Dt(0) = 0)                 '零が2種類ある！
End Function

'多倍長整数_合同式 x mod m
Private Function LngMod(x As Lng, m As Lng) As Lng
Dim a As Lng
    a = LngDiv(x, m)                                                            '余りをMdSvに保存
    LngMod = MdSv                                                               '余りを返却
End Function

'多倍長整数_最大公約数 GCD(x,y)
Private Function LngGcd(x As Lng, y As Lng) As Lng
Dim a As Lng, xx As Lng, yy As Lng
    If LngSub(x, y).Pm = 0 Then xx = x: yy = y _
                           Else xx = y: yy = x                                  'xx >= yy
    For i% = 1 To 1000                                                          '無限ループ回避
        a = LngMod(xx, yy)                                                      'Euclidの互除法(再帰呼び出し)
        If LngZeroChk(a) Then Exit For                                          '余りが0まで繰り返す
        xx = yy: yy = a
    Next i%
    LngGcd = yy
End Function

'多倍長整数_最小公倍数 LCM(x,y)
Private Function LngLcm(x As Lng, y As Lng) As Lng
    LngLcm = LngDiv(LngMul(x, y), LngGcd(x, y))                                 'LCM = x*y/GCD(x,y)
End Function

'オイラー関数φ(n)  nと互いに素の数
Private Function LngFai(n As Lng) As Long
    an# = LngDbl(n): s$ = LngSoinsuu(n, False)
    If an# < 2 Then LngFai = 1: Exit Function
    For i% = 1 To 1023
        a1% = InStr(s$, " * ")
        If a1% = 0 Then
                an# = an# * GetP(s$)
                Exit For
        Else:   an# = an# * GetP(Left(s$, a1% - 1))
                s$ = Mid(s$, a1% + 3)
        End If
    Next i%
    LngFai = CLng(an#)
End Function

'素因数取得
Private Function GetP(s As String) As Double
    a% = InStr(s, "^")
    If a% = 0 Then p# = Val(s) _
              Else p# = Val(Left(s, a% - 1))
    GetP = (1 - 1 / p#)
End Function

'多倍長整数_累乗 x^n
Private Function LngPow(x As Lng, n As Long) As Lng
    If n = 0 Then LngPow = LngLongSet(1): Exit Function                         'x ^ 0 = 1
    If n = 1 Then LngPow = x: Exit Function                                     'x ^ 1 = x
    If (n Mod 2) = 0 Then LngPow = LngMul(LngPow(x, n \ 2), _
                                          LngPow(x, n \ 2)) _
                     Else LngPow = LngMul(LngPow(x, n \ 2), _
                                          LngPow(x, n \ 2 + 1))                 '再帰呼び出し
End Function

'多倍長整数_合同累乗 x^p mod m
Private Function LngModPow(x As Lng, p As Lng, m As Lng) As Lng                 'スタックオバー回避の為再帰呼び出し不可
Dim n1 As Lng, n2 As Lng, n3 As Lng, m1 As Lng, c1 As Lng, c2 As Lng
Dim m2(1 To 1000) As Lng                                                        'm2^2の再利用テーブル
    c1 = LngLongSet(1): c2 = LngLongSet(2): m1 = c1                             'm1=1 : m2=x (mod m)
    n1 = LngLongSet(0): n2 = c1: m2(1) = LngMod(x, m): d2% = 1                  'n1=0 : n2=1
    For i& = 1 To 501000                                                        '無限ループ回避 x=2^p → p(p-1)/2 p<1,000
        n3 = LngAdd(n2, n2): d2% = d2% + 1                                      'n3=2*n2
        If LngSub(p, LngAdd(n1, n3)).Pm = 0 Then                                'if p => n1+n3 then
                If LngZeroChk(m2(d2%)) Then _
                    m2(d2%) = LngMod(LngMul(m2(d2% - 1), m2(d2% - 1)), m)       'm2^2の計算を再利用判定
                n2 = n3                                                         'n2=2*n2
        Else:   n1 = LngAdd(n1, n2)                                             'n1=n1+n2/2 : n2=1
                m1 = LngMod(LngMul(m1, m2(d2% - 1)), m): n2 = c1: d2% = 1       'm1=m1*m2   : m2=x (mod m)
        End If
        If LngZeroChk(LngSub(p, n1)) Then Exit For                              'p=n1なら終了
    Next i&: If i& > 501000 Then Error 1
    LngModPow = m1
End Function

'多倍長整数_階乗 n!
Private Function LngFact(n As Long) As Lng
Dim a As Lng
    a = LngLongSet(1)
    For i& = 2 To n
        a = LngMul(a, LngLongSet(i&))
    Next i&
    LngFact = a
End Function

'多倍長整数_順列 nPk=n!/(n-k)!
Private Function LngPermut(n As Long, k As Long) As Lng
Dim a As Lng
    a = LngLongSet(1)
    For i& = n To n - k + 1 Step -1
        a = LngMul(a, LngLongSet(i&))
    Next i&
    LngPermut = a
End Function

'多倍長整数_組合 nCk=nPk/k!=n!/(n-k)!/k!
Private Function LngCombin(n As Long, k As Long) As Lng
    LngCombin = LngDiv(LngPermut(n, k), LngFact(k))
End Function

'多倍長整数_素因数分解
Private Function LngSoinsuu(n As Lng, Optional f As Boolean = True) As String
Dim m As Lng, m1 As Lng, s As Lng, c2 As Lng: c2 = LngLongSet(2)
    If LngSub(n, c2).Pm = 1 Then Exit Function                                  '2未満に素数無し
    m = n: d$ = ""
    For k& = 1 To Nsosuu(0)                                                     '素数リスト分繰返し
        s = LngLongSet(Nsosuu(k&))                                              'k&番目の素数
        For Ss% = 0 To 2000
            m1 = LngDiv(m, s)                                                   'm\素数
            If LngZeroChk(MdSv) Then m = m1 Else Exit For                       'mがSで割り切れなくなるまで繰返し
        Next Ss%
        Select Case Ss%                                                         '同じ素数の数で分岐
            Case 0                                                              '素数無し
            Case 1:     d$ = d$ & " * " & Lng2Str(s, f)                         '同じ素数が1個
            Case Else:  d$ = d$ & " * " & Lng2Str(s, f) & "^" & Ss%             '同じ素数がSs%個
        End Select
        If LngSub(m, LngMul(s, s)).Pm = 1 Then Exit For                         'mが素数^2より大きいと終了
    Next k&
    If LngSub(m, c2).Pm = 0 Then d$ = d$ & " * " & Lng2Str(m, f)                '最後の素数
    LngSoinsuu = Mid(d$, 4)                                                     '先頭の" * "を削除
End Function

'フェルマーテストによる素数検査   m:検査回数
Private Function LngSosuuTest(n As Lng, Optional m As Long = 7) As Boolean
Dim s As Lng, n1 As Lng, p As Lng, c1 As Lng
    c1 = LngLongSet(1): n1 = LngSub(n, c1): SosuuTest = False
    For i& = 1 To m
        s = LngLongSet(Nsosuu(i&))                                              's:i&番目の素数
        If LngSub(s, n1).Pm = 0 Then Exit For                                   's >= n-1 then 素数
        If Not LngZeroChk(LngSub(LngGcd(n, s), c1)) Then Exit Function          'nとsが互いに素出なければ合成数
        If Not LngZeroChk(LngSub(LngModPow(s, n1, n), c1)) Then Exit Function   's^(n-1)<>1 Mod n → 合成数
    Next i&
    LngSosuuTest = True                                                         '→ 素数の可能性が高い
End Function