/*
   Copyright(C) 1999 松元隆二

   本プログラムは「GNU General Public License」です．

1999/05/19 初版作成．

別ファイルのrama_table.m 「Ramanujan型級数一覧」 の検査用のプログラム
です．mathematicaのプログラムでも検査可能ですが，mathematicaはたまに計
算ミスします．

本プログラムをコンパイルするには，ファイル「rama_check_init.cc」が必要
です．これは配列の初期化用のデータです．関数int main()内部で直接include
しています．ファイルrama_table.mにあるmathematicaの関数printnum[]を引数
1000で実行し，出力されたデータをファイルにセーブして生成してください．

本プログラムは，数値計算ライブラリLiDIA Version 1.3.2に依存しています．
LiDIAは次の所から入手してください．

          http://www.informatik.tu-darmstadt.de/TI/LiDIA/Welcome.html

コンパイル引数:
        g++ rama_check.cc 環境依存の引数...

*/
#include <LiDIA/bigfloat.h>
#include <LiDIA/polynomial.h>
#include <LiDIA/bigint_lattice.h>

#define EXP_LOG 1000
//#define EXP_LOG 100

bigfloat LOG_E_10;
bigfloat PI;
bigfloat INV_PI;
bigfloat EXP;

double bigfloat2double(const bigfloat &x)
{
  char buff[EXP_LOG * 2];
  bigfloat_to_string(x, buff);

  double y = atof(buff);

  return y;
}

int calc(bigfloat X, bigfloat Y, bigfloat Z2, int is, int flag, bigfloat *ans)
{
  bigfloat tmp;

  bigfloat s[3];

  s[0] = 1/bigfloat(2);

  if(is > 0)
    {
      tmp = 1/bigfloat(is);
      s[1] = s[0] + tmp;
      s[2] = s[0] - tmp;
    }

  bigfloat sum = 0;
  bigfloat diff = 1;

  bigfloat R = 1;
  bigfloat ZZ = 1;
  bigfloat XY = X;
  unsigned long n;
  bool final = false;

  for(n = 0;;)
    {
      /* cacl term n*/
      diff = R * XY * ZZ;
      
      if(flag == -1 && n%2 == 1)
	sum -= diff;
      else
	sum += diff;

      ZZ *= Z2;
      XY += Y;

      /*cout.form("R[%d]=",n);
      cout << R << endl;*/

      /* init. next n */


      if(is > 0)
	R *= ((s[0] + n) * (s[1] + n) *(s[2] + n));
      else
	{
	  tmp = s[0] + n;
	  R *= (tmp*tmp*tmp);
	}

      n++;

      if(n < 1000)
	R /= (n*n*n);
      else
	R /= n, R /= n, R /= n;

      
      if(final)
	break;

      if(diff < EXP)
	final = true;/* あと一回で終り */
    }

  *ans = sum;

  return n;
}

void check(bigfloat X[], bigfloat Y[], bigfloat Z2[], 
	   int size, int type, char symbol[])
{
  int loop_last;
  bigfloat ans, diff;
  bigfloat isqrt3 = sqrt(1/bigfloat(3));

  for(int i = 1; i <= size; i++)
    {
      if(X[i] > 0 && Y[i] > 0 && abs(Z2[i]) > 0)
	{
	  switch(type)
	    {
	    case 1: /*type A*/
	      loop_last  = calc(X[i], Y[i], Z2[i], 0, +1, &ans);
	      break;
	    case 2: /*type B*/
	      loop_last  = calc(X[i], Y[i], Z2[i], 0, -1, &ans);
	      break;
	    case 3: /*type D*/
	      loop_last  = calc(X[i], Y[i], Z2[i], 4, +1, &ans);
	      ans *= sqrt(Z2[i]);
	      break;
	    case 4: /*type E*/
	      loop_last  = calc(X[i], Y[i], Z2[i], 4, -1, &ans);
	      ans *= sqrt(Z2[i]);
	      break;
	    case 5: /*type F-1*/
	      loop_last  = calc(X[i], Y[i], Z2[i], 3, +1, &ans);
	      ans *= (sqrt(Z2[i])*isqrt3);
	      break;
	    case 6: /*type F-2*/
	      loop_last  = calc(X[i], Y[i], -Z2[i], 3, -1, &ans);
	      ans *= (sqrt(Z2[i]/bigfloat(-1728)));
	      break;
	    case 7: /*type F-3*/
	      loop_last  = calc(X[i], Y[i], -Z2[i], 3, -1, &ans);
	      ans *= (sqrt(Z2[i]/bigfloat(-12)));
	      break;
	    case 8: /*type G s=1/3 */
	      loop_last  = calc(X[i], Y[i], Z2[i], 3, +1, &ans);
	      break;
	    case 9: /*type G s=1/4 */
	      loop_last  = calc(X[i], Y[i], Z2[i], 4, +1, &ans);
	      break;
	    case 10: /*type G s=1/6 */
	      loop_last  = calc(X[i], Y[i], Z2[i], 6, +1, &ans);
	      break;
	    default:
	      cerr<<"Bad Argument!\n";
	      exit(1);
	    }
	  
	  diff = ans - INV_PI;
	  diff = abs(diff);

	  double diff_log = bigfloat2double(log(diff)/LOG_E_10);
	  double p = - diff_log / loop_last;
	  if(diff > EXP)
	    {
	      cout.form("%s, N=%d: Bad", symbol, i);
	    }
	  else
	    {
	      cout.form("%s, N=%d: OK", symbol, i);
	    }
	  cout.form(", looplast=%d, diff=%.15f, P(%.3fn)\n", 
		    loop_last, diff_log, p);
	}
    }
}

int main()
{
  bigfloat::set_precision(EXP_LOG);
  if(EXP_LOG < 30)
    power(EXP, 10, -EXP_LOG+3);
  else if(EXP_LOG < 100)
    power(EXP, 10, -EXP_LOG+ 5);
  else
    power(EXP, 10, -EXP_LOG + int(EXP_LOG * 0.01));

  LOG_E_10 = log(bigfloat(10));  
  constant_Pi(PI);
  INV_PI = 1 / PI;
  
#define SIZE 163+1
  bigfloat AX[SIZE],  AY[SIZE],  AZ2[SIZE];
#undef SIZE
#define SIZE 232+1
  bigfloat BX[SIZE],  BY[SIZE],  BZ2[SIZE];

#undef SIZE
#define SIZE 190+1
  bigfloat DX[SIZE],  DY[SIZE],  DZ2[SIZE];

#undef SIZE
#define SIZE 793+1
  bigfloat EX[SIZE],  EY[SIZE],  EZ2[SIZE];

#undef SIZE
#define SIZE 190+1
  bigfloat F1X[SIZE],  F1Y[SIZE],  F1Z2[SIZE];

#undef SIZE
#define SIZE 1555+1
  bigfloat F2DX[SIZE],  F2DY[SIZE],  F2DZ2[SIZE];

#undef SIZE
#define SIZE 58+1
  bigfloat G3X[SIZE],  G3Y[SIZE],  G3Z2[SIZE];

#undef SIZE
#define SIZE 116+1
  bigfloat G4X[SIZE],  G4Y[SIZE],  G4Z2[SIZE];

#undef SIZE
#define SIZE 6+1
  bigfloat G6X[SIZE],  G6Y[SIZE],  G6Z2[SIZE];

#include "rama_check_init.cc"

  /* 計算開始 */

  check(AX, AY, AZ2, 163, 1, "Type A");
  BZ2[2]=0;/*N=2の時はとても収束が遅いので無視*/
  check(BX, BY, BZ2, 232, 2, "Type B");
  check(DX, DY, DZ2, 190, 3, "Type D");
  check(EX, EY, EZ2, 793, 4, "Type E");
  check(F1X, F1Y, F1Z2, 190, 5, "Type F-1");
  check(F2DX, F2DY, F2DZ2, 1555, 6, "Type F-2'");
  check(G3X, G3Y, G3Z2, 58, 8, "Type G-3");
  check(G4X, G4Y, G4Z2, 116, 9, "Type G-4");
  check(G6X, G6Y, G6Z2, 6, 10, "Type G-6");

  return 0;
}


/* EOF */