みなさまこんにちは。Nzt3です。

SageMathはPythonベースの言語でさまざまな数学ソフトウェアを使用できるオープンソースのソフトウェアです。計算に必要な機能がPythonに似た記述で利用できます。

こちらからダウンロードできます。https://www.sagemath.org/
SageMath Cell を使えばオンラインで実行することもできます。

対話型インタプリタが使えますが、jupyter notebookでの利用もできます。

今回はそのSageMathのさまざまな機能を使ってフィボナッチ数列のn項目を計算します。

定義 フィボナッチ数列

次の漸化式で定義します。

$$\begin{gathered} a_0=a_1=1 \\ {a}_n=a_{n-1}+a_{n-2}(n\ge 2) \end{gathered}$$

今回の目的は$f(n)=a_n$となる関数$f$を実装することです。

愚直実装

漸化式を素直に実装してみましょう。

def fib(n):
  ans=[0]*(n+2)
  ans[0],ans[1]=1,1
  for i in range(2,n+1):
    ans[i]=ans[i-1]+ans[i-2]
  return ans[n]

実行してみると
$1,1,2,3,5,8,13,21,34,55\dots$ という列が実際に計算できていることがわかります。

SageMathではPythonの記述は大体はそのまま動きます。上記のコードもSageMathで動作します。

途中の計算結果はいらないので変数を使い回すこともできます。

def fib(n):
  a,b=1,1
  for _ in range(n):
    a,b=b,a+b
  return a

これもSageMathで問題なく動作します。

一般項の利用

フィボナッチ数列は隣接3項間漸化式で定義されます。これの一般項の導出は高校数学の定番問題です。

$x^2=x+1$ の2つの解を$\alpha ,\beta$ とすれば定数$A,B$を用いて
$a_n=A{\alpha }^n+B{\beta }^n$ と表すことができます。

ところで、$\alpha ,\beta ,A,B$を計算するのが面倒ですね。
SageMathに計算してもらいましょう。

var('x')
solve([x^2==x+1],x)

出力: [x == -1/2*sqrt(5) + 1/2, x == 1/2*sqrt(5) + 1/2]
変数$x$を宣言して方程式$x^2=x+1$を$x$について解いてもらいました。
解が得られたので、$\alpha =\frac{1-\sqrt{5}}{2},\beta =\frac{1+\sqrt{5}}{2}$としましょう。
あとは最初の3項に合うように$A,B$を計算してもらいます。

alpha = -1/2*sqrt(5) + 1/2
beta = 1/2*sqrt(5) + 1/2
var('A','B')
solve([A*alpha^0+B*beta^0==1,A*alpha^1+B*beta^1==1,A*alpha^2+B*beta^2==2],A,B)

出力: [[A == -1/10*sqrt(5) + 1/2, B == 1/10*sqrt(5) + 1/2]]
解が得られました。
しれっとalpha^2みたいに書いていますが、SageMathではこれが累乗の記法になります。xor演算子は^^となるのでPythonのプログラムを移植する際には注意しましょう。

def fib(n):
  alpha = -1/2*sqrt(5) + 1/2
  beta = 1/2*sqrt(5) + 1/2
  A = -1/10*sqrt(5) + 1/2
  B = 1/10*sqrt(5) + 1/2
  return A*alpha^n+B*beta^n

これで良さそうなので試しにfib(10)を計算してみましょう。

fib(10)

出力: 1/10240*(sqrt(5) + 5)*(sqrt(5) + 1)^10 - 1/10240*(sqrt(5) - 1)^10*(sqrt(5) - 5)
数式がそのまま出てきてしまいました。これでは使えないので項を展開して計算してもらいましょう。

expandを使うと展開してくれます。

fib(10).expand()

出力:89
計算できましたね。

mod P での計算

大きな$n$についてfib(n)を計算したいです。
しかし、$n$を大きくするとfib(n)も巨大になって計算に時間がかかります。
値は非常に大きくなる可能性があるので$998244353$で割ったあまりを計算することにしましょう。

一般項を手に入れましたが、これを$998244353$で割ったあまりにするにはどうしたら良いでしょうか。
整数の計算は問題ありませんが、$\sqrt{5}$という数が厄介です。ですが、最終的な計算結果では$\sqrt{5}$は出てこないので何かしら$x^2=5$を満たす$x$であれば何を使っても計算結果は正しくなるはずです。

$x^2\equiv 5(\mathrm{mod}998244353)$ の解を調べましょう。といってもどうやって解けばいいのかわかりません。SageMathに調べてもらいましょう。

$x^2-5$の根を調べてもらいます。

var('x')
R.<x>=PolynomialRing(Zmod(998244353))
(x^2-5).roots()

出力: []
根はないらしいです。終わりだ……
そう言えば$\mathbb{Z}/p\mathbb{Z}$って体になりますね。それに$\sqrt{5}$を添加して体を拡大してしまいましょう。
一般項を ${\mathbb{F}}_p[x]/(x^2-5)$ 上の計算にします。

def fib(n):
  K.<s>=GF(998244353^2,modulus=x^2 -5)
  alpha = -K(1/2)*s + K(1/2)
  beta = K(1/2)*s + K(1/2)
  A = -K(1/10)*s + K(1/2)
  B = K(1/10)*s + K(1/2)
  return A*alpha^n+B*beta^n

fib(100000000000000)

出力: 273740723
計算できていそうです！よかった〜

mod m での計算

素数以外の$m$についてもfib(n)を$m$で割ったあまりを知りたいです。

一般項には $1/2,1/5$ と $\sqrt{5}$ が出てくるのでこれをどうにかしましょう。

$m=3^{10}$とすれば$1/2,1/5$は存在します。これに$\sqrt{5}$を追加して計算してみます。

def fib(n):
  R.<x> = PolynomialRing(Zmod(3^10))
  K.<s> = R.quotient(x^2 - 5)
  alpha = -K(1/2)*s + K(1/2)
  beta = K(1/2)*s + K(1/2)
  A = -K(1/10)*s + K(1/2)
  B = K(1/10)*s + K(1/2)
  return A*alpha^n+B*beta^n

fib(30)

出力:47191

$R=\mathbb{Z}/m\mathbb{Z}[x],K=R/(x^2-5)$です。
実際に$a_{30}=1346269$を$3^{10}$で割ったあまりが計算できました！

$1/2,1/5$が存在すれば問題なく一般項を計算できます。

行列の利用

線形漸化式は行列積で表現できます。
フィボナッチ数列の場合は次のようになります。

$$\left(\begin{array}{c}{a}_n\\ a_{n-1}\end{array}\right)=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}{a}_{n-1}\\ a_{n-2}\end{array}\right)$$

これを初項まで繰り返していくことで行列の累乗で一般項を表記できます。

$$\left(\begin{array}{c}{a}_n\\ a_{n-1}\end{array}\right)={\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)}^n\left(\begin{array}{c}1\\ 0\end{array}\right)(n\ge 0)$$

$a_{-1}$という謎の要素が出ていますが、$a_n(n\ge 0)$は正しく計算できているので問題ありません。

これをSageMathで実装しましょう。

def fib(n):
  M = matrix([[1,1],[1,0]])
  A = matrix([[1],[0]])
  return (M^n*A)[0,0]

$M=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right),A=\left(\begin{array}{c}1\\ 0\end{array}\right)$ として$M^{n}A$の第一成分を取り出します。

fib(10)

出力: 89
計算できていますね！

同じ方法で$a_n$を$998244353$で割ったあまりを求めることにします。
今回は全部整数の計算なので困る要素がありません。
Zmodを使って割ったあまりの計算にします。

def fib(n):
  R = Zmod(998244353)
  M = matrix(R,[[1,1],[1,0]])
  A = matrix(R,[[1],[0]])
  return (M^n*A)[0,0]

fib(100000000000000)

出力: 273740723
簡単に変更することができました。

$M^n=M^{100000000000000}$を計算しているので非常に長い時間がかかりそうなものですが、実際には$O(\log (n))$回くらいの計算で済むように賢い計算をしてくれるので変な数を入れてもちゃんと計算できます。

母関数の利用

数列$\{a_n\}$について$A(x)={\sum }_{n=0}^{\infty }{a}_n{x}^n$ とした関数を考えることができます。これがその数列の母関数です。
この関数が値を持つかは知りませんが、係数列だけ考えます。

$a_n$がフィボナッチ数列のとき母関数はどうなるでしょうか。答えは$\frac{1}{1-x-x^2}$です。

この母関数の係数列を求められればフィボナッチ数列が計算できます。

SageMathで冪級数環を使って計算してみましょう。

def fib(n):
  R.<x> = PowerSeriesRing(ZZ,n+1)
  return (1/(1-x-x^2))[n]

fib(10)

出力: 89

PowerSeriesRing(ZZ,n+1)は整数係数の冪級数の係数を最初の$n+1$項($x^0,x^1,\dots ,x^n$)計算することを意味します。
あとは係数列を計算させて$x^n$の項の係数を見ればフィボナッチ数が計算できます。

記述が随分とシンプルになりましたね！

おまけに ${\sum }_{k=0}^n{a}_k$ も計算してみましょう。

def fibsum(n):
  R.<x> = PowerSeriesRing(ZZ,n+1)
  return (1/(1-x-x^2)/(1-x))[n]

fibsum(10)

出力: 232
$1-x$ で割るだけで総和になるので変更もシンプルです。

おわり

SageMathを使うと$\sqrt{5}$だったり${\mathbb{F}}_p^2$だったりといったちょっと面倒なものを手軽に使うことができます。
適当に計算をするとき用のPythonに加えて、適当に数学的操作をするとき用にSageMathも使えると便利です。

みなさまも急に適当な計算をして家族や友達を驚かせましょう！

アドベントカレンダーの明日の担当は @mumumu です！
明日もお楽しみに！