3 GMP

GMP (GNU Multiple Precision arithmetic library) は、 名前が示しているように GNU Project でサポートされている多倍長計算パッケージで、 C言語とアセンブリ言語で記述されている。

付属するマニュアル (WWWサイトでも入手可能) が十分詳しい。 アルゴリズムについても1章を使って説明されている。

普段使っている Mac には、気がついたときには既に GMP が入っていた (いつインストールしたのか覚えていない、何かの依存関係?)。 何でも GCC のインストールには GMP が必要になっているそうで、 UNIX 風の環境ではあるのが当たり前に近いライブラリィなのかもしれない。

mpz_t 整数型
mpq_t 有理数型
mpf_t 浮動小数点数型 (Version 2から)
mpfr_t 浮動小数点数型 (Version 4から, IEEE 754 互換)

typedef  struct  {
   mpfr_prec_t    _mpfr_prec;
   mpfr_sign_t    _mpfr_sign;
   mp_exp_t          _mpfr_exp;
   mp_limb_t      *_mpfr_d;
}  __mpfr_struct;

Verion 4 から C++ インターフェイスが導入された。

pi.c -- 100万桁

/*
The Gauss AGM algorithm using Schonhage variation.
[1] Jorg Arndt, Christoph Haenel, "Pi Unleashed", pp. 93, Springer, 2001.
(C) 2003 Hironobu SUZUKI, licensed by GPL2 
$Id: pi.c,v 1.1 2003/09/16 01:20:41 hironobu Exp hironobu $
*/
#include <stdio.h>
#include <gmp.h>
#include <sys/timeb.h>
#include <stdlib.h>

#define LOG_TEN_TWO  3.32192809488736234789
#define bprec(n) (int)(((n+10)*LOG_TEN_TWO)+2)

int
main (int ac, char *av[])
{
  long int k, loopmax;
  mpf_t A, B, a, b, s, s_1, t, t1, t2, t3, c2;
  long int prec, dprec;
  struct timeb start_t, end_t;

  dprec = 1000000L;		/* decimal precision */
  prec = bprec(dprec);		/* binary precision (plus alpha) */
  mpf_set_default_prec (prec);
  loopmax = 21;

  mpf_init (A);			/* big A */
  mpf_init (B);			/* big B */
  mpf_init (a);			/* a */
  mpf_init (b);			/* b */
  mpf_init (s);			/* s(n) */
  mpf_init (s_1);		/* s(n-1) */
  mpf_init (t);			/* temporary */
  mpf_init (t1);		/* temporary */
  mpf_init (t2);		/* temporary */
  mpf_init (t3);		/* temporary */
  mpf_init (c2);		/* 2 constant */
  mpf_set_ui (A, 1);
  mpf_set_ui (a, 1);
  mpf_set_ui (t1, 1);
  mpf_div_ui (B, t1, 2);
  mpf_div_ui (s, t1, 2);
  mpf_set_ui (t1, 10);
  mpf_set_ui (c2, 2);
  for (k = 1; k <= loopmax; k++)
    {
      mpf_add (t1, A, B);	/* (A+B) */
      mpf_div_ui (t, t1, 4);	/*  t = (A+B)/4 */
      mpf_sqrt (b, B);		/* b = sqrt(B) */
      mpf_add (t1, a, b);	/* (a+b) */
      mpf_div_ui (a, t1, 2);	/* a = (a+b)/2 */
      mpf_mul (A, a, a);	/* A = a * a */
      mpf_sub (t1, A, t);	/*  (A-t) */
      mpf_mul_ui (B, t1, 2);	/* B = (A - t) * 2 */
      mpf_sub (t1, B, A);	/* (B-A) */
      mpf_pow_ui (t2, c2, k);	/* 2^k */
      mpf_mul (t3, t1, t2);	/* (B-A) * 2^k  */
      mpf_set (s_1, s);
      mpf_add (s, s_1, t3);	/* s = s + (B-A) * 2^k */
    }
  mpf_add (t1, a, b);
  mpf_pow_ui (t2, t1, 2);
  mpf_mul_ui (t1, s, 2);
  mpf_div (A, t2, t1);
  mpf_out_str (stdout, 10, dprec + 10, A);

  exit(0);
}
MacBook Airでのコンパイルと実行
[katsurada:~/study-c-c++-2010] mk% gcc pi-gmp.c -lgmp
[katsurada:~/study-c-c++-2010] mk% time ./a.out>pi.out
5.446u 0.093s 0:05.55 99.6%     0+0k 0+3io 0pf+0w
[katsurada:~/study-c-c++-2010] mk%
5秒強で100万桁の数値を求めることが出来ている。

桂田 祐史
2017-09-13