8n+1型の素数を y^2-2z^2 という形であらわす

概要

タイトルのとおりだが、素数を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
ed9ktsxznLcGi.png

参考文献

森 正武, and 杉原 正顯. 2003. 複素関数論. 単行本. 岩波書店. https://lead.to/amazon/jp/?op=bt&la=ja&key=4000059505.

ironcladでハッシュ関数MD5を使ってみる

概要

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 []

参考文献

数学セミナー編集部, ed. 2024. 数学セミナー2024年4月号 通巻750号【特集】数学とのつきあい方2024. Kindle版. 日本評論社. https://lead.to/amazon/jp/?op=bt&la=ja&key=B0CWDVX4KN.

GMP_ECMによる素因数分解

概要

(岩波書店『科学』編集部 2024) にあった”RSAと素因数分解”に GMP-ECM が紹介されていたので、実際に動かしてみた。

GMP-ECMの準備

ZIMMERMANN Paul / ecm · GitLab を利用した。 

git clone https://gitlab.inria.fr/zimmerma/ecm.git

INSTALL.dev にある記載のように実行して、コンパイルした。 

cd ecm
autoreconf -i
./configure
make -j4
make check

コンパイル後に実行ファイル ecm がフォルダに生成されている。

素数の生成

2つの素数を生成して、素因数分解のための合成数を用意する。 

import sympy

N = 100
p1 = sympy.randprime(2**(N+1), 2**(N+2))
p2 = sympy.randprime(2**(N+1), 2**(N+2))

print(f'{p1=}')
print(f'{p2=}')
print(f'{p1*p2=}')

p1=2869069113607969685726929882709
p2=3871722605512470843352197003061 
p1*p2=11108239743933603608670116823884965417896735085597594043972249

2つの素数p1, p2から合成数: 11108239743933603608670116823884965417896735085597594043972249 を得た。

実行

GMP-ECM の使い方 の USAGE を参考にコマンドのオプションの引数を決定した。未知の素因数の桁数が30桁なので表から以下の値を採用した。

  • -c(反復): 500
  • B1: 250000

結果の最後部をみると、合成数が元の2つの素数に分解されていることが分る。 

echo 11108239743933603608670116823884965417896735085597594043972249 | ./ecm -c 500 250000


GMP-ECM 7.0.6-dev [configured with GMP 6.3.0, --enable-asm-redc, --enable-assert] [ECM]
Input number is 11108239743933603608670116823884965417896735085597594043972249 (62 digits)
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:4010257250
Step 1 took 77ms
Step 2 took 70ms
Run 2 out of 500:
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:2451171661
Step 1 took 70ms
Step 2 took 69ms
Run 3 out of 500:
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:2557109986
Step 1 took 70ms
Step 2 took 69ms

...

Run 18 out of 500:
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:504643144
Step 1 took 70ms
Step 2 took 70ms
Run 19 out of 500:
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:2975852467
Step 1 took 70ms
Step 2 took 71ms
********** Factor found in step 2: 2869069113607969685726929882709
Found prime factor of 31 digits: 2869069113607969685726929882709
Prime cofactor 3871722605512470843352197003061 has 31 digits

参考文献

岩波書店『科学』編集部, ed. 2024. 科学2024年3月号［雑誌］. Kindle版. 岩波書店. https://lead.to/amazon/jp/?op=bt&la=ja&key=B0CW1DZS1W.

ラゲール法による多項式の求解

概要

(長田 2024) の第4章代数方程式:1ラゲール法をもとに実装する。以下のサイトも参考にした。

  1. ラゲール法による求根アルゴリズム～Pythonコード付き～ – LASCODE

方針

根 1, 2, 3, 4, 5 を持つ5次の多項式を解くことにする. 

import numpy as np

p5deg = np.polynomial.polynomial.Polynomial.fromroots([1, 2, 3, 4, 5])
print(p5deg)

-120.0 + 274.0·x - 225.0·x² + 85.0·x³ - 15.0·x⁴ + 1.0·x⁵

実装

def laguerre(f, x0, delta=0.01, max_iter=100):
    ans = []

    f = f.copy()
    while 1 < f.degree():
        n = f.degree()
        df = f.deriv()
        ddf = df.deriv()

        x = x0
        iter = 0
        while not (np.abs(f(x)) < delta):
            numerator = n * f(x)
            denom_plus = df(x) + np.sqrt((n-1)**2 * df(x)**2 - n*(n-1)*f(x)*ddf(x))
            denom_minus = df(x) - np.sqrt((n-1)**2 * df(x)**2 - n*(n-1)*f(x)*ddf(x))

            denominator = denom_plus if np.abs(denom_minus) < np.abs(denom_plus) else denom_minus
            x -= numerator / denominator

            if np.abs(f(x)) < delta:
                break
            if max_iter < iter:
                break
            iter += 1

        ans.append(x)
        f = f // np.polynomial.polynomial.Polynomial.fromroots([x])

    if f.degree() == 1:
        ans.append(-f.coef[0])
    return ans

検証

5つの根の近似値が得られた。 

roots = laguerre(p5deg, x0=0, delta=0.00001)
print(roots)

[1.0000000000000002, 1.9999999999999907, 3.0000000000000107, 4.000000000000017, 4.999999999999981]

参考文献

長田 直樹. 2024. 数値解析 非線形方程式と数値積分. 単行本. 現代数学社. https://lead.to/amazon/jp/?op=bt&la=ja&key=4768706266.

ラグランジュ基底多項式とSympy

以下の多項式をラグランジェ基底多項式と呼ぶ。

\begin{equation}
L_k(x) := \frac{\Pi_{j \neq k} (x-x_j)}{\Pi_{j \neq k} (x_k-x_j)}
\end{equation}

手計算では大変そうなので、sympyを使ってみた。 Polynomials Manipulation Module Reference – SymPy 1.12 documentation を参考にすると、 3次のラグランジェ基底多項式 \(L_3(x)\) は下記のようになる。 

import sympy
from sympy.core import Add, Mul, symbols
from sympy.abc import x, y

n = 3
X = symbols("%s1:%s" % ('x', n+1))
# Y = symbols("%s1:%s" % ('Y', n+1))

numert = Mul(*[x - X[i] for i in range(n)])

i = 2
numer = numert/(x - X[i])
denom = Mul(*[(X[i] - X[j]) for j in range(n) if i != j])
sympy.latex(numer / denom)

\frac{\left(x - x_{1}\right) \left(x - x_{2}\right)}{\left(- x_{1} + x_{3}\right) \left(- x_{2} + x_{3}\right)}

\[ \frac{\left(x – x_{1}\right) \left(x – x_{2}\right)}{\left(- x_{1} + x_{3}\right) \left(- x_{2} + x_{3}\right)} \]

matplotlibによるアニメーション

下記を参考にmatplotlibのアニメーションを試してみた。

  1. matplotlib.animationでグラフをアニメーションさせる – Qiita
  2. [Python/matplotlib] FuncAnimationを理解して使う – Qiita

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig = plt.figure()

MAX_FRAME = 100

def plot(frame):
    """."""
    plt.cla()
    xs1 = np.arange(0, 10, .01)
    ys1 = np.sin(xs1)
    xs2 = np.arange(0, 10.0 / MAX_FRAME * frame, .01)
    ys2 = np.sin(xs2)
    _ = plt.scatter(xs1, ys1)
    _ = plt.scatter(xs2, ys2)

ani = animation.FuncAnimation(fig, plot, frames=range(MAX_FRAME), interval=100)
ani.save(f, writer="imagemagick")
plt.close()
f
matplotlib_anim.gif

udev

サスペンドからの復帰可能デバイスの設定

下記の2つのページを参考に設定する。

  1. 20.04 – Mouse movement wakes system up from sleep if the USB adapter is not reconnected – Ask Ubuntu
  2. 復帰トリガー – ArchWiki

具体的にはマウスによるサスペンドからの復帰を阻止したい。まず、マウスのデバイス情報を lsusb で取得する。 

lsusb


Bus 002 Device 001: ID 1d6b:0003 Linux Foundation 3.0 root hub
Bus 001 Device 007: ID 25a7:fa10 Areson Technology Corp 2.4G Wireless Receiver
Bus 001 Device 004: ID 2be8:0002 ARCHISS PTR87 ARCHISS PTR87
Bus 001 Device 002: ID 05e3:0608 Genesys Logic, Inc. Hub
Bus 001 Device 006: ID 0b05:19af ASUSTek Computer, Inc. AURA LED Controller
Bus 001 Device 001: ID 1d6b:0002 Linux Foundation 2.0 root hub

25a7:fa10 がマウスに対応する情報である。ワイヤレス・マウスであるため、レシーバがUSB機器として認識されている。 25a7 がベンダーIDで fa10 が製品IDとなる。

/etc/udev/rules.d に wireless_mouse.rules を作成する。 

ls /etc/udev/rules.d


keyboard.rules
wireless_mouse.rules

wireless_mouse.rules に下記を記載する。 idVendor にベンダーID、 idProduct に製品IDを割り当てる。 ATTR{power/wakeup} がサスペンド時に考慮される機器の設定でこのマウスでは復帰させたくないため diabled を設定する。 

ACTION=="add", SUBSYSTEM=="usb", DRIVERS=="usb", ATTRS{idVendor}=="25a7", ATTRS{idProduct}=="fa10", ATTR{power/wakeup}="disabled"

PCを終了し、起動すれば設定が反映されサスペンド時にマウスでの復帰は行なわれなくなる。