GMPの使い方

GMPの使い方	2007-12-19(Wed) 23:28

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

Cygwinのインストール
GMPのインストール
使用例1（CでGMPを使う）
使用例2（PerlでGMPを使う）

1. Cygwinのインストール

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

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

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

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

1.2 setup.exeをダウンロードする

Cygwinは数多くのパッケージで構成されており、Cygwin Information and Installationにあるsetup.exeが専用のインストーラになっています。これをローカルパッケージディレクトリD:\cygwinpackages\にダウンロードします。

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

1.3 setup.exeを実行する

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

1.4 作業内容を選択する

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

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

Select Root Install DirectoryではRoot Directoryにルートディレクトリ「D:\cygwin\」を入力し、Install Forは「All Users (RECOMMENDED)」、Default Text File Typeは「Unix / binary (RECOMMENDED)」を選択します。

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

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

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

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

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

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

ローカルパッケージディレクトリD:\cygwinpackages\の下にftp%3a%2f%2fring.nict.go.jp%2farchives%2fpc%2fgnu-win32のような名前のディレクトリができます。

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

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

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

Mathにはコンパイル済みのGMPのパッケージが含まれていますが、使用環境に適した最適化が行われていないのでそのまま使うと本来の性能を発揮できません。そのため、ここではGMPのパッケージを外しておいて、後で自分でGMPをインストールすることにします。

注意：CygwinのGMPのパッケージをインストールしたまま後で自分でGMPをインストールしてしまうと厄介なことになるので、どちらを入れるにしても必ず他方をアンインストールしてから行ってください。

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

注意：以下の(1)(2)(4)(5)(6)はこの順序で行わないと失敗することがあります。

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

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

(4) libgmp-develの「○～」の部分を何度かクリックして「○Skip」（既にインストールしてしまっているときは「○Uninstall」）に変えます。

(5) gmpの「○～」の部分を何度かクリックして「○Skip」（既にインストールしてしまっているときは「○Uninstall」）に変えます。

(6) libgmp3の「○～」の部分を何度かクリックして「○Skip」（既にインストールしてしまっているときは「○Uninstall」）に変えます。

Warning! Unmet Dependencies Foundという画面が表示されます。この警告が出るのはGMPに依存しているパッケージがあるのにGMPだけ外したためです。しかし、依存しているパッケージを外そうとするとさらにそれに依存しているパッケージについて警告が出てしまい、すべて解決するのは面倒です。ここでは警告を無視してGMPなしのインストールを強行します。

(8) 下のほうにある「Install these packages to meet dependencies (RECOMMENDED)」の左側のチェックボックスのチェックを消します。

WARNING - Required Packages Not Selectedというダイアログが表示されます。

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

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

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

注意：パッケージをアンインストールするときはSelect Packagesで「○Uninstall」を選択します。これ以外の方法で削除してはいけません。

1.10 setup.exeを終了する

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

1.11 環境変数を設定する

Windows 98の場合はAUTOEXEC.BATに次のコマンドを追加してから再起動します。ルートディレクトリD:\cygwin\の部分は各自の環境に合わせます。usernameの部分はユーザ名を半角英数字で記述します。

set PATH=%PATH%;D:\cygwin\usr\local\bin;D:\cygwin\bin
set USER=username
set HOME=/home/username

Windows XPのときはシステムのプロパティの詳細設定にある環境変数で同様の設定を行います。

パスの追加はCygwinのウインドウを開かなくてもLinuxライクなコマンドを実行できるようにするためのものです。例えば、xyzzyの*Shell*でコンパイラが使えるようになります。usernameを指定しなかったときはWindowsのユーザ名が使用されますが、ユーザ名に半角英数字以外の文字が使われていると何かと不便だと思います。

2. GMPのインストール

注意：GMPの4.2.1までのバージョンにはGGNFSを誤動作させるバグ（[1],[2],[3],[4]を参照）があります。4.2.1までのバージョンをインストールしている人は4.2.2以上に更新してください。

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

~> cd gmp-4.2.1
~/gmp-4.2.1> make uninstall
~/gmp-4.2.1> cd ..
~>

GMPはソースコードの形で配布されています。The GNU MP Bignum Libraryからソースコードのアーカイブgmp-4.2.2.tar.bz2をダウンロードしてきて展開します。

~> tar jxf gmp-4.2.2.tar.bz2
~> cd gmp-4.2.2
~/gmp-4.2.2>

~/gmp-4.2.2> ./configure

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

~/gmp-4.2.1> make check

~/gmp-4.2.2> make install

ついでにPerlのモジュールもインストールしておきましょう。Perlで多倍長計算を手軽に利用できるようになります。

~/gmp-4.2.2> cd demos/perl
~/gmp-4.2.2/demos/perl> perl Makefile.PL
~/gmp-4.2.2/demos/perl> make
~/gmp-4.2.2/demos/perl> make test
~/gmp-4.2.2/demos/perl> make install

3. 使用例1（CでGMPを使う）

mpf_pi.c

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

  Reference:
    http://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 -Wall -O6 -ffast-math -fomit-frame-pointer test_mpf_pi.c mpf_pi.c -lgmp -o test_mpf_pi
*/

#include <stdio.h>
#include <time.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;
  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 = clock ();
  while ((start = clock ()) == stop);
  mpf_pi (pi);
  stop = clock ();
  printf ("%.3f sec.\n", ((double)(stop - start)) / CLOCKS_PER_SEC);
  if (output) {
    gmp_printf ("%.*Ff\n", prec10, pi);
    stop = clock ();
    printf ("%.3f sec.\n", ((double)(stop - start)) / CLOCKS_PER_SEC);
  }

  mpf_clear (pi);

  return 0;
}

> gcc -Wall -O6 -ffast-math -fomit-frame-pointer test_mpf_pi.c mpf_pi.c -lgmp -o test_mpf_pi

実行します。コマンドラインに求めたい円周率の桁数を書きます。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
15.485 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
17.109 sec.

手元のPentium 4 3.06GHzの環境では100万桁の円周率が15.5秒で求まりました。出力は1.6秒で完了しており、算術演算だけでなく10進数への変換もかなり速いことがわかります。参考までに、同じ環境でスーパーπは104万桁の円周率を求めるのに56秒かかりました。一方、大浦さんのpi_caは11秒台で140万桁を超えており、GMPのライブラリをアルゴリズム通りに呼び出すだけのプログラムでは大浦さんのプログラムには及びませんでした。なお、GNU MP Pi computationではチュドノフスキーの公式を用いる円周率計算プログラムが公開されており、同じ環境で3.4秒で100万桁に達しました。もちろん、いずれの結果も小数点以下100万桁目まで一致していました。

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^53-1)/9を分解してみましょう。