Scala で『Java によるアルゴリズム事典』のコードを書き換えてみようシリーズ(目次)。 今回はメルセンヌ数(Mersenne Number)が素数であるかどうかを判定するリュカ・テスト (Lucas–Lehmer primality test) を見ていきます。
の形の数のことです。 ただし、この記事では が素数でなくても
をメルセンヌ数ということにします。
参考
アルゴリズム
リュカ・テストは、奇素数
で定義される数列 に対して
が成り立つことを使って、メルセンヌ数 が素数かどうかを判定します。 証明は 『wikipedia:リュカ–レーマー・テストの証明』 に日本語で書いてくれているので、興味のある方はどうぞ。
補足
リュカ・テストが成り立つ証明は上記のリンク先に載ってますが、その記事では
で定義される数列 に対して
を漸化式の各評価ごとにとるか、最後にとるかの違いですが、結果には影響しません。 プログラム中では漸化式の途中でとっていかないと数列の各項が大きくなりすぎるので、
の方の定義を用います。
また、リンク先の記事の「必要性」の箇所に
() とありますが、対偶をとると
となります。 これを見ると、もし が2(偶素数)もしくは合成数なら
が素数でなくてもよさそうに見えます。 しかし、
が合成数なら
とおいて
と変形でき、どちらの因子も3以上なので素因数分解できることが分かります。 つまり が合成数なら
も合成数となります。
が2の場合は
よりテストが通りませんが2は素数なので、この場合はリュカ・テストとは別に素数判定する必要があります。
の一般項
リュカ・テストの証明はちょっと込み入ってるようですが、
で定義すると、 の漸化式から
となります。 これと を見比べると、
の漸化式を
と定義すればよいことが分かります。 また、初項 は
複号のうちどちらをとっても には影響しないので
としましょう。 このとき
となって、 の一般項が求まりました*1。
実装1
では実装してみましょう。 メルセンヌ数はそんなに大きくない整数import spire.implicits._ case class MersenneNumber(p: Int){ require(p >= 0) lazy val toBigInt: BigInt = BigInt(2)**p - 1 def isPrime: Boolean = p match { case 0 | 1 => false case 2 => true // p = 2 の場合 case _ => val m = toBigInt val seq = Stream.iterate[BigInt](4)(x => (x * x - 2) % m) seq(p - 2) == 0 } }
- 初項と漸化式が与えられている数列の Seq を生成したい場合、Stream.iterate メソッドを使うのがいいでしょう。
- 数学の mod は、コード中では言わずと知れた「%」で。
まぁ別にクラスにする必要はなかったかも知れませんが。 使い方は以下のようになります:
// p = 3 の場合 val mn3 = MersenneNumber(3) assert( mn3.toBigInt == BigInt(7) ) assert( mn3.isPrime == true ) // p = 11 の場合(メルセンヌ数は素数ではない) val mn11 = MersenneNumber(11) assert( mn11.toBigInt == BigInt(2047) ) assert( mn11.isPrime == false ) // p が 640 まででメルセンヌ数が素数になるもの assert( (0 to 640).filter(MersenneNumber(_).isPrime) == Seq(2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607) )
p = 607 のときのメルセンヌ数ってかなり大きい数になるでしょうけど、このコードは大した時間かからずに処理が帰ってきます。 おぉすごい!
実装2
『Java によるアルゴリズム事典』には、BigInt (BigInteger) を使った実装以外に、配列を使った実装も載ってあるのですが、そのアルゴリズムを Scala の Seq 等で書いてみるとかなりパフォーマンスがクソなことになったので、ここではメルセンヌ数を法とする和・差・積を定義した MersenneModuloUInt というクラスを作って、計算をどのようにするかの簡単な例としたいと思います。ぶっちゃけ、BigInt の実装で何の問題もないので、これ以降は特に読む必要はないかと思います。
メルセンヌ数を法とする代数
メルセンヌ数
を同一視します。
最大の数 は
となります。
を法とした場合、ビット数が
より大きい数を mod (プログラムでは %)によって
ビットに収まる数に変換するのが、比較的簡単なビット演算で計算できます。 上で見たように
個のビットが全て1ならその数は0となるので、その次の数を考えて
となりますが、これは下位 個の0を取り除くことに対応しています。 ビット演算で言えば、符号なしで
だけ右シフトすることに相当します。
負の数(というか和に関する逆元)は単に全ビットを反転すればいいだけです。 実際、全ビットが逆の2数を足せば全ビットが1の数になりますが、これは0に等しいんでしたね。 ただし、0の場合はそのまま0です。 全ビットを反転して全て1になってもそれは0です。
ちょっと具体例
もう少し理解を深めるために、具体的な数字で計算してみましょう。
なので確かに合ってますね。 上記の計算では少々助長に途中式を書いてますが、要は右から
個ごとにビットを区切って足せば OK です:
負の数は、例えば に対して、全ビットを反転した
となります。 実際にこの2数を加えてみると
となって、確かに和に関して逆元になっています。
コード
せっかく代数の軽い説明をしたので、その代数の演算が定義されたクラスを作りましょう。 何桁でも扱えるようにするのは非常に大変なので、Int で表せる32ビットまでを扱うことにします(import spire.math.UInt import spire.implicits._ import scala.annotation.tailrec case class MersenneModulus(p: Int){ require(2 <= p && p <= 32) val mask64 = 2L**p-1L val mask32 = UInt(mask64.toInt) // For p = 3 : // mask = 00...0111 // ~mask = 11...1000 /** UInt から MersenneModularUInt オブジェクトを生成するメソッド */ def fromUInt(n: UInt): MersenneModularUInt = { @tailrec def mod(n: UInt): UInt = n & ~mask32 match { case x if x == UInt(0) => if(n == mask32) UInt(0) else n case x => mod((n & mask32) + (n >>> p)) } MersenneModularUInt(mod(n), this) } /** Int から MersenneModularUInt オブジェクトを生成する apply ファクトリ・メソッド */ def apply(n: Int): MersenneModularUInt = if(n < 0) -fromUInt(UInt(-n)) else fromUInt(UInt(n)) /** Long から MersenneModularUInt オブジェクトを生成する apply ファクトリ・メソッド */ def apply(n: Long): MersenneModularUInt = n match { case _ if n < 0L => -apply(-n) case _ => @tailrec def mod(n: Long): Long = n & ~mask64 match { case 0L => if(n == mask64) 0L else n case x => mod((n & mask64) + (n >>> p)) } MersenneModularUInt(UInt(mod(n)), this) } }
- Int 値や Long 値をとって、後で定義する MersenneModularUInt オブジェクトを生成する apply() ファクトリ・メソッドを定義しています。
- UInt 値をとる apply() ファクトリ・メソッドも定義しようとしたんですが、Int をとるものが定義されているとダメなようなので、fromUInt() メソッドとして定義しています。
- p ビットを超える数や負の数を、ビットシフトを使って p ビットで表せる数に変換しています。
さて、次は各種演算が定義された MersenneModularUInt クラス:
import spire.math.UInt import spire.implicits._ case class MersenneModularUInt(digit: UInt, modulus: MersenneModulus){ def p: Int = modulus.p def isZero: Boolean = this.digit == UInt(0) def isOne : Boolean = this.digit == UInt(1) def unary_- : MersenneModularUInt = if(this.isZero) this else MersenneModularUInt(modulus.mask32 & ~digit, modulus) def +(that: MersenneModularUInt): MersenneModularUInt = if(this.isZero) that else if(that.isZero) this else calculateSum(that) /** 足し算実行 */ private def calculateSum(that: MersenneModularUInt): MersenneModularUInt = this.modulus(this.digit.toLong + that.digit.toLong) def -(that: MersenneModularUInt): MersenneModularUInt = if(this.isZero) -that else if(that.isZero) this else this + (-that) def *(that: MersenneModularUInt): MersenneModularUInt = if(this.isZero | that.isOne) this else if(this.isOne | that.isZero) that else calculateProduct(that) /** 掛け算実行 */ private def calculateProduct(that: MersenneModularUInt): MersenneModularUInt = this.modulus(this.digit.toLong * that.digit.toLong) def toInt: Int = this.digit.toInt def toLong: Long = this.digit.toLong def toBigInt: BigInt = this.digit.toBigInt override def toString: String = s"$digit (mod 2^$p-1)" }
- 演算に0や1が現れた場合には計算せずに済むことが多いので、いろいろ条件分岐していますが、和や積の実質的な計算は calculateXxxx という名前の private メソッドで行っています
- 数値として保持しているのは UInt ですが、和や積の計算途中では結果の桁が溢れると困るので Long に変換してから計算しています。
- p ビット以下への変換は先に定義した MersenneModulus クラスの apply() ファクトリ・メソッドに任せてます。
さて、この MersenneModularUInt クラスを使ってリュカ・テストを実装すると以下のようになります(まぁ、ほとんど BigInt を使った実装と同じですが):
object UIntMersenneNumber{ def isPrime(p: Int): Boolean = p match { case 0 | 1 => false case 2 => true case _ => val mmod = MersenneModulus(p) val v2 = mmod(2) val xs = Stream.iterate(mmod(4))(x => x*x - v2) xs(p-2).isZero } }
この isPrime メソッドを使って、32までのメルセンヌ素数を全て求めてみると
assert( (0 to 32).filter(UIntMersenneNumber.isPrime(_)) == Seq(2, 3, 5, 7, 13, 17, 19, 31) )
となってキチンと答えが出ました。
ここで実装した MersenneModularUInt クラスでは 32 以下の p しか扱えませんが、spire の Natural の実装などを参考にすると任意の p に対して計算できるようにすることもできるかと思います。 ただ、Natural のソースを見ると「パフォーマンスを気にする場合は BigInt を使え」とかコメントされてるので、あんまり作る気が起きてきませんがw
【修正】
- BigInt の実装でのサンプルコード部分が間違っていたので修正しました。
- 0の和に関する逆元はそのままであることを追記しました。
- MersenneModulo を MersenneModulus に、MersenneModuloUInt を MersenneModularUInt に名前を変更しました。

- 作者: Martin Odersky,Lex Spoon,Bill Venners,羽生田栄一,水島宏太,長尾高弘
- 出版社/メーカー: インプレス
- 発売日: 2016/09/20
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (1件) を見る

- 作者: 奥村晴彦,杉浦方紀,津留和生,首藤一幸,土村展之
- 出版社/メーカー: 技術評論社
- 発売日: 2003/05
- メディア: 単行本
- 購入: 2人 クリック: 61回
- この商品を含むブログ (60件) を見る
*1:どちらかと言うと、この一般項がありきで、コードに落とし込むために漸化式に書き換えたのではないかという気がしますがどうなんでしょうかね。 というか、証明では漸化式ではなく一般項を使ってるのでそういうことなんでしょう。