月: 2024年10月

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

概要

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

計算

正弦関数（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
6NWdAb_jDAk_M.png

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

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（アリーナ・アロケーター）」 というものも使われている

“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)

参考文献

有理算術演算: 高精度数値計算のためのアルゴリズムとプログラミング. 2023. 森北出版. https://books.google.co.jp/books?id=PVAP0AEACAAJ.

“Jujutsu in practice”を読んだ

Jujutsu in practice を読んだ。

  • Jujutsu は新しいバージョン管理ツールであるである
  • 既存のgitレポジトリにも適用可能である
    • .jjというフォルダが作成される
  • Jujutsu はブランチでなく、リビジョンを主に扱う
  • Jujutsu には ステージング はない、リビジョンが全てである

ワークフローの全容はこの記事のみでは掴めなかった。次に作るプログラムで試してみようか？

“Everything I Know About The Fast Inverse Square Root Algorithm” を読んだ

githublog/2024/5/29/fast-inverse-sqrt.md at main · francisrstokes/githublog を読んだ。

  • \(\frac{1}{\sqrt x}\) を高速に計算するC言語のコードを説明する記事である
  • float(単精度浮動小数点、32bit)のbitの表現を利用している
    • 仮数(23bit)、指数(8bit)、符号(1bit)
  • float型をlong型に変更して1bitシフトしてfloat型に変更すると、 log2 相当の処理になる
  • さらにニュートン法を利用して、近似の精度を向上させている