3 GMP

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

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

普段使っている 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++ インターフェイスが導入された。

• 初期化 mpz_init(変数名) または mpz_init2(変数名, ビット数)

    mpz_t i,j;
  
    mpz_init(i);
    mpz_init2(j, 1024);
  

• 定数の代入 mpz_set_str(変数名, 文字列, 基数)

    mpz_t i;
  
    mpz_init(i);
    mpz_set_str(i, "1234567890123456", 10);
  

• 演算 mpz_add(), mpz_sub(), mpz_mul()
• 出力 mpz_out_str(stdout, 基数, mpz型の式)

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%

５秒強で100万桁の数値を求めることが出来ている。