自作Lispで円周率1万桁計算に挑戦した話

# はじめに
自作Lisp,Easy-ISLisp(以下[EISL」）のBIGNUMのテストも兼ねて円周率1万桁計算をしていました。気が付いたことなど後日のために書きとどめます。

# アルゴリズム
Brent-Salamin の平方根を用いるアルゴリズムによりました。Lispコードは下記の通りです。整数平方根を計算するISQRTの計算が実行速度に大きく影響します。あれこれと改良を加え、現在では円周率1面桁計算におよそ1秒程度になりました。当初は5秒ほどかかっていました。

```lisp

(defun square-root (x)
    (isqrt x))

(defun square (x)
    (* x x))


#|
;;For Common-Lisp
(defun div (x y)
    (floor x y))
|#

; Compute pi using the 'brent-salamin' method.
(defun pi-brent-salamin (nb-digits)
    (let ((one (expt 10 nb-digits)))
        (pi-brent-salamin1 one
                           one
                           (square-root (div (square one) 2))
                           (div one 4)
                           1)))

(defun pi-brent-salamin1 (one a b tt x)
    (if (= a b)
        (div (square (+ a b)) (* 4 tt))
        (let ((new-a (div (+ a b) 2)))
            (pi-brent-salamin1 one
                               new-a 
                               (square-root (* a b))
                               (- tt (div (* x (square (- new-a a))) one))
                               (* 2 x)))))


Easy-ISLisp Ver2.54
> (load "tests/pi.lsp")
T
Easy-ISLisp Ver2.54
> (time (pi-brent-salamin 10000))
Elapsed Time(second)=1.076235
<undef>
> 

```

# isqrt　整数平方根
「ハッカーの楽しみ　本物のプログラマはいかにして問題を解くか」に整数平方根の話題があります。ニュートン法が最善の方法のように思えます。1万桁の計算ともなりますと、その積の整数平方根をとりこととなり、かなり大きな桁数になります。効率よく計算しなければなりません。当初の実装は極めて素朴でありニュートン法の初期値を求めようとする値そのものにしていました。計算はできますが時間とメモリを大食いします。

初期値を工夫して真の値に近づける工夫をしました。https://qiita.com/sym_num/items/a2adfaf4e1c2ec39da98
これにより2秒程度にまで短縮できました。しかし、ジュリアンさんのOpenLispは1万桁を0.3秒程度で計算します。何か工夫の余地はないものだろうか？と考えていましたら、初期値もニュートン法で求めればいいことを思いつきました。100桁の整数の整数平方根は５０桁程度です。これを求めるのに25桁（4分の１）以上の桁数で初期値を求めて、それによって得られた近似値を初期値にしてさらにニュートン法で計算するという方法です。多段式のロケットのようなものです。コードは下記のとおりです。isqrt3は普通にニュートン法で計算するものです。

```c
//newton method for bignum
int
isqrt2 (int n, int init)
{
  int s;

  s = init;
  while (greaterp (mult (s, s), n))
    {
      //print(s);printf("\n");
      s = divide (plus (s, divide (n, s)), makeint (2));
    }
  return (s);
}

// basic isqrt (for small bignum)
int
isqrt3 (int x)
{
  int len, msb, pointer, init, i;
  long long int lmsb;
  len = get_length (x);
  pointer = get_pointer (x);
  if (len % 2 == 0)
    {
      lmsb =
	(long long int) bigcell[pointer] * BIGNUM_BASE +
	(long long int) bigcell[pointer - 1];
      init = gen_big ();
      len = (len - 2) / 2;
      for (i = 0; i < len; i++)
	{
	  bigcell[big_pt0++] = 0;
	}
      bigcell[big_pt0++] = (int) sqrt (lmsb) + 10000;
      set_sign (init, 1);
      set_pointer (init, big_pt0 - 1);
      set_length (init, len + 1);
    }
  else
    {
      msb = bigcell[pointer];
      init = gen_big ();
      len = (len - 1) / 2;
      for (i = 0; i < len; i++)
	{
	  bigcell[big_pt0++] = 0;
	}
      bigcell[big_pt0++] = (int) sqrt (msb) + 1;
      set_sign (init, 1);
      set_pointer (init, big_pt0 - 1);
      set_length (init, len + 1);
    }

  return (isqrt2 (x, init));


}

/*
 * fast isqrt
 * Calculate the initial value of Newton's method by Newton's method.
 * The value of ISQRT (X) is half the number of digits in X. 
 * Find a value that is close to the true value with half the number of digits, 
 * and use that as the initial value.
*/
int
isqrt4 (int x)
{
  int len, p, init;

  len = get_length (x);

  p = 2;
  while ((len - p) > len / 4)
    {
      p = 2 * p;
    }
  p = p / 2;

  init = bigx_shift_left (x, p);
  init = plus (isqrt3 (init), makeint (1));
  init = bigx_shift_right (init, p / 2);

  return (isqrt2 (x, init));


}

int
isqrt (int x)
{
  if (integerp (x))
    return (makeint (floor (sqrt (GET_INT (x)))));
  else if (floatp (x))
    return (makeint (floor (sqrt (GET_FLT (x)))));
  else if (longnump (x))
    {
      long long int n, init;
      n = GET_LONG (x);
      init = (long long int) sqrt (n) + 1;
      return (makeint ((int) isqrt1 (n, init)));
    }
  else
    {
      if (get_length (x) < 100)
	return (isqrt3 (x));
      else
	return (isqrt4 (x));
    }
}



```

# 円周率計算はバグ発見に役に立つ
上記の工夫により2秒程度だったのが1秒程度にまで短縮されました。しかし、奇妙な値がでてくるケースがありました。それらは細部の見落としによるバグでした。ニュートン法で利用する除算にはクヌース先生のDアルゴリズムを使用しています。この時レアなケースで判定を間違えていました。他にも加算のキャリーがある場合の桁数のバグなどいくつか補正をしました。普通に加減乗除の計算をしているだけではバグは顕在化しませんでした。円周率計算はバグ発見に役に立ちました。


# より良いアルゴリズムを求めて
ジュリアンさんのものとはまだ3倍の速度の差があります。また、GMPを利用しているSBCLなどでは0.1秒程度です。Reddit-Lispに問いかけると誰か応えてくれるのでは？と期待してい投稿、質問してみました。コメントが付いたもののSBCLのコードを読めというものでした。GMPの解説によればISQRTはカラツバsquare-rootというアルゴリズムを用いているということでした。これについてはまだ理解できていません。なかなか世界レベルに達するのはたいへんなことです、あれこれと考え、工夫を、勉強をするのは楽しいものです。

EISLの実装はGithubにおいてあります。
https://github.com/sasagawa888/eisl