月: 2024年8月

ベルヌーイ数のアルゴリズム

概要

正則素数を調べていると、ベルヌーイ数により正則素数を判定する記事を幾つか見つけた。

特に 非正則素数チェッカー #日曜数学会 | PPT では論文 (Harvey 2008) による高速なベルヌーイ数の計算方法が紹介されていた。

この論文のアルゴリズムを実際に実装してみたいと考えたため、C言語で実装されたアルゴリズムを論文と比較して理解に努めた。

アルゴリズムの概要

(Harvey 2008) が提案しているアルゴリムズには以下のような特徴がある。

  • ゼータ関数に依存しない
  • 中国剰余定理 (chinese remainder theorem) を活用する
  • 並列演算が可能である

論文の疑似コード

論文の疑似コードを引用すると以下のようになる。

TtlpkHrFVzo.jpg
Figure 1: Algorithm1: Compute \(B_k\) modulo \(p\); (Harvey 2008) から引用

ソースコード

論文の著者の C++ 実装は以下のベージに存在するが、非常に古くコンパイルもできなかった。

上記のページに SageMath – Open-Source Mathematical Software System に導入されたとあるので、レポジトリを見ると以下の実装を発見した。こちらはコンパイル、実行まで可能であった。

コンパイルは公式のmakefileを参考に以下のようなmakefikeを自作して対応した。 

CFLAGS=-O2 -DNDEBUG -DUSE_THREADS -DTHREAD_STACK_SIZE=4096 -std=c++11

bernmm-test: bernmm-test.cpp bern_modp.cpp bern_modp.h bern_rat.cpp bern_rat.h bern_modp_util.cpp bern_modp_util.h
    g++ $(CFLAGS) -o bernmm-test \
            bernmm-test.cpp bern_modp.cpp bern_rat.cpp bern_modp_util.cpp \
            -lntl -lgmp -lpthread

clean:
    rm bernmm-test

実際にベルヌーイ数の100番目を計算して、オンライン整数列大辞典 でベルヌーイ数の100番目を確認したが一致していた。

Algorithm1 の対応するコードは以下になるようだ。 

/******************************************************************************

   Computing the main sum (general case)

******************************************************************************/

/*
   Returns (1 - g^k) B_k / 2k mod p.

   PRECONDITIONS:
      5 <= p < NTL_SP_BOUND, p prime
      2 <= k <= p-3, k even
      pinv = PrepMulMod(p)
      g = a multiplicative generator of GF(p), in [0, p)
*/
long bernsum_powg(long p, mulmod_t pinv, long k, long g)
{
   long half_gm1 = (g + ((g & 1) ? 0 : p) - 1) / 2;    // (g-1)/2 mod p
   long g_to_jm1 = 1;
   long g_to_km1 = PowerMod(g, k-1, p, pinv);
   long g_to_km1_to_j = g_to_km1;
   long sum = 0;
   muldivrem_t g_pinv = PrepMulDivRem(g, p);
   mulmod_precon_t g_to_km1_pinv = PrepMulModPrecon(g_to_km1, p, pinv);

   for (long j = 1; j <= (p-1)/2; j++)
   {
      // at this point,
      //    g_to_jm1 holds g^(j-1) mod p
      //    g_to_km1_to_j holds (g^(k-1))^j mod p

      // update g_to_jm1 and compute q = (g*(g^(j-1) mod p) - (g^j mod p)) / p
      long q;
      g_to_jm1 = MulDivRem(q, g_to_jm1, g, p, g_pinv);

      // compute h = -h_g(g^j) = q - (g-1)/2
      long h = SubMod(q, half_gm1, p);

      // add h_g(g^j) * (g^(k-1))^j to running total
      sum = SubMod(sum, MulMod(h, g_to_km1_to_j, p, pinv), p);

      // update g_to_km1_to_j
      g_to_km1_to_j = MulModPrecon(g_to_km1_to_j, g_to_km1, p, g_to_km1_pinv);
   }

   return sum;
}

参考文献

Harvey, David. 2008. “A Multimodular Algorithm for Computing Bernoulli Numbers.” https://arxiv.org/abs/0807.1347.