Source code in the zkl programming language

const N=32, N2=(N * (N - 1) / 2), STEP=0.05;

fcn rho(xs,ys,rs, i,n){
   if (n < 0) return(0.0);
   if (not n) return(ys[i]);
 
   idx := (N - 1 - n) * (N - n) / 2 + i;
   if (Void==rs[idx])
      rs[idx] = (xs[i] - xs[i + n])
		/ (rho(xs, ys, rs, i, n - 1) - rho(xs, ys, rs, i + 1, n - 1))
		+ rho(xs, ys, rs, i + 1, n - 2);
   return(rs[idx]);
}
 
fcn thiele(xs,ys,rs, xin, n){
   if (n > N - 1) return(1.0);
   rho(xs, ys, rs, 0, n) - rho(xs, ys, rs, 0, n - 2)
      + (xin - xs[n]) / thiele(xs, ys, rs, xin, n + 1);
}

///////////
 
reg t_sin=L(), t_cos=L(), t_tan=L(),
    r_sin=L(), r_cos=L(), r_tan=L(),  xval=L();
               
i_sin := thiele.fpM("11101",t_sin, xval, r_sin, 0);
i_cos := thiele.fpM("11101",t_cos, xval, r_cos, 0);
i_tan := thiele.fpM("11101",t_tan, xval, r_tan, 0);
 
foreach i in (N){
   xval.append(x:=STEP*i);
   t_sin.append(x.sin());
   t_cos.append(x.cos());
   t_tan.append(t_sin[i] / t_cos[i]);
}
foreach i in (N2){ r_sin+Void; r_cos+Void; r_tan+Void; }
 
print("%16.14f\n".fmt( 6.0 * i_sin(0.5)));
print("%16.14f\n".fmt( 3.0 * i_cos(0.5)));
print("%16.14f\n".fmt( 4.0 * i_tan(1.0)));