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

double random_number(double *rmin, double *rmax){
  return (*rmin) + (float)rand()/(float)RAND_MAX*((*rmax)-(*rmin));
}

int main() {
  double *a, *b, *c, *d, *x, *a_k, *b_k, *c_k, *d_k;
  double rmin,rmax;
  int i,j,k,n,of,im,ip;

  n = 7;

  k = (int) ceil( log(n) / log(2) );
  of = 1<<(k-1) ;     /* bit shift: equivalent 2^(k-1) */

  a = (double *)malloc( (n+2*of) * sizeof(double) );
  b = (double *)malloc( (n+2*of) * sizeof(double) );
  c = (double *)malloc( (n+2*of) * sizeof(double) );
  d = (double *)malloc( (n+2*of) * sizeof(double) );

  srand( (unsigned)time( (long *)0 ));

  for(i = 0      ; i < of+1   ; i++)  a[i] = 0.0;
  for(i = of+1   ; i < of+n   ; i++)  a[i] = random_number(&rmin,&rmax);
  for(i = of+n   ; i < 2*of+n ; i++)  a[i] = 0.0;

  for(i = 0      ; i < of     ; i++)  b[i] = 1.0;
  for(i = of     ; i < of+n   ; i++)  b[i] = random_number(&rmin,&rmax);
  for(i = of+n   ; i < 2*of+n ; i++)  b[i] = 1.0;

  for(i = 0      ; i < of     ; i++)  c[i] = 0.0;
  for(i = of     ; i < of+n-1 ; i++)  c[i] = random_number(&rmin,&rmax);
  for(i = of+n-1 ; i < 2*of+n ; i++)  c[i] = 0.0;

  for(i = 0      ; i < of     ; i++)  d[i] = 0.0;
  for(i = of     ; i < of+n   ; i++)  d[i] = random_number(&rmin,&rmax);
  for(i = of+n   ; i < 2*of+n ; i++)  d[i] = 0.0;

  a_k = (double *)malloc( n * sizeof(double) );
  b_k = (double *)malloc( n * sizeof(double) );
  c_k = (double *)malloc( n * sizeof(double) );
  d_k = (double *)malloc( n * sizeof(double) );

  for(j=1; j<k+1; j++) {
    for(i=0; i<n; i++) {
      im = i + of - ( 1<<(j-1) );
      ip = i + of + ( 1<<(j-1) );
      a_k[i] = -a[i+of] / b[im] * a[im];
      c_k[i] = -c[i+of] / b[ip] * c[ip];
      b_k[i] = -a[i+of] / b[im] * c[im] + b[i+of] - c[i+of] / b[ip] * a[ip];
      d_k[i] = -a[i+of] / b[im] * d[im] + d[i+of] - c[i+of] / b[ip] * d[ip];
    }
    memcpy( &(a[of]) , a_k , n * sizeof(double) );
    memcpy( &(b[of]) , b_k , n * sizeof(double) );
    memcpy( &(c[of]) , c_k , n * sizeof(double) );
    memcpy( &(d[of]) , d_k , n * sizeof(double) );
 }

 x = (double *)malloc( n * sizeof(double) );

 for(i=0; i<n; i++) {    
   x[i] = d[i+of] / b[i+of]; 
   printf("x[%d] = %6.2f\n",i,x[i]);
 }
}