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

#define DIMENSION 7

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

   n = DIMENSION;

   k = (int) ceil( log(n) / log(2) );
   of = 1<<(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 ));
   randmax=2147483647.0 ; 

   for(i = 0      ; i < of+1   ; i++)  a[i] = 0.0;
   for(i = of+1   ; i < of+n   ; i++)  a[i] = rand()/randmax;
   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] = rand()/randmax;
   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] = rand()/randmax;
   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] = rand()/randmax;
   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]);
   }
}