pthread 対応の PARI/GP を入れる

# WSL2 の Fedora に PARI/GP を入れる

[WSL で Fedora を使う](http://var.blog.jp/archives/81457333.html) を見て Fedora をインストールし、WSL2 で動かし、[馬鹿の一つ覚え: Fedoraをアップグレード](https://qiita.com/glires/items/6bc71c0392c638661e13) で Fedora 32にアップグレードした。

しかし、することが思いつかなかったので PARI/GP を入れてみた。

```bash
% gp
Reading GPRC: /etc/gprc ...Done.

                                      GP/PARI CALCULATOR Version 2.11.4 (released)
                              amd64 running linux (x86-64/GMP-6.1.2 kernel) 64-bit version
                      compiled: Jun 29 2020, gcc version 10.1.1 20200507 (Red Hat 10.1.1-1) (GCC)
                                                threading engine: single
                                     (readline v8.0 enabled, extended help enabled)
```

なるほど、Debian だと 2.11.1 だけど、Fedora は最新の 2.11.4 が入るのか。さすが Fedora。(Debian はリリース後に大きな変更を加えないよう管理しているので、パッケージのパッチは当たってもバージョンアップはあまりされない。)
そして気づく。**threading engine: single** とはどういう意味だ？マルチスレッドで動くのか？ということで野良ビルドしてみた。

## コンパイルに必要なパッケージを入れる。(これは野良ビルド対象外)
gcc 系 (C, C++, Object-C, Fortran) m4, readline のライブラリをパッケージでインストール

```bash
% sudo yum install gcc.x86_64 gcc-c++.x86_64 gcc-objc++.x86_64 gcc-gfortran.i686 m4.x86_64 readline-devel.x86_64
```

## 必要なソースをダウンロード。

```bash
% mkdir ~/opt ; cd ~/opt
% wget http://ftp.jaist.ac.jp/pub/GNU/gmp/gmp-6.2.0.tar.bz2
% wget http://ftp.jaist.ac.jp/pub/GNU/mpfr/mpfr-4.1.0.tar.bz2
% wget https://pari.math.u-bordeaux.fr/pub/pari/unix/pari-2.11.4.tar.gz
```
PARI/GP の追加パッケージが必要な場合は [CentOS 8.1におけるPARI/GP 2.11.4のインストール手順](https://yoshiiz.blog.fc2.com/blog-entry-937.html) を参照。

## gmp をコンパイル
依存関係から gmp を最初にビルド
CPU が Ryzen5 3500U なので最適化を入れておく CFLAGS="-m64 -mtune=znver1 -march=znver1"
最適化をしたいけど、何を指定すればよいか分からない場合は -mtune=native -march=native で、自動的にコンパイラが最適と思われるものを選んでくれます。(他のマシンで動かす場合には指定してはいけません。)

```bash
% export CFLAGS="-m64 -mtune=znver1 -march=znver1"
% tar jxf gmp-6.2.0.tar.bz2 ; cd gmp-6.2.0
% ./configure --prefix=/usr/local --enable-cxx
% make -j8 && make -j8 check
% sudo make install
```
どうやら、gmp にはマルチスレッドは関係なかった模様。

## mpfr をコンパイル
次は mpfr をビルド

```bash
% tar jxf mpfr-4.1.0.tar.bz2 && cd mpfr-4.1.0
% ./configure --prefix=/usr/local --with-gmp=/usr/local --enable-thread-safe
% make -j8 && make -j8 check && make check-gmp-symbols
% sudo make install
```
こいつは意味があった模様。

## PARI/GP をコンパイル
そして、目的の PARO/GP をビルド & インストール。
設定ファイルの場所を /usr/local/etc/gprc に変更。(--prefix が効くようにするのが良いのでしょうが、決め打ちで)
Configure の --time=ftime オプションは、経過時間ではなくて CPU 使用時間を報告するようにする。

```bash
% tar jxf pari-2.11.4.tar.gz && cd pari-2.11.4
% vi src/language/gplib.c
% diff -u src/language/gplib.c.orig src/language/gplib.c
--- src/language/gplib.c.orig   2020-07-19 06:50:12.830000000 +0900
+++ src/language/gplib.c        2020-07-19 06:51:16.580000000 +0900
@@ -804,7 +804,7 @@
     f = gprc_chk(str); /* in $HOME */
     if (!f) f = gprc_chk(s); /* in . */
 #ifndef _WIN32
-    if (!f) f = gprc_chk("/etc/gprc");
+    if (!f) f = gprc_chk("/usr/local/etc/gprc");
 #else
     if (!f)  /* in basedir */
     {

% ./Configure --prefix=/usr/local --with-gmp=/usr/local --mt=pthread --time=ftime
% make -j8 all

tex が無いとこけるが、ドキュメントは web などを見ることにしてここでは無視。

% make test-all
...
+++ Total bench for gp-sta is 332867
+++ Total bench for gp-dyn is 436899
The following tests were skipped: ellglobalred ellmodulareqn galois galoischartable galpol lfunartin member programming
...
%  make test-parallel
Making test-parallel in Olinux-x86_64
make[1]: Entering directory '/root/opt/pari-2.11.4/Olinux-x86_64'
* Testing parallel         gp-sta..TIME=31620   gp-dyn..TIME=40839
+++ Total bench for gp-sta is 31620
+++ Total bench for gp-dyn is 40839
make[1]: Leaving directory '/root/opt/pari-2.11.4/Olinux-x86_64'

% make test-rnfkummer
Making test-rnfkummer in Olinux-x86_64
make[1]: Entering directory '/root/opt/pari-2.11.4/Olinux-x86_64'
* Testing rnfkummer        gp-sta..TIME= 5745   gp-dyn..TIME= 7685
+++ Total bench for gp-sta is 5745
+++ Total bench for gp-dyn is 7685
make[1]: Leaving directory '/root/opt/pari-2.11.4/Olinux-x86_64'

% sudo make install
% sudo cp misc/gprc.dft /usr/local/etc/gprc
% sudo cp Olinux-x86_64/gp-sta /usr/local/bin/gp-sta-2.11
% sudo strip /usr/local/bin/gp-sta-2.11
```
parigmp へスタティックリンクしている gp もインストールしておく。

## ベンチマークを実行
make ファイルに、ベンチマークがあるらしいので走らせてみる。(5度実行した最速の値)

```bash
% make bench
Making bench in Olinux-x86_64
make[1]: Entering directory '/root/opt/pari-2.11.4/Olinux-x86_64'
* Testing objets           gp-sta..TIME=    1   gp-dyn..TIME=    1
* Testing analyz           gp-sta..TIME=   10   gp-dyn..TIME=   20
* Testing number           gp-sta..TIME=   11   gp-dyn..TIME=   14
* Testing polyser          gp-sta..TIME=    3   gp-dyn..TIME=    4
* Testing linear           gp-sta..TIME=    4   gp-dyn..TIME=    9
* Testing elliptic         gp-sta..TIME=   13   gp-dyn..TIME=   22
* Testing sumiter          gp-sta..TIME=    7   gp-dyn..TIME=    9
* Testing graph            gp-sta..TIME=    4   gp-dyn..TIME=    6
* Testing program          gp-sta..TIME=    3   gp-dyn..TIME=    4
* Testing trans            gp-sta..TIME=   10   gp-dyn..TIME=   13
* Testing nfields          gp-sta..TIME=   96   gp-dyn..TIME=  121
+++ Total bench for gp-sta is 85
+++ Total bench for gp-dyn is 126
make[1]: Leaving directory '/root/opt/pari-2.11.4/Olinux-x86_64'
%
```

試しに -pthred なしにコンパイルしてみると

```bash
% ./Configure --prefix=/usr/local --with-gmp=/usr/local --time=ftime
% make -j88 all
% make bench
Making bench in Olinux-x86_64
make[1]: Entering directory '/root/opt/pari-2.11.4/Olinux-x86_64'
* Testing objets           gp-sta..TIME=    1   gp-dyn..TIME=    1
* Testing analyz           gp-sta..TIME=   13   gp-dyn..TIME=   14
* Testing number           gp-sta..TIME=   12   gp-dyn..TIME=   14
* Testing polyser          gp-sta..TIME=    4   gp-dyn..TIME=    4
* Testing linear           gp-sta..TIME=    4   gp-dyn..TIME=    5
* Testing elliptic         gp-sta..TIME=   14   gp-dyn..TIME=   13
* Testing sumiter          gp-sta..TIME=    5   gp-dyn..TIME=    6
* Testing graph            gp-sta..TIME=    3   gp-dyn..TIME=    4
* Testing program          gp-sta..TIME=    3   gp-dyn..TIME=    3
* Testing trans            gp-sta..TIME=   11   gp-dyn..TIME=   12
* Testing nfields          gp-sta..TIME=   50   gp-dyn..TIME=   47
+++ Total bench for gp-sta is 80
+++ Total bench for gp-dyn is 85
```
うーん、-pthred  無しの方が若干早いかな？
実行時間は cpu 時間より実時間の方が良いかな。ということで --time=ftime なしでコンパイルして、インストールしなおし。

# 実行時間の比較
せっかくビルドしたので使ってみる。最初に -pthread 付きでコンパイルした gp、次に yum で入れた gp で測定した。

- $ \displaystyle\sum_{i=1}^{100000}\frac1n $ を 1000桁の精度で実行し、実行時間を測定

- 時間が短かったので [Sums, products, integrals and similar functions](https://pari.math.u-bordeaux.fr/dochtml/html/Sums__products__integrals_and_similar_functions.html) から少し時間がかかりそうな式
? \pb 4000
? exponent(sumpos(i = 1, 1 / i^2) - zeta(2))
の実行時間を測定

両方の時間を比較してみた。

```bash
% echo 'parisizemax = 1024M' >> ~/.gprc
% gp
Reading GPRC: /root/.gprc ...Done.

                                      GP/PARI CALCULATOR Version 2.11.4 (released)
                              amd64 running linux (x86-64/GMP-6.2.0 kernel) 64-bit version
                      compiled: Jul 19 2020, gcc version 10.1.1 20200507 (Red Hat 10.1.1-1) (GCC)
                                               threading engine: pthread
                                     (readline v8.0 enabled, extended help enabled)

                                         Copyright (C) 2000-2018 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 ?17 for how to get moral (and possibly technical) support.

parisizemax = 1024000000, primelimit = 500000, nbthreads = 8
? \p 1000
   realprecision = 1001 significant digits (1000 digits displayed)
? default(timer,1)
? sum(i=1, 100000, 1/i);
time = 1,178 ms.
? sum(i=1, 100000, 1/i);
time = 1,152 ms.
? sum(i=1, 100000, 1/i);
time = 1,200 ms.
? sum(i=1, 100000, 1/i);
time = 1,155 ms.
? sum(i=1, 100000, 1/i);
time = 1,139 ms.

? \p 10
   realprecision = 19 significant digits (10 digits displayed)
? ( 1178 + 1152 + 1200 + 1155 +1139 ) / 5.0
%7 = 1164.800000
?
? \pb 4000
   realbitprecision = 4000 significant bits (1204 decimal digits displayed)
? exponent(sumpos(i = 1, 1 / i^2) - zeta(2))
  ***   Warning: increasing stack size to 16000000.
time = 14,560 ms.
%8 = -4031

$ /usr/bin/gp
Reading GPRC: /root/.gprc ...Done.

                                      GP/PARI CALCULATOR Version 2.11.4 (released)
                              amd64 running linux (x86-64/GMP-6.1.2 kernel) 64-bit version
                      compiled: Jun 29 2020, gcc version 10.1.1 20200507 (Red Hat 10.1.1-1) (GCC)
                                                threading engine: single
                                     (readline v8.0 enabled, extended help enabled)

                                         Copyright (C) 2000-2018 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 ?17 for how to get moral (and possibly technical) support.

parisizemax = 1024000000, primelimit = 500000
? \p 1000
   realprecision = 1001 significant digits (1000 digits displayed)
? default(timer,1)
? sum(i=1, 100000, 1/i);
time = 1,397 ms.
? sum(i=1, 100000, 1/i);
time = 1,336 ms.
? sum(i=1, 100000, 1/i);
time = 1,336 ms.
? sum(i=1, 100000, 1/i);
time = 1,370 ms.
? sum(i=1, 100000, 1/i);
time = 1,379 ms.

それぞれ 5回の平均の計算時間は？ (1000桁は長いので10桁ぐらいで)
? \p 10
   realprecision = 19 significant digits (10 digits displayed)
? ( 1178 + 1152 + 1200 + 1155 + 1139 ) / 5.0
%7 = 1164.800000
? ( 1397 + 1336 + 1336 + 1370 + 1379) / 5.0
%8 = 1363.600000
? 100 * (1363.6 - 1164.8) / 1363.6
%9 = 14.57905544

? \pb 4000
   realbitprecision = 4000 significant bits (1204 decimal digits displayed)
? exponent(sumpos(i = 1, 1 / i^2) - zeta(2))
  ***   Warning: increasing stack size to 16000000.
time = 15,918 ms.
%10 = -4031

? 100.0 * ( 15918 - 14560 ) / 15918
%11 = 8.531222515
```
間違えてシングルスレッドの gp も作ってしまったのでこちらでも測定。

```bash
$ /usr/local/bin/gp
Reading GPRC: /root/.gprc ...Done.

                                      GP/PARI CALCULATOR Version 2.11.4 (released)
                              amd64 running linux (x86-64/GMP-6.2.0 kernel) 64-bit version
                      compiled: Jul 19 2020, gcc version 10.1.1 20200507 (Red Hat 10.1.1-1) (GCC)
                                                threading engine: single
                                     (readline v8.0 enabled, extended help enabled)

                                         Copyright (C) 2000-2018 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 ?17 for how to get moral (and possibly technical) support.

parisizemax = 1024000000, primelimit = 500000
? \p 1000
   realprecision = 1001 significant digits (1000 digits displayed)
? default(timer,1)
? sum(i=1, 100000, 1/i);
time = 1,128 ms.
? sum(i=1, 100000, 1/i);
time = 1,129 ms.
? sum(i=1, 100000, 1/i);
time = 1,139 ms.
? sum(i=1, 100000, 1/i);
time = 1,143 ms.
? sum(i=1, 100000, 1/i);
time = 1,143 ms.
?
? \p 10
   realprecision = 19 significant digits (10 digits displayed)
? ( 1128 + 1129 + 1139 + 1143 + 1143 ) / 5.0
%7 = 1136.400000

? \p 10
   realprecision = 19 significant digits (10 digits displayed)
? ( 1128 + 1129 + 1139 + 1143 + 1143 ) / 5.0
%7 = 1136.400000
? \pb 4000
   realbitprecision = 4000 significant bits (1204 decimal digits displayed)
? exponent(sumpos(i = 1, 1 / i^2) - zeta(2))
  ***   Warning: increasing stack size to 16000000.
time = 13,701 ms.
%8 = -4031
```

実行時間

|                                             コマンド | yum パッケージ | 野良 pthred | 野良 single |
|:----------------------------------------------------|:----------:|:--------:|:-----------:|
| ? \p 1000<br>? sum(i=1, 100000, 1/i)<br>5回の平均        |  1363.6 ms | 1164.8 ms<br>(-14.6%)| 1136.4 ms<br>(-16.6%)|
| ? \pb 4000<br>? exponent(sumpos(i = 1, 1 / i^2)-zeta(2))|  15,918 ms | 14,560 ms<br>(-8.5%)| 13,701 ms<br>(-13.9%)|
どうやらマルチスレッドは効いてないけど CPU の最適化で 8～14% ぐらい高速化したようです。マルチスレッドで動くのって何かあるのかな？

PARI/GP のドキュメント [Introduction to parallel GP](https://pari.math.u-bordeaux.fr/pub/pari/manuals/2.12.1/parallel.pdf) によると、--mt=pthread を指定するとダイナミックリンクで 25% の速度低下、スタティックリンクで 5% の速度低下を起こすそうです。ドキュメント通り pthred より、single thred の方が早いですね。

isprime(N) N が素数かどうかを判定する関数がマルチスレッドで動くみたいですね。
prime(N) (N より大きな一番近い素数を返す) は、シングルスレッドですねぇ。

```gp:$i&nbsp;でスレッド数を指定して実行
default(timer,1)
default(nbthreads,$i)
now_nbthreads
default(nbthreads)
isprime(2^600+187);
```
を printf 文に組み込んで、スレッド数を変えながらループさせてみる。

```bash
$ for i in 8 4 2 1
do
  printf "default(timer,1)\ndefault(nbthreads,$i)\nnow_nbthreads\ndefault(nbthreads)\nisprime(2^600+187);\n" | ( time /usr/local/bin/gp -q )
  echo
done
now_nbthreads
8
time = 1,223 ms.

real    0m0.200s
user    0m1.230s
sys     0m0.011s

now_nbthreads
4
time = 930 ms.

real    0m0.254s
user    0m0.936s
sys     0m0.010s

now_nbthreads
2
time = 804 ms.

real    0m0.411s
user    0m0.809s
sys     0m0.000s

now_nbthreads
1
time = 659 ms.

real    0m0.670s
user    0m0.660s
sys     0m0.010s

%
```

|スレッド数 | 経過時間 | CPU使用時間 | システム時間 |
|---------:|--------:|-----------:|--------------:|
|         8|  0.200s |  1.230s | 0.011s |
|         4|  0.254s |  0.936s | 0.010s |
|         2|  0.411s |  0.809s | 0.000s |
|         1|  0.670s |  0.660s | 0.010s |

 - 経過時間 (real) : 実際に時計の針が進んだ時間。
 - CPU使用時間 (user) : 各スレッドが使った時間を合計したもの。 4スレッドで 1秒間、フルで CPU を使っていたら 1秒× 4スレッド = 4秒
 - システム時間 (sys) : I/O待ち、タスク切り替えの時間などのために OS が使った時間

スレッドが増えると CPU 時間も増えるのが、マルチスレッドの効率が悪そうだねぇとなるところ。でも、経過時間はちゃんと減ってるので効いてますね。


最終的に全部入れてしまった😅

```bash:
% ls -l /usr/local/bin/gp*
lrwxrwxrwx 1 root root       21 Jul 20 11:41 /usr/local/bin/gp -> gp-2.11-single-static
-rwxr-xr-x 1 root root    78240 Jul 19 20:17 /usr/local/bin/gp-2.11-pthread-dyn
-rwxr-xr-x 1 root root 10171624 Jul 19 20:26 /usr/local/bin/gp-2.11-pthread-static
-rwxr-xr-x 1 root root  9602272 Jul 20 11:41 /usr/local/bin/gp-2.11-single-static
-rwxr-xr-x 1 root root    37114 Jul 19 20:17 /usr/local/bin/gphelp
```

# マチンの式で円周率を求めてみる

$\displaystyle \frac {10000}{log_{10}5}=14306.7655807...$ と $\displaystyle \frac {10000}{log_{10}239}= 4204.51094243...$ で1万桁ぐらいまで円周率を求めて比較してみる。

```pg:
% /usr/local/bin/gp-2.11-pthread-static
? default(timer,1)
? \p 50000
? 4*(atan(1) - (parsum(n=1,14306,(4/(4*n-3) - 4/((4*n-1)*5^2))/5^(4*n-3)) - parsum(n=1,4204,(1/(4*n-3) - 1/((4*n-1)*239^2))/239^(4*n-3))))
time = 47,955 ms.
%5 = -6.48140576723279.... E-40002

? 4*(atan(1) - (sum(n=1,14306,(4/(4*n-3) - 4/((4*n-1)*5^2))/5^(4*n-3)) - sum(n=1,4204,(1/(4*n-3) - 1/((4*n-1)*239^2))/239^(4*n-3))))
time = 46,807 ms.
%5 = -6.48140576723279.... E-40002

? 4*(atan(1) - (sum(n=1,14306,(4/(4*n-3) - 4/((4*n-1)*5^2))/5^(4*n-3) - (1/(4*n-3) - 1/((4*n-1)*239^2))/239^(4*n-3))))
time = 5min, 15,332 ms.
%7 = 7.4302477881593534352... E-40003
```
あら？

任意精度のライブラリとして Apfloat なんてのが使いやすいと書いてあるところがありましたが、さて？ pthred も使えるそうな。