GMP の使い方

1. GMP とは

The GNU MP Bignum Library 略して GMP は多倍長計算を非常に高速に行うライブラリです。アセンブリ言語を少しかじったことがある程度では GMP の計算速度には太刀打ちできないでしょう。あえて欠点をあげれば Linux ライクな開発環境がないと make できないということでしょうか。GMP は GMP-ECM や GGNFS などの素因数分解ソフトウェアで利用されています。もちろん、このライブラリを呼び出すプログラムを自分で書くこともできます。

2. Cygwin のインストール

Cygwin は Windows に Linux ライクな環境を構築するためのソフトウェア群です。GMP をインストールするために Linux ライクな環境が必要なので、先に Cygwin をインストールしておきます。

メモ : Cygwin のインストールには時間がかかります。休日などに余裕を持って行いましょう。

2.1. オペレーティングシステムの種類を確認する

「コントロールパネル」→「システム」で「システムの種類」が「32 ビットオペレーティングシステム」（以下、32bit OS）または「64 ビットオペレーティングシステム」（以下、64bit OS）のどちらであるかを確認します。

2.2. Cygwin をインストールする場所を決める

十分な空き容量のあるドライブを選び、パッケージをダウンロードするためのディレクトリ (ローカルパッケージディレクトリ) とインストール先のディレクトリ (ルートディレクトリ) を決めます。ここではローカルパッケージディレクトリ C:\cygwinpack64\ にダウンロードしてルートディレクトリ C:\cygwin64\ にインストールすることにします。自分の選んだ場所に応じて適宜読み替えてください。

2.3. setup-x86.exe または setup-x86_64.exe をダウンロードする

Cygwin は数多くのパッケージで構成されており、32bit OS のときは setup-x86.exe、64bit OS のときは setup-x86_64.exe が専用のインストーラになっています。オペレーティングシステムの種類に合っているほう（以下、setup-*.exe）をローカルパッケージディレクトリ C:\cygwinpack64\ にダウンロードします。

補足 : 64bit OS 用のインストーラが公開される以前のインストーラは setup.exe という名前でした。

注意 : セキュリティ関連のソフトウェアを使用している場合は setup-*.exe がインターネットにアクセスできるように設定しておきます。

2.4. setup-*.exe を実行する

ローカルパッケージディレクトリ C:\cygwinpack64\ にダウンロードした setup-*.exe を実行します。

「次へ」。

2.5. 作業内容を選択する

Choose A Download Source の選択肢は次の 3 つです。

「Install from Internet」 … ダウンロードしてインストールする
「Download Without Installing」 … 後でインストールできるようにダウンロードだけ行う
「Install from Local Directory」 … ダウンロードしておいたパッケージをインストールする

ダウンロードしてインストールするので「Install from Internet」を選択します。

「次へ」。

2.6. ルートディレクトリを設定する

Select Root Install Directory では Root Directory にルートディレクトリ C:\cygwin64\ を入力し、Install For は「All Users (RECOMMENDED)」を選択します。

「次へ」。

2.7. ローカルパッケージディレクトリを設定する

Select Local Package Directory では Local Package Directory にローカルパッケージディレクトリ C:\cygwinpack64\ を入力します。

「次へ」。

2.8. ダウンロードに用いるインターネットへの接続方法を設定する

Select Your Internet Connection では「Use System Proxy Settings」または「Direct Connection」を選択します。環境によっては他の接続方法を選択する必要があるかも知れません。

「次へ」。

ミラーサイトの一覧がダウンロードされます。

2.9. パッケージのダウンロード元のサイトを選択する

Choose A Download Site では Available Download Sites でパッケージのダウンロード元のサイトを選択します。利用できることがわかっているミラーサイトがリストになければ User URL に入力して「Add」ボタンでリストに追加します。私は「https://ftp.iij.ad.jp」を利用しています。

注意 : これを書いている時点で「ftp://ftp.iij.ad.jp」は Cygwin のミラーサイトとして利用できません。setup-*.exe を攻撃と誤認するようになったのか、ダウンロードが中途半端に失敗して復旧が困難になります。

「次へ」。

ローカルパッケージディレクトリ C:\cygwinpack64\ の下に https%3a%2f%2fftp.iij.ad.jp%2fpub%2fcygwin%2f のような名前のディレクトリができます。

選択したサイトからパッケージの一覧がダウンロードされます。

2.10. パッケージを選択してインストールする

Select Packages でインストールするパッケージを選択します。デフォルトでインストールされるもの以外に GMP を利用するにあたって必要なパッケージは開発関係の Devel と数学関係の Math に分類されているものです。また、GMP のアーカイブを展開するために xz を使用します。

メモ : ウインドウの右下をつまんでドラッグしてウインドウを広くしておくと作業しやすいでしょう。

View を Category にします。

All は「○Default」のままにしておきます。

(1) Devel の「○ ～」の部分を何度かクリックして「○Install」に変えます。

(2) Math の「○ ～」の部分を何度かクリックして「○Install」に変えます。

(4)「次へ」。

Resolving Dependencies の画面が出るときは RECOMMENDED に従います。

パッケージのダウンロードとインストールが始まります。時間がかかるので終わるまで放置しておきます。

注意 : Progress の画面のプログレスバーはかなり大雑把です。特に TeX のメタフォントの処理中は止まってしまったように見えますが中断しないでください。

注意 : インストールが完了する前にインストールされた新しいコマンドがネットワークを探します。セキュリティ関連のソフトウェアが原因でネットワークへの接続が正常に処理されないとインストールを完了できなくなる場合があるかも知れません。

注意 : パッケージをアンインストールするときは Select Packages で「○Uninstall」を選択します。初心者はこれ以外の方法で削除しないほうが無難です。

補足 : All を「○Install」に変えてフルインストールする場合も Debug だけは「○Default」または「○Uninstall」にすることをお勧めします。

2.11. setup-*.exe を終了する

Cygwin のアイコンをデスクトップに置くときは「Create icon on Desktop」をチェックします。また、スタートメニューに登録するときは「Add icon to Start Menu」をチェックします。

「完了」。

setup-*.exe が終了します。

2.12. 環境変数 Path に登録する

私はテキストエディタとして xyzzy を利用しているので、Cygwin64 Termilal を起動しなくても *Shell* バッファで自作のプログラムをコンパイルできるように環境変数 Path に Cygwin の bin ディレクトリを追加しています。具体的には、「コントロールパネル」→「システム」→「システムの詳細設定」→「環境変数」（OS によって多少経路が異なります）にあるシステム環境変数の Path の行を選択して「編集」をクリックし、変数値の末尾に以下を追加します。

;C:\cygwin64\usr\local\bin;C:\cygwin64\usr\local\lib;C:\cygwin64\bin

注意 : Path の変数値は表示領域の何倍もの長さがあり、編集を開始した時点ですべて選択されているので、うっかり間違ったキーを押すと他のアプリケーションで必要な部分に何が書いてあったのかわからないまま全部消してしまう恐れがあります。間違えたと思ったときは「キャンセル」ボタンをクリックして編集をやり直します。

補足 : 変数値をすべて選択・コピーして xyzzy の *scratch* バッファなどの文字列の全体が見えるところに貼り付け、そこで Cygwin の bin ディレクトリを追加してから再びコピーして変数値にすべて選択・貼り付け、という手順で書き換えると少し安心できます。気休め程度ですが、余計な緊張はしないに越したことはありません。

補足 : 多少不便とはいえ、環境変数を編集する方法が標準で提供されている以上、この目的のためだけに信用してよいかわからないツールを新たに導入することはお勧めしません。

環境変数を登録した後、xyzzy を開き直せばすぐに *Shell* バッファでコンパイルできるようになります。

2.13. Cygwin の更新

ときどき (月に 1 回くらい？) setup-*.exe を動かして Cygwin を最新版に更新します。手順はインストールと同じで、更新されたパッケージだけが自動的に選択されます。setup-*.exe 自体も更新されることがあるので setup-x86.exe または setup-x86_64.exe をダウンロードするところから行います。

3. GMP のインストール

これを書いている時点で GMP の最新版は 6.3.0 です。Cygwin にはもともとコンパイル済みの GMP が含まれていますが、バージョンが古かったり、特に 32bit OS では実行環境に合った最適化が施されていないため十分なパフォーマンスを発揮できない場合があるので、最新版をインストールします。

注意 : GMP の 4.2.1 までのバージョンには GGNFS を誤動作させるバグ (`mpz_set_d' causes segmentation faults on parameters less than 2^-32.) があります。4.2.1 までのバージョンをインストールしている人は最新版に更新してください。

以下の作業は Cygwin64 Terminal で行います。

3.1. GMP の古いバージョンをアンインストールする

古いバージョンがインストールされているときは新しいバージョンをインストールする前に古いバージョンをアンインストールしておきましょう。

~> cd ~/gmp-6.2.1
~/gmp-6.2.1> make uninstall
~/gmp-6.2.1> cd ~
~>

3.2. GMP のソースコードをダウンロードする

GMP はソースコードの形で配布されています。The GNU MP Bignum Library からソースコードのアーカイブ gmp-6.3.0.tar.xz を C:\cygwin64\home\ユーザ名 にダウンロードします。

3.3. GMP をインストールする

アーカイブを展開します。アーカイブの拡張子が .tar.xz なので tar コマンドにオプション J を指定します。このとき xz コマンドがインストールされている必要があります。

~> tar Jxf gmp-6.3.0.tar.xz
~> cd ~/gmp-6.3.0
~/gmp-6.3.0>

./configure を実行します。

~/gmp-6.3.0> ./configure

make します。大きいライブラリなので少し時間がかかります。

~/gmp-6.3.0> make

make check してエラーが出ないことを確かめます。少し時間がかかります。

~/gmp-6.3.0> make check

make install します。

~/gmp-6.3.0> make install

これで GMP 本体のインストールは完了です。

3.4. Perl のモジュールをインストールする

ついでに Perl のモジュールもインストールしておきましょう。Perl で多倍長計算を手軽に利用できるようになります。標準の Math::BigInt モジュールよりも圧倒的に高速です。

~/gmp-6.3.0> cd demos/perl
~/gmp-6.3.0/demos/perl> perl Makefile.PL GMP_BUILDDIR=../..
~/gmp-6.3.0/demos/perl> make
~/gmp-6.3.0/demos/perl> make test
~/gmp-6.3.0/demos/perl> make install
~/gmp-6.3.0/demos/perl> cd ~
~>

3.5. GMP-ECM を再インストールする

GMP と GMP-ECM がインストールされている状態で GMP を更新したときは、GMP-ECM の使い方を参照して GMP-ECM を再インストールします。

メモ : バージョンの組み合わせにもよるかも知れませんが、GMP と GMP-ECM がインストールされている状態で GMP だけを更新して libecm を使用するプログラムをコンパイルすると ecm_factor () が正常に動作しなくなる場合があるようです。手元の環境では、gmp-5.0.5 + ecm-4.6.3 から gmp-5.1.1 に更新したときに gmp-5.0.5 でコンパイルしたプログラムではすぐに見つかる素因数が gmp-5.1.1 でコンパイルしたプログラムでは何度繰り返しても発見できなくなり、ecm-4.6.4 をインストールしてプログラムを再コンパイルしたところ再び発見できるようになりました。

4. 使用例

4.1. C で GMP を使う

GNU MP 6.3.0: Top を参照して C でプログラムを書いてみましょう。ここでは大浦拓哉さんが FFT と AGM による円周率計算プログラムで示されている円周率計算のアルゴリズムをそのまま GMP で実装してみました。大浦さんに感謝します。

mpf_pi.c

/*
  Function: void mpf_pi (mpf_t rop)
    Set the value of pi.

  Reference:
    https://www.kurims.kyoto-u.ac.jp/~ooura/pi_fft-j.html

  ---- a formula based on the AGM (Arithmetic-Geometric Mean) ----
    c = sqrt (0.125);
    a = 1 + 3 * c;
    b = sqrt (a);
    e = b - 0.625;
    b = 2 * b;
    c = e - c;
    a = a + e;
    npow = 4;
    do {
        npow = 2 * npow;
        e = (a + b) / 2;
        b = sqrt (a * b);
        e = e - b;
        b = 2 * b;
        c = c - e;
        a = e + b;
    } while (e > SQRT_SQRT_EPSILON);
    e = e * e / 4;
    a = a + b;
    pi = (a * a - e - e / 2) / (a * c - e) / npow;
  ---- modification ----
    This is a modified version of Gauss-Legendre formula
    (by T.Ooura). It is faster than original version.
*/

#include <gmp.h>

void mpf_pi (mpf_t rop);

void mpf_pi (mpf_t a) {
  unsigned long int prec;
  mpf_t SQRT_SQRT_EPSILON, c, b, e;
  unsigned long int log2_npow;

  prec = mpf_get_prec (a);
  mpf_init2 (SQRT_SQRT_EPSILON, prec);
  mpf_init2 (c, prec);
  mpf_init2 (b, prec);
  mpf_init2 (e, prec);

  mpf_set_ui (SQRT_SQRT_EPSILON, 1);
  mpf_div_2exp (SQRT_SQRT_EPSILON, SQRT_SQRT_EPSILON, prec >> 2);

  mpf_set_d (c, 0.125);     /* c = sqrt (0.125); */
  mpf_sqrt (c, c);
  mpf_mul_ui (a, c, 3);     /* a = 1 + 3 * c; */
  mpf_add_ui (a, a, 1);
  mpf_sqrt (b, a);          /* b = sqrt (a); */
  mpf_set_d (e, 0.625);     /* e = b - 0.625; */
  mpf_sub (e, b, e);
  mpf_add (b, b, b);        /* b = 2 * b; */
  mpf_sub (c, e, c);        /* c = e - c; */
  mpf_add (a, a, e);        /* a = a + e; */
  log2_npow = 2;            /* npow = 4; */
  do {
    log2_npow++;            /* npow = 2 * npow; */
    mpf_add (e, a, b);      /* e = (a + b) / 2; */
    mpf_div_2exp (e, e, 1);
    mpf_mul (b, a, b);      /* b = sqrt (a * b); */
    mpf_sqrt (b, b);
    mpf_sub (e, e, b);      /* e = e - b; */
    mpf_add (b, b, b);      /* b = 2 * b; */
    mpf_sub (c, c, e);      /* c = c - e; */
    mpf_add (a, e, b);      /* a = e + b; */
  } while (mpf_cmp (e, SQRT_SQRT_EPSILON) > 0);  /* e > SQRT_SQRT_EPSILON */
  mpf_mul (e, e, e);        /* e = e * e / 4; */
  mpf_div_2exp (e, e, 2);
  mpf_add (a, a, b);        /* a = a + b; */
  mpf_mul (c, c, a);        /* pi = (a * a - e - e / 2) / (a * c - e) / npow; */
  mpf_sub (c, c, e);
  mpf_mul (a, a, a);
  mpf_sub (a, a, e);
  mpf_div_2exp (e, e, 1);
  mpf_sub (a, a, e);
  mpf_div (a, a, c);
  mpf_div_2exp (a, a, log2_npow);

  mpf_clear (e);
  mpf_clear (b);
  mpf_clear (c);
  mpf_clear (SQRT_SQRT_EPSILON);
}

test_mpf_pi.c

/*
  test_mpf_pi precision(10-) output-flag(0-1)

  Compile:
    gcc -I/usr/local/include -L/usr/local/lib -Wall -O3 -o test_mpf_pi test_mpf_pi.c mpf_pi.c -lgmp
*/

#include <stdio.h>
#include <sys/times.h>
#include <unistd.h>
#include <gmp.h>

#define LOG_2_10 3.32192809488736234787031942949

void mpf_pi (mpf_t rop);

int main (int argc, char *argv[]) {
  int error = 0, output = 1;
  unsigned long int prec10 = 100;
  mpf_t pi;
  struct tms t_buf;
  clock_t start, stop;

  switch (argc) {
  case 3:
    error |= sscanf (argv[2], "%d", &output) != 1;
  case 2:
    error |= sscanf (argv[1], "%lu", &prec10) != 1;
    if (prec10 < 10) {
      prec10 = 10;
    }
    break;
  default:
    error = 1;
  }

  if (error) {
    printf ("usage: %s precision(10-) output-flag(0-1)\n", argv[0]);
    return error;
  }

  mpf_init2 (pi, (int) ((double) (2 + prec10) * LOG_2_10) + 1);

  stop = times (&t_buf);
  while ((start = times (&t_buf)) == stop);
  mpf_pi (pi);
  stop = times (&t_buf);
  printf ("%.3f sec.\n", (double) (stop - start) / (double) sysconf (_SC_CLK_TCK));
  if (output) {
    gmp_printf ("%.*Ff\n", prec10, pi);
    stop = times (&t_buf);
    printf ("%.3f sec.\n", (double) (stop - start) / (double) sysconf (_SC_CLK_TCK));
  }

  mpf_clear (pi);

  return 0;
}

コンパイルします。このとき必ず -I/usr/local/include -L/usr/local/lib を指定します。

注意 : -I/usr/local/include の -I は「小文字のエル」ではなくて「大文字のアイ」です。

~> gcc -I/usr/local/include -L/usr/local/lib -Wall -O3 -o test_mpf_pi test_mpf_pi.c mpf_pi.c -lgmp

実行します。コマンドラインに求めたい円周率の桁数を書きます。2 番目の引数は円周率の値を出力するとき 1、しないとき 0 です。とりあえず 1000 桁まで計算してみましょう。

~> ./test_mpf_pi 1000 1
0.000 sec.
3.141592653589793238462643383279502884197169399375105820974944592307816406286208
99862803482534211706798214808651328230664709384460955058223172535940812848111745
02841027019385211055596446229489549303819644288109756659334461284756482337867831
65271201909145648566923460348610454326648213393607260249141273724587006606315588
17488152092096282925409171536436789259036001133053054882046652138414695194151160
94330572703657595919530921861173819326117931051185480744623799627495673518857527
24891227938183011949129833673362440656643086021394946395224737190702179860943702
77053921717629317675238467481846766940513200056812714526356082778577134275778960
91736371787214684409012249534301465495853710507922796892589235420199561121290219
60864034418159813629774771309960518707211349999998372978049951059731732816096318
59502445945534690830264252230825334468503526193118817101000313783875288658753320
83814206171776691473035982534904287554687311595628638823537875937519577818577805
321712268066130019278766111959092164201989
0.000 sec.

円周率の値の前後に計算開始から計算終了までの所要時間と計算開始から出力終了までの所要時間が表示されています。1000 桁のときは一瞬で終わってしまいました。これではよくわからないので次は 100 万桁まで計算してみます。

~> ./test_mpf_pi 1000000 1 > test_mpf_pi.out
~> cat test_mpf_pi.out
0.592 sec.
3.141592653589793238462643383279502884197169399375105820974944592307816406286208
99862803482534211706798214808651328230664709384460955058223172535940812848111745
02841027019385211055596446229489549303819644288109756659334461284756482337867831
65271201909145648566923460348610454326648213393607260249141273724587006606315588
17488152092096282925409171536436789259036001133053054882046652138414695194151160
94330572703657595919530921861173819326117931051185480744623799627495673518857527
24891227938183011949129833673362440656643086021394946395224737190702179860943702
77053921717629317675238467481846766940513200056812714526356082778577134275778960
91736371787214684409012249534301465495853710507922796892589235420199561121290219
60864034418159813629774771309960518707211349999998372978049951059731732816096318
59502445945534690830264252230825334468503526193118817101000313783875288658753320
83814206171776691473035982534904287554687311595628638823537875937519577818577805
32171226806613001927876611195909216420198938095257201065485863278865936153381827
─── 中略 ───
70770999399922886212243248802062634850888530360107234368901360642758142528398785
94917997961121963797576519245218670960880921371119775000878159304307293448839309
57574159241375285977797291893453850508038319867745900251865791723708085741642971
53807884060713068680361982419715774763895072534684045691927595319372237022290155
80065607604738547359904477996748749969769427137668695533195125337764098587096683
86326392616494560868414037456842071940595070174303546918215090046649399855174138
93851975731215682616228622318810967297476060130283311937161140874727067625585677
75119956667486151964912970193318084994109618139296492789360902125354433273750642
60624299412032736255824417498345094730945343661590728416319368307571979806823153
57371555718161221567879364250138871170232755557793022667858031999308108305763076
52332050740013939095807901637717629259283764874790177274125678190555562180504876
74699114083997791937654232062337471732470336976335792589151526031561403332127284
91944184371506965520875424505989567879613033116462839963464604220901061057794581
51
0.660 sec.

手元の AMD Ryzen 9 3900X の環境では 100 万桁の円周率が 0.660 秒で求まりました。出力は 0.068 秒で完了しており、算術演算だけでなく 10 進数への変換もかなり速いことがわかります。

4.2. Perl で GMP を使う

Perl を利用するとプログラムをコンパイルする手間も省けて手軽に多倍長計算を利用できます。サンプルとして P-1 法で素因数分解を行うスクリプトを書いてみました。

pm1.pl

use GMP::Mpz qw(:all);
$| = 1;
$expr = $ARGV[0];  #expression
$expr =~ s/\^/**/g;  #translate power operators
$n0 = eval('use GMP::Mpz qw(:constants);' . $expr);
print "$n0\n";
$t = time();
@p = ();  #prime factors
@c = ();  #composite factors
probab_prime_p($n0, 10) ? push(@p, $n0) : push(@c, $n0);
while (@c) {
  print '= ' . (@p ? join(' * ', sort { $a <=> $b } @p) . ' * ' : '') .
    '[' . join('] * [', sort { $a <=> $b } @c) . "]\n";
  $n = shift(@c);
  $m = mpz(3);
  for ($p = mpz(2); !divisible_p($n, $f = $p); $p = nextprime($p)) {
    $f *= $f while ($f < 1000);
    $m = powm($m, $f, $n);
    $f = gcd($n, $m - 1);
    $f > 1 && $f < $n and last;
  }
  probab_prime_p($f, 10) ? push(@p, $f) : push(@c, $f);
  if ($f < $n) {
    $f = $n / $f;
    probab_prime_p($f, 10) ? push(@p, $f) : push(@c, $f);
  }
}
print '= ' . join(' * ', sort { $a <=> $b } @p) . "\n";
print ((time() - $t) . " sec.\n");

運がよければ短い時間で素因数が見つかります。試しに 53 桁のレピュニット R(53)=(10⁵³-1)/9 を分解してみましょう。

~> perl pm1.pl "(10^53-1)/9"
11111111111111111111111111111111111111111111111111111
= [11111111111111111111111111111111111111111111111111111]
= 107 * [103842159916926272066458982346832814122533748701973]
= 107 * 1659431 * [62576967597282605945326429569432422392093283]
= 107 * 1659431 * 1325815267337711173 * 47198858799491425660200071
0 sec.

目次