Machin型の円周率公式を格子上のCVPを利用して導く方法について

はじめに

本記事では，ある格子上の最近ベクトル問題（CVP）を用いてSchnorrSch21風の方法によりMachin型の円周率公式を導出するアルゴリズムについて書いていこうと考えている．具体的には，次の形をした円周率公式をCVPへ帰着して導出する．

格子について

本節では，格子について簡単にまとめる．詳しくは格子19を参照せよ．

格子を利用した円周率公式の導出方法について

次の行列$\bm{B}_{c, a_i}$を基底行列に持つ格子$L_{c, a_i}$と目標ベクトル$\bm{t}_c$を考える．
$$ \bm{B}_{c, a_i}\coloneqq \begin{bmatrix} 1 & 0 & \cdots & 0 & 10^c\arctan \frac{1}{a_1}\\ 0 & 1 & \cdots & 0 & 10^c\arctan \frac{1}{a_2}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & 10^c\arctan \frac{1}{a_n} \end{bmatrix} $$
$$ \bm{t}_c\coloneqq \left(0,\cdots,0,\frac{10^c}{4}\pi\right)\in \mathbb{R}^{n+1} $$
このとき，$\bm{v}\in L_{c, a_i}$が$\operatorname{CVP}(\bm{B}_{c,a_i}, \bm{t}_c)$の解だとすると
$$ \begin{aligned} \bm{v}&={}^{\exists}\bm{x}\bm{B}_{c, a_i}\\ &=\sum_{i=1}^n \left[{}^{\exists}x_i\bm{e}_i+10^c\arctan\frac{1}{a_i}\cdot{}^{\exists} x_i \bm{e}_{n+1}\right]\\ &=\left(x_1,\ldots, x_n, 10^c\sum_{i=1}^n x_i\arctan\frac{1}{a_i}\right) \end{aligned} $$
であるため，$\bm{v}\approx \bm{t}$より
$$ \left(x_1,\ldots,x_n, 10^c \sum_{i=1}^n x_i\arctan \frac{1}{a_i}\right)\approx \left(0,\cdots,0,\frac{10^c}{4}\pi\right) $$
なので，
$$ 10^c \sum_{i=1}^n x_i\arctan \frac{1}{a_i}\approx \frac{10^c}{4}\pi\iff \frac{\pi}{4}\approx \sum_{i=1}^n x_i\arctan \frac{1}{a_i} $$
であることが期待される．特に，これは，パラメタ$c$が大きいほど精度が高くなり，より適切な$x_i$が発見できると思われる．

Babaiのアルゴリズムで十分か

CVPの（近似）解を考える上で，例えば数え上げ法やKannanの埋め込み法を利用すればCVPの厳密解を得ることができる．しかし，計算量が非常に大きいため，多項式時間で停止するBabaiの最近平面アルゴリズムで十分な係数が得られるかが非常に興味のある点である．そこで，本小節において，少しその辺りについて議論を行う．

まず，Babaiの最近平面アルゴリズムについて確認する．

上の命題を基に開発されたアルゴリズムがBabaiの最近平面アルゴリズムBabai86である．Babaiのアルゴリズムは，最近ベクトルを出力するとは限らないが，ある程度近いベクトルを多項式時間で出力するという性格をしており，且つ実装も非常に簡単なためBabiのアルゴリズムで十分であるのならば，こんなに素晴らしいことはない．

擬似コードをMathlog上で書くのはやや骨が折れるので，擬似コードの代わりとして$\texttt{SageMath}$での実装例を以下に示す．
実際には，$\texttt{C++}$などで実装するのが良い．

      def babai(B, t):
    bs, _ = B.gram_schmidt()
    b = copy(t)
    for i in xsrange(n - 1, -1, -1):
        q = round(b.inner_product(bs[i])) / (bs[i].norm() ^ 2)
        b -= q * B[i]
    return t - b
    

命題2から
$$ \begin{aligned} &\norm{\left(x_1,\ldots,x_n, 10^c \sum_{i=1}^n x_i\arctan \frac{1}{a_i}\right)- \left(0,\cdots,0,\frac{10^c}{4}\pi\right)}^2\\ =&\norm{\left(x_1,\ldots,x_n, 10^c\left(\sum_{i=1}^n x_i\arctan \frac{1}{a_i}-\frac{\pi}{4}\right)\right)}^2\\ =&\sum_{i=1}^n x_i^2+10^{2c}\left(\sum_{i=1}^n x_i\arctan \frac{1}{a_i}-\frac{\pi}{4}\right)^2\le \frac{2^n-1}{4}\|\bm{b}_n^\star\|^2 \end{aligned} $$
更に，GSAを信じると
$$ \sum_{i=1}^n x_i^2+10^{2c}\left(\sum_{i=1}^n x_i\arctan \frac{1}{a_i}-\frac{\pi}{4}\right)^2\lessapprox 1.02\cdot \frac{2^n-1}{4}\|\bm{b}_1\|^2 $$
となるNguyen09．また，各$x_i$は$x_i\approx \mathcal{O}(1)$だろうと期待できるので，
$$ 10^{2c}\left(\sum_{i=1}^n x_i\arctan \frac{1}{a_i}-\frac{\pi}{4}\right)^2\lessapprox 1.02\cdot \frac{2^n-1}{4}\|\bm{b}_1\|^2+\mathcal{O}(n) $$
つまり
$$ \left(\sum_{i=1}^n x_i\arctan \frac{1}{a_i}-\frac{\pi}{4}\right)^2\lessapprox \frac{1}{10^{2c}}\left\{1.02\cdot \frac{2^n-1}{4}\|\bm{b}_1\|^2+\mathcal{O}(n)\right\} $$
が期待されるため，特に$c$を大きくとればBabaiでも十分であると思われる．

Machin型円周率公式の探索報告

ここでは，実装上の簡単のため，整数格子を考える．即ち
$$ \round{\bm{B}_{c, a_i}}\coloneqq \begin{bmatrix} 1 & 0 & \cdots & 0 & \round{10^c\arctan \frac{1}{a_1}}\\ 0 & 1 & \cdots & 0 & \round{10^c\arctan \frac{1}{a_2}}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & \round{10^c\arctan \frac{1}{a_n}} \end{bmatrix} $$
を基底行列にもち，
$$ \round{\bm{t}_c}\coloneqq \left(0,\cdots,0,\round{\frac{10^c}{4}\pi}\right)\in \mathbb{R}^{n+1} $$
を目標ベクトルに持つCVPを考える．実装コードは
https://github.com/satoshin-des/MachinLattice
で確認可能である．
まずは，実際にMachinの公式が導出できるのかを確かめてみよう．

      $ ./a.exe 5 239
[4 -1 7853981632]
    

より，$a_1=5, a_2=239$のときは，$p_1=4, p_2=-1$であると言っている．即ち，
$$ \frac{\pi}{4}=4\arctan \frac{1}{5}-\arctan \frac{1}{239} $$
ということである．これはMachinの公式そのものである．
実際に幾つか検証を行ってみると以下のようなMachin型の円周率公式が得られた．

$$ \begin{aligned} \frac{\pi}{4}&=\arctan\left(\frac{1}{30}\right)+1 \arctan\left(\frac{1}{2}\right)-1\arctan\left(\frac{1}{242}\right)+2\arctan\left(\frac{1}{4}\right)-1\arctan\left(\frac{1}{5}\right) \end{aligned} $$

また，以下の近似式も得られた

$$ \begin{aligned} \frac{\pi}{4}&\approx -2 \arctan\left(\frac{1}{360}\right)+29 \arctan\left(\frac{1}{43}\right)-5 \arctan\left(\frac{1}{440}\right)-1 \arctan\left(\frac{1}{751}\right)+13 \arctan\left(\frac{1}{99}\right)-2 \arctan\left(\frac{1}{420}\right)+2 \arctan\left(\frac{1}{713}\right) \end{aligned} $$

但し，松元00によれば具体的な構成方法は既にあるらしく，少なくとも数百万個見つかっているらしい．