Коллекция алгоритмов от Johna Smith - C++ - Ответ 4579895

//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (substraction)
//  (c) Johna Smith, 1996
//
//  Given: A,B - matrixes M x N
//  This algorithm substract these matrixes and display the result.
//
//////////////////////////////////////////////////////////////////////////////
#define M   5
#define N   6
#include <stdio.h>
int a[N][M]={{4,4,4,4,4},{5,4,3,2,1},{1,2,3,4,5},{5,4,3,2,1},{1,2,3,4,5},{0,1,0,1,0}};
int b[N][M]={{3,4,4,4,4},{5,3,3,2,1},{1,2,2,4,5},{5,4,3,1,1},{1,2,3,4,4},{0,1,0,1,0}};
int c[N][M];
void show_matrix(int *matrix) // this functions displays matrix
{
 for(int i=0;i<N;i++)
 {
    for(int j=0;j<M;j++)
      printf("%d ",*(matrix+i*M+j));
    printf("\n");
 }
 printf("\n");
}
void main(void)
{
  // display matrixes A and B
  show_matrix(a[0]);
  show_matrix(b[0]);
  // substract these matrixes
  for(int i=0;i<N;i++)
    for(int j=0;j<M;j++)
      c[i][j]=a[i][j]-b[i][j];
  // display the difference
  show_matrix(c[0]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (transponating)
//  (c) Johna Smith, 1996
//
//  Given: A (N x N),                                 T
//  This algorithm transponates matrix A and returns A  .
//          T
//   A(ji)=A (ij)
//
//////////////////////////////////////////////////////////////////////////////
#define N  3
#include <stdio.h>
float a[N][N]={{4,8,0},{8,8,8},{2,0,1}};
void Transponate(float *matrix)
{
  float swp;
  for (int i=0;i<N;i++)
    for (int j=i+1;j<N;j++)
    {
      swp=*(matrix+i*N+j);
      *(matrix+i*N+j)=*(matrix+j*N+i);
      *(matrix+j*N+i)=swp;
    }
}
void show_matrix(float *matrix) // this functions displays matrix
{
 for(int i=0;i<N;i++)
 {
   for(int j=0;j<N;j++)
     printf("%f ",*(matrix+i*N+j));
   printf("\n");
 }
 printf("\n");
}
void main(void)
{
  // display matrix A
  show_matrix(a[0]);
  // Transponate it
  Transponate(a[0]);
  // display transponated matrix
  show_matrix(a[0]);
  // Transponate it again
  Transponate(a[0]);
  // display the transponation of the transponated matrix
  show_matrix(a[0]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (addition and substraction in a+bi form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=a1+i*b1, z2=a2+i*b2
//  z1+z2=(a1+a2)+i*(b1+b2)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
struct complex
{
  float re;
  float im;
};
void show_complex(complex c) // this functions displays complex value
{
  printf("%f%+fi",c.re,c.im);
}
complex Add(complex a, complex b)
{
  complex c;
  c.re=a.re+b.re;
  c.im=a.im+b.im;
  return c;
}
complex Substract(complex a, complex b)
{
  complex c;
  c.re=a.re-b.re;
  c.im=a.im-b.im;
  return c;
}
complex a={2,3},b={4,6};
void main(void)
{
  // addition
  show_complex(a);
  printf(" + ");
  show_complex(b);
  printf(" = ");
  show_complex(Add(a,b));
  // substraction
  printf("\n");
  show_complex(b);
  printf(" - ");
  show_complex(a);
  printf(" = ");
  show_complex(Substract(b,a));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (division in a+bi form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=a1+i*b1, z2=a2+i*b2
//  z1/z2=(a1*a2-b1*b2)/(a2^2+b2^2)+i*(b1*a2-a1*b2)/(a2^2+b2^2)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
struct complex
{
  float re;
  float im;
};
void show_complex(complex c) // this functions displays complex value
{
  printf("%f%+fi",c.re,c.im);
}
complex Div(complex a, complex b)
{
  complex c;
  c.re=(a.re*b.re+a.im*b.im)/(b.re*b.re+b.im*b.im);
  c.im=(a.im*b.re-a.re*b.im)/(b.re*b.re+b.im*b.im);
  return c;
}
complex a={2,3},b={-1,2};
void main(void)
{
  show_complex(a);
  printf(" / ");
  show_complex(b);
  printf(" = ");
  show_complex(Div(a,b));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (multiplication in a+bi form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=a1+i*b1, z2=a2+i*b2
//  z1*z2=(a1*a2-b1*b2)+i*(a1*b2+a2*b1)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
struct complex
{
  float re;
  float im;
};
void show_complex(complex c) // this functions displays complex value
{
  printf("%f%+fi",c.re,c.im);
}
complex Mul(complex a, complex b)
{
  complex c;
  c.re=a.re*b.re-a.im*b.im;
  c.im=a.im*b.re+a.re*b.im;
  return c;
}
complex a={2,3},b={-1,2};
void main(void)
{
  show_complex(a);
  printf(" * ");
  show_complex(b);
  printf(" = ");
  show_complex(Mul(a,b));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (division in M*exp(i*phi) form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=M1*exp(i*phi1), z2=M2*exp(i*phi2)
//  z1/z2=M1/M2*exp(i*(phi1-phi2))
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define Pi      3.1415926536
struct exp_complex
{
  float M;
  float phi;
};
void show_exp_complex(exp_complex c) // this functions displays complex value
{
  printf("%f*exp(%f*i)",c.M,c.phi);
}
exp_complex Div(exp_complex a,exp_complex b)
{
  exp_complex c;
  c.M=a.M/b.M;
  c.phi=a.phi-b.phi;
  return c;
}
exp_complex a={1,0}, b={-1,-Pi/2};
void main(void)
{
  show_exp_complex(a);
  printf(" / ");
  show_exp_complex(b);
  printf(" = ");
  show_exp_complex(Div(a,b));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (multiplication in M*exp(i*phi) form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=M1*exp(i*phi1), z2=M2*exp(i*phi2)
//  z1*z2=M1*M2*exp(i*(phi1+phi2))
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define Pi      3.1415926536
struct exp_complex
{
  float M;
  float phi;
};
void show_exp_complex(exp_complex c) // this functions displays complex value
{
  printf("%f*exp(%f*i)",c.M,c.phi);
}
exp_complex Mul(exp_complex a,exp_complex b)
{
  exp_complex c;
  c.M=a.M*b.M;
  c.phi=a.phi+b.phi;
  return c;
}
exp_complex a={-1,Pi/2}, b={-1,-Pi/2};
void main(void)
{
  show_exp_complex(a);
  printf(" * ");
  show_exp_complex(b);
  printf(" = ");
  show_exp_complex(Mul(a,b));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating <x>, (x-dx,x+dx)
//  (c) Johna Smith, 1996
//
//  Method description:
//    dx=t(alpha,n)*Sn/sqrt(n)
//    <x>=Sum x/n
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N   6 // number of experiments
float x[N]={3,5,4,5,3,4}; // experimental data
double Student095(unsigned int n)
{
  double result;
  result=1.96+2.4/(n-1)+5.901610*pow(n-1,-2.372);
  return result;
}
double Student099(unsigned int n)
{
  double result;
  result=2.576+5.0/(n-1)+29.12178*pow(n-1,-2.591843);
  return result;
}
// this function calculates average x
double Average(void)
{
  double sum=0;
  for (int i=0;i<N;i++) sum+=x[i];
  return sum/N;
}
// this function calculates delta x
double DeltaX(void)
{
  double ave,sum=0,Sn;
  ave=Average();
  for (int i=0;i<N;i++) sum+=x[i]*x[i];
  Sn=sqrt(((double)N/(N-1))*((double)sum/N-ave*ave));
  return Student095(N)*Sn/sqrt(N); // you can use Student099 if you want
}
void main(void)
{
  printf("Data: ");
  for (int i=0;i<N;i++) printf("%f ",x[i]);
  printf("\nx=%f +/- %f",Average(),DeltaX());
}

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating Student's coefficients t(alpha,n), if alpha=0.95 or 0.99
//  (c) Johna Smith, 1996
//
//  Method description:
//   To calculate Student's coefficients we use the following
//  approximation:
//    t(n)=C0+C1/(n-1)+C2(n-1)^(-C3)
//  Coefficients Ci for alpha=0.95 and alpha=0.99 is given
//  n - number of experiments
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
double Student095(unsigned int n)
{
  double result;
  result=1.96+2.4/(n-1)+5.901610*pow(n-1,-2.372);
  return result;
}
double Student099(unsigned int n)
{
  double result;
  result=2.576+5.0/(n-1)+29.12178*pow(n-1,-2.591843);
  return result;
}
void main(void)
{
  printf("alpha=0.95 n=5 t5= %g\n",Student095(5));
  printf("alpha=0.99 n=15 t15= %g",Student099(15));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating theoretical corellation coefficient for two values
//  (c) Johna Smith, 1996
//
//  Method description:
//    Corellation coefficient shows dependencies between two values
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N       3 // number of experiments
double x[N]={1,2,3};
double y[N]={2,4,5};
double TheoreticalCorellation(void)
{
  double sxy=0,sx2=0,sy2=0;
  for (int i=0;i<N;i++)
  {
    sxy+=x[i]*y[i];
    sx2+=x[i]*x[i];
    sy2+=y[i]*y[i];
  }
  return sxy/sqrt(sx2*sy2);
}
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%f %f\n",x[i],y[i]);
  printf("\nCorellation coefficient r= %f\n",TheoreticalCorellation());
}

//////////////////////////////////////////////////////////////////////////////
//
//  Smoothing experimental data by 3 points
//  (c) Johna Smith, 1996
//
//  Method description:
//    We smooth experimental data to enclose them to real curve.
//    All formulas use smallest sqares method to approximate data.
//    This program uses smoothing using 3-points groups
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N 10      // number of experiments
float data[N]={0.9,2.12,2.92,4.15,4.9,6.1,6.92,8.15,9.05,9.8};
float smoothed[N];
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%.4f ",data[i]);
  //smoothing
  smoothed[0]=(5*data[0]+2*data[1]-data[2])/6;
  for (i=1;i<N-1;i++) smoothed[i]=(data[i-1]+data[i]+data[i+1])/3;
  smoothed[N-1]=(5*data[N-1]+2*data[N-2]-data[N-3])/6;
  printf("\nSmoothed data:\n");
  for (i=0;i<N;i++) printf("%.4f ",smoothed[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Smoothing experimental data by 5 points
//  (c) Johna Smith, 1996
//
//  Method description:
//    We smooth experimental data to enclose them to real curve.
//    All formulas use smallest sqares method to approximate data.
//    This program uses smoothing using 5-points groups
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N  10      // number of experiments
float data[N]={0.9,2.12,2.92,4.15,4.9,6.1,6.92,8.15,9.05,9.8};
float smoothed[N];
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%.4f ",data[i]);
  //smoothing
  smoothed[0]=(3*data[0]+2*data[1]+data[2]-data[4])/5;
  smoothed[1]=(4*data[0]+3*data[1]+2*data[2]+data[3])/10;
  for (i=2;i<N-2;i++) smoothed[i]=(data[i-2]+data[i-1]+data[i]+
                                   data[i+1]+data[i+2])/5;
  smoothed[N-2]=(data[N-4]+2*data[N-3]+3*data[N-2]+4*data[N-1])/10;
  smoothed[N-1]=(3*data[N-1]+2*data[N-2]+data[N-3]-data[N-5])/5;
  printf("\nSmoothed data:\n");
  for (i=0;i<N;i++) printf("%.4f ",smoothed[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Smoothing experimental data by 7 points
//  (c) Johna Smith, 1996
//
//  Method description:
//    We smooth experimental data to enclose them to real curve.
//    All formulas use smallest sqares method to approximate data.
//    This program uses smoothing using 7-points groups
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N 10      // number of experiments
float data[N]={0.7,1.08,1.39,1.64,1.76,1.99,2.04,2.22,2.28,2.42};
float smoothed[N];
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%.4f ",data[i]);
  //smoothing
  smoothed[0]=(39*data[0]+8*data[1]-4*(data[2]+data[3]-data[5]) + data[4]-2*data[6])/42;
  smoothed[1]=(8*data[0]+19*data[1]+16*data[2]+6*data[3]-4*data[4] - 7*data[5]+4*data[7])/42;
  smoothed[2]=(-4*data[0]+16*data[1]+19*data[2]+12*data[3]+2*data[4] - 4*data[5]+data[6])/42;
  for (i=3;i<N-3;i++)
    smoothed[i]=(7*data[i]+6*(data[i-1]+data[i+1]) + 
          3*(data[i-2]+data[i+2])-2*(data[i-3]+data[i+3]))/21;
  smoothed[N-3]=(data[N-7]-4*data[N-6]+2*data[N-5]+12*data[N-4] + 
          9*data[N-3]+16*data[N-2]-4*data[N-1])/42;
  smoothed[N-2]=(4*data[N-7]-7*data[N-6]-4*data[N-5]+6*data[N-4] + 
          16*data[N-3]+19*data[N-2]+8*data[N-1])/42;
  smoothed[N-1]=(-2*data[N-7]+4*data[N-6]+data[N-5]-4*data[N-4] - 
          4*data[N-3]+8*data[N-2]+39*data[N-1])/42;
  printf("\nSmoothed data:\n");
  for (i=0;i<N;i++) printf("%.4f ",smoothed[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating empirical corellation coefficient for two values
//  (c) Johna Smith, 1996
//
//  Method description:
//    Corellation coefficient shows dependencies between two values
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N 3 // number of experiments
double x[N]={1,2,3};
double y[N]={2,4,5};
double EmpiricalCorellation(void)
{
  double sxy=0,sx2=0,sy2=0,sx=0,sy=0;
  for (int i=0;i<N;i++)
  {
    sxy+=x[i]*y[i];
    sx2+=x[i]*x[i];
    sy2+=y[i]*y[i];
    sx+=x[i];
    sy+=y[i];
  }
  return (sxy-sx*sy/N)/(sqrt(sx2-sx*sx/N)*sqrt(sy2-sy*sy/N));
}
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%f %f\n",x[i],y[i]);
  printf("\nCorellation coefficient r= %f\n",EmpiricalCorellation());
}

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating corellation coefficient for three values
//  (c) Johna Smith, 1996
//
//  Method description:
//    Corellation coefficient shows dependencies between three values
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N   3 // number of experiments
double x[N]={1,2,3};
double y[N]={2,4,5};
double z[N]={1,7,0};
double EmpiricalCorellation2D(double *i, double *j)
{
  double sij=0,si2=0,sj2=0,si=0,sj=0;
  for (int a=0;a<N;a++)
  {
    sij+=*(i+a)**(j+a);
    si2+=*(i+a)**(i+a);
    sj2+=*(j+a)**(j+a);
    si+=*(i+a);
    sj+=*(j+a);
  }
  return (sij-si*sj/N)/(sqrt(si2-si*si/N)*sqrt(sj2-sj*sj/N));
}
double Corellation3D(void)
{
  double rxz,ryz,rxy;
  // calculating partial corellation coefficients
  rxy=EmpiricalCorellation2D(x,y);
  rxz=EmpiricalCorellation2D(x,z);
  ryz=EmpiricalCorellation2D(y,z);
  return sqrt((rxz*rxz-ryz*ryz-2*rxy*rxz*ryz)/(1-rxy*rxy));
}
void main(void)
{
  printf("Data:\n");
  for (int i=0;i<N;i++) printf("%f %f %f\n",x[i],y[i],z[i]);
  printf("\nCorellation coefficient R= %f\n",Corellation3D());
}

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating partial corellation coefficient for three values
//  (c) Johna Smith, 1996
//
//  Method description:
//   Partial corellation coefficient shows dependencies between two values
//   using set of three values (we exclude influence of one of these values)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
// excluding influence of z, recieving corellation
// coefficients between x,y  y,z and x,z
double PartialCorellation3D(double rxy,double rxz,double ryz)
{
  return (rxy-rxz*ryz)/sqrt((1-rxz*rxz)*(1-ryz*ryz));
}
void main(void)
{
  double rxy=0.998753,rxz=0.69753,ryz=0.7473;
  printf("Data:\n");
  printf("%f %f %f\n",rxy,rxz,ryz);
  printf("\nPartial corellation coefficient r= %.10f\n",
  PartialCorellation3D(rxy,rxz,ryz));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (trapecias method)
//  (c) Johna Smith, 1996
//
//  Method description:
//   We approximate area under function with trapecias, using
//   the same intervals by the X axe.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
double f(double x)
{
  // !!! you should specify your own function here
  return x;
}
double Integrate(double x1, double x2, double step)
{
  double x=x1;
  double I=0; // integral value
  while (x<x2-step)
  {
    // calculating trapecia's area as f(x)+f(x+step) * halfheight of it
    I += (f(x)+f(x+step))*step/2;
    x += step;
  }
  return I;
}
void main(void)
{
  printf("I=%f",Integrate(0,1,1e-4));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (Simpson's method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We approximate area under function with parabolas, using
//    three points. We repeat calculation with less step if
//    the diffference between I(h) and I(h/2) is greater than eps
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
double f(double x)
{
  // you should specify your own function here
  return x*x*sin(x);
}
double Integrate(double x1, double x2, double step, double precision)
{
  double x=x1;
  double I1=0,I2=0; // integral value
  double s;
  do
  {
    I1=I2;
    s=(f(x1)-f(x2))/2;
    x=x1+step;
    while (x<x2)
    {
      s+=2*f(x)+f(x+step);
      x+=2*step;
    }
    I2=2*step*s/3;
    step/=2.0;  // try once more using less step
  } while (fabs(I1-I2)>precision);
  return I2;
}
void main(void)
{
  printf("I=%f",Integrate(0,1,1e-4,1e-6));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (rectanges method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We approximate area under function with rectangles, using
//    the same intervals by the X axe.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
double f(double x)
{
  // you should specify your own function here
  return x;
}
double Integrate(double x1, double x2, double step)
{
  double x=x1;
  double I=0; // integral value
  
  while (x<x2-step)
  {
    // calculating function at the halfpoint of interval
    I+=step*f(x+step/2);
    x+=step;
  }
  
  return I;
}
void main(void)
{
  printf("I=%f",Integrate(0,1,1e-4));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (Monte-Karlo method)
//  (c) Johna Smith, 1996
//
//  Method description:
//   We generate random numbers from interval [a;b] where we want to
//   integrate our function and sum function values in this random points.
//   Here we use approximation of the integral:
//   I=<y>*sigma, where <y> is average value of the function f on [a;b]
//   and sigma is volume of integrating area (if we integrate by interval it
//   is its length)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
double f(double x)
{
  // you should specify your own function here
  return x*x*sin(x);
}
double Integrate(double x1, double x2, long int steps)
{
  double I=0; // integral value
  randomize();  // initializing random numbers generator
  for (long int i=0;i<steps;i++)
    I+=f((random(10000)/10000.0)*(x2-x1)+x1);
  I*=(x2-x1)/steps;
  return I;
}
void main(void)
{
  printf("I=%f",Integrate(0,1,1e5));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating number of days between two dates
//  (c) Johna Smith, 1996
//
//  Method description:
//    We calculte two Julian dates and substract them
//    Result will be number of days between two dates
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
enum Month {January=1,February,March,April,May,June,July,August,
            September,October,November,December};
struct Date
{
  float day;
  Month month;
  int year;
};
double JulianDate(Date date)
{
  long int A,B,Y,M;
  if (date.month==January || date.month==February)
  {
    Y=date.year-1;
    M=date.month+12;
  } else
  {
    Y=date.year;
    M=date.month;
  }
  if ((date.year>1581) ||
      (date.year==1581 && date.month>October) ||
      (date.year==1581 && date.month==October && date.day>=15))
  {
    A=Y/100;
    B=2-A+A/4;
  } else B=0;
  return B+(long int)(365.25*Y)+(long int)(30.6001*(M+1))+date.day+1720994.5;
}
Date Revolution={17.5,October,1917};
Date Today={11,October,1996};
void main(void)
{
  printf("Days between %f/%d/%d and %f/%d/%d = %f",
         Revolution.day,Revolution.month,Revolution.year,
         Today.day,Today.month,Today.year,
         JulianDate(Today)-JulianDate(Revolution));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating day of week by date
//  (c) Johna Smith, 1996
//
//  Method description:
//  Here we use special formulas different for January & February
//  and other monthes. Use this algorithm only for dates after
//  October, 15, 1528
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
enum Month {January=1,February,March,April,May,June,July,August,
            September,October,November,December};
struct Date
{
  float day;
  Month month;
  int year;
};
int DayOfWeek(Date date)
{
  float F;
  if (date.month<March)
    F=365*date.year+date.day+31*(date.month-1)+
      (int)((date.year-1)/4)-(int)(3*(int)((date.year-1)/100+1)/4);
  else
    F=365*date.year+date.day+31*(date.month-1)-(int)(0.4*date.month+2.3)+
      (int)(date.year/4)-(int)(3*(int)(date.year/100+1)/4);
  return  F-7*(int)(F/7)-1;
}
Date Today={11,October,1996};
void main(void)
{
  printf("%f/%d/%d is ",Today.day,Today.month,Today.year);
  switch (DayOfWeek(Today))
  {
    case -1:printf("Sunday\n");break;
    case 0:printf("Monday\n");break;
    case 1:printf("Tuesday\n");break;
    case 2:printf("Wednesday\n");break;
    case 3:printf("Thursday\n");break;
    case 4:printf("Friday\n");break;
    case 5:printf("Saturday\n");break;
  }
}

//////////////////////////////////////////////////////////////////////////////
//
//  Julian date calculating
//  (c) Johna Smith, 1996
//
//  Method description:
//    Julian date is the number of days since afternoon (GMT) of
//    January, 1, 4713 BC.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
enum Month {January=1,February,March,April,May,June,July,August,
            September,October,November,December};
struct Date
{
  float day;
  Month month;
  int year;
};
double JulianDate(Date date)
{
  long int A,B,Y,M;
  if (date.month==January || date.month==February)
  {
    Y=date.year-1;
    M=date.month+12;
  } else
  {
    Y=date.year;
    M=date.month;
  }
  if ((date.year>1581) ||
      (date.year==1581 && date.month>October) ||
      (date.year==1581 && date.month==October && date.day>=15))
  {
    A=Y/100;
    B=2-A+A/4;
  } 
  else 
    B=0;
  return B+(long int)(365.25*Y)+(long int)(30.6001*(M+1))+date.day+1720994.5;
}
Date Revolution={17.5,October,1917};
void main(void)
{
  printf("%f/%d/%d = %f",Revolution.day,Revolution.month,
                         Revolution.year,JulianDate(Revolution));
}

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (square root method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    From given matrix A we build two auxulary triangular matrixes S & D
//    Then solving equations Gy=b (where g[I][j]=s[j][I]*d[j]) and Sx=y
//    we obtain vector x.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N    4     // size of matrix
#define N1   N+1
float matrix[N][N1]=
    {{10.9,1.2,2.1,0.9,23.2},
  {1.2,11.2,1.5,2.5,38.1},
  {2.7,1.5,9.8,1.3,40.3},
  {0.9,2.5,1.3,12.1,58.2}
    };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  register char i,j,k;
  float s[N][N],x[N],y[N],d[N],sum;
  // Printing given matrix
  ShowMatrix();
  // Building matrixes S and D
  for (j=0;j<N;j++)
    for (i=0;i<=j;i++)
    {
      sum=matrix[i][j];
      for (k=0;k<i;k++)
        sum-=s[k][i]*d[k]*s[k][j];
      if (i==j)
      {
        d[i]=(sum>0?1:-1);
        s[i][j]=sqrt(fabs(sum));
      }
      else s[i][j]=sum/(d[i]*s[i][i]);
    }
  // Solving equation Gy=b (G: g[I][j]=s[j][I]*d[j])
  for (i=0;i<N;i++)
  {
    y[i]=matrix[i][N];
    for (j=0;j<i;j++) y[i]-=s[j][i]*d[j]*y[j];
    y[i]/=s[i][i]*d[i];
  }
  // Solving equation Sx=y
  for (i=N-1;i>=0;i--)
  {
    x[i]=y[i];
    for (j=i+1;j<N;j++) x[i]-=s[i][j]*x[j];
    x[i]/=s[i][i];
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
  printf("x%d=%f\n",i+1,x[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Nekrasov method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This is an iterative method of solving systems of linear equations.
//    Criteria for end of iterations is ||xn||-||x(n-1)||<epsilon, where
//    ||x|| is norm of vector x.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N     4     // size of matrix
#define N1    N+1
float matrix[N][N1]=
      {{10.9,1.2,2.1,0.9,23.2},
       {1.2,11.2,1.5,2.5,38.1},
       {2.7,1.5,9.8,1.3,40.3},
       {0.9,2.5,1.3,12.1,58.2}
     };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
// this function calculates ||x||
float norma(float *x)
{
  float sum=0;
  for (unsigned i=0;i<N;i++) sum+=fabs(*(x+i));
  return sum;
}
void main(void)
{
  // Variables declaration
  register char i,j;
  float x_current[N],x_next[N],sum1,sum2,epsilon;
  epsilon=1e-7;
  // Printing given matrix
  ShowMatrix();
  // Solving
  do
  {
    for (i=0;i<N;i++) x_current[i]=x_next[i];
    for (i=0;i<N;i++)
    {
      sum1=0;sum2=0;
      for(j=0;j<i;j++) sum1+=matrix[i][j]*x_next[j];
      for(j=i+1;j<N;j++) sum2+=matrix[i][j]*x_current[j];
      x_next[i]=(matrix[i][N]-sum1-sum2)/matrix[i][i];
    }
  } while (fabs(norma(x_current)-norma(x_next))>epsilon);
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x_current[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (3-diagonal matrix)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This method is developed for 3-diagonal matrixes:
//    (a11 a12 0 ..............)
//    (a21 a22 a23 0 ..........)
//    (0   a32 a33 a34 0 ......) = M
//    (........................)
//    (.......... 0 an(n-1) ann)
//    Cause all elements except three diagonals is 0 we'll store
//    only non zero elements in the A,B,C arrays. Array D stores
//    vector B (M*X=B). We build two auxulary vectors P and Q and
//    then calculate vector X using the following recurrent formula:
//      X(k-1)=P(k)*X(k)+Q(k)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#define N    20
#define N1   N+1
void main(void)
{
  // Variables declaration
  register int i,i1,k;  // using processor registers for better performance
  float a[N]={0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
  float b[N]={2,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,2};
  float c[N]={-1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0};
  float d[N]={0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0};
  float x[N1],p[N1],q[N1],tmp;
  // Some initializations
  p[0]=0;
  q[0]=0;
  x[N]=0;
  c[N-1]=0; // to avoid overflow
  // Building vectors P and Q
  for (i=0;i<N;i++)
  {
    i1=i+1;
    if ((tmp=b[i]+a[i]*p[i])==0)
    {
      printf("Wrong method for such system...\n");
      return;
    };
    p[i1]=-c[i]/tmp;
    q[i1]=(d[i]-a[i]*q[i])/tmp;
  }
  // Building vector X
  for (k=N;k>0;k--) x[k-1]=p[k]*x[k]+q[k];
  // Printing solution
  printf("Solution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Gauss method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This methos based on excluding variables from equations.
//    Using first equation we exclude x1 from all other equations
//    Using second equation equation we exclude x2 from all equations
//   (excluding first) and so on. So we have triangular matrix from which
//    we can easily get vector X. (AX=B)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N     4     // size of matrix
#define N1    N+1
float matrix[N][N1]=
                  {{2,5,4,1,28},
                   {1,3,2,1,17},
                   {2,10,9,7,77},
                   {3,8,9,2,54}
                  };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float tmp,x[N1];
  register short int i,j,k;
  // Printing given matrix
  ShowMatrix();
  // Main loop
  for (i=0;i<N;i++)
  {
    // Excluding variable x[i] from equations
    tmp=matrix[i][i];
    for (j=N;j>=i;j--) matrix[i][j]/=tmp;
    for (j=i+1;j<N;j++)
    {
      tmp=matrix[j][i];
      for (k=N;k>=i;k--)
        matrix[j][k]-=tmp*matrix[i][k];
    }
  }
  // Calculating vector x
  x[N-1]=matrix[N-1][N];
  for (i=N-2;i>=0;i--)
  {
    x[i]=matrix[i][N];
    for (j=i+1;j<N;j++) x[i]-=matrix[i][j]*x[j];
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Gauss method with max element selecting)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This methos based on excluding variables from equations.
//    Using first equation we exclude x1 from all other equations
//    Using second equation equation we exclude x2 from all equations
//    (excluding first) and so on. So we have triangular matrix from which
//    we can easily get vector X. (AX=B)
//    'Selecting max element' means that on each step we select equation
//    with max element in current column and use it as next equation to exclude
//    variable.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N      4     // size of matrix
#define N1     N+1
float matrix[N][N1]=
    {{2,5,4,1,28},
  {1,3,2,1,17},
   {2,10,9,7,77},
  {3,8,9,2,54}
        };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float max,swp,tmp,x[N1];
  register short int row_with_max_element,i,j,k;
  // Printing given matrix
  ShowMatrix();
  // Main loop
  for (i=0;i<N;i++)
  {
    // Searching for max element in the current column (i)
    max=fabs(matrix[i][i]);
    row_with_max_element=i;
    for (j=i+1;j<N;j++)
      if (fabs(matrix[j][i])>max)
      {
        row_with_max_element=j;
        max=fabs(matrix[j][i]);
      }
    // Swapping 2 lines of matrix - row_with_max_element & i
    if (row_with_max_element!=i)
    for (j=0;j<N1;j++)
    {
      swp=matrix[row_with_max_element][j];
      matrix[row_with_max_element][j]=matrix[i][j];
      matrix[i][j]=swp;
    }
    // Excluding variable x[i] from equations
    tmp=matrix[i][i];
    for (j=N;j>=i;j--) matrix[i][j]/=tmp;
    for (j=i+1;j<N;j++)
    {
      tmp=matrix[j][i];
      for (k=N;k>=i;k--)
        matrix[j][k]-=tmp*matrix[i][k];
    }
  }
  // Calculating vector x
  x[N-1]=matrix[N-1][N];
  for (i=N-2;i>=0;i--)
  {
    x[i]=matrix[i][N];
    for (j=i+1;j<N;j++) x[i]-=matrix[i][j]*x[j];
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (iterations method)
//  (c) Johna Smith, 1996
//
//  Method description:
//  This program automatically transforms system of equations to
//  iterative form and solves it.
//  Iterative form:
//     x1=a11*x1+a12*x2+... + a1n*xn+b1
//     x2=a21*x1+a22*x2+... + a2n*xn+b2
//     ...
//     xn=an1*x1+an2*x2+... + ann*xn+bn
//  Setting first iteration of vector X we can get next iteration from
//  such form of system of equations. Iterations will be finished when
//  difference between two consequitive iterations will be less than epsilon
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N      3     // size of matrix
#define N1     N+1
float matrix[N][N1]=
    {{14.38,-2.41,1.39,5.86},
  {1.84,25.36,-3.31,-2.28},
  {2.46,-3.49,16.37,4.47}
    };
int maxiterations=10;  // maximum number of iterations
float epsilon=0.0001;  // required accuracy
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float x[N],y[N],t;
  register short int i,j,k;
  int iterations=0;
  // Printing given matrix
  ShowMatrix();
  // setting first iteration of vector X
  for (i=0;i<N;i++) x[i]=matrix[i][N];
  do
  {
    for (i=0;i<N;i++)
    {
      t=-matrix[i][N];
      for (j=0;j<N;j++)
        t+=matrix[i][j]*x[j];
        y[i]=(-t+matrix[i][i]*x[i])/matrix[i][i];
    }
    iterations++;
    k=0;
    // checking solution
    while (fabs(x[k]-y[k])<epsilon && k<N) k++;
    // new iteration becomes old
    for(i=0;i<N;i++) x[i]=y[i];
  } while (k!=N && iterations<maxiterations);
  if (iterations==maxiterations)
  {
    printf("Iterations are very slow...");
  } else
  {
    // Printing solution
    printf("\nSolution:\n");
    for (i=0;i<N;i++)
      printf("x%d=%f\n",i+1,x[i]);
    printf("%d iterations were made",iterations);
  }
}

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Monte-Karlo method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We should translate given system of linear equations to special form:
//      x1=a(11)x1+a(12)x2+...+a(1n)xn+a(1n+1)
//      x2=a(21)x2+a(22)x2+...+a(2n)xn+a(2n+1)
//      ....
//      xn=a(n1)x1+a(n2)x2+...+a(nn)xn+a(nn+1)
//            n
//     where SUM |aij|<1 (i=1, 2, ... ,n)
//           j=1
//    Lets divide interval [0;1] into N+1 smaller intervals.
//    Imagine a particle that is moving randomly at [0;1] interval.
//    Remember its moves until it reachs one of the bounds of the interval.
//    So we'll get a trajectory of the particle s(i),s(j),s(k),s(l),s(m),...
//    If particle moves from interval s(i) to s(j) we'll write v(ij)
//    y[i]=v(ij)*v(jk)*...*v(tm)*w(m), v(ij)=sign(aij), v(jk)=sign(ajk)
//                    n
//   w(m)=a(mn+1)/(1-SUM a(mj)),   x[i]=y/M, where M is number of particle runs
//                   j=1
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N     3     // size of matrix
#define N1    N+1
#define M     1000 // number of tries
float matrix[N][N1]=
    {{3,1,-1,-2},
  {1,4,-2,-3},
  {2,-1,5,-15}
    };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    printf("x%d=",i+1);
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],j+1);
    printf("%+f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float b[N][N],w[N],y,c,x[N],l;
  register short int i,j,v;
  unsigned char finish;
  int row,   // currently analysed row
      tries; // number of tries to calculate true root by generating randoms
  // Transforming matrix
  for (i=0;i<N;i++)
  {
    l=0;
    for (j=0;j<N1;j++) l+=fabs(matrix[i][j]);
    v=(matrix[i][i]>0?-1:1);
    for (j=0;j<N1;j++) matrix[i][j]=v*(matrix[i][j])/l+(i==j?1:0);
  }
  // Printing transformed matrix
  ShowMatrix();
  // filling matrix B
  for (i=0;i<N;i++)
  {
    b[i][0]=fabs(matrix[i][0]);
    for(j=1;j<N;j++)
      b[i][j]=b[i][j-1]+fabs(matrix[i][j]);
  }
  // filling matrix w
  for(i=0;i<N;i++) w[i]=matrix[i][N]/(1-b[i][N-1]);
  // for each equation in the system
  for(i=0;i<N;i++)
  {
    tries=0;
    row=i;
    y=0;
    v=1;
    while (tries<M)
    {
      // generating random number
      c=random(32767)/32767.0;
      j=N-1;
      finish=0;
      while(j>=0 && !finish)
      {
        if (c>b[row][j])
        {
          if (j==N-1)
          {
            tries++;
            y+=v*w[row];
            row=i;v=1;
          }
          else
          {
            v*=(matrix[row][j+1]>=0?1:-1);
            row=j+1;
          }
          finish=1;
        }
        j--;
      }
      if (!finish)
      {
        v*=(matrix[row][j+1]>=0?1:-1);
        row=0;
      }
    }
    // calculating root
    x[i]=y/M;
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (orthogonalization method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given system of linear equations:
//      a11 x1 + a12 x2 + ... + a1n xn + a(1n+1) = 0
//      a21 x1 + a22 x2 + ... + a2n xn + a(2n+1) = 0
//
//      an1 x1 + an2 x2 + ... + ann xn + a(nn+1) = 0
//
//    Left part of each equation is a result of scalar multiplication of
//    two vectors: ai=(ai1,ai2,...,ain,a(in+1)) and x=(x1,x2,..,xn,1)
//    So to solve this system we need only to build vector x which will be
//    orthogonal to each of ai vectors.
//                            n+1
//    Let u1=a1,   b1=u1/sqrt(SUM u(1j)^2)
//                            j=1
//
//    Other ui we can get from the following iterative formula:
//                      k
//       u(k+1)=u(k+1)-SUM {u(k+1), bj}bj, where {..} means scalar multiplication
//                     j=1
//                          n+1
//       b(k+1)=u(k+1)/sqrt(SUM u(k+1,j)^2)
//                          j=1
//    And finally we can obtain roots of this system from this formula:
//
//    xi=b(n+1,i)/b(n+1n+1)
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N     3     // size of matrix
#define N1    N+1
double matrix[N1][N1]=
    {{3,1,-1,-2},
  {1,4,-2,-3},
  {2,-1,5,-15}
    };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("%+f=0\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  double c[N],x[N],s;
  int i,j,k,l,m;
  ShowMatrix();
  // expanding system by vector (0,0,0,...,0,1)
  for (i=0;i<N;i++) matrix[N][i]=0;
  matrix[N][N]=1;
  // for all equations
  for (i=0;i<N1;i++)
  {
    if (i>0)
    {
      k=0;
      // make some iterations to increase accuracy of calculations
      while (k<=3)
      {
        for (l=0;l<i;l++)
        {
          c[l]=0;
          for(m=0;m<N1;m++) c[l]+=matrix[l][m]*matrix[i][m];
        }
        for (j=0;j<N1;j++)
        {
          s=0;
          for(l=0;l<i;l++) s+=c[l]*matrix[l][j];
          matrix[i][j]-=s;
        }
        k++;
      }
    }
     // normalizing vector
     s=0;
     for (k=0;k<N1;k++) s+=matrix[i][k]*matrix[i][k];
     s=sqrt(s);
     for (j=0;j<N1;j++) matrix[i][j]/=s;
  }
  s=matrix[N][N];
  for (j=0;j<N;j++) x[j]=matrix[N][j]/s;
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
}

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Zaidel method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: system of linear equations
//       a11 x1 + a12 x2 + ... + a1n xn = b1
//       a21 x1 + a22 x2 + ... + a2n xn = b2
//       ...
//       an1 x1 + an2 x2 + ... + ann xn = bn
//
//    We use following iterative formula to solve this system:
//
//                     1    i-1                N
//     xi(j+1)=xi(j)- --- ( SUM aik xk(j+1) + SUM aik xk(j) -bi)
//                    aii   k=1               k=i
//
//     where xi(j+1) is j+1-th iteration, xi(j) is j-th iteration
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N    3     // size of matrix
#define N1   N+1
float matrix[N][N1]=
    {{10,1,1,12},
  {2,10,1,13},
  {2,2,10,14}
    };
float z[N]={0,0,0};  // first iteration
float epsilon=0.0001;  // required accuracy
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
void main(void)
{
  // Variables declaration
  float x[N];
  short int i,j;
  int iterations=0,finish=0;
  // Printing given matrix
  ShowMatrix();
  while (!finish)
  {
    finish=1;
    for (i=0;i<N;i++)
    {
      x[i]=-matrix[i][N];
      for (j=0;j<N;j++) x[i]+=matrix[i][j]*z[j];
      // don't stop iterations until required accuracy is reached
      if (fabs(x[i]/matrix[i][i])>=epsilon) finish=0;
      x[i]=z[i]-x[i]/matrix[i][i];
      // next iteration
      z[i]=x[i];
    }
    iterations++;
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
    printf("x%d=%f\n",i+1,x[i]);
  printf("%d iterations were made",iterations);
}

//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (chords method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    given: f(x)=0, [a;b], f(a)*f(b)<0
//    find x0: f(x0)=0
//    We approximate curve by chord between its borders and get
//    point where the chord crosses X axe as a new border of
//    interval with the root. When interval becomes small enough
//    we can take one of its borders as a solution.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
// this function returns value of f(x)
double f(double x)
{
  // sin 2x - ln x = 0
  return sin(2*x)-log(x);
}
double Solve(double a, double b, double epsilon)
{
  double u,v,x,y;
  u=f(a); v=f(b);
  do
  {
    // x - new border
    x=(a*v-b*u)/(v-u);
    y=f(x);
    // determine whether x is left or right border
    // f(a)*f(b)<0 condition must be true
    if (y*u<0)
    {
      b=x;
      v=y;
    } else
    {
      a=x;
      u=y;
    }
  } while (b-a>epsilon);
  return x;
}
void main(void)
{
  printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
  printf("\nx0=%f",Solve(1.3,1.5,1e-9));
}

//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (dichotomia method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    given: f(x)=0, [a;b], f(a)*f(b)<0
//    find x0: f(x0)=0
//    We split interval into two parts and select that part that
//    contains root of the equation by checking condition
//    f(a)*f(b)<0 We take middle point of the interval as a new border of
//    interval with the root. When interval becomes small enough
//    we can take one of its borders as a solution.
//
//////////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
// this function returns value of f(x)
double f(double x)
{
  // sin 2x - ln x = 0
  return sin(2*x)-log(x);
}
double Solve(double a, double b, double epsilon)
{
  double x;
  do
  {
    // x - new border
    x=(b+a)/2;
    // determine whether x is left or right border
    // f(a)*f(b)<0 condition must be true
    if (f(x)*f(a)<0) b=x;
    else a=x;
  } while (b-a>epsilon);
  return x;
}
void main(void)
{
  printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
  printf("\nx0=%f",Solve(1.3,1.5,1e-9));
}