/*
 * testx-5.c --- 菊地文雄・山本昌宏『微分方程式と計算機演習』山海堂
 *              pp.125--127 のプログラムを C 言語に翻訳したものを
 *              -(p u')'+r u=λu, u(0)=u(1)=0 
 *              にした.  
 */

#include <stdio.h>
#include <math.h>

#define NDIM (1000)

#include <stdlib.h>
#include <limits.h>

int number_p = 0;
int number_r = 0;

double sqr(double x) { return x * x; }

/* [0,1] 上の一様乱数 */
double drandom()
{
  return random() / (double) INT_MAX;
}

/* プロトタイプ宣言 (定義はずっと後に出てくる) */
void TRID(double A[], double B[], double C[], double F[], int N, int ID);
void INPR(double *X, double *Y, int N, double *S);
void shinv(double a[], double b[], double c[], double x[], double y[],
           double eps, double h, double SHIFT, int N, int M, double *Eval,
           int *NITER, int *NPIVOT);

void show_vector(int n, char *name, double x[])
{
  int i;
  for (i = 1; i <= n; i++)
    printf("%s[%d]=%f\n", name, i, x[i]);
}

double Pa, Pb, Ra, Rb, Rc, Oa, Ob, Oc, Od, Qa, Qb, Qc, Qd, Qe;
double Ha, Hb, Ia, Ib, Ic, Ja, Jb, Jc, Jd;

double P(double x){
  if (number_p == 0)
    return 1.0;
  else if (number_p == 1)
    return Pa * x + Pb;   
  else if (number_p == 2)
    return Ra * x * x + Rb * x + Rc;
  else if (number_p == 3)
    return Oa * x * x * x + Ob * x * x + Oc * x + Od;
  else if (number_p == 4)
    return Qa * x * x * x * x + Qb * x * x * x + Qc * x * x + Qd * x + Qe;
}
double R(double x){
  if (number_r == 0)
    return 0.0;
  else if (number_r == 1)
    return Ha*x+Hb;
  else if (number_r == 2)
    return Ia * x * x + Ib * x + Ic;
  else if (number_r == 3)
    return Ja * x * x * x - Jb * x * x + Jc * x + Jd;
}

int main(void)
{
  int N, M, i, k, di;
  double h, h2, h3, x, SHIFT, EPS, EVAL, MINEIGEN;
  int NITER, NPIVOT; 
  double *A, *B, *C, *U, *V;
  char fname[128];
  FILE *fp;
  double pi = 4 * atan(1.0);

  printf("INPUT : N=\n");
  scanf("%d",&N);
  if(N<=1)  exit(0);
  A = malloc((N+1) * sizeof(double));
  B = malloc((N+1) * sizeof(double));
  C = malloc((N+1) * sizeof(double));
  U = malloc((N+1) * sizeof(double));
  V = malloc((N+1) * sizeof(double));
  EPS = 1e-15; M = 100;

  printf("0: P(x)=1.0,\n"); 
  printf("1: P(x)=ax+b,\n");
  printf("2: P(x)=ax^2+bx+c,\n");
  printf("3: P(x)=ax^3+bx^2+cx+d, \n"); 
  printf("4: P(x)=ax^4+bx^3+cx^2+dx+e \n"); 
 printf("P(x)=");
  scanf("%d", &number_p);
  if (number_p == 1) {
    printf("a: "); scanf("%lf", &Pa);
    printf("b: "); scanf("%lf", &Pb);
  }
  else if (number_p == 2){
    printf("a: "); scanf("%lf", &Ra);
    printf("b: "); scanf("%lf", &Rb);
    printf("c: "); scanf("%lf", &Rc);
  }
  else if (number_p == 3){
    printf("a: "); scanf("%lf", &Oa);
    printf("b: "); scanf("%lf", &Ob);
    printf("c: "); scanf("%lf", &Oc);
    printf("d: "); scanf("%lf", &Od);
  }
  else if (number_p == 4){
    printf("a: "); scanf("%lf", &Qa);
    printf("b: "); scanf("%lf", &Qb);
    printf("c: "); scanf("%lf", &Qc);
    printf("d: "); scanf("%lf", &Qd);
    printf("e: "); scanf("%lf", &Qe);
  }

  printf("0: R(x)=0,\n");
  printf("1: R(x)=ax+b,\n");
  printf("2: R(x)=ax^2+bx+c,\n");
  printf("3: R(x)=ax^3+bx^2+cx+d \n");
  printf("R(x)="); 
  scanf("%d", &number_r);
  if (number_r == 1) {
    printf("a: "); scanf("%lf", &Ha);
    printf("b: "); scanf("%lf", &Hb);
  }
  else if (number_r == 2){
    printf("a: "); scanf("%lf", &Ia);
    printf("b: "); scanf("%lf", &Ib);
    printf("c: "); scanf("%lf", &Ic);
  }
  else if (number_r == 3){
    printf("a: "); scanf("%lf", &Ja);
    printf("b: "); scanf("%lf", &Jb);
    printf("c: "); scanf("%lf", &Jc);
    printf("d: "); scanf("%lf", &Jd);
  }

  printf("number p=%d\n", number_p);
  printf("number_r=%d\n", number_r);

  printf("n= "); scanf("%d", &di);
 
  h = 1.0/N;
  h2 = h*h;
  h3 = h/2;

  /* 計算結果をファイルに書き出す */
  printf("file name: ");
  scanf("%s", fname);
  fp = fopen(fname, "w");
  fprintf(fp, "# N, EPS, NITER\n");
  fprintf(fp, "# %d %g %g %d\n", N, EPS, M);
  fprintf(fp, "# number_p=%d, number_r=%d\n",number_p,number_r);
  if (number_p == 9) {
    fprintf(fp, "# a=%f, b=%f\n", Pa, Pb);
  }
  for (k = 0; k < di; k++) {
    if (k == 0) {
      SHIFT = 0;
    }
    else {
      SHIFT = sqr((k + 1.0) / k) * EVAL;
    }
    for(i = 1; i < N; i++){
      x = i * h;
      A[i] = -P(x-h3)/h2;
      B[i] = (P(x+h3)+P(x-h3)+h2*R(x))/h2-SHIFT;
      C[i] = -P(x+h3)/h2;
    } 
    shinv(A,B,C,U,V,EPS,h,SHIFT,N-1,M,&EVAL,&NITER,&NPIVOT);
    printf("反復回数 = %3d  EPS = %12.4e SHIFT = %12.4e NPIVOT = %3d\n",
	   NITER, EPS, SHIFT, NPIVOT);
    if (k == 0) {
      MINEIGEN = EVAL;
      printf("固有値 = %20.15f\n",EVAL);
      fprintf(fp, "# 固有値 = %20.15f\n",EVAL);
    }
    else {
      printf("固有値 = %20.15f (%f)\n",EVAL, sqrt(EVAL / MINEIGEN));
      fprintf(fp, "# 固有値 = %20.15f (%f)\n",EVAL, sqrt(EVAL / MINEIGEN));
    }
#ifdef OLD
    printf("** 節点での固有関数の値 **\n");
    for(i=0; i<5; i++)
      printf("  I     U(I)    ");
    printf("\n");
    for(i=1; i<N; i++){
      printf("%3d%12.4e ", i, U[i]);
      if(i%5==0)  printf("\n");
    }
    if((N-1)%5!=0)  printf("\n");
#endif
    for (i = 0; i <= N; i++)
      fprintf(fp, "%f %g\n", i * h, U[i]);
  }
  fclose(fp);
  return(0);
}

void shinv(double a[], double b[], double c[], double x[], double y[],
           double eps, double h, double SHIFT, int N, int M,
	   double *Eval, int *NITER, int *NPIVOT)
{
  int i;
  double s,t,ENEW,ERROR;
  /* J=1;
     これは RND() に渡す種(seed)だったが、RND() の代りに drandom() を
     使うことになったので不要になった */
  for(i=1; i<=N; i++){
    x[i]=drandom();
  }
  *NITER = 0;
  *Eval = SHIFT;
  /* ITERATION */
  do {
    for (i=1; i<=N; i++){
      y[i]=x[i];
    }
    TRID(a, b, c, y, N, *NITER);
    INPR(x, y, N, &s);
    INPR(y, y, N, &t);
    ENEW=s/t+SHIFT;
    t=sqrt(h*t);
    for(i=1; i<=N; i++)
      x[i]=y[i]/t;
    ERROR=fabs((ENEW-*Eval)/ENEW);
    *Eval=ENEW;
    *NITER=*NITER+1;
#ifdef OLD
    printf("NITER=%d, Eval=%f, ERROR=%f\n", *NITER, *Eval, ERROR);
#endif
  }
  while (ERROR > eps && *NITER < M);
  *NPIVOT=0;
  for (i=1; i<=N; i++) {
    if (b[i]<0.0) *NPIVOT=*NPIVOT+1;
  }
}

void TRID(double A[], double B[], double C[], double F[], int N, int ID)
{
  int i;
  double AA;
  for(i=1; i< N; i++){
    AA=A[i+1]/B[i];
    if(ID==0) B[i+1]=B[i+1]-AA*C[i];
    F[i+1]=F[i+1]-AA*F[i];
  }
  F[N]=F[N]/B[N];
  for(i=N-1; i>=1; i--)
    F[i]=(F[i]-C[i]*F[i+1])/B[i];
}

void INPR(double *X, double *Y, int N, double *S)
{
  int i;
  *S=0.0;
  for (i=1; i<=N; i++){
    *S=*S+X[i]*Y[i];
  }
}