k-sunakoの技術ログ – ページ 5 – Emacs + 暗号 + 数学 + プログラム

“Can a Rubik’s Cube be brute-forced?” を読んだ

Can a Rubik’s Cube be brute-forced? を読んだ。

ルービックキューブをどのように解けるかを解説した記事である
ルービックキューブを解く有名なアルゴリズムは以下がある
- Korf’s algorithm
- Thistlethwaite’s algorithm
ルービックキューブを順列パズルと考え、総当たりを考える
パターンは43京ほどあるが、10億まで減らせる方法がある
順列を表わすデータ構造として、 “Permutaion Trie” がある
作者がCommon Lispで実装したものが stylewarning/cl-permutation: Permutations and permutation groups in Common Lisp. にある

ルービックキューブと群論は (Joyner 2010) で興味をもったことを思いだした。Common Lispも会得しようと考えていたが、いまだ覚つかない。これらを高レベルに融合するのが私の人生の目標の一つであったが、それは果せるだろうか？

参考文献

Joyner, David. 2010. 群論の味わいー置換群で解き明かすルービックキューブと15パズルー. Translated by 川辺治之. 単行本. 共立出版. https://lead.to/amazon/jp/?op=bt&la=ja&key=4320019415.

正弦関数のマクローリン展開の計算

概要

解析学の級数の理解のために、正弦関数のマクローリン展開の数値例を使おうと考えた。

計算

正弦関数（sin関数）の高階微分とテイラー展開（マクローリン展開）| 関数の微分 | 微分積分 | 数学 | ワイズを参考に以下のようにコードにした。

import math
import numpy as np

def calc_sin_Maclaurin(x, n):
    ans = []
    for k in range(n):
        term = (-1)**k * x**(2*k+1) / math.factorial(2*k+1)
        ans.append(term)
    return sum(ans), ans
np.sin(np.sqrt(2)), calc_sin_Maclaurin(np.sqrt(2), 10)[0]

確認のためにプロットしてみる。

import matplotlib.pyplot as plt
import numpy as np

vfunc = np.vectorize(lambda x: calc_sin_Maclaurin(x, 10)[0])

fig, ax = plt.subplots(1, 1)
x = np.linspace(0, 10, 1000)
ax.plot(x, vfunc(x), label="Maclaurin(x, 10)")
ax.plot(x, np.sin(x), label="sin(x)")
ax.set_ylim(-2, 2)
plt.legend()
fig.savefig(outfile)
outfile

これで以下の級数を利用しよう。

calc_sin_Maclaurin(np.sqrt(2), 5)[1]

1.4142135623730951

-0.4714045207910318

0.04714045207910319

-0.0022447834323382474

6.23550953427291e-05

参考文献

“Unsafe Rust Is Harder Than C” を読んだ

Unsafe Rust Is Harder Than C を読んだ。

Rustのunsafeなコードの説明である
侵入型リンクリスト(intrusive linked lists)のようなデータ構造を実装した際の経験を記事にしている
Photohashというソフトウェアを作成するために、効率的な並列処理が必要になった
そのため、独自のチャネルを実装することにしたようだ

ペル方程式を解くコード

連分数の計算を元にペル方程式を計算する。ペル方程式の連分数を用いた魔法の解法 – tsujimotterのノートブックを参考にした。

ペル方程式は以下のように表わされる。特にDが素数の時に興味がある。

\[ x^{2} – Dy^{2} = 1 \]

まず、D=61で計算してみる。計算手順は以下のようになる。

\(\sqrt{61}\) の連分数展開を求める
連分数の周期を確認し、奇数なら2周期とする
近似分数 \(\frac{P}{Q}\) を求める
\(P^{2} – 61Q^{2} = 1\) となる

cfrac = get_continued_fraction(61)
if 1 == (len(cfrac[1:]) % 2):
    cfrac += cfrac[1:]
get_approx_fraction(cfrac[:-1])

つぎに115以下の素数をペル方程式 – Wikipedia にある一覧表で確認する。

import sympy

lst_pq = [['D(prime)', 'P', 'Q'], None]
for i, p in enumerate(sympy.primerange(115)):
    cfrac = get_continued_fraction(p)
    if 1 == (len(cfrac[1:]) % 2):
        cfrac += cfrac[1:]
    pq = get_approx_fraction(cfrac[:-1])
    lst_pq.append([p, pq.numerator, pq.denominator])
lst_pq

D(prime)	P	Q
2	3	2
3	2	1
5	9	4
7	8	3
11	10	3
13	649	180
17	33	8
19	170	39
23	24	5
29	9801	1820
31	1520	273
37	73	12
41	2049	320
43	3482	531
47	48	7
53	66249	9100
59	530	69
61	1766319049	226153980
67	48842	5967
71	3480	413
73	2281249	267000
79	80	9
83	82	9
89	500001	53000
97	62809633	6377352
101	201	20
103	227528	22419
107	962	93
109	158070671986249	15140424455100
113	1204353	113296

EmacsのTypescript環境

久々にTypeScriptを書く必要があったので、環境を構築した。

まず、 nvm で 22.2.0 をインストール。(バージョンは適当。)

nvm install 22.2.0

ボイラープレートを利用して雛形を作成する。

git clone https://github.com/jsynowiec/node-typescript-boilerplate.git
cd node-typescript-boilerplate

nvm で Node.js のバージョンを 22.2.0 に指定する

nvm use 22.2.0

依存パッケージをインストールする。

npm install

以下のようにemacsの設定ファイルに記載する。

;; 拡張子.tsのファイルを開いたら、typescript-ts-modeにする
(add-to-list 'auto-mode-alist '("\\.ts" . typescript-ts-mode))

;; eslintなどはプロジェクトのnode_module内のものを使いたいため
(use-package add-node-modules-path
  :straight t
  :hook ((typescript-ts-mode . add-node-modules-path))
  :config
  ;; https://github.com/codesuki/add-node-modules-path/issues/23 より
  (setq add-node-modules-path-command '("echo \"$(npm root)/.bin\"")))

;; 必要か？
(use-package nvm
  :straight (:host github :repo "rejeep/nvm.el")
  :config
  ;; Optionally set a default node version
  (nvm-use "22.2.0"))

Emacsのコマンドで M-x lsp-install-server を実行して、 ts-ls をインストールする。

これで環境が構築できた。

“My NumPy Year: Creating a DType for the Next Generation of Scientific Computing” を読んだ

My NumPy Year: Creating a DType for the Next Generation of Scientific Computing を読んだ。

Numpy2.xの文字列のDTypeについて記載された記事である
Python3はUnicode文字列を扱うため、Numpy1.xでは文字列の型は dtype='<U5' のようになる
- アスキー文字でも4バイト使うため、メモリ効率が悪い
- 文字列の操作が遅い
Numpy1.xでは文字列の型は dtype=object としても表現できる
- これはpandasの性能に大きな影響があった
- また、astropyという天文学(?)関連のライブラリにも
オブジェクトの配列でなくポインターの配列で実装した？
Numpy2.xで文字列の型 np.dtype.StringDtype が導入された
メモリ管理に「Arena Allocator（アリーナ・アロケーター）」というものも使われている
- [Zig] Arena Allocatorを使った “ふつうの” Zigプログラミング

“Smolderingly fast b-trees”を読んだ

Smolderingly fast b-trees を読んだ。

順序があるデータ構造(B-Tree)とハッシュマップの性能差を述べた記事である
WASMではいくつかの制約が追加される
順序があるデータ構造とハッシュマップの性能差には人々の間で認識に隔たりがある
rust と zig を使って、性能調査を実施している
- ランダムな整数と文字列でのアクセス、順序のある整数と文字列のアクセスなど
B-Treeの改善を低水準(CPU)で記載されている

連分数の計算

(有理算術演算: 高精度数値計算のためのアルゴリズムとプログラミング 2023) を参考にして、連分数を計算するコードを書いた。

関数: get_continued_fraction で連分数展開 [a0, a1, …] を取得する
関数: get_approx_fraction で連分数展開から近似分数を計算する

import math
from fractions import Fraction

def get_continued_fraction(n):
    """入力は素数と仮定して、循環の最初の周期を返却する"""
    r = 0
    s = 1
    a = int(math.floor(math.sqrt(n)))
    a0 = a
    lst_a = [a]
    cnt = 100
    while 0 < cnt:
        _s = (n-(r-a*s)**2) / s
        _r = -(r-a*s)
        a = int(math.floor( (math.sqrt(n) + _r) / _s))
        r = _r
        s = _s
        cnt -= 1
        lst_a.append(a)
        if a == 2*a0:
            break
    return lst_a

def get_approx_fraction(cfrac):
    assert 0 < len(cfrac)

    if 1 == len(cfrac):
        return Fraction(cfrac[0], 1)
    elif 2 == len(cfrac):
        return Fraction(cfrac[0] * cfrac[1] + 1, cfrac[1])
    else:
        pi_1 = cfrac[0] * cfrac[1] + 1
        pi = cfrac[2] * pi_1 + cfrac[0]
        qi_1 = cfrac[1]
        qi = cfrac[2] * qi_1 + 1
        for ai in cfrac[3:]:
            pi, pi_1 = ai * pi + pi_1, pi
            qi, qi_1 = ai * qi + qi_1, qi
        return Fraction(pi, qi)