数値実験 de PARI/GP

PARI/GP

$\hspace{3pt}$ PARI/GP とは，数論特化の計算機代数アプリケーションで，完全無料です．ど素人ながら，今回はこの PARI/GP を用いて数値実験をしてみたいと思います．自分は級数と積分しか知らないので，初等数論などを嗜まれる御方は直ちにブラウザバックなさることが強く推奨されると思います．

始め方：青字の裏のリンク先で "GP in your browser" をクリック
マニュアル："documentation" をクリックすると User's Guide をPDF形式でダウンロードするページになる

$\hspace{5pt}$細々とした操作方法や関数の記法などは最後に「コマンドまとめ」で自分なりにまとめているので，もしエラー出て進まない時などに読んでいただければと思います．

lindep 関数

$\hspace{5pt}$数値実験にあたり，lindep 関数さんに活躍してもらおうと思います．このlindep 関数，簡単に言うと「和が$0$になる係数の組を返す」です．しかしながら，あくまで観測できる数値の範囲でというもので，許容観測範囲の外で誤差が生じている可能性も考慮する必要があり，「厳密にこれが成り立つ」とは言えないと思います．たぶん．
$\hspace{5pt}$実際のコマンドは，r 個の数値 i1,i2,..,ir に対して "lindep([i1,i2,...,ir])" と入力すると，設定した精度のもとで係数の組 "[o1,o2,...,or]" が出力されます．小さい係数が返ってくれば，（これは成り立ちそうだ…），なん十桁もの係数が返ってくれば，（これは成り立たなそうだ…）という感覚です．

具体例 1　級数

$\hspace{5pt}$級数（とくに無限級数）の数値実験をします．使用する関数は sumalt 関数です．

sumalt 関数：sumalt(n=a,f(n)) で$\quad\D\sum_{n=a}^\infty f(n)\quad$を表します．

次の等式があります．ソースは無いです．

これを数値で見ていきます．実際のコマンドを以下に示します．

      ? \p300
   realprecision = 308 significant digits (300 digits displayed)
? B(x,y)=x!/((x-y)!*y!)
%2 = (x,y)->x!/((x-y)!*y!)
? A=sumalt(n=1,(-1)^(n-1)/(2*n)^2*B(2*n,n)^2/(B(3*n,n)*B(6*n,3*n)))
%3 = 0.0165069937265260166676913437968416377675540126881728385374890903187901826225326669644032884891609429020281857824156775146802567001136858734676123528440741379439409947492965160087696240702377038183569546337224613912692692996316598411617397400149051841712195780928871324056100093586452554123651029121737
? B=sumalt(n=0,(-1)^n/(2*n+1)^2*B(2*n,n)^2/(B(3*n+1,n)*B(6*n+1,3*n)))
%4 = 0.996864636344851471884063906265720096191902348336982765763828391813278591815440124478819313465351788923619502846851720916585388832271216196841629179504964248865404582470235472662412378565426655169101108057167693328728218731683449778223221785390968936560901088107291484274347892942190460874086172183199
? lindep([zeta(2),A,B])
%5 = [-3, -3, 5]~
    

"?" が付いている行が入力（実際の入力は "?" より後ろ部分）で，"%数字" が出力です．
だいたいの流れは

$\hspace{20pt}$１．精度の設定

$\hspace{20pt}$２．関数 B(x,y) の定義$\qquad\D\small B(x,y)=\frac{x!}{(x-y)!y!}$

$\hspace{20pt}$３．数 A の定義$\qquad\D\small A=\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n)^2}\frac{\binom{2n}{n}^2}{\binom{3n}{n}\binom{6n}{3n}}$

$\hspace{20pt}$４．数 B の定義$\qquad\D\small B=\sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)^2}\frac{\binom{2n}{n}^2}{\binom{3n+1}{n}\binom{6n+1}{3n}}$

$\hspace{20pt}$５．$\zeta(2)$と$A$と$B$の和が$0$になる係数の出力

です．結果をみると，

を示唆していることになります．
$\hspace{5pt}$また，少し乱暴な主張ですが，こういう式はたいてい weight を１上げても成り立つのでそれも数値で見てしまいます．

      ? \p300
   realprecision = 308 significant digits (300 digits displayed)
? B(x,y)=x!/((x-y)!*y!)
%2 = (x,y)->x!/((x-y)!*y!)
? C=sumalt(n=1,(-1)^(n-1)/(2*n)^3*B(2*n,n)^2/(B(3*n,n)*B(6*n,3*n)))
%3 = 0.00829319533331464343025861754414363981361156213223294552475033315238564646326796460009190579067969346408578036836708702147502921058709882869816135554142615132391631560916655672716392607297744976003823399042811643637263159348554128414055681993736245647846645744959400137647769392005565737941547789307612
? D=sumalt(n=0,(-1)^n/(2*n+1)^3*B(2*n,n)^2/(B(3*n+1,n)*B(6*n+1,3*n)))
%4 = 0.998949687521890356689695595411493779032951389573004752985309518400736622667504920288682592744219483967137371373261570508251653595365369735247610245826108084103646047254744536690271670910876261949955225261244047263841702619642748448019746397608975227830524045108286251241030379165637783808808994820359
? lindep([zeta(3),C,D])
%5 = [-5, 2, 6]~
    

ということで

が示唆されます．

具体例 2　積分

$\hspace{5pt}$単なる定積分では，intnum 関数を用います．

intnum 関数：intnum(x=a,b,f(x)) で$\quad\D\int_a^b f(x)\,dx\quad$を表します．

第一種完全楕円積分$K(x)$に関する式を調べてみました．

      ? \p300
   realprecision = 308 significant digits (300 digits displayed)
? A(n)=intnum(x=0,1,x^n*ellK(sqrt(1-x^2))^3)
%2 = (n)->intnum(x=0,1,x^n*ellK(sqrt(1-x^2))^3)
? lindep([Pi^4,A(0),A(2)])
%3 = [-1, 8, -40]~
    

ということで，

が示唆されます．

コマンドまとめ

出力

$\hspace{5pt}$"Shift + Enter" で出力します．

精度

$\hspace{5pt}$"\p200" で200桁の精度で計算します．数字を上げるほど精度が上がりますが，同時に処理速度は落ちます．

定義

数値を文字でおく場合は，そのまま "文字=・・・" と入力すればよいです．変数をもつ関数の場合は "f(x,y,...)=・・・"とすればよいです．f 以外の文字も使えますが，I（大文字のアイ）は使えません．

乗算

乗算する場合は必ず "$\ast$" を挟みます．でないとエラーが生じます．

代入

変数に数値を代入する場合は頭に "x=a;y=b;..." のように入力します．

特殊関数

三角関数：sin(x), asin(x), sinh(x), asinh(x) など
リーマンゼータ関数：zeta(x)
多重ゼータ関数：zetamult([a1,a2,...,ar])　　注 インデックスの向きは右から左です．
ガンマ関数：gamma(x)
第一種完全楕円積分：ellK(x)
超幾何級数：hypergeom([a1,a2,...,ap], [b1,b2,...,bq], x)
第一種ルジャンドル多項式：pollegendre(n,x)

数学定数他

無限大：-oo, +oo　（小文字のオーが二つ）
虚数単位：I（大文字のアイ）
円周率：Pi
オイラーの定数：Euler
カタランの定数：Catalan

級数

sumalt(n=a,f(n)) で，$\D\small \sum_{n=a}^\infty f(n)$

sum(n=a,b,f(n),0.) で，$\D\small \sum_{n=a}^b f(n)$　　注 f(n) の後ろに ",0." を付けます．

積分

intnum(x=a,b,f(x)) で，$\D\small \int_a^b f(x)\,dx$

出力結果の利用

"%数字" で出力された数値は，そのラベル "%数字" で使うことができます．