素数の Lucas テスト - 倭算数理研究所

Scala で『Java によるアルゴリズム事典』のコードを書き換えてみようシリーズ（目次）。　今回はメルセンヌ数（Mersenne Number）が素数であるかどうかを判定するリュカ・テスト (Lucas–Lehmer primality test) を見ていきます。

メルセンヌ数とは ${ p }$ を素数として

　　 ${ \displaystyle\begin{align*} M_p = 2^p - 1 \end{align*}}$

の形の数のことです。　ただし、この記事では ${ p }$ が素数でなくても ${ M_p }$ をメルセンヌ数ということにします。

参考

アルゴリズム

リュカ・テストは、奇素数 ${ p }$ に対して

　　 ${ \displaystyle\begin{align*} \begin{cases} x_0 \;\!\! &= 4 \\ x_{i+1} \!\!&= (x_i^2 -2) \mod M_p \qquad (i \geqq 0) \end{cases} \end{align*}}$

で定義される数列 ${ x_i }$ に対して

　　 ${ \displaystyle\begin{align*} x_{p-2} = 0 &\quad\Longrightarrow\quad M_p\,\textrm{は素数} \\ x_{p-2} \ne 0 &\quad\Longrightarrow\quad M_p\,\textrm{は合成数} \end{align*}}$

が成り立つことを使って、メルセンヌ数 ${ M_p }$ が素数かどうかを判定します。　証明は『wikipedia:リュカ–レーマー・テストの証明』に日本語で書いてくれているので、興味のある方はどうぞ。

補足

リュカ・テストが成り立つ証明は上記のリンク先に載ってますが、その記事では ${ x_i }$ の代わりに以下で定義される ${ S_i }$ を使っています：

　　 ${ \displaystyle\begin{align*} \begin{cases} S_0 \;\!\! &= 4 \\ S_{i+1} \!\!&= S_i^2 -2 \qquad (i \geqq 0) \end{cases} \end{align*}}$

で定義される数列 ${ S_i }$ に対して

　　 ${ \displaystyle\begin{align*} S_{p-2} \equiv 0 \mod M_p &\quad\Longrightarrow\quad M_p\,\textrm{は素数} \\ S_{p-2} \not\equiv 0 \mod M_p &\quad\Longrightarrow\quad M_p\,\textrm{は合成数} \end{align*}}$

${ \mod M_p }$ を漸化式の各評価ごとにとるか、最後にとるかの違いですが、結果には影響しません。　プログラム中では漸化式の途中でとっていかないと数列の各項が大きくなりすぎるので、 ${ x_i }$ の方の定義を用います。

また、リンク先の記事の「必要性」の箇所に

${ p }$ が奇素数であり、かつ ${ Q }$ が素数 ⇒ ${ S_{p-2} \equiv 0 \mod Q }$

( ${ Q = M_p }$ ) とありますが、対偶をとると

${ S_{p-2} \not\equiv 0 \mod Q }$ ⇒ ${ p }$ が奇素数でない、または ${ Q }$ が素数でない

となります。　これを見ると、もし ${ p }$ が2（偶素数）もしくは合成数なら ${ Q }$ が素数でなくてもよさそうに見えます。　しかし、 ${ p }$ が合成数なら ${ p = mn\quad (m,\,n \geqq 2) }$ とおいて

　　 ${ \displaystyle\begin{align*} M_p &= 2^{mn} - 1 \\ &= (2^n-1)\left(2^{(m-1)n} + 2^{(m-2)n} + 2^{(m-3)n} + \cdots + 2^{n} + 1\right) \end{align*}}$

と変形でき、どちらの因子も3以上なので素因数分解できることが分かります。　つまり ${ p }$ が合成数なら ${ Q }$ も合成数となります。　 ${ p }$ が2の場合は ${ S_{p-2} = S_0 = 4 \ne 0 }$ よりテストが通りませんが2は素数なので、この場合はリュカ・テストとは別に素数判定する必要があります。

${ S_i }$ の一般項

リュカ・テストの証明はちょっと込み入ってるようですが、 ${ S_i }$ の一般項はそんなに難しくなく出せるのでちょっとやってみましょう。　数列 ${ \omega_i }$ を

　　 ${ \displaystyle\begin{align*} S_i = \omega_i + \frac{1}{\omega_i} \end{align*}}$

で定義すると、 ${ S_i }$ の漸化式から

　　 ${ \displaystyle\begin{align*} S_{i+1} &= S_i^2 - 2 \\ &= \left(\omega_i + \frac{1}{\omega_i}\right)^2 - 2 \\ &= \omega_i^2 + \frac{1}{\omega_i^2} \end{align*}}$

となります。　これと ${ S_{i+1} = \omega_{i+1} + \dfrac{1}{\omega_{i+1}} }$ を見比べると、 ${ \omega_i }$ の漸化式を

　　 ${ \displaystyle\begin{align*} \omega_{i+1} = \omega_i^2 \end{align*}}$

と定義すればよいことが分かります。　また、初項 ${ \omega_0 }$ は

　　 ${ \displaystyle\begin{align*} &4 = S_0 = \omega_0 + \frac{1}{\omega_0} \\ &\omega_0^2 - 4\omega_0 + 1 = 0 \\ \therefore \, &\omega_0 = 2 \pm \sqrt{3} \end{align*}}$

複号のうちどちらをとっても ${ S_i }$ には影響しないので ${ \omega_0 = 2 + \sqrt{3} }$ としましょう。　このとき

　　 ${ \displaystyle\begin{align*} \omega_i &= (2+\sqrt{3})^{2^i} \\[2mm] \therefore\, S_i &= (2+\sqrt{3})^{2^i} + (2-\sqrt{3})^{2^i} \end{align*}}$

となって、 ${ S_i }$ の一般項が求まりました*1。

実装1

では実装してみましょう。　メルセンヌ数はそんなに大きくない整数 ${ p }$ に対してもすぐに Int や Long の最大値を超えるので（それぞれ 31, 63 まで）BigInt を使うようにしましょう：

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 の実装で何の問題もないので、これ以降は特に読む必要はないかと思います。

メルセンヌ数を法とする代数

メルセンヌ数 ${ M_p = 2^p-1 }$ （ ${ p }$ は素数でなくてもよい）を法とした代数の要素は0から ${ 2^p-2 }$ までの ${ 2^p-1 }$ 個の数で表せるので、 ${ p }$ ビットで表すことができます。　 ${ p }$ ビットでは ${ 2^p }$ 個の数を表せるので、実際には1つ余計で（右下の添字2で2進数を表し、合同 ${ \equiv }$ も = を使うことにして）

　　 ${ \displaystyle\begin{align*} \underbrace{111\cdots1}_{p\textrm{個}}{}_2 = \underbrace{000\cdots0}_{p\textrm{個}}{}_2 = 0 \end{align*}}$

を同一視します。

最大の数 ${ 2^p-2 }$ は

　　 ${ \displaystyle\begin{align*} 2^p-2 = \underbrace{111\cdots1}_{(p-1)\textrm{個}}0_2 \end{align*}}$

となります。

${ M_p }$ を法とした場合、ビット数が ${ p }$ より大きい数を mod （プログラムでは %）によって ${ p }$ ビットに収まる数に変換するのが、比較的簡単なビット演算で計算できます。　上で見たように ${ p }$ 個のビットが全て1ならその数は0となるので、その次の数を考えて

　　 ${ \displaystyle\begin{align*} 1\underbrace{000\cdots0}_{p\textrm{個}}{}_2 &= \underbrace{000\cdots0}_{(p-1)\textrm{個}}{1}_2 = 1 \end{align*}}$

となりますが、これは下位 ${ p }$ 個の0を取り除くことに対応しています。　ビット演算で言えば、符号なしで ${ p }$ だけ右シフトすることに相当します。

負の数（というか和に関する逆元）は単に全ビットを反転すればいいだけです。　実際、全ビットが逆の2数を足せば全ビットが1の数になりますが、これは0に等しいんでしたね。　ただし、0の場合はそのまま0です。　全ビットを反転して全て1になってもそれは0です。

ちょっと具体例

もう少し理解を深めるために、具体的な数字で計算してみましょう。　 ${ p = 3 }$ のとき ${ M_3 = 7 }$ ですね。　このとき、最大の数は ${ 110_2 = 6 }$ 、また ${ 111_2 = 7 \equiv 0 \!\!\mod 7 }$ となります。　 ${ p }$ ビットより大きい数に対しては

　　 ${ \displaystyle\begin{align*} 100110_2 &= 100\underline{000}_2 + 110_2 \\ &= 100_2 + 110_2 \\ &= 1010_2 \\ &= 1\underline{000}_2 + 10_2 \\ &= 1_2 + 10_2 \\ &= 11_2 \end{align*}}$

${ 100110_2 = 38 \equiv 3 \!\!\mod 7 }$ なので確かに合ってますね。　上記の計算では少々助長に途中式を書いてますが、要は右から ${ p }$ 個ごとにビットを区切って足せば OK です：

　　 ${ \displaystyle\begin{align*} 100110_2 &= 100_2 + 110_2 \\ &= 1010_2 \\ &= 1_2 + 10_2 \\ &= 11_2 \end{align*}}$

負の数は、例えば ${ 101_2 (= 5) }$ に対して、全ビットを反転した ${ 010_2 = 2 }$ となります。　実際にこの2数を加えてみると

　　 ${ \displaystyle\begin{align*} 101_2 + 010_2 &= 111_2 \\ &= 000_2 \end{align*}}$

となって、確かに和に関して逆元になっています。

コード

せっかく代数の軽い説明をしたので、その代数の演算が定義されたクラスを作りましょう。　何桁でも扱えるようにするのは非常に大変なので、Int で表せる32ビットまでを扱うことにします（ ${ p \leqq 32 }$ ）。　実質的に非負の整数しか出てこないので、プロパティとしては Int ではなく spire に定義されている UInt を使う方がよさそうです（ ${ p = 32 }$ も扱えるので）。　そのクラスを MersenneModularUInt として各種演算を定義しますが、まずはこのオブジェクトを生成する MersenneModulus クラスを定義しましょう：

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 に名前を変更しました。