フェルマーテスト(確率的素数判定アルゴリズム)｜Excel VBAで学ぶ数学とアルゴリズム

今回は確率的素数判定アルゴリズムの「フェルマーテスト」について解説していきます。
フェルマーというと『フェルマーの最終定理』で聞いたことがあるという人も多いのではないでしょうか。最終定理という名前で有名なこの定理ですが別名『フェルマーの大定理』とも呼ばれています。

そしてこの大定理と区別するためにあえて”小”と付けられた『フェルマーの小定理』という素数の性質に関した定理も存在しますが、フェルマーテストはこの『フェルマーの小定理』の考えを利用した素数判定アルゴリズムです。(※大定理と小定理は全く持って別のものです)

本ページではこの「フェルマーテスト」のアルゴリズムがどういった仕組みなのかを解説し、最終的にExcel VBAのコードとして実装していきます。

1. フェルマーテストとは
2. サンプルコード
3. まとめ

フェルマーテストとは

フェルマーテストとは入力された値が素数であるかを判定する、素数判定アルゴリズムです。素数判定アルゴリズムは大きく分けて「決定的素数判定」と「確率的素数判定」の2種類ありますが、このフェルマーテストは『フェルマーの小定理』に基づいた、”確率的”素数判定に属されます。
詳しい説明は後にして、まずはフェルマーテストのアルゴリズムを紹介します。

アルゴリズム自体は非常に単純ですが何故このアルゴリズムでできるのかは少しだけ勉強が必要です。

① パラメータとして、2以上 $n$ 未満の自然数 $a$ を決める
② $a$ と $n$ が「互いに素」(2つの数の最大公約数が1)の場合、 $n$ は合成数である ⇒ 終了
③「$a^{n-1} ÷ n$」の余りが「1」でない場合、 $n$ は合成数である ⇒ 終了
④「$a^{n-1} ÷ n$」の余りが「1」である場合、 $n$ は素数である”可能性”がある
⑤ $a$ の値を変更して②～④の判定を指定した回数繰り返す
⑥ 指定回数繰り返してもなお合成数判定が出なければ素数と判定する

フェルマーテストは素数判定のアルゴリズムなので、上記の通り「合成数」とわかった時点で処理は終わります。ただこのフェルマーテスト、合成数の判定は確証を持ってできますが素数の判定は100％の確証がありません。これが「確率的素数判定」といわれる所以です。(※詳しくは以降で解説します)

『フェルマーの小定理』とは

フェルマーテストのアルゴリズムを理解するには『フェルマーの小定理』を理解する必要があります。
フェルマーの小定理は下記の通りの内容で、素数の性質についての定理です。

―『フェルマーの小定理』とは ―

$$a^{p-1} \equiv 1　(mod \ p)$$

$p$ が素数で $a$ と $p$ が「互いに素」の自然数の場合、上記の式が成立する。
式の意味は「$a^{p-1}$ を $p$ で割った時の余りは１である」

「互いに素」とは、2つの数を共に割り切ることのできる整数が「1」のみである、つまりは2つの数の最大公約数が「1」であるもののことをいいます。

たとえば「3」と「10」の場合、共に割り切ることのできる整数は「1」だけなので、この2つの数は互いに素であるといえます。逆に「5」と「10」の場合、共に「5」で割り切ることができるのでこの2つの数は互いに素ではありません。

　
『フェルマーの小定理』の対偶を利用したアルゴリズム

フェルマーテストは上記の『フェルマーの小定理』の対偶を利用して素数を判定するアルゴリズムです。対偶とは「AならばB」を「BでないならばAでない」というように前提と結論を否定にして両者を逆転させたもののことをいいます。

フェルマーの小定理の対偶は下記のような内容になります。

―『フェルマーの小定理』の対偶 ―

$$ a^{n-1} ≢ 1　(mod \ n) $$

$a$ と $n$ が「互いに素」の自然数の場合、上記の式を満たせば $n$ は素数ではない。
式の意味は「$a^{n-1}$ を $n$ で割った時の余りは１ではない」

このようにフェルマーの小定理の対偶を考えることで $n$ の値が”素数ではない”ということを判定することができます。

ただ、ここで注意しないといけない点は「$a^{n-1} ÷ n$」の余りが「1」であるからといって素数であるとは必ずしもいえないということです。フェルマーの小定理の対偶はあくまでも”素数ではない”ということがわかるだけです。つまり、「$a^{n-1} ÷ n$」の余りが「1」であったとしてもそれは素数であるかもしれないし素数ではないかもしれないという訳です。

ここでフェルマーテストではもう1つの自然数 $a$ の値に注目します。
$a$ の値は $n$ と互いに素となる自然数という縛り以外は何もありません。そこでこの条件を満たす範囲内で $a$ の値をランダムに変化させていきます。このとき、指定の回数(たとえば100回)だけ値を変更したとき常に「$a^{n-1} ÷ n$」の余りが「1」であればそれは素数であると判定しようという確率的な判定を行います。(※ $a$ の値の変更回数が多いほど素数の判定精度はよくなります)

このような確率的に判定された素数のことを「確率的素数」といい、確率的素数を使って素数判定を行うアルゴリズムのことを「確率的素数判定アルゴリズム」といいます。これに対して100%素数を判定するようなアルゴリズムのことは決定的素数判定アルゴリズムといいます。

必ず素数が判定できるなら決定的アルゴリズムだけ使えばいいと思う人も多いと思いますが、決定的素数判定は確率的素数判定に比べてかなりの処理時間を必要とします。そのため、処理の高速な確率的素数判定の判定精度を上げて決定的素数判定の代わりに使われるケースも多いです。

ちなみに、判定精度に関わる $a$ はランダム性を出すために乱数扱いとなりますが、このようにアルゴリズム内で乱数を使うようなもののことを「乱択アルゴリズム」ともいいます。乱択アルゴリズムについては「モンテカルロ法」も合わせて参照下さい。

　　
カーマイケル数(絶対擬素数)

フェルマーテストによって求められた素数は「確率的素数」といいますが、その中でも誤って素数として判定されてしまった合成数のことは「擬素数」といいます。

フェルマーテストは判定する回数を増やせば増やすほど精度を上げて素数を判定することができます。
しかし、ある一定の数に対しては何回判定しても素数と判定されてしまう合成数が存在します。

このような、「$a^{n-1} ÷ n$」の「$a$」の数値をいくら変えても素数と判定されてしまう合成数のことを「カーマイケル数」もしくは「絶対擬素数」といい、無限に存在していることがわかっています。

つまりこのカーマイケル数が存在するとわかった時点で、同時にフェルマーテストでの完全な素数判定はできないということもわかっています。

サンプルコード

フェルマーテストをVBAで書くと下記の通りです。

Sub main()

    Dim p As Long
    p = 137
    
    'フェルマーテスト
    If FermatTest(p) = True Then
        Debug.Print p & "は素数"
    Else
        Debug.Print p & "は合成数"
    End If

End Sub
'------------------------------------------------------------------
'   フェルマーテスト
'       N       ：判定する数値
'       戻り値  ：[True]素数, [False]合成数
'------------------------------------------------------------------
Function FermatTest(N As Long) As Boolean

    Const NUMBER_OF_TRIALS As Long = 100    '試行回数
    
    Dim a As Long
    Dim r As Long
    Dim i As Long
    
    '2以下の数字は指定して返す
    If N = 1 Then
        FermatTest = False
        Exit Function
    ElseIf N = 2 Then
        FermatTest = True
        Exit Function
    End If
    
    '確率的素数判定
    For i = 1 To NUMBER_OF_TRIALS
    
        '2以上n未満の範囲で乱数作成
        a = Int((N - 2) * Rnd + 2)
    
        'nとaが「互いに素」でない場合は合成数
        If EuclideanAlgorithm(N, a) <> 1 Then
            FermatTest = False
            Exit Function
        End If
        
        'フェルマーの小定理の対偶
        r = modBignum.Bignum_ModExp(a, N - 1, N) '※ｵｰﾊﾞｰﾌﾛｰ回避用　[r = a ^(n-1) Mod n]と同じ
        If r <> 1 Then
            FermatTest = False
            Exit Function
        End If
        
    Next i
    
    
    FermatTest = True

End Function
'------------------------------------------------------------------
'   ユークリッドの互除法
'       iNum1   ：入力値①
'       iNum2   ：入力値②
'       戻り値  ：最大公約数
'------------------------------------------------------------------
Function EuclideanAlgorithm(iNum1 As Long, iNum2 As Long) As Long

    Dim a As Long
    Dim b As Long
    Dim r As Long
    
    '数値の大きい方の入力値を[a],小さい方の入力値を[b]に代入
    If (iNum1 > iNum2) Then
        a = iNum1
        b = iNum2
    Else
        a = iNum2
        b = iNum1
    End If
    
    '余りが0になるまで除算を繰り返す
    r = a Mod b
    Do While r <> 0
        a = b
        b = r
        r = a Mod b
    Loop
    
    '最大公約数を返す
    EuclideanAlgorithm = b

End Function

Sub main()

Dim p As Long

p = 137

'フェルマーテスト

If FermatTest(p) = True Then

Debug.Print p & "は素数"

Else

Debug.Print p & "は合成数"

End If

End Sub

'------------------------------------------------------------------

' フェルマーテスト

' N ：判定する数値

' 戻り値：[True]素数, [False]合成数

'------------------------------------------------------------------

Function FermatTest(N As Long) As Boolean

Const NUMBER_OF_TRIALS As Long = 100 '試行回数

Dim a As Long

Dim r As Long

Dim i As Long

'2以下の数字は指定して返す

If N = 1 Then

FermatTest = False

Exit Function

ElseIf N = 2 Then

FermatTest = True

Exit Function

End If

'確率的素数判定

For i = 1 To NUMBER_OF_TRIALS

'2以上n未満の範囲で乱数作成

a = Int((N - 2) * Rnd + 2)

'nとaが「互いに素」でない場合は合成数

If EuclideanAlgorithm(N, a) <> 1 Then

FermatTest = False

Exit Function

End If

'フェルマーの小定理の対偶

r = modBignum.Bignum_ModExp(a, N - 1, N) '※ｵｰﾊﾞｰﾌﾛｰ回避用　[r = a ^(n-1) Mod n]と同じ

If r <> 1 Then

FermatTest = False

Exit Function

End If

Next i

FermatTest = True

End Function

'------------------------------------------------------------------

' ユークリッドの互除法

' iNum1 ：入力値①

' iNum2 ：入力値②

' 戻り値：最大公約数

'------------------------------------------------------------------

Function EuclideanAlgorithm(iNum1 As Long, iNum2 As Long) As Long

Dim a As Long

Dim b As Long

Dim r As Long

'数値の大きい方の入力値を[a],小さい方の入力値を[b]に代入

If (iNum1 > iNum2) Then

a = iNum1

b = iNum2

Else

a = iNum2

b = iNum1

End If

'余りが0になるまで除算を繰り返す

r = a Mod b

Do While r <> 0

a = b

b = r

r = a Mod b

Loop

'最大公約数を返す

EuclideanAlgorithm = b

End Function

本コードでは「互いに素」なのかを確認するため、「ユークリッドの互除法」という2つの数の最大公約数を求めるアルゴリズムを利用しています。

また、フェルマーテストのアルゴリズムは上記の内容で完結していますが、フェルマーの小定理の対偶を確認するコードとして「modBignum.Bignum_ModExp(a, n-1, n)」というよくわからない関数が存在していると思います。この関数は引数として「基数[a],指数[n-1],除数[n]」の順に数値を入れたら「$a^{n-1} ÷ n$」の余りを返す関数ということだけわかっていれば大丈夫です。

なぜこのような関数を挟むかというと、単純にVBAで「a ^ (n-1) Mod n」と書くと入力された数値が2桁になるだけでオーバーフローを起こすためです。そのため、オバーフローを起こさないような仕組みで計算が行えるBignum_ModExp関数を使い「$a^{n-1} ÷ n$ 」の余りを求めています。
もしオーバーフローを起こさないのであれば「a ^ (n-1) Mod n」という書き方で問題ありません。

※上記コードをコピペして実行する際はBignum_ModExp関数のコードも必要になります。
「modBignum」という名称の新規標準モジュールを作成し下項目のコードをコピペしてください。
　

巨大桁数演算 (Bignum_ModExpについて)

VBAにおいてべき乗の計算は非常に厄介です。
というのも「$a^n$」のような簡単な数式に対して「$a$=10」「$n$=10」のような小さい数値を入れたとしてもその計算結果は「10,000,000,000」と非常に膨大な数になりすぐにオーバーフローを起こしてしまうためです。

そんなオーバーフローを回避する方法として”文字列”を使うといった手法があります。
下記リンクのサイトではそんな仕組みを使った巨大桁数の計算のサンプルコードが公開されています。

巨大桁整数演算 - 仕事に役立つエクセル実践問題集

巨大桁整数演算 - 仕事に役立つエクセル実践問題集このページを表示するには、フレームをサポートして ...

https://mt-soft.sakura.ne.jp/kyozai/excel_vba/310_vba_chu/30...

　
下記のコードは全て上記リンクのサイト様より拝借したものです。
よってコードの詳細はリンク先を参照してください。
非常に長いコードですがフェルマーテストのアルゴリズムだけ理解できればいいという場合は、このコードで何をしているのかを理解する必要は一切ありません。

Option Explicit

'********************************************************************************************
'
'   巨大桁数演算 関数
'
'   標準モジュール「modBignum」
'   出典：https://mt-soft.sakura.ne.jp/kyozai/excel_vba/310_vba_chu/30_bignum.htm#bignum_exp
'
'********************************************************************************************

'--------------------------------------------------------------
' 巨大桁数の演算
'--------------------------------------------------------------
Public Function Bignum_Sgn(ByVal 数値 As String) As Integer
'
' 符号判定
'
' 　-1:負 0:０ 1:正
'
    数値 = ZeroTrim(数値)
    Select Case Val(Left(数値, 2))
        Case Is < 0:    Bignum_Sgn = -1
        Case 0:         Bignum_Sgn = 0
        Case Else:      Bignum_Sgn = 1
    End Select

End Function

Public Function Bignum_Comp(ByVal 数値1 As String, ByVal 数値2 As String) As Integer
'
' ２つの数値を比較する
'
' 数値1と数値2を比較 数値1>=数値2 → True
'
' 　-1:数値1<数値2 0:数値1=数値2 1:数値1>数値2
'
    Dim i As Long, L As Long
    Dim Sgn1 As Integer, Sgn2 As Integer '符号

    数値1 = ZeroTrim(数値1)
    数値2 = ZeroTrim(数値2)

    Sgn1 = Bignum_Sgn(数値1)    '符号
    Sgn2 = Bignum_Sgn(数値2)

    Bignum_Comp = 0

    If Sgn1 > Sgn2 Then
        Bignum_Comp = 1  '数値1 > 数値2
    ElseIf Sgn1 < Sgn2 Then
        Bignum_Comp = -1 '数値1 < 数値2
    ElseIf Sgn1 = 0 And Sgn2 = 0 Then
        Bignum_Comp = 0  '数値1 = 数値2 = 0

    End If

    If Bignum_Comp <> 0 Then Exit Function  '比較決定済み

    '------------------------------------------------
    ' 符号が同じ場合 Sgn1=1,Sgn2=1 or Sgn1=-1,Sgn2=-1
    '------------------------------------------------
    数値1 = Bignum_Abs(数値1)   '絶対値
    数値2 = Bignum_Abs(数値2)

    L = Max(Len(数値1), Len(数値2))

    数値1 = String(L - Len(数値1), "0") & 数値1
    数値2 = String(L - Len(数値2), "0") & 数値2


    For i = 1 To L
        Select Case StrComp(Mid(数値1, i, 1), Mid(数値2, i, 1))
            Case 0
            Case 1
                Bignum_Comp = 1 * Sgn1
                Exit For
            Case -1
                Bignum_Comp = -1 * Sgn1
                Exit For
        End Select
    Next i

End Function

Public Function Bignum_Add(ByVal 数値1 As String, ByVal 数値2 As String) As String
'
' 数値文字列同士の足し算
'
' 　返却値：計算結果
' 　　　　　エラー時には空白
'
    Bignum_Add = ""

    If Not (Bignum_IsNumeric(数値1) And Bignum_IsNumeric(数値2)) Then Exit Function '無効文字列

    If Bignum_Sgn(数値1) >= 0 Then
        If Bignum_Sgn(数値2) >= 0 Then
            Bignum_Add = Bignum_Add_S(数値1, 数値2)
        Else
            Bignum_Add = Bignum_Subt(数値1, Bignum_Abs(数値2))
        End If
    Else
        If Bignum_Comp(数値2, 0) >= 0 Then
            Bignum_Add = Bignum_Subt(数値2, 数値1)
        Else
            Bignum_Add = "-" & Bignum_Add_S(Bignum_Abs(数値2), Bignum_Abs(数値1))
        End If
    End If

    If Bignum_Sgn(Bignum_Add) = 0 Then Bignum_Add = "0"

End Function

Private Function Bignum_Add_S(ByVal 数値1 As String, ByVal 数値2 As String) As String
'
' 同、サブ
'

    Dim strResult As String

    Dim i As Integer
    Dim Tmp As Integer
    Dim L As Integer

    Bignum_Add_S = ""

    strResult = ""

    数値1 = ZeroTrim(数値1)
    数値2 = ZeroTrim(数値2)
    L = Max(Len(数値1), Len(数値2))
    数値1 = String(L - Len(数値1), " ") & 数値1
    数値2 = String(L - Len(数値2), " ") & 数値2

    Dim r As Integer    '繰り上がり

    r = 0
    strResult = ""
    L = Len(数値1)

    For i = L To 1 Step -1
        Tmp = Val(Mid(数値1, i, 1)) + Val(Mid(数値2, i, 1)) + r
        If Tmp > 9 Then
            r = Tmp \ 10
            strResult = (Tmp Mod 10) & strResult
        Else
            r = 0
            strResult = Tmp & strResult
        End If
    Next i
    strResult = r & strResult
    Bignum_Add_S = ZeroTrim(strResult)


End Function

Public Function Bignum_Subt(ByVal 数値1 As String, ByVal 数値2 As String) As String
'
' 数値文字列同士の引き算
'
'
    If Bignum_Sgn(数値1) >= 0 Then
        If Bignum_Sgn(数値2) >= 0 Then
            If Bignum_Comp(数値1, 数値2) >= 0 Then
                Bignum_Subt = Bignum_Subt_S(数値1, 数値2)
            Else
                Bignum_Subt = "-" & Bignum_Subt_S(数値2, 数値1)
            End If
        Else
            Bignum_Subt = Bignum_Add(数値1, Bignum_Abs(数値2))
        End If
    Else
        If Bignum_Sgn(数値2) >= 0 Then
            Bignum_Subt = "-" & Bignum_Add(Bignum_Abs(数値2), 数値1)
        Else
            Bignum_Subt = Bignum_Subt_S(Bignum_Abs(数値2), Bignum_Abs(数値1))
        End If
    End If

    If Bignum_Sgn(Bignum_Subt) = 0 Then Bignum_Subt = "0"

End Function

Private Function Bignum_Subt_S(ByVal 数値1 As String, ByVal 数値2 As String) As String
'
' 数値文字列同士の引き算
'
' 　本ルーチンでは、S1>0,S2>0,S1>=S2 が前提となる
'
    Dim strResult As String

    Dim i As Integer
    Dim L As Integer
    Dim tmp1 As Integer, tmp2 As Integer

    strResult = ""


    数値1 = ZeroTrim(数値1)
    数値2 = ZeroTrim(数値2)
    L = Max(Len(数値1), Len(数値2))
    数値1 = String(L - Len(数値1), "0") & 数値1
    数値2 = String(L - Len(数値2), "0") & 数値2

    Dim r As Integer


    For i = Len(数値1) To 1 Step -1
        tmp1 = Val(Mid(数値1, i, 1)) - r
        tmp2 = Val(Mid(数値2, i, 1))
        If tmp1 >= tmp2 Then
            r = 0
            strResult = (tmp1 - tmp2) & strResult
        Else
            r = 1
            strResult = (10 + tmp1 - tmp2) & strResult
        End If
    Next i

    Bignum_Subt_S = ZeroTrim(strResult)


End Function


Public Function Bignum_Mul(ByVal 数値1 As String, ByVal 数値2 As String) As String
'
' 数値文字列の掛け算
'
    Dim strResult As String

    Dim i As Integer
    Dim Tmp As Long
    Dim strTmp As String
    Dim r As Long
    Dim strR As String

    Dim intSgn As Integer    '符号

    '--------------------------
    ' 整形
    '--------------------------
    数値1 = ZeroTrim(数値1)
    数値2 = ZeroTrim(数値2)

    '--------------------------
    ' 符号判定＆絶対値化
    '--------------------------
    intSgn = 1
    If Bignum_Sgn(数値1) < 0 Then
        数値1 = Bignum_Abs(数値1)
        intSgn = -intSgn
    End If
    If Bignum_Sgn(数値2) < 0 Then
        数値2 = Bignum_Abs(数値2)
        intSgn = -intSgn
    End If


    r = 0
    strResult = ""

    For i = Len(数値1) To 1 Step -1


        strTmp = Bignum_StrMuls(数値2, Val(Mid(数値1, i, 1)))
        strTmp = Bignum_Add(strTmp, strR)
        strResult = Right(strTmp, 1) & strResult
        strR = Left(strTmp, Len(strTmp) - 1)

        If strR = "" Then strR = "0"

    Next i

    Bignum_Mul = ZeroTrim(strR & strResult)
    If intSgn = -1 And Bignum_Sgn(Bignum_Mul) <> 0 Then Bignum_Mul = "-" & Bignum_Mul '符号付加

End Function

Private Function Bignum_StrMuls(ByVal strS As String, ByVal M As Integer) As String
'
' 数値文字列に１桁掛ける
'
'
    Dim strResult As String

    Dim i As Integer
    Dim Tmp As Long
    Dim r As Long
    Dim S1 As Long, S2 As Long

    strS = ZeroTrim(strS)

    r = 0
    strResult = ""

    For i = Len(strS) To 1 Step -1
        Tmp = Fix(Mid(strS, i, 1)) * M + r
        r = Tmp \ 10
        strResult = (Tmp Mod 10) & strResult
    Next i
    strResult = ZeroTrim(r & strResult)
    If strResult = "" Then strResult = "0"

    Bignum_StrMuls = strResult

End Function

'Public Function Bignum_Div_Quotient(ByVal 被序数 As String, ByVal 除数 As String)
''
'' 商を求める
''
'    Dim strRemainder As String
'    Bignum_Div_Quotient = Bignum_Div(被序数, 除数, strRemainder)
'End Function
'
Public Function Bignum_Div_Remainder(ByVal 被序数 As String, ByVal 除数 As String)
'
' 余りを求める
'
    Dim rc, strRemainder As String
    rc = Bignum_Div(被序数, 除数, strRemainder)
    If IsError(rc) Then
        Bignum_Div_Remainder = rc
    Else
        Bignum_Div_Remainder = strRemainder
    End If
End Function

Public Function Bignum_Div(ByVal strDivident As String, ByVal strDivisor As String, ByRef strRemainder As String)

'
' 割り算
'
'
'strDivident: 被除数 strDivisor: 除数 strQuotient: 商 strRemainder: 余り

    Dim strQuotient As String

    If Bignum_Abs(strDivisor) = "0" Then
        Bignum_Div = CVErr(xlErrDiv0)   ' #DIV/0!
        Exit Function
    End If

    strQuotient = Bignum_Div_S(Bignum_Abs(strDivident), Bignum_Abs(strDivisor), strRemainder)

    If Bignum_Sgn(strRemainder) = 0 Then    '余りが０

        If Bignum_Sgn(strDivident) * Bignum_Sgn(strDivisor) < 0 Then
            strQuotient = "-" & strQuotient
        End If

    Else

        If Bignum_Sgn(strDivident) < 0 And Bignum_Sgn(strDivisor) < 0 Then  '共に負
            strRemainder = "-" & strRemainder
        ElseIf Bignum_Sgn(strDivident) * Bignum_Sgn(strDivisor) < 0 Then    '片方負
            strQuotient = "-" & Bignum_Add(strQuotient, 1)
            strRemainder = Bignum_Subt(strDivident, Bignum_Mul(strQuotient, strDivisor))
        End If

    End If

    Bignum_Div = strQuotient

    If Bignum_Sgn(strRemainder) = 0 Then strRemainder = "0"
    If Bignum_Sgn(Bignum_Div) = 0 Then Bignum_Div = "0"

End Function

Private Function Bignum_Div_S(ByVal strDivident As String, ByVal strDivisor As String, ByRef strRemainder As String)
'
' 割り算　前提：strDivident>=0　　strDivisor>0
'
' strDivident:被除数 strDivisor:除数 strQuotient:商 strRemainder:余り
'
    Dim strQuotient As String

    Dim i As Integer, k As Integer
    Dim S As Integer, L As Integer
    Dim KK As Integer
    Dim strTmp As String
    Dim strTmpS As String

    Dim strF As String    '符号

    strDivident = ZeroTrim(strDivident) '被除数
    strDivisor = ZeroTrim(strDivisor)   '除数
    strQuotient = ""                    '商

    strF = ""                           '符号
    If Bignum_Sgn(strDivident) * Bignum_Sgn(strDivisor) = -1 Then strF = "-"

    L = Len(strDivident)


    For i = 1 To L

        strTmp = Left(strDivident, i)
        k = 0

        Do While Bignum_Comp(strTmp, strDivisor) >= 0
            k = k + 1
            strTmp = Bignum_Subt(strTmp, strDivisor)
        Loop

        If k > 0 Then
            strDivident = strTmp & Mid(strDivident, i + 1, L)
            strDivident = String(L - Len(strDivident), "0") & strDivident
        End If

        strQuotient = Bignum_Add(strQuotient, k & String(L - i, "0"))

        If Bignum_Comp(strDivident, strDivisor) < 0 Then Exit For

    Next i

    strRemainder = ZeroTrim(strDivident)

    Bignum_Div_S = ZeroTrim(strQuotient)


End Function

Function Bignum_ModExp(基数, 指数, 除数)
'
' 基数を指数乗した数値を除数で割った余りを求める
'
    Dim strResult As String
    Dim strB As String, strE As String, strM As String

    Dim strQ As String, strR As String
    Dim strTmp As String

    strB = 基数
    strE = 指数
    strM = 除数

    If 指数 = 0 Then
        Bignum_ModExp = "1"
        Exit Function
    End If

    strQ = Bignum_Div(指数, 2, strR)
    strResult = Bignum_ModExp(基数, strQ, 除数)
    strResult = Bignum_Mul(strResult, strResult)
    strResult = Bignum_Div_Remainder(strResult, strM)

    If Bignum_Div_Remainder(strE, "2") = "1" Then
        strResult = Bignum_Mul(strResult, strB)
        strResult = Bignum_Div_Remainder(strResult, strM)
    End If

    Bignum_ModExp = strResult
    
End Function
'----------------------------------------------------------
' 巨大桁数数値　判定
'----------------------------------------------------------

Public Function Bignum_IsNumeric(ByVal strNo As String) As Boolean
'
' 文字列の有効性判定 0～9 または　-0～9 の文字列
'
' 　返却値　True:有効 False:無効
'
    Dim i As Integer

    Bignum_IsNumeric = True

    strNo = ZeroTrim(strNo)

    For i = 1 To Len(strNo)
        Select Case Mid(strNo, i, 1)
            Case "0" To "9", "-"
            Case Else:  Bignum_IsNumeric = False
        End Select
    Next i

End Function

'-------------------------------------------------------
' 巨大桁数数値　整形
'-------------------------------------------------------

Public Function ZeroTrim(strSrc As String) As String
'
' 左側から空白、０が現れるまで無視する
'
    Dim i As Integer

    ZeroTrim = "0"
    For i = 1 To Len(strSrc)
        Select Case Mid(strSrc, i, 1)
            Case " ", "0"
            Case Else
                ZeroTrim = Trim(Mid(strSrc, i, Len(strSrc)))
                Exit Function
        End Select
    Next i

End Function

Public Function Bignum_Abs(strNo As String) As String
'
' 絶対値を返す
'
    Bignum_Abs = strNo
    If Left(strNo, 1) = "-" Then Bignum_Abs = Right(strNo, Len(strNo) - 1)

End Function



'-------------------------------------------------------
' 内部関数
'-------------------------------------------------------
Public Function Min(n1, n2) '大きい方を返す
    Min = n1
    If n1 > n2 Then Min = n2
End Function

Public Function Max(n1, n2) '小さい方を返す
    Max = n1
    If n1 < n2 Then Max = n2
End Function

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

Option Explicit

'********************************************************************************************

' 巨大桁数演算関数

' 標準モジュール「modBignum」

' 出典：https://mt-soft.sakura.ne.jp/kyozai/excel_vba/310_vba_chu/30_bignum.htm#bignum_exp

'********************************************************************************************

'--------------------------------------------------------------

' 巨大桁数の演算

'--------------------------------------------------------------

Public Function Bignum_Sgn(ByVal 数値 As String) As Integer

' 符号判定

' 　-1:負 0:０ 1:正

数値 = ZeroTrim(数値)

Select Case Val(Left(数値, 2))

Case Is < 0: Bignum_Sgn = -1

Case 0: Bignum_Sgn = 0

Case Else: Bignum_Sgn = 1

End Select

End Function

Public Function Bignum_Comp(ByVal 数値1 As String, ByVal 数値2 As String) As Integer

' ２つの数値を比較する

' 数値1と数値2を比較数値1>=数値2 → True

' 　-1:数値1<数値2 0:数値1=数値2 1:数値1>数値2

Dim i As Long, L As Long

Dim Sgn1 As Integer, Sgn2 As Integer '符号

数値1 = ZeroTrim(数値1)

数値2 = ZeroTrim(数値2)

Sgn1 = Bignum_Sgn(数値1) '符号

Sgn2 = Bignum_Sgn(数値2)

Bignum_Comp = 0

If Sgn1 > Sgn2 Then

Bignum_Comp = 1 '数値1 > 数値2

ElseIf Sgn1 < Sgn2 Then

Bignum_Comp = -1 '数値1 < 数値2

ElseIf Sgn1 = 0 And Sgn2 = 0 Then

Bignum_Comp = 0 '数値1 = 数値2 = 0

End If

If Bignum_Comp <> 0 Then Exit Function '比較決定済み

'------------------------------------------------

' 符号が同じ場合 Sgn1=1,Sgn2=1 or Sgn1=-1,Sgn2=-1

'------------------------------------------------

数値1 = Bignum_Abs(数値1) '絶対値

数値2 = Bignum_Abs(数値2)

L = Max(Len(数値1), Len(数値2))

数値1 = String(L - Len(数値1), "0") & 数値1

数値2 = String(L - Len(数値2), "0") & 数値2

For i = 1 To L

Select Case StrComp(Mid(数値1, i, 1), Mid(数値2, i, 1))

Case 0

Case 1

Bignum_Comp = 1 * Sgn1

Exit For

Case -1

Bignum_Comp = -1 * Sgn1

Exit For

End Select

Next i

End Function

Public Function Bignum_Add(ByVal 数値1 As String, ByVal 数値2 As String) As String

' 数値文字列同士の足し算

' 　返却値：計算結果

' 　　　　　エラー時には空白

Bignum_Add = ""

If Not (Bignum_IsNumeric(数値1) And Bignum_IsNumeric(数値2)) Then Exit Function '無効文字列

If Bignum_Sgn(数値1) >= 0 Then

If Bignum_Sgn(数値2) >= 0 Then

Bignum_Add = Bignum_Add_S(数値1, 数値2)

Else

Bignum_Add = Bignum_Subt(数値1, Bignum_Abs(数値2))

End If

Else

If Bignum_Comp(数値2, 0) >= 0 Then

Bignum_Add = Bignum_Subt(数値2, 数値1)

Else

Bignum_Add = "-" & Bignum_Add_S(Bignum_Abs(数値2), Bignum_Abs(数値1))

End If

If Bignum_Sgn(Bignum_Add) = 0 Then Bignum_Add = "0"

End Function

Private Function Bignum_Add_S(ByVal 数値1 As String, ByVal 数値2 As String) As String

' 同、サブ

Dim strResult As String

Dim i As Integer

Dim Tmp As Integer

Dim L As Integer

Bignum_Add_S = ""

strResult = ""

数値1 = ZeroTrim(数値1)

数値2 = ZeroTrim(数値2)

L = Max(Len(数値1), Len(数値2))

数値1 = String(L - Len(数値1), " ") & 数値1

数値2 = String(L - Len(数値2), " ") & 数値2

Dim r As Integer '繰り上がり

r = 0

strResult = ""

L = Len(数値1)

For i = L To 1 Step -1

Tmp = Val(Mid(数値1, i, 1)) + Val(Mid(数値2, i, 1)) + r

If Tmp > 9 Then

r = Tmp \ 10

strResult = (Tmp Mod 10) & strResult

Else

r = 0

strResult = Tmp & strResult

End If

Next i

strResult = r & strResult

Bignum_Add_S = ZeroTrim(strResult)

End Function

Public Function Bignum_Subt(ByVal 数値1 As String, ByVal 数値2 As String) As String

' 数値文字列同士の引き算

If Bignum_Sgn(数値1) >= 0 Then

If Bignum_Sgn(数値2) >= 0 Then

If Bignum_Comp(数値1, 数値2) >= 0 Then

Bignum_Subt = Bignum_Subt_S(数値1, 数値2)

Else

Bignum_Subt = "-" & Bignum_Subt_S(数値2, 数値1)

End If

Else

Bignum_Subt = Bignum_Add(数値1, Bignum_Abs(数値2))

End If

Else

If Bignum_Sgn(数値2) >= 0 Then

Bignum_Subt = "-" & Bignum_Add(Bignum_Abs(数値2), 数値1)

Else

Bignum_Subt = Bignum_Subt_S(Bignum_Abs(数値2), Bignum_Abs(数値1))

End If

If Bignum_Sgn(Bignum_Subt) = 0 Then Bignum_Subt = "0"

End Function

Private Function Bignum_Subt_S(ByVal 数値1 As String, ByVal 数値2 As String) As String

' 数値文字列同士の引き算

' 　本ルーチンでは、S1>0,S2>0,S1>=S2 が前提となる

Dim strResult As String

Dim i As Integer

Dim L As Integer

Dim tmp1 As Integer, tmp2 As Integer

strResult = ""

数値1 = ZeroTrim(数値1)

数値2 = ZeroTrim(数値2)

L = Max(Len(数値1), Len(数値2))

数値1 = String(L - Len(数値1), "0") & 数値1

数値2 = String(L - Len(数値2), "0") & 数値2

Dim r As Integer

For i = Len(数値1) To 1 Step -1

tmp1 = Val(Mid(数値1, i, 1)) - r

tmp2 = Val(Mid(数値2, i, 1))

If tmp1 >= tmp2 Then

r = 0

strResult = (tmp1 - tmp2) & strResult

Else

r = 1

strResult = (10 + tmp1 - tmp2) & strResult

End If

Next i

Bignum_Subt_S = ZeroTrim(strResult)

End Function

Public Function Bignum_Mul(ByVal 数値1 As String, ByVal 数値2 As String) As String

' 数値文字列の掛け算

Dim strResult As String

Dim i As Integer

Dim Tmp As Long

Dim strTmp As String

Dim r As Long

Dim strR As String

Dim intSgn As Integer '符号

'--------------------------

' 整形

'--------------------------

数値1 = ZeroTrim(数値1)

数値2 = ZeroTrim(数値2)

'--------------------------

' 符号判定＆絶対値化

'--------------------------

intSgn = 1

If Bignum_Sgn(数値1) < 0 Then

数値1 = Bignum_Abs(数値1)

intSgn = -intSgn

End If

If Bignum_Sgn(数値2) < 0 Then

数値2 = Bignum_Abs(数値2)

intSgn = -intSgn

End If

r = 0

strResult = ""

For i = Len(数値1) To 1 Step -1

strTmp = Bignum_StrMuls(数値2, Val(Mid(数値1, i, 1)))

strTmp = Bignum_Add(strTmp, strR)

strResult = Right(strTmp, 1) & strResult

strR = Left(strTmp, Len(strTmp) - 1)

If strR = "" Then strR = "0"

Next i

Bignum_Mul = ZeroTrim(strR & strResult)

If intSgn = -1 And Bignum_Sgn(Bignum_Mul) <> 0 Then Bignum_Mul = "-" & Bignum_Mul '符号付加

End Function

Private Function Bignum_StrMuls(ByVal strS As String, ByVal M As Integer) As String

' 数値文字列に１桁掛ける

Dim strResult As String

Dim i As Integer

Dim Tmp As Long

Dim r As Long

Dim S1 As Long, S2 As Long

strS = ZeroTrim(strS)

r = 0

strResult = ""

For i = Len(strS) To 1 Step -1

Tmp = Fix(Mid(strS, i, 1)) * M + r

r = Tmp \ 10

strResult = (Tmp Mod 10) & strResult

Next i

strResult = ZeroTrim(r & strResult)

If strResult = "" Then strResult = "0"

Bignum_StrMuls = strResult

End Function

'Public Function Bignum_Div_Quotient(ByVal 被序数 As String, ByVal 除数 As String)

'' 商を求める

' Dim strRemainder As String

' Bignum_Div_Quotient = Bignum_Div(被序数, 除数, strRemainder)

'End Function

Public Function Bignum_Div_Remainder(ByVal 被序数 As String, ByVal 除数 As String)

' 余りを求める

Dim rc, strRemainder As String

rc = Bignum_Div(被序数, 除数, strRemainder)

If IsError(rc) Then

Bignum_Div_Remainder = rc

Else

Bignum_Div_Remainder = strRemainder

End If

End Function

Public Function Bignum_Div(ByVal strDivident As String, ByVal strDivisor As String, ByRef strRemainder As String)

' 割り算

'strDivident: 被除数 strDivisor: 除数 strQuotient: 商 strRemainder: 余り

Dim strQuotient As String

If Bignum_Abs(strDivisor) = "0" Then

Bignum_Div = CVErr(xlErrDiv0) ' #DIV/0!

Exit Function

End If

strQuotient = Bignum_Div_S(Bignum_Abs(strDivident), Bignum_Abs(strDivisor), strRemainder)

If Bignum_Sgn(strRemainder) = 0 Then '余りが０

If Bignum_Sgn(strDivident) * Bignum_Sgn(strDivisor) < 0 Then

strQuotient = "-" & strQuotient

End If

Else

If Bignum_Sgn(strDivident) < 0 And Bignum_Sgn(strDivisor) < 0 Then '共に負

strRemainder = "-" & strRemainder

ElseIf Bignum_Sgn(strDivident) * Bignum_Sgn(strDivisor) < 0 Then '片方負

strQuotient = "-" & Bignum_Add(strQuotient, 1)

strRemainder = Bignum_Subt(strDivident, Bignum_Mul(strQuotient, strDivisor))

End If

Bignum_Div = strQuotient

If Bignum_Sgn(strRemainder) = 0 Then strRemainder = "0"

If Bignum_Sgn(Bignum_Div) = 0 Then Bignum_Div = "0"

End Function

Private Function Bignum_Div_S(ByVal strDivident As String, ByVal strDivisor As String, ByRef strRemainder As String)

' 割り算　前提：strDivident>=0　　strDivisor>0

' strDivident:被除数 strDivisor:除数 strQuotient:商 strRemainder:余り

Dim strQuotient As String

Dim i As Integer, k As Integer

Dim S As Integer, L As Integer

Dim KK As Integer

Dim strTmp As String

Dim strTmpS As String

Dim strF As String '符号

strDivident = ZeroTrim(strDivident) '被除数

strDivisor = ZeroTrim(strDivisor) '除数

strQuotient = "" '商

strF = "" '符号

If Bignum_Sgn(strDivident) * Bignum_Sgn(strDivisor) = -1 Then strF = "-"

L = Len(strDivident)

For i = 1 To L

strTmp = Left(strDivident, i)

k = 0

Do While Bignum_Comp(strTmp, strDivisor) >= 0

k = k + 1

strTmp = Bignum_Subt(strTmp, strDivisor)

Loop

If k > 0 Then

strDivident = strTmp & Mid(strDivident, i + 1, L)

strDivident = String(L - Len(strDivident), "0") & strDivident

End If

strQuotient = Bignum_Add(strQuotient, k & String(L - i, "0"))

If Bignum_Comp(strDivident, strDivisor) < 0 Then Exit For

Next i

strRemainder = ZeroTrim(strDivident)

Bignum_Div_S = ZeroTrim(strQuotient)

End Function

Function Bignum_ModExp(基数, 指数, 除数)

' 基数を指数乗した数値を除数で割った余りを求める

Dim strResult As String

Dim strB As String, strE As String, strM As String

Dim strQ As String, strR As String

Dim strTmp As String

strB = 基数

strE = 指数

strM = 除数

If 指数 = 0 Then

Bignum_ModExp = "1"

Exit Function

End If

strQ = Bignum_Div(指数, 2, strR)

strResult = Bignum_ModExp(基数, strQ, 除数)

strResult = Bignum_Mul(strResult, strResult)

strResult = Bignum_Div_Remainder(strResult, strM)

If Bignum_Div_Remainder(strE, "2") = "1" Then

strResult = Bignum_Mul(strResult, strB)

strResult = Bignum_Div_Remainder(strResult, strM)

End If

Bignum_ModExp = strResult

End Function

'----------------------------------------------------------

' 巨大桁数数値　判定

'----------------------------------------------------------

Public Function Bignum_IsNumeric(ByVal strNo As String) As Boolean

' 文字列の有効性判定 0～9 または　-0～9 の文字列

' 　返却値　True:有効 False:無効

Dim i As Integer

Bignum_IsNumeric = True

strNo = ZeroTrim(strNo)

For i = 1 To Len(strNo)

Select Case Mid(strNo, i, 1)

Case "0" To "9", "-"

Case Else: Bignum_IsNumeric = False

End Select

Next i

End Function

'-------------------------------------------------------

' 巨大桁数数値　整形

'-------------------------------------------------------

Public Function ZeroTrim(strSrc As String) As String

' 左側から空白、０が現れるまで無視する

Dim i As Integer

ZeroTrim = "0"

For i = 1 To Len(strSrc)

Select Case Mid(strSrc, i, 1)

Case " ", "0"

Case Else

ZeroTrim = Trim(Mid(strSrc, i, Len(strSrc)))

Exit Function

End Select

Next i

End Function

Public Function Bignum_Abs(strNo As String) As String

' 絶対値を返す

Bignum_Abs = strNo

If Left(strNo, 1) = "-" Then Bignum_Abs = Right(strNo, Len(strNo) - 1)

End Function

'-------------------------------------------------------

' 内部関数

'-------------------------------------------------------

Public Function Min(n1, n2) '大きい方を返す

Min = n1

If n1 > n2 Then Min = n2

End Function

Public Function Max(n1, n2) '小さい方を返す

Max = n1

If n1 < n2 Then Max = n2

End Function

まとめ

今回は素数判定アルゴリズムの「フェルマーテスト」の解説でした。
フェルマーテストのアルゴリズムは下記の通りです。

① パラメータとして、2以上 $n$ 未満の自然数 $a$ を決める。
② $a$ と $n$ が「互いに素」(2つの数の最大公約数が1)の場合、 $n$ は合成数である。 ⇒ 終了
③「$a^{n-1} ÷ n$」の余りが「1」でない場合、 $n$ は合成数である。⇒ 終了
④「$a^{n-1} ÷ n$」の余りが「1」である場合、 $n$ は素数である”可能性”がある。
⑤ $a$ の値を変更して②～④の判定を指定した回数繰り返す。
⑥ 指定回数繰り返してもなお合成数判定が出なければ素数と判定する。

フェルマーテストはいわゆる確率的素数を求めるアルゴリズムなので、運が悪いと合成数でも素数として判定されてしまうことがあります。

このフェルマーテストは素数判定アルゴリズムのはしりであり、以降に改良版の素数判定アルゴリズムとして「ミラーラビン素数判定法」や「AKS素数判定法」といったものが誕生しています。

特に「AKS素数判定法」は”確率的”素数ではなく決定的に素数を判定することが可能になった、世界初の素数判定アルゴリズムといわれています。興味があればぜひこれらのアルゴリズムも勉強してみて下さい。