PARI/GPに入門する

こんにちは｜こんばんは。カエルのアイコンで活動しております @kyamaz です。

はじめに

名前のある素数を集めたエントリ『名前のついた素数コレクション 〜Haskellを添えて〜』を寄稿した際に、オンライン整数列大辞典（ https://oeis.org/ ）の参照リンクをつけしましたが、そのオンライン整数列大辞典には各項目に PROG という欄があり、その数列を生成するプログラムの記載があります。そのプログラムには幾つかのプログラミング言語での紹介がされています。特に、（PARI）が紹介されているケースが多く、 が好きな（Haskell）も紹介されていますが、ある項目には（Haskell）のプログラム例がないものもあります。一方で（PARI）のプログラム例が紹介されていないケースはほとんどありません。
そこで本稿では、PARIに入門してみようと思います。

PARI

PARIは、ボルドー大学（ https://www.u-bordeaux.fr/ ）のアンリ・コーエンのチームにより、数論に関する様々な演算を行うために開発された計算機代数アプリケーションです。現在は同大学の多くのボランティア開発者が管理、開発しています。

PARI は、高速な演算を行う関数を提供するC言語のライブラリです。PARIの機能を使うための対話型のコマンドラインインターフェイスである gp コマンドが付属しており、gp コマンドが解釈できるスクリプト言語を GP と呼びます。この基本ライブラリとスクリプト言語を合わせて、PARI/GP と表すようです。PARI は現在 GPLv3 でライセンスされており、マルチプラットフォームで動作します。
また、GPで書かれたスクリプトをC言語に変換し、コンパイルされた関数をgpコマンドで効率的に読み込むようなコンパイラl gp2c も提供されております。

PARI/GPの各種情報はボルドー大学のサイトから提供されており、日本語ページ（下記参照URL）も準備されております。

PARI/GP のインストール

PARI/GPはマルチプラットフォームで動作します。ソースコードからビルドする方法もありますが、既に動作可能なプログラムをDownloadサイト https://pari.math.u-bordeaux.fr/download.html からダウンロードしてインストールすることができます。なお、WSL含むLinux環境では apt を使い、macOS では brew を使って次のようにインストールも可能です。

apt でインストール

$ apt install pari-gp
$ apt install pari-gp2c # gp2c はPARI/GPとは別でインストールします

brew でインストール

$ brew install pari
# gp2c はインストールされません

PARI/GP を試してみる

上記でインストールした環境では、gp コマンドを使ってコマンドラインインターフェイスで計算を行うことができます。チュートリアルにある例を試してみましょう。

チュートリアル

% gp
                  GP/PARI CALCULATOR Version 2.17.0 (released)
         arm64 running darwin (aarch64/GMP-6.3.0 kernel) 64-bit version
      compiled: Sep 28 2024, Apple clang version 15.0.0 (clang-1500.3.9.4)
                    threading engine: pthread, nbthreads = 8
                 (readline v8.2 enabled, extended help enabled)
                     Copyright (C) 2000-2024 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?18 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 1048576, factorlimit = 1048576
? \\ -*-*-*- Tutorial Start -*-*-*-
? 57!
%1 = 40526919504877216755680601905432322134980384796226602145184481280000000000000
? 2 / 6
%2 = 1/3
? (1+I)^2
%3 = 2*I
? (x+1)^(-2)
%4 = 1/(x^2 + 2*x + 1)
? Mod(2,5)^3
%5 = Mod(3, 5)
? Mod(x, x^2+x+1)^3
%6 = Mod(1, x^2 + x + 1)
? w = ffgen([3,5],'w); w^12 \\ in F_3^5
%7 = 2*w^4 + 2*w^3 + 2
? Pi
%8 = 3.1415926535897932384626433832795028842
? log(2)
%9 = 0.69314718055994530941723212145817656807
? \p100
   realprecision = 115 significant digits (100 digits displayed)
? log(2)
%10 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875
? exp(%)
%11 = 2.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
? log(1+x)
%12 = x - 1/2*x^2 + 1/3*x^3 - 1/4*x^4 + 1/5*x^5 - 1/6*x^6 + 1/7*x^7 - 1/8*x^8 + 1/9*x^9 - 1/10*x^10 + 1/11*x^11 - 1/12*x^12 + 1/13*x^13 - 1/14*x^14 + 1/15*x^15 + O(x^16)
? exp(%12)
%13 = 1 + x + O(x^16)
? V = [1,2,3];
? W = [4,5,6]~;
? M = [1,2,3;4,5,6]
%16 =
[1 2 3]

[4 5 6]

? V*W
%17 = 32
? M*W
%18 = [32, 77]~
? \o0
   output = 0 (raw)
? factor(91)
%19 = [7,1;13,1]
? factor(x^4+4)
%20 = [x^2-2*x+2,1;x^2+2*x+2,1]
? factor((x^4+1)*Mod(1,a^2-2))
%21 = [x^2+Mod(-a,a^2-2)*x+1,1;x^2+Mod(a,a^2-2)*x+1,1]
? factor((x^4+4)*Mod(1,13))
%22 = [Mod(1,13)*x+Mod(4,13),1;Mod(1,13)*x+Mod(6,13),1;Mod(1,13)*x+Mod(7,13),1;Mod(1,13)*x+Mod(9,13),1]
? factor(x^4+1,Mod(1,a^2-2))
%23 = [x^2+Mod(-a,a^2-2)*x+1,1;x^2+Mod(a,a^2-2)*x+1,1]
? factor(x^4+1,Mod(1,13))
%24 = [Mod(1,13)*x^2+Mod(5,13),1;Mod(1,13)*x^2+Mod(8,13),1]
? intnum(x=0,1,1/(1+x^2))/Pi
%25 = 0.2500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
? sumnum(n=1,1/n^2)/Pi^2
%26 = 0.1666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667
? sumalt(n=0,(-1)^n/(2*n+1))*4
%27 = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
? sumalt(n=1,(-1)^n*log(n)) \\ diverging!
%28 = 0.2257913526447274323630976149474410717858973392775281586964715309893720739575656820888799716395355101
? 2*exp(2*%)
%29 = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
? [n^2|n<-[1..10]]
%30 = [1,4,9,16,25,36,49,64,81,100]
? [n^2|n<-[1..10],isprime(n)]
%31 = [4,9,25,49]
? [n^2|n<-primes([1,10])]
%32 = [4,9,25,49]
? [a,b] = [1,2];
? print("a=",a," b=",b)
a=1 b=2

さすが数学的な計算が得意そうです。数論に関する様々な演算を行うプログラム言語という感じがします。もう少し前出のオンライン整数列大辞典のプログラム例を実行してみましょう。

ピタゴラス素数（ https://oeis.org/A002144 ）

? select(p->p%4==1, primes(1000))
%1 = [5, 13, 17, 29, 37, 41, 53, 61, 73, 89, 97, 101, 109, 113, 137, 149, 157, 173, 181, 193, 197, 229, 233, 241, 257, 269, 277, 281, 293, 313, 317, 337, 349, 353, 373, 389, 397, 401, 409, 421, 433, 449, 457, 461, 509, 521, 541, 557, 569, 577, 593, 601, 613, 617, 641, 653, 661, 673, 677, 701, 709, 733, 757, 761, 769, 773, 797, 809, 821, 829, 853, 857, 877, 881, 929, 937, 941, 953, 977, 997, 1009, 1013, 1021, 1033, 1049, 1061, 1069, 1093, 1097, 1109, 1117, 1129, 1153, 1181, 1193, 1201, 1213, 1217, 1229, 1237, 1249, 1277, 1289, 1297, 1301, 1321, 1361, 1373, 1381, 1409, 1429, 1433, 1453, 1481, 1489, 1493, 1549, 1553, 1597, 1601, 1609, 1613, 1621, 1637, 1657, 1669, 1693, 1697, 1709, 1721, 1733, 1741, 1753, 1777, 1789, 1801, 1861, 1873, 1877, 1889, 1901, 1913, 1933, 1949, 1973, 1993, 1997, 2017, 2029, 2053, 2069, 2081, 2089, 2113, 2129, 2137, 2141, 2153, 2161, 2213, 2221, 2237, 2269, 2273, 2281, 2293, 2297, 2309, 2333, 2341, 2357, 2377, 2381, 2389, 2393, 2417, 2437, 2441, 2473, 2477, 2521, 2549, 2557, 2593, 2609, 2617, 2621, 2633, 2657, 2677, 2689, 2693, 2713, 2729, 2741, 2749, 2753, 2777, 2789, 2797, 2801, 2833, 2837, 2857, 2861, 2897, 2909, 2917, 2953, 2957, 2969, 3001, 3037, 3041, 3049, 3061, 3089, 3109, 3121, 3137, 3169, 3181, 3209, 3217, 3221, 3229, 3253, 3257, 3301, 3313, 3329, 3361, 3373, 3389, 3413, 3433, 3449, 3457, 3461, 3469, 3517, 3529, 3533, 3541, 3557, 3581, 3593, 3613, 3617, 3637, 3673, 3677, 3697, 3701, 3709, 3733, 3761, 3769, 3793, 3797, 3821, 3833, 3853, 3877, 3881, 3889, 3917, 3929, 3989, 4001, 4013, 4021, 4049, 4057, 4073, 4093, 4129, 4133, 4153, 4157, 4177, 4201, 4217, 4229, 4241, 4253, 4261, 4273, 4289, 4297, 4337, 4349, 4357, 4373, 4397, 4409, 4421, 4441, 4457, 4481, 4493, 4513, 4517, 4549, 4561, 4597, 4621, 4637, 4649, 4657, 4673, 4721, 4729, 4733, 4789, 4793, 4801, 4813, 4817, 4861, 4877, 4889, 4909, 4933, 4937, 4957, 4969, 4973, 4993, 5009, 5021, 5077, 5081, 5101, 5113, 5153, 5189, 5197, 5209, 5233, 5237, 5261, 5273, 5281, 5297, 5309, 5333, 5381, 5393, 5413, 5417, 5437, 5441, 5449, 5477, 5501, 5521, 5557, 5569, 5573, 5581, 5641, 5653, 5657, 5669, 5689, 5693, 5701, 5717, 5737, 5741, 5749, 5801, 5813, 5821, 5849, 5857, 5861, 5869, 5881, 5897, 5953, 5981, 6029, 6037, 6053, 6073, 6089, 6101, 6113, 6121, 6133, 6173, 6197, 6217, 6221, 6229, 6257, 6269, 6277, 6301, 6317, 6329, 6337, 6353, 6361, 6373, 6389, 6397, 6421, 6449, 6469, 6473, 6481, 6521, 6529, 6553, 6569, 6577, 6581, 6637, 6653, 6661, 6673, 6689, 6701, 6709, 6733, 6737, 6761, 6781, 6793, 6829, 6833, 6841, 6857, 6869, 6917, 6949, 6961, 6977, 6997, 7001, 7013, 7057, 7069, 7109, 7121, 7129, 7177, 7193, 7213, 7229, 7237, 7253, 7297, 7309, 7321, 7333, 7349, 7369, 7393, 7417, 7433, 7457, 7477, 7481, 7489, 7517, 7529, 7537, 7541, 7549, 7561, 7573, 7577, 7589, 7621, 7649, 7669, 7673, 7681, 7717, 7741, 7753, 7757, 7789, 7793, 7817, 7829, 7841, 7853, 7873, 7877, 7901]

ヴィーフェリッヒ素数（ https://oeis.org/A001220 ）

? N=10^4; default(primelimit, N); forprime(n=2, N, if(Mod(2, n^2)^(n-1)==1, print1(n, ", ")));
1093, 3511,

素数だけを対象としたループ処理（forprime）という記述があるのも数学的（数論的）な計算に特化したスクリプトとして適した処理系だと思われます。

PARI/GP は、多倍長精度の数値計算が可能です。GP のプロンプトで \p 10000 のような指定をすると、10000桁まで計算されて結果も10000桁出力されます。例えば、Piの出力が桁数を指定した前後で変わることが次のように確かめられます。

π の値の出力例

? Pi
%1 = 3.1415926535897932384626433832795028842
? \p100
   realprecision = 115 significant digits (100 digits displayed)
? Pi
%2 = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
? \p10000
   realprecision = 10018 significant digits (10000 digits displayed)
? Pi
%3 = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303820 ... 省略 ... 5525637568

Inside PARI （ソースコードを読む）

様々な計算が可能なライブラリですのでやはり実装が気になります。PARI/GP のソースコードは、ボルドー大学で運営されているリポジトリで管理されています。ソースコードを読むためには、リポジトリから git clone してください。

ソースコードの取得

$ git clone http://pari.math.u-bordeaux.fr/git/pari.git

ビルド方法は、INSTALLファイルに記載があります。autogen.shして ./Configureのあとmake all; make bench すれば、その環境でモジュールは作れます。いまは実装を確認したいのでソースコードを掘ってみます。
ソースコードは src ディレクトリの下にありますが、src/languageの下には言語仕様に関わる処理があり、src/basemath の下に基本とばる数学的な処理があります。少し高度な数学的処理は src/modules の下にあります。

Pi の実装を確認してみると、src/basemath/trans1.c に以下のようなコメントと共に実装があります。

Piの実装にあるコメント

/*                         ----
 *  53360 (640320)^(1/2)   \    (6n)! (545140134 n + 13591409)
 *  -------------------- = /    ------------------------------
 *        Pi               ----   (n!)^3 (3n)! (-640320)^(3n)
 *                         n>=0
 *
 * Ramanujan's formula + binary splitting */

以前の のエントリ『円周率を計算する（Haskell実装）』で紹介した手法である Binary Splitting Method（BS法） が使われていることが分かります。

おわりに

$\pi$ の実装の他にも、数学的な計算として参考になりそうな実装も多く、C言語の実装としてコードリーディングしてみると勉強になると思います。ご興味が湧いた方は、お時間のあるときにでもご覧頂ければと思います。

最後に、本稿を記載するために検証した環境を記しておきます。お手元の環境で検証する際に、動作が異なるときには参考になるかもしれません。

本稿の環境

本稿のために使用した環境は以下となります。

• macOS: Sonoma 14.6.1 (chip: Apple M1)
• Homebrew 4.4.15
• Homebrew LLVM version 12.0.1
• Homebrew clang version 12.0.1
• Homebrew bison (GNU Bison) 3.8.2

それでは、ご一読いただきまして有り難うございます。

（●）（●） Happy Hacking!
/""　＿＿""＼