#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include "str01.h"

using namespace qrk;


static matrix ret_mat;
typedef double (*daryp)[Max][Max];
typedef double (dary)[Max][Max];

double abs1(double);


/*--- 行列とベクトルの積 */
vector matvec(matrix inmat, vector v1)
{
  int i, j;
  vector ret_v;

  for (i=0; i < Max; i++) {
    ret_v.a[i]=0.0;
  }
  ret_v.err = 0;
  if((inmat.n > Max) || (inmat.n<0) || (inmat.m>Max) || (inmat.m<0)) {
    ret_v.err=1;
  }
  if ((v1.n > Max) || (v1.n<0)) {
    ret_v.err=1;
  }
  if ((v1.n) != (inmat.m)) {
    ret_v.err=1;
  }
  if (ret_v.err == 1) {
    ret_v.n=0;
    return ret_v;
  }
  ret_v.n=inmat.n;
  for (i=0;i<inmat.n;i++) {
    for (j=0;j<inmat.m;j++) {
      ret_v.a[i] = ret_v.a[i] + inmat.a[i][j] * v1.a[j];
    }
  }
  return ret_v;
}


/*--- 行列の積 */
matrix matmat(matrix *m1, matrix *m2)
{
  int i, j, k, mcount;

  ret_mat.err = 0;
  if ((m1->n>Max) || (m1->n<0) || (m1->m>Max) || (m1->m<0)) {
    ret_mat.err=1;
  }
  if ((m2->n>Max) || (m2->n<0) || (m2->m>Max) || (m2->m<0)) {
    ret_mat.err=1;
  }

  if ((m1->m) != (m2->n)) {
    ret_mat.err=1;
  }
  mcount = m1->m;
  if (ret_mat.err == 1) {
    strcpy(ret_mat.message,"Something went wrong.");
    printf("Something went wrong!\n");
    ret_mat.n = 0;
    ret_mat.m = 0;
    return ret_mat;
  }

  ret_mat.n = m1->n;
  ret_mat.m = m2->m;

  for (i=0;i<ret_mat.n;i++)  {
    for (j=0;j<ret_mat.m;j++) {
      ret_mat.a[i][j]=0.0;
    }
  }

  for (i=0;i<ret_mat.n;i++){
    for (j=0;j<ret_mat.m;j++){
      for (k=0;k<mcount;k++){
        ret_mat.a[i][j] += m1->a[i][k] * m2->a[k][j] ;
      }
    }
  }
  return ret_mat;
}


/*--- 転置行列を求める */
int transpose(matrix *dest, matrix *src)
{
  int i, j;

  dest->err = 0;
  if ((src->n>Max) || (src->n<0)) {
    dest->err = 1;
  }
  if ((src->m>Max) || (src->m<0)) {
    dest->err = 1;
  }
  if (dest->err == 1) {
    strcpy(dest->message, "Define the size of your matrix, stupid.");
    return 0 ;
  }
  strcpy(dest->message,"Done.");
  dest->n = src->m;
  dest->m = src->n;
  for (i=0;i<src->n;i++) {
    for (j=0;j<src->m;j++) {
      dest->a[j][i] = src->a[i][j];
    }
  }
  return 1 ;
}


double nfunc(double x, double mk, double e)
{
  return(x - e * sin(x) - mk);  /*  Kepler's equation  */
}


/*--- ニュートン法 */
double newton(double x0, double mk, double e)
{
  int i;
  double x1, fx0, dfdx0, aeps;
  double soln;
  int iterations, err;

  int imax=20;
  double eps=1e-13;

  err=1;
  aeps=x0 * eps / 10;

  for (i=0;i<imax;i++) {
    fx0 = nfunc(x0, mk, e);
    dfdx0 = (nfunc(x0 + aeps, mk, e)-fx0) / aeps;
    x1= x0 - fx0 / dfdx0;
    if (fabs(x1-x0) < 10*aeps) {
      err=0;
      break;
    }
    x0=x1;
  }
  iterations = i + 1;
  soln = x1;
  if (err) {
    printf("Error: not convergence in newton()\n");
  }
  return(soln);
}


void matdisp(matrix mm)
{
  int i,j;
  printf("(rows, columns) = (%d, %d)\n",mm.n, mm.m);
  for (i=0;i<mm.n;i++) {
    for (j=0;j<mm.m;j++) printf("%6.2f ",mm.a[i][j]);
    printf("\n");
  }
  printf("\n");
}
