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

#define A 16807
#define M 2147483647
#define Q 127773
#define R 2836
#define S 260572

#define N 10

float  random(int *p);
double chi_0(float *x,float *y,double *coef,int n);
double chi_min(float *x,float *y,double *coef,int n);

int main()
{
  int i,n=N,semilla=S;
  int *p=&semilla;
  float e,s2,data_x[N],data_y[N];
  float *x=&data_x[0],*y=&data_y[0];
  double salida[4],*coef=&salida[0];
  
  for(i=0;i<N;i++) 
    {
      *(x+i)=1.0*i;
      e=random(p);
      e=e-0.5;
      *(y+i)=6.0*i+0.1*e;
      /*printf("%f\t%f\n",*(x+i),*(y+i));*/
    }

  s2=chi_0(x,y,coef,n);
  s2=chi_min(x,y,coef,n);

  return 1;

}

double chi_min(float *x,float *y,double *coef,int n)
{
  /*
  Devuelve el coeficiente minimo cuadratico con ordenada.
  */

  int i;
  float s_x=0.0,s_y=0.0,s_xx=0.0,s_xy=0.0;
  double delta,s,s2=0.0;


  for(i=0;i<n;i++)
    {
      s_x+=*(x+i);
      s_y+=*(y+i);
      s_xx+=*(x+i)*(*(x+i));
      s_xy+=*(x+i)*(*(y+i));
    }

  delta=(double)(s_xx*n-s_x*s_x);
  *coef=(double)(s_xx*s_y-s_x*s_xy)/delta;
  *(coef+2)=(double)(s_xy*n-s_x*s_y)/delta;
  
  for(i=0;i<n;i++)
    {
      s=*(y+i)-(*coef)-(*(coef+2))*(*(x+i));
      s2+=s*s;
    }

  s2=s2/(n-2.0);

  *(coef+1)=s2*s_xx/delta;
  *(coef+4)=s2*(1.0*n)/delta;

  printf("a=\t%f\n",*coef);
  printf("var_a=\t%f\n",*(coef+1));
  printf("b=\t%f\n",*(coef+2));
  printf("var_b=\t%f\n",*(coef+4));
  printf("s2=\t%f\n",s2);

  return s2;

}


double chi_0(float *x,float *y,double *coef,int n)
{
  /*
  Devuelve el coeficiente minimo cuadratico sin ordenada.
  */

  int i;
  float s_xx=0.0,s_xy=0.0;
  double s,s2=0.0;

  for(i=0;i<n;i++)
    {
      s_xx+=*(x+i)*(*(x+i));
      s_xy+=*(x+i)*(*(y+i));
    }

  *coef=(double)(s_xy)/s_xx;
  
  for(i=0;i<n;i++)
    {
      s=*(y+i)-(*coef)*(*(x+i));
      s2+=s*s;
    }

  s2=s2/(n-2.0);

  *(coef+1)=s2/s_xx;

  printf("b=\t%f\n",*coef);
  printf("var_b=\t%f\n",*(coef+1));
  printf("s2=\t%f\n",s2);

  return s2;

}

float random(int *semilla)
{

/* 
     Devuelve un numero aleatorio entre 0.0 y 1.0
     a partir de una semilla dada. Se usa el metodo
     de Park y Miller
*/

int k;
float x;

k=(*semilla)/Q;
*semilla=A*(*semilla-k*Q)-R*k;
if (*semilla<0) *semilla+=M;

x=(*semilla)*(1.0/M);

return x;
}