// ----------------------------------------------------------------------
//  式(1.8)の百万番目までの和と同精度のπ/4の値をエイトケン加速で計算
//
// ＜変数の説明＞
//  S: 式(1.19)の計算結果 S_n^(k) (p.10グレー枠)の値を保存する2次元配列
//  m: 配列Sの縦方向の添字 (m = n + 1) n = 0,1,2,…, m = 1,2,3,…
//  h: 配列Sの横方向の添字 (h = k + 1) k = 0,1,2,…, h = 1,2,3,…
//
// ※ Scilabでは，配列の添字が1から始まる(C言語では0から始まる)．
//   式(1.19)で，n,k は0から始まるが，0は配列Sの添字として使えない．
//   そこで添字用の変数 m,h を用意し，S_n^(k)の結果を配列S(m,h)に保存
// ----------------------------------------------------------------------
clear;                  // 定義済の変数の値を全て消去 (なくてもOK)
stacksize(10000000);    // 使用可能なメモリのサイズを増やす (繰り返し回数が非常多いため必要)

itr = 1000000;               // 式(1.8)で何番目まで和を計算するか(百万番目)
S_real = %pi/4;              // 式(1.8)の理論値 π/4=0.7853981634…
exec('pi4.sci');             // 式(1.8)の関数のファイルpi4.sciの読み込み
S_slow = pi4(itr);           // 式(1.8)の百万番目までの和を計算
err1 = abs(S_real - S_slow); // 式(1.8)の百万番目まで和と理論値との絶対誤差
                             // (↑エイトケンの繰返しの停止条件として利用)
for n=0:itr
   m = n + 1;                // mは配列Sの縦方向の添字(nは0から，mは1から始まる)
   S(m,1) = pi4(n);          // p.10グレー枠の一番左の列のS_n^(0)の値の計算
   err2 = abs(S(m,1) - S_real);
   for k = 0:((n/2)-1)       // p.10グレー枠内を左下から右上へ向かって計算
      if err2 <= err1        // 誤差が小さくなったら終了(breakで内側のfor文から脱出)
         break
      end
      h = k + 1;             // hは配列Sの横方向の添字(kは0から，hは1から始まる)
      m = m - 1;             // 計算を行う位置を一段上の行へ移動する
      S(m,h+1) = (S(m+1,h)*S(m-1,h)-S(m,h)^2)/(S(m+1,h)+S(m-1,h)-2*S(m,h));
      err2 = abs(S(m,h+1) - S_real);  // ↑ エイトケン加速の計算式(1.19)そのまま
   end
   if err2 <= err1           // 誤差が小さければ終了(breakで外側のfor文から脱出)
      break
   end
end

format(11);                  // 9桁まで数値を表示するように変更 (9=11-2)
disp(S);                     // 配列 S の中身を表示

printf( "error1 = %.16f\n", err1 );              // 式(1.8)の結果の理論値との誤差
printf( "error2 = %.16f\n", err2 );              // エイトケンの結果の理論値との誤差
printf( "S_real = %.16f\n", S_real );            // π/4の理論値
printf( "S_slow = %.16f\n", S_slow );            // 式(1.8)の百万番目までの和
printf( "S(%d,%d) = %.16f\n", m-1, h, S(m,h+1)); // エイトケン加速の結果