// To solve the equations described at
// http://cap-lore.com/MathPhys/Implicit/1D.html
// In routine solve argument n is the number of unknowns, excludng the ends.
typedef double D;
typedef struct{D x; D y; D z;} mr;
void solve(int n, mr L[n], D fh, D lh, D ans[n+1]){
  void qry(D A, D B, int k){
    if(k==n) return; 
    else {
      D r = 1./(1. - L[k].x*B), a = (L[k].x*A + L[k].y)*r, b=L[k].z*r;
      qry(a, b, k+1);
      ans[k] = a + b*ans[k+1];}}
  qry(fh, 0, 0);}

#if 1
#include <stdio.h>
#include <math.h>
int main(){
  D fh = 4., lh = 5.;     
  void tx(int n, mr t[n]){
    D ans[n+2]; ans[0]=fh; ans[n+1]=lh;
    solve(n, t, fh, lh, ans+1);
    {int j=n+2; while(j--) {
       printf("%15.10f", ans[j]);
       if(j&&j<=n) printf(" %15.10f %15.10f %15.10f",
         t[j-1].x, t[j-1].y, t[j-1].z);
       printf("\n");}}
     {int j = n; while(j--) {
       printf ("%e ", t[j].x*ans[j] + t[j].y + t[j].z*ans[j+2] - ans[j+1]);}}
     printf("\n");}
  {mr t[] = {{1,3,2}, {2,7,3}}; tx(2, t);}
  {mr t[] = {{0.03, 3, 0.02}, {0.034, 3.6, 0.28}, {0.022, 4, 0.018}};
      tx(3, t);}
  
  return 0;}

/*
   5.0000000000
 -12.0000000000    2.0000000000    7.0000000000    3.0000000000
 -17.0000000000    1.0000000000    3.0000000000    2.0000000000
   4.0000000000
0.000000e+00 0.000000e+00 
   5.0000000000
   4.1974632083    0.0220000000    4.0000000000    0.0180000000
   4.8846912884    0.0340000000    3.6000000000    0.2800000000
   3.2176938258    0.0300000000    3.0000000000    0.0200000000
   4.0000000000
-8.881784e-16 0.000000e+00 0.000000e+00 
*/
#endif