投稿者: k-sunako

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

• 素数pがベルヌーイ数を割り切るなら非正則素数である | Mathlog
• 非正則素数チェッカー #日曜数学会 | PPT

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

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

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

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

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

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

• bernmm

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

• sage/src/sage/rings/bernmm at develop · sagemath/sage

コンパイルは公式の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番目を確認したが一致していた。

• A027641 – OEIS
• A027642 – OEIS

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;
}

ライプニッツは1689年にイソクロナ・パラケントリカを求めよと書いた。イソクロナ・パラケントリカは最速降下曲線(ブラッキストクロン, brachistochrone)に類似した曲線である。

最速降下曲線は2点間A, Bを結ぶ曲線で、その曲線に沿って質点が重力による影響のみで移動した場合に最短時間でAからBに移動できるという条件を満す曲線である。これはサイクロイドであることが知られている。

イソクロナ・パラケントリカは質点がその曲線に沿って質点が重力による影響のみで移動した場合に開始点から遠ざかる速度と近づく速度が一定になるという条件を満す曲線である。

イソクロナ・パラケントリカの研究からヨハン・ベルヌーイ/ヤコブ・ベルヌーイがレムニスケート曲線を得たという史実がある。

この時以下のような方程式が得られ、レムニスケート曲線に至ったようだ。

\[ (x dx + y dy)\sqrt{y} = (x dy – y dx) \sqrt{a} \]

イソクロナ・パラケントリカの導出から上述の式までを関連づけられないかと考えたため、この文章を書いた。

Paracentric curve を参考に以下の2式を使用する.

\[ \frac{d\rho}{dt} = c t e = v \tag{1}\label{eq:cond1} \]

\[ (\frac{ds}{dt})^{2} = 2 g (y-y_{0}) \tag{2}\label{eq:cond2} \]

数式 \eqref{eq:cond1} はイソクロナ・パラケントリカの特有の条件を表わしている。すなわち、開始点から遠ざかる速度と近づく速度が一定というものである。

数式 \eqref{eq:cond2} はエネルギー保存則を表現しているものである。

数式 \eqref{eq:cond2} を数式 \eqref{eq:cond1} を二乗したもので割ると

\[ (\frac{ds}{d\rho})^{2} = \frac{2g(y-y_{0})}{v^{2}} \tag{3}\label{eq:polar1} \]

極座標の線素を以下のように表示する。

\[ ds = \sqrt{(d\rho)^{2}+(\rho d \theta)^{2}} = \sqrt{1 + \rho^{2}(\frac{d\theta}{d\rho})} d\rho \tag{4}\label{eq:polar2} \]

数式 \eqref{eq:polar1} の左辺を数式 \eqref{eq:polar1} を用いて変形する。

\[ (\frac{ds}{d\rho})^{2} = \frac{1 + \rho^{2}(\frac{d\theta}{d\rho})^{2} d\rho^{2}}{d\rho^{2}} = 1 + \rho^{2} (\frac{d\theta}{d\rho})^{2} \tag{5}\label{eq:polar3} \]

数式 \eqref{eq:polar1} の右辺を \(y_{0}=\frac{-v^{2}}{2g} \Leftrightarrow \frac{2g}{v^{2}}=-\frac{1}{y_{0}}\) を用いて変形する。

\[ \frac{2g(y-y_{0})}{v^{2}} = – \frac{1}{y_{0}}(y-y+0) = 1 – \frac{y}{y_{0}} \tag{6}\label{eq:polar4} \]

数式 \eqref{eq:polar3} と数式 \eqref{eq:polar4} を等値して、変形する。

\begin{align*}
1 + \rho^{2} (\frac{d\theta}{d\rho})^{2} &= 1 – \frac{y}{y_{0}} \\
\rho^{2} (\frac{d\theta}{d\rho})^{2} &= – \frac{y}{y_{0}} \\
\rho^{2} &= – \frac{y}{y_{0}} (\frac{d\rho}{d\theta})^{2} \\
– y_{0} \rho^{2} &= \frac{y}{\rho} (\frac{d\rho}{d\theta})^{2} \\
– y_{0} \rho^{2} &= \sin \theta (\frac{d\rho}{d\theta})^{2} \tag{7}\label{eq:polar5} \\
\end{align*}

極座標の微分方程式を直交座標の微分方程式に変換する。

\[ \frac{dy}{dx} = \frac{\frac{d\rho}{d\theta} \sin \theta + r \sin \theta}{\frac{d\rho}{d\theta} \cos \theta – r \sin \theta} \tag{8}\label{eq:orth1} \]

数式 \eqref{eq:orth1} を \(\frac{d\rho}{d\theta}\) に関する式に直す。

\[ \frac{d\rho}{d\theta} = \sqrt{x^{2}+y^{2}} \frac{\frac{dy}{dx}+x}{\frac{dy}{dx}-y} \tag{9}\label{eq:orth2} \]

数式 \eqref{eq:polar5} を数式 \eqref{eq:orth2} を使って変形してゆく。

\begin{align*}
-y_{0}\sqrt{x^{2}+y^{2}} &= \frac{y}{\sqrt{x^{2}+y^{2}}}(\frac{d\rho}{d\theta})^{2} \\
\frac{-y_{0}(x^{2}+y^{2})}{y} &= (\frac{d\rho}{d\theta})^{2} \\
\sqrt{\frac{-y_{0}(x^{2}+y^{2})}{y}} &= \sqrt{x^{2}+y^{2}} \frac{\frac{dy}{dx}y+x}{\frac{dy}{dx}x-y} \\
\sqrt{\frac{-y_{0}}{y}} &= \frac{\frac{dy}{dx}y+x}{\frac{dy}{dx}x-y} \\
\sqrt{\frac{-y_{0}}{y}} (\frac{dy}{dx}x-y) &= \frac{dy}{dx}y+x \\
\sqrt{-y_{0}} (\frac{dy}{dx}x-y) &= \sqrt{y} \frac{dy}{dx}y+x \\
\sqrt{-y_{0}} (x dy – y dx) &= \sqrt{y} (y dy + x dx) \\
\sqrt{y} (y dy + x dx) &= \sqrt{-y_{0}} (x dy – y dx) \\
(x dx + y dy) \sqrt{y} &= (x dy – y dx) \sqrt{-y_{0}}
\end{align*}

\(a=-y_{0}\) とすれば、冒頭の式が得られる。