Emacs + 暗号 + 数学 + プログラム
Avyならなんでもできる | Emacs JP を読んだ。Avy can do anything | Karthinks の翻訳となる。前に一回読もうとした記事だったので和訳は大変ありがたい。
さっそく以下のようによう設定してみたが、アクションが動作しなかった。
(global-set-key (kbd "M-j") 'avy-goto-char-timer)
The magic kernel には画像サイズを変更する際に使用されるアルゴリズムを紹介している記事である。
Gleam is Pragmatic | software is fun を読んだ。
use を備えているCursed Rust — binarycat を読んだ。Rustの奇妙な箇所を指摘しているらしいが、Rust勉強中の身なのでなにが奇妙なのか分かりかねた。
正則素数を調べていると、ベルヌーイ数により正則素数を判定する記事を幾つか見つけた。
特に 非正則素数チェッカー #日曜数学会 | PPT では論文 (Harvey 2008) による高速なベルヌーイ数の計算方法が紹介されていた。
この論文のアルゴリズムを実際に実装してみたいと考えたため、C言語で実装されたアルゴリズムを論文と比較して理解に努めた。
(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; }
ライプニッツは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}\) とすれば、冒頭の式が得られる。
タイトルのとおりだが、素数を2次形式で表現する手法がある。これをスクリプトで求めてみる。
まず1000以下の素数を収集する。
import sympy primes = [] for p in sympy.primerange(1000): if p % 8 == 1: primes.append(p) print(primes)
[17, 41, 73, 89, 97, 113, 137, 193, 233, 241, 257, 281, 313, 337, 353, 401, 409, 433, 449, 457, 521, 569, 577, 593, 601, 617, 641, 673, 761, 769, 809, 857, 881, 929, 937, 953, 977]
前記であつめた素数を走査して、その素数以下のzを走査して \(y^2 = p+2z^2\) を求め、完全平方数かを判定する。
for p in primes: for z in range(p): y2 = p + 2 * z * z if y2 <= 0: break y, flg = sympy.integer_nthroot(y2, 2) if flg: assert p == y * y - 2 * z * z print(f"{p:3} = {y:2}^2 - 2 * {z:2}^2") # 複数あるが、最初のひとつ以外は除外する break
17 = 5^2 - 2 * 2^2 41 = 7^2 - 2 * 2^2 73 = 9^2 - 2 * 2^2 89 = 11^2 - 2 * 4^2 97 = 13^2 - 2 * 6^2 113 = 11^2 - 2 * 2^2 137 = 13^2 - 2 * 4^2 193 = 15^2 - 2 * 4^2 233 = 19^2 - 2 * 8^2 241 = 21^2 - 2 * 10^2 257 = 17^2 - 2 * 4^2 281 = 17^2 - 2 * 2^2 313 = 21^2 - 2 * 8^2 337 = 25^2 - 2 * 12^2 353 = 19^2 - 2 * 2^2 401 = 23^2 - 2 * 8^2 409 = 21^2 - 2 * 4^2 433 = 21^2 - 2 * 2^2 449 = 29^2 - 2 * 14^2 457 = 23^2 - 2 * 6^2 521 = 23^2 - 2 * 2^2 569 = 31^2 - 2 * 14^2 577 = 33^2 - 2 * 16^2 593 = 25^2 - 2 * 4^2 601 = 27^2 - 2 * 8^2 617 = 25^2 - 2 * 2^2 641 = 29^2 - 2 * 10^2 673 = 31^2 - 2 * 12^2 761 = 31^2 - 2 * 10^2 769 = 29^2 - 2 * 6^2 809 = 29^2 - 2 * 4^2 857 = 37^2 - 2 * 16^2 881 = 41^2 - 2 * 20^2 929 = 31^2 - 2 * 4^2 937 = 35^2 - 2 * 12^2 953 = 31^2 - 2 * 2^2 977 = 37^2 - 2 * 14^2
(森 and 杉原 2003) の 第3章の複素変数の初等関数のプロットを実際にやってみる。
import numpy as np import matplotlib.pyplot as plt def plot(in_curves, func): f, (ax1, ax2) = plt.subplots(1, 2) ax1.axis('equal') ax2.axis('equal') for crv in in_curves: ax1.plot(crv.real, crv.imag) y = func(crv) ax2.plot(y.real, y.imag) in_curves = [] N = 10 for i in range(N): crv = np.linspace(-2, 2, 100) + 1j * (-3*np.pi + i * (3*np.pi-(-3*np.pi)) / N) in_curves.append(crv) for i in range(N): crv = -2 + i * (2 - (-2)) / N + 1j * np.linspace(-3*np.pi, 3*np.pi, 100) in_curves.append(crv) print(in_curves) plt.cla() plot(in_curves, np.exp) plt.savefig(ofile) ofile
sharplispers/ironclad: A cryptographic toolkit written in Common Lisp でハッシュ関数MD5を使ってみた。
quicklisp で ironclad をインストールする。
(ql:quickload "ironclad")
文字列をバイト配列に変換する関数とバイト配列を16進数表記にする関数を定義する。
(defun plain2bytes (message) (sb-ext:string-to-octets message :external-format :utf-8)) (defun bytes2hexstr (bytes) (reduce (lambda (x y) (format nil "~x~x" x y)) bytes :initial-value ""))
ironclad の doc/ironclad.html を参考ハッシュ関数MD5を呼び出す。
(let* ((message "hello!")) (bytes2hexstr (ironclad:digest-sequence :md5 (plain2bytes message))))
5A8DD3AD756A93DED72B823B19DD877
RFC 1321 – The MD5 Message-Digest Algorithm の”A.5 Test suite”を参照して、正しく呼び出せているか確認する。
MD5 (“message digest”) = f96b697d7cb7938d525a2f31aaf161d0
(let* ((message "message digest")) (bytes2hexstr (ironclad:digest-sequence :md5 (plain2bytes message))))
F96B697D7CB7938D525A2F31AAF161D0
(数学セミナー編集部 2024) でミュンヒハウゼン数を知ったので計算してみた。
import itertools import functools @lru_cache(maxsize=1000) def _calc(i, ind): return i * (10**ind) def Munchhausen(n_digits=4): ret = [] values = [(i, i**i) for i in range(1, 10)] for x in itertools.product(*[values]*n_digits): a = 0 b = 0 for ind, (i, iv) in enumerate(x): a += _calc(i, ind) b += iv if a == b: ret.append(a) return ret for i in range(1, 7): print(i, Munchhausen(i))
1 [1] 2 [] 3 [] 4 [3435] 5 [] 6 []