中国の剰余定理のアルゴリズム(Algorithm.10.14)

概要

036.crt.html に引き続き、 (Gathen and Gerhard 2013) を参考に中国の剰余アルゴリズムを実装している。

“10.3 Fast Chinese remaindering” に6つの疑似コードがある。

  • 済: Algorithm.10.3
  • Algorithm.10.14
  • Algorithm.10.16
  • Algorithm.10.18
  • Algorithm.10.20
  • Algorithm.10.22

本記事では “Algorithm.10.14” を調査、実装する。

Algorithm.10.14

このアルゴリズムは一つの整数(被乗数)に対して、複数の剰余を求めるものである。

入力を以下とする。( \(r=2^{k}\) とする.)

  • \(ms=m_{1},\cdots,m_{r-1}\) : 互いに素な整数(除数)のリスト
  • \(M\) : alg10_3 に \(ms\) を入力して得られる行列
  • \(f\) : 整数(被乗数)

出力を以下とする。

  • \(f \mod m_{1}, \cdots ,f \mod m_{r-1}\) : 複数の剰余

このアルゴリズムは再帰呼び出しで定義されている。 ms=[2, 3, 5, 7, 11, 13, 17, 19], f=9699791 とすると以下のような呼び出しとなる。

  • alg10_14(ms=[2, 3, 5, 7, 11, 13, 17, 19], M=algo10_3(ms), f=9699791);
    • alg10_14(ms=[2, 3, 5, 7], M=algo10_3(ms), f=101)
      • alg10_14(ms=[2, 3], M=algo10_3(ms), f=5)
        • alg10_14(ms=[2], M=algo10_3(ms), f=1)
        • alg10_14(ms=[3], M=algo10_3(ms), f=2)
      • alg10_14(ms=[5, 7], M=algo10_3(ms), f=31)
        • alg10_14(ms=[5], M=algo10_3(ms), f=1)
        • alg10_14(ms=[7], M=algo10_3(ms), f=3)
    • alg10_14(ms=[11, 13, 17, 19], M=algo10_3(ms), f=101)
      • alg10_14(ms=[11, 13], M=algo10_3(ms), f=101)
        • alg10_14(ms=[11], M=algo10_3(ms), f=2)
        • alg10_14(ms=[13], M=algo10_3(ms), f=10)
      • alg10_14(ms=[17, 19], M=algo10_3(ms), f=101)
        • alg10_14(ms=[17], M=algo10_3(ms), f=16)
        • alg10_14(ms=[19], M=algo10_3(ms), f=6)

実装

実装は疑似コードを忠実に再現したものと、効率を考えたもので2つ作成した。

def alg10_14(ms, f):
    if len(ms) == 1:
        return [f]
    else:
        M = alg10_3(ms)
        f0 = f % M[-2][0]
        f1 = f % M[-2][1]
        return alg10_14(ms[:len(ms)//2], f0) + alg10_14(ms[len(ms)//2:], f1)

def alg10_14x(i_ms_s, i_ms_e, f, M, i_f0, j_f0):
    if i_ms_e - i_ms_s == 1:
        return [f]
    else:
        f0 = f % M[i_f0][j_f0]
        f1 = f % M[i_f0][j_f0+1]
        return alg10_14x(i_ms_s, i_ms_s+(i_ms_e-i_ms_s)//2, f0, M, i_f0-1, 2*j_f0) \
             + alg10_14x(i_ms_s+(i_ms_e-i_ms_s)//2, i_ms_e, f1, M, i_f0-1, 2*j_f0+2)

ms = [2, 3, 5, 7, 11, 13, 17, 19]
f = 9699791

res0 = [f % m for m in ms]

res1 = alg10_14(ms=ms, f=f)

M = alg10_3(ms)
res2 = alg10_14x(0, len(ms), f, M, -2, 0)

assert all([x==y for x, y in zip(res0, res1)])
assert all([x==y for x, y in zip(res1, res2)])
res1
1 2 1 3 2 10 16 6

参考文献

Gathen, Joachim von zur, and Gerhard Jürgen. 2013. Modern computer algebra. 3rd ed. Cambridge: Cambridge University Press.

中国の剰余定理のアルゴリズム(Algorithm.10.14)

概要

036.crt.html に引き続き、 (Gathen and Gerhard 2013) を参考に中国の剰余アルゴリズムを実装している。

“10.3 Fast Chinese remaindering” に6つの疑似コードがある。

  • 済: Algorithm.10.3
  • Algorithm.10.14
  • Algorithm.10.16
  • Algorithm.10.18
  • Algorithm.10.20
  • Algorithm.10.22

本記事では “Algorithm.10.14” を調査、実装する。

Algorithm.10.14

このアルゴリズムは一つの整数(被乗数)に対して、複数の剰余を求めるものである。

入力を以下とする。( \(r=2^{k}\) とする.)

  • \(ms=m_{1},\cdots,m_{r-1}\) : 互いに素な整数(除数)のリスト
  • \(M\) : alg10_3 に \(ms\) を入力して得られる行列
  • \(f\) : 整数(被乗数)

出力を以下とする。

  • \(f \mod m_{1}, \cdots ,f \mod m_{r-1}\) : 複数の剰余

このアルゴリズムは再帰呼び出しで定義されている。 ms=[2, 3, 5, 7, 11, 13, 17, 19], f=9699791 とすると以下のような呼び出しとなる。

  • alg10_14(ms=[2, 3, 5, 7, 11, 13, 17, 19], M=algo10_3(ms), f=9699791);
    • alg10_14(ms=[2, 3, 5, 7], M=algo10_3(ms), f=101)
      • alg10_14(ms=[2, 3], M=algo10_3(ms), f=5)
        • alg10_14(ms=[2], M=algo10_3(ms), f=1)
        • alg10_14(ms=[3], M=algo10_3(ms), f=2)
      • alg10_14(ms=[5, 7], M=algo10_3(ms), f=31)
        • alg10_14(ms=[5], M=algo10_3(ms), f=1)
        • alg10_14(ms=[7], M=algo10_3(ms), f=3)
    • alg10_14(ms=[11, 13, 17, 19], M=algo10_3(ms), f=101)
      • alg10_14(ms=[11, 13], M=algo10_3(ms), f=101)
        • alg10_14(ms=[11], M=algo10_3(ms), f=2)
        • alg10_14(ms=[13], M=algo10_3(ms), f=10)
      • alg10_14(ms=[17, 19], M=algo10_3(ms), f=101)
        • alg10_14(ms=[17], M=algo10_3(ms), f=16)
        • alg10_14(ms=[19], M=algo10_3(ms), f=6)

実装

実装は疑似コードを忠実に再現したものと、効率を考えたもので2つ作成した。

def alg10_14(ms, f):
    if len(ms) == 1:
        return [f]
    else:
        M = alg10_3(ms)
        f0 = f % M[-2][0]
        f1 = f % M[-2][1]
        return alg10_14(ms[:len(ms)//2], f0) + alg10_14(ms[len(ms)//2:], f1)

def alg10_14x(i_ms_s, i_ms_e, f, M, i_f0, j_f0):
    if i_ms_e - i_ms_s == 1:
        return [f]
    else:
        f0 = f % M[i_f0][j_f0]
        f1 = f % M[i_f0][j_f0+1]
        return alg10_14x(i_ms_s, i_ms_s+(i_ms_e-i_ms_s)//2, f0, M, i_f0-1, 2*j_f0) \
             + alg10_14x(i_ms_s+(i_ms_e-i_ms_s)//2, i_ms_e, f1, M, i_f0-1, 2*j_f0+2)

ms = [2, 3, 5, 7, 11, 13, 17, 19]
f = 9699791

res0 = [f % m for m in ms]

res1 = alg10_14(ms=ms, f=f)

M = alg10_3(ms)
res2 = alg10_14x(0, len(ms), f, M, -2, 0)

assert all([x==y for x, y in zip(res0, res1)])
assert all([x==y for x, y in zip(res1, res2)])
res1
1 2 1 3 2 10 16 6

参考文献

Gathen, Joachim von zur, and Gerhard Jürgen. 2013. Modern computer algebra. 3rd ed. Cambridge: Cambridge University Press.

中国の剰余定理のアルゴリズム(Algorithm.10.3)

概要

(Gathen and Gerhard 2013) を参考に中国の剰余アルゴリズムを実装している。

“10.3 Fast Chinese remaindering” に6つの疑似コードがある。

  • Algorithm.10.3
  • Algorithm.10.14
  • Algorithm.10.16
  • Algorithm.10.18
  • Algorithm.10.20
  • Algorithm.10.22

上記をPythonで実装しようと考えた。

Algorithm.10.3 の解説

このアルゴリズムは他のアルゴリズムから利用される補助的なものである。本文では多項式環を扱えるように一般化しているが、今回は \(\mathbb{Z}/p\mathbb{Z}\) だけを扱うことにする。

入力は互いに素な数である。入力のサイズはパラメータ \(k\) で決定する。 \(r=2^{k}\) 個の値を入力とすることになる。出力は \(r \times r\) の行列となる。

k=2として具体的な計算手順を実行する。入力は \(r=2^{k}=4\) の数 \(m_{1},m_{2},m_{3},m_{4}\) とする。出力となる行列 \(M\) を \(r\times r=4×4\) の大きさで初期化する。

\begin{pmatrix}
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
\end{pmatrix}

初期化した行列の1行目に入力の数を配置する。

\begin{pmatrix}
m_{1} & m_{2} & m_{3} & m_{4}\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
\end{pmatrix}

行列の2行目に1行目の隣同士の積を左詰めで設定する。

\begin{pmatrix}
m_{1} & m_{2} & m_{3} & m_{4}\\
m_{1}m_{2} & m_{3}m_{4} & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
\end{pmatrix}

行列の3行目に2行目の隣同士の積を左詰めで設定する。

\begin{pmatrix}
m_{1} & m_{2} & m_{3} & m_{4}\\
m_{1}m_{2} & m_{3}m_{4} & 0 & 0\\
m_{1}m_{2}m_{3}m_{4} & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
\end{pmatrix}

上記が出力となる。

Python実装

参考文献の疑似コードをPythonで実装した。

def alg10_3(ms):
    k=0
    while 2**k < len(ms):
        k += 1
    r = 2**k
    ms = ms + ([0] * (r-len(ms)))
    M = [[0]*r for _ in range(k+1)]

    for i in range(r):
        M[0][i] = ms[i]

    for i in range(1, k+1):
        for j in range(0, 2**(k-i)):
            M[i][j] = M[i-1][2*j] * M[i-1][2*j+1]

    return M

出力は以下のようなる。

alg10_3([2, 3, 5, 7, 11, 13, 17, 19])
2 3 5 7 11 13 17 19
6 35 143 323 0 0 0 0
210 46189 0 0 0 0 0 0
9699690 0 0 0 0 0 0 0

参考文献

Gathen, Joachim von zur, and Gerhard Jürgen. 2013. Modern computer algebra. 3rd ed. Cambridge: Cambridge University Press.

“You could have designed state of the art positional encoding”を読んだ

You could have designed state of the art positional encoding を読んだ。

  • LLMに利用される位置エンコーディングについて解説した記事である
  • 整数位置エンコーディング、バイナリー位置エンコーディング、正弦位置エンコーディングが紹介されている
  • 文章内の単語の位置は絶対位置よりも相対位置が重要である
  • RoFormer: Enhanced Transformer with Rotary Position Embedding | Abstract が紹介されている
  • 4つの位置エンコーディングにアニメーションが利用されている

TikzをEmacsのOrgモードで使う

概要

TikZによるLaTeXグラフィックス | Stefan Kottwitz, 黒川 利明 を購入した。 TikZはやろうやろうと考えていたのでこれを気にある程度勉強したい。

書籍中のコードは PacktPublishing/LaTeX-graphics-with-TikZ: LaTeX graphics with TikZ, by Packt Publishing にある。

Emacsの環境

書籍のP.008に掲載されているコードは以下にある。

上記をEmacsのOrgモードのコードブロックで使えるようにする。ちょっと複雑になってしまうが、あとで実用的な方法を考えよう。

#+header: :file "images/v0ud5YVim1jv4.png"
#+header: :results raw :fit yes :border 0cm
#+header: :imagemagick t
#+header: :iminoptions -density 400
#+header: :imoutoptions -geometry 200 -flatten
#+begin_src latex :results file
  \begin{tikzpicture}
    \draw[thin,dotted] (-3,-3) grid (3,3);
    \draw[->] (-3,0) -- (3,0);
    \draw[->] (0,-3) -- (0,3);
  \end{tikzpicture}
#+end_src
v0ud5YVim1jv4.png

参考文献

Mergiraf: Gitのマージドライバー?

概要

MergirafTree-sitter を利用したgitのマージを補助するツールのようだ。

gitで他の開発者が実装したブランチをマージする際に開発箇所が被ると衝突の回避をする必要がある。gitがうまく解決してくれることがあるが、手動は気を使う作業になる。

Margiral はその衝突を回避するために構文(syntax)を考慮してマージ作業をしてくれるようだ。

gitやバージョン管理のツールはこれ以外にもいくつかあるので、併せて試してみたい。

参考文献

“Understanding Machine Learning” Lemma19.2 について

概要

(Shalev-Shwartz and Ben-David 2009) を読んでいる。ふと定理の数式を実験できないかと考えた。

lemma 19.2

詳細は省くが、以下のような数式がある。 これを実験できないか?

\[ \underset{S \sim \mathcal{D}^m}{\mathbb{E}}\left[\sum_{i: C_i \cap S=\emptyset} \mathbb{P}\left[C_i\right]\right] \leq \frac{r}{m e} \]

コード

\(\mathcal{D}^m\) を 0 から 1000000 の値を一様に出力する確率分布だとして実験する。

import numpy as np

def D(size, low=0, high=1000000):
    return np.random.randint(X_low, X_high, size)
def P(n, low=0, high=1000000):
    return n / (X_high - X_low)

m = 2000
r = 10
C = D((r, 10))

N = 100
lst_sum_P = []
for _ in range(N):
    sum_P = 0
    S = D(m)
    for i in range(r):
        if np.intersect1d(S, C[i]).size == 0:
            pass
        else:
            sum_P += P(C[i].size)
    lst_sum_P.append(sum_P)

np.mean(lst_sum_P) <= r / (m * np.e)
True

参考文献

Shalev-Shwartz, Shai, and Shai Ben-David. 2009. Understanding Machine Learning. Cambridge University Press. https://doi.org/10.1017/cbo9781107298019.

“Breaking CityHash64, MurmurHash2/3, wyhash, and more…” を読んだ

Breaking CityHash64, MurmurHash2/3, wyhash, and more… を読んだ。

  • 非暗号学的ハッシュ関数(non-criptographic hash function)をセキュリティ目的に使った場合の脆弱性について
  • 暗号学的ハッシュ関数の安全性は Pre-image resistance, Second pre-image resistance, Collision resistance で満足すべき
  • ハッシュ関数の実装調べて、逆演算を用いることで脆弱性を突いている

卒論でMD6のコンポーネントに対して高階差分解析をやったことを思いだした。また、ハッシュ関数を実装してみよう。

winstonでログ出力

winston を使ったスタックトレースの出しかた。

logger.ts

import { createLogger, format, transports } from 'winston';

const { combine, timestamp, errors, printf, cli, metadata } = format;


const myFormat = printf(({ level, message, metadata }) => {
  if (metadata && metadata.stack) {
    return `${metadata.timestamp} ${level}: ${metadata.stack}`;
  } else {
    return `${metadata.timestamp} ${level}: ${message}`;
  }
});

export const logger = createLogger({
  format: combine(
    errors({ stack: true }),
    timestamp(),
    cli(),
    metadata(),
    myFormat,
  ),
  transports: [new transports.Console()],
});

logger.test.ts

import {logger} from '../../src/logger.js'
import { test } from 'vitest'


test("winstonの通常ログ", () => {
  logger.info("これはINFOです。");
});


test("winstonの例外時にスタックトレースを表示する", () => {
  try {
    throw new RangeError("これは範囲エラーです.");
  }catch(err) {
    logger.error(err);
  }
});

テスト出力

 RUN  v2.1.0 /home/hoge/projects/hello_typescript

stdout | __tests__/unit/logger.test.ts > winstonの通常ログ
2024-11-04T10:58:55.786Z info:     これはINFOです。

stdout | __tests__/unit/logger.test.ts > winstonの例外時にスタックトレースを表示する
2024-11-04T10:58:55.789Z error: RangeError: これは範囲エラーです.
    at /home/hoge/projects/hello_typescript/__tests__/unit/logger.test.ts:12:11
    at file:///home/hoge/projects/hello_typescript/node_modules/@vitest/runner/dist/index.js:146:14
    at file:///home/hoge/projects/hello_typescript/node_modules/@vitest/runner/dist/index.js:529:11
    at runWithTimeout (file:///home/hoge/projects/hello_typescript/node_modules/@vitest/runner/dist/index.js:61:7)
    at runTest (file:///home/hoge/projects/hello_typescript/node_modules/@vitest/runner/dist/index.js:982:17)
    at processTicksAndRejections (node:internal/process/task_queues:95:5)
    at runSuite (file:///home/hoge/projects/hello_typescript/node_modules/@vitest/runner/dist/index.js:1131:15)
    at runFiles (file:///home/hoge/projects/hello_typescript/node_modules/@vitest/runner/dist/index.js:1188:5)
    at startTests (file:///home/hoge/projects/hello_typescript/node_modules/@vitest/runner/dist/index.js:1197:3)
    at file:///home/hoge/projects/hello_typescript/node_modules/vitest/dist/chunks/runBaseTests.9YDrdSI4.js:130:11

 ✓ __tests__/unit/logger.test.ts (2)
 ✓ __tests__/unit/main.test.ts (3)
 ✓ build/__tests__/unit/main.test.js (2)

 Test Files  3 passed (3)
      Tests  7 passed (7)
   Start at  19:58:55
   Duration  209ms (transform 52ms, setup 0ms, collect 95ms, tests 12ms, environment 0ms, prepare 141ms)

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