B..2.1 C++プログラム

/*
 * 竹内・安藤に載っている計算法を明確化して実装した。
 */

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
#include "easy_elliptic.hpp"

int main(void)
{
  int n, m;
  double k, kn, kp, a, an, bn, phi, newphi, newa, M, phin, pi, exact, estimate;
  pi = 4.0 * atan(1.0);

  cout << setprecision(16);

  k=0.99; phi=1;
  exact = F(k, phi);
  kn=k; kp=sqrt(1-kn*kn); phin=phi;
  a=1; an=a; bn=kp;
  for (n = 1; n < 100; n++) {
    cout << endl << "n=" << n << endl;
    newa = (an+bn)/2;
    bn = sqrt(an*bn);
    an = newa;
    newphi = phin + atan(kp*tan(phin));
    cout << "とりあえずnewphi=" << newphi << ", 2*phi" << (n-1)
	 << "=" << 2*phin << endl;
    if (abs(newphi - 2 * phin) <= pi / 2)
      cout << "|newphi-2*phi|=" << abs(newphi - 2 * phin) << "<=pi/2 なのでOK"
	   << endl;
    else {
      cout << "少し修正" << endl;
      m = rint((2 * phin - newphi) / pi);
      newphi += m * pi;
    }
    // 多分不要
    while (abs(newphi - 2 * phin) > pi/2) {
      cout << "小さすぎるので足す" << endl;
      newphi += pi;
      cout << "修正newphi=" << newphi << ", 2*phi=" << 2*phi << endl;
    }
    kn = (1-kp)/(1+kp);
    kp = sqrt(1-kn*kn);
    phin = newphi;
    estimate = a / an * phin / (1 << n);
    cout << n << " ステップの近似値=" << estimate <<
      ", 誤差=" << abs(exact-estimate) << endl;
    if (abs(an-bn) < 1e-15) {
      cout << "n=" << n << endl;
      break;
    }
  }
  if (n == 100)
    cout << "n==" << n << endl;
  M = an;
  estimate = a / M * phin / (1 << n);
  cout << "k=" << k << ", phi=" << phi << endl;
  cout << "M=" << M << ", phi" << n << "=" << phin << endl;
  cout << "評価=" << estimate << endl;
  cout << "boost によると F(k,phi)=" << exact << endl;
  return 0;
}