https://naoyat.hatenablog.jp/entry/2017/06/17/030000

こないだのJune Long ChallengeのEuler Sumを考えてる時に
Mark Jason Dominus氏のスライド
Arithmetic with Continued Fractions
を読んだ時のメモ。

連分数とかいじるのって多分昔SICPを読んだ時以来。

$\displaystyle x=[a_0;a_1,a_2,a_3,\cdots]=a_0+\frac{1}{a_1+\frac{1}{a_2+\frac{1}{a_3+\cdots}}}$
のような記法を用いる。

$\sqrt2=[1;2,2,2,\cdots$ ]
$e=[2;1,2,1,1,4,1,1,6,1,1,8,1,\cdots$ ]
のような規則的な形に展開される。

連分数 $x$ に対し、 $\displaystyle\frac{a+bx}{c+dx}$ に相当する連分数を求めるのは
$\begin{pmatrix} a&b \\ c&d \end{pmatrix}$
に、例えば $e=[2;1,2,1,1,4,1,1,6,1,1,8,1,\cdots]$ なら
$\begin{pmatrix} 0&1 \\ 1&\bf2 \end{pmatrix}\begin{pmatrix} 0&1 \\ 1&\bf1 \end{pmatrix}\begin{pmatrix} 0&1 \\ 1&\bf2 \end{pmatrix}\begin{pmatrix} 0&1 \\ 1&\bf1 \end{pmatrix}\begin{pmatrix} 0&1 \\ 1&\bf1 \end{pmatrix}\begin{pmatrix} 0&1 \\ 1&\bf4 \end{pmatrix}\begin{pmatrix} 0&1 \\ 1&\bf1 \end{pmatrix}\cdots$ を右から食わせて
$\displaystyle\lfloor\frac{a}{c}\rfloor=\lfloor\frac{b}{d}\rfloor$
になった時点で随時 $\displaystyle\frac{a}{c}$ を取り出して $\begin{pmatrix} c&d \\ a\%c&b\%d \end{pmatrix}$ に差し替えていくことで求められる。
10進表記が得たいなら、1つ取り出した後の残りに $\displaystyle\frac{10}{x}$ すなわち $\begin{pmatrix} 10&0 \\ 0&1 \end{pmatrix}$ を適用していけば良い。

連分数 $x, y$ に対し、 $\displaystyle\frac{a+bx+cy+dxy}{e+fx+gy+hxy}$ に相当する連分数を求めるのは
$\begin{pmatrix} a&b&c&d \\ e&f&g&h \end{pmatrix}$ に
$x$ 側から $\begin{pmatrix} 0&1&0&0 \\ 1&\bf x_i&0&0 \\ 0&0&0&1 \\ 0&0&1&\bf x_i \end{pmatrix}$
$y$ 側から $\begin{pmatrix} 0&0&1&0 \\ 0&0&0&1 \\ 1&0&\bf y_j&0 \\ 0&1&0&\bf y_j \end{pmatrix}$
を右から食わせて
$\displaystyle\lfloor\frac{a}{e}\rfloor=\lfloor\frac{b}{f}\rfloor=\lfloor\frac{c}{g}\rfloor=\lfloor\frac{d}{h}\rfloor$
になった時点で随時 $\displaystyle\frac{a}{e}$ を取り出して $\begin{pmatrix} e&f&g&h \\ a\%e&b\%f&c\%g&d\%h \end{pmatrix}$ に差し替えていくことで求められる。

次にどちらを食わせるかは
$\displaystyle|\frac{b}{f}-\frac{a}{e}|>|\frac{c}{g}-\frac{a}{e}|$ なら $x$ から、そうでなければ $y$ から。

$\displaystyle\frac{(x+3)(y+4)}{x-y}$ みたいな演算は、分子・分母をそれぞれ計算してから…でなくても、 $\displaystyle\frac{12+4x+3y+xy}{x-y}$ なので
$\begin{pmatrix} 12&4&3&1 \\ 0&1&-1&0 \end{pmatrix}$ で一発で求めることができる。

以下、pythonで実装してみたやつ

$[2;1,2,1,1,4,1,1,6,1,\cdots$ ] みたいな規則的な数列を生成するジェネレータを用意する

10進展開とかしようとすると500桁も行かないうちに

RuntimeError: maximum recursion depth exceeded

とか言われるので sys.setrecursionlimit() の出番。

$ python renbunsuu e 1000
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117301238197068416140397019837679320683282376464804295311802328782509819455815301756717361332069811250996181881593041690351598888519345807273866738589422879228499892086805825749279610484198444363463244968487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016768396424378140592714563549061303107208510383750510115747704171898610687396965521267154688957035035

これで所要時間は5秒ぐらい。

分子・分母ともに巨大な整数になるので、np.dot()でポンみたいなわけにはいかない...

#!/usr/bin/env python
# -*- encoding: utf-8 -*-
#
# > python renbunsuu.py e 100
#
import sys
import re
import click
import math


class ContUnary:
    def __init__(self,a=0,b=1,c=1,d=0):
        self.args = (a,b,c,d)

    def input(self, p):
        a,b,c,d = self.args
        # ab 01
        # cd 1p
        self.args = (b, a+b*p, d, c+d*p)

    def ready(self):
        a,b,c,d = self.args
        if c>0 and d>0:
            return a/c == b/d
        else:
            return False

    def output(self):
        # 01 ab
        # 1-q cd
        assert self.ready()
        a,b,c,d = self.args
        q = a/c
        self.args = (c, d, a%c, b%d)
        return q

    def apply(self, g, stop=-1):
        # self.g = g
        cnt = 0
        while True:
            while not self.ready():
                self.input(g.next())
            yield self.output()
            cnt += 1
            if cnt == stop: break

    def mul(self,e,f,g,h):
        a,b,c,d = self.args
        self.args = (a*e+b*g, a*f+b*h, c*e+d*g, c*f+d*h)


def e_gen():
    yield 2
    k = 2
    while True:
        yield 1
        yield k
        k += 2
        yield 1


def phi_gen():
    while True:
        yield 1


def root_gen(n):
    ren_dic = {
        # 2: [1, 2],
        3: [1, 1, 2],
        # 5: [2, 4],
        6: [2, 2, 4],
        7: [2, 1, 1, 1, 4],
        # 10: [3, 6],
        19: [4, 2, 1, 3, 1, 2, 8],
        29: [5, 2, 1, 1, 2, 10],
        43: [6, 1, 1, 3, 1, 5, 1, 3, 1, 1, 12],
        54: [7, 2, 1, 6, 1, 2, 14],
        76: [8, 1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1, 16],
        94: [9, 1, 2, 3, 1, 1, 5, 1, 8, 1, 5, 1, 1, 3, 2, 1, 18],
        1000: [31, 1, 1, 1, 1, 1, 6, 2, 2, 15, 2, 2, 6, 1, 1, 1, 1, 1, 62],
    }

    q = int(math.sqrt(n - 1))
    if q*q+1 == n:
        yield q
        q *= 2
        while True:
            yield q
    elif n in ren_dic:
        a0, rest = ren_dic[n][0], ren_dic[n][1:]
        assert rest[-1] == a0 * 2
        yield a0
        while True:
            for a in rest:
                yield a


def render(src, base=10):
    U = ContUnary(b=1)
    gen = U.apply(src)
    yield gen.next()
    while True:
        gen = ContUnary(a=base,b=0,c=0,d=1).apply(gen)
        yield gen.next()


@click.command()
@click.argument('type', type=str)
@click.argument('nterms', type=int, default=10)
@click.argument('base', type=int, default=10)
def main(type, nterms, base):
    if type == 'e':
        gen = e_gen()
    elif type == 'phi':
        gen = phi_gen()
    elif re.match(r'\d+', type):
        gen = root_gen(int(type))
    else:
        return

    sys.setrecursionlimit(1000000)

    gen = render(gen, base)
    sys.stdout.write('%d.' % gen.next())

    for i, d in enumerate(gen):
        if i == nterms-1: break
        sys.stdout.write('%d' % d)
        sys.stdout.flush()

    sys.stdout.write('\n')
    return


if __name__ == '__main__':
    main()