// gcc a.c -Wall -std=c99 -lm
#include <stdio.h>
typedef double R;
typedef struct {R l; R d; R u;} tr;
void solve(int n, tr A[n], R b[n], R x[n]){
  void s(R p, R q, int j) {
    printf("I am %e + you*%e\n", p, q);
    R r = 1./(A[j].l*q + A[j].d), mec = b[j] - A[j].l*p;
    if(j==n-1) {x[j] = r*mec; return;}
    s(r*mec, - A[j].u*r, j+1);
    x[j] = r*(mec - A[j].u*x[j+1]);
    printf("I am %e\n", x[j]);
    return;}
  R r = 1./A[0].d, p = r * b[0], q = -r * A[0].u;
  s(p, q, 1); x[0] = p + x[1]*q;}

void ml(int n, tr A[n], R x[n], R b[n]){
  b[0] = A[0].d*x[0] + A[0].u*x[1];
  for(int j = 1; j<n-1; ++j)
    b[j] = A[j].l*x[j-1] + A[j].d*x[j] + A[j].u*x[j+1];
  b[n-1] = A[n-1].l*x[n-2] + A[n-1].d*x[n-1];}

#include <math.h>
double fmax(double, double); // should not be necessary
#define z 4
int main(){
  tr t[] = {
    {0, 1.4, 0.21},
    {0.21, 0.3, 0.18},
    {0.1, 1.6, 0.14},
    {0.16, -1.44, 0}};
  R x[] = {2, 6, 5, 7};
  R b[z], xz[z];
  ml   (z, t, x, b);
  solve(z, t, b, xz);
  {R mx = 0;
    for(int j=0; j<z; ++j) {
      printf("%e %e %e\n", x[j], b[j], xz[j]);
      mx = fmax(mx, fabs(x[j] - xz[j]));}
  printf("max error = %e\n", mx);}
  return 0;}