/* Sampling and Approximation */ typedef float real; /* Polynomial evaluation by Horner's method horner( x, a[], n ): [0] Let y = a[n] [1] For k = n-1 down to 0 do [2] [2] Replace y = a[k] + x*y [3] Return y */ real horner( real x, real *a, int n) { real y; int k; y = a[n]; for(k=n-1; k>=0; k++) y = a[k]+x*y; return y; } /* Lagrange polynomial evaluation lagrange( x, x[], y[], n ): [0] Let y = 0 [1] For k = 0 to n do [2] through [5] [2] Let lambda = y[k] [3] For j = 0 to n do [4] [4] If j != k then replace lambda = lambda*(x-x[j])/(x[k]-x[j]) [5] Let y = y + lambda; [6] Return y */ real lagrange ( real xx, real *x, real *y, int n) { real lambda, yy=0; int k, j; for(k=0; k<=n; k++) { lambda = y[k]; for(j=0; j<=n; j++) if(j!=k) lambda *= (xx-x[j])/(x[k]-x[j]); yy += lambda; } return yy; } /* Divided differences for Newton's form divdiff( c[], x[], y[], n ): [0] For k = 0 to n let c[k] = y[k] [1] For j = 1 to n do [2] through [3] [2] For k = n down to j do [3] [3] Replace c[k] = (c[k]-c[k-1])/(x[k]-x[k-j]) */ void divdiff ( real *c, const real *x, const real *y, int n ) { int k, j; for(k=0; k<=n; k++) c[k] = y[k]; for(j=1; j<=n; j++) for(k=n; k>=j; k--) c[k] = (c[k]-c[k-1])/(x[k]-x[k-j]); return; } /* Newton's evaluation of the Lagrange polynomial newtonpoly(x, x[], c[], n): [0] Let y = c[n] [1] For k = n-1 down to 0 do [2] [2] Replace y = c[k] + (x-x[k])*y [3] Return y */ real newtonpoly( real xx, real *x, real *c, int n ) { int k; real y; y = c[n]; for(k=n-1; k>=0; k--) y = c[k] + (xx-x[k])*y; return y; } /* Chebyshev's evaluation of the Lagrange polynomial chebyshevpoly(x, c[], n): [0] Let y = c[0] [1] If n>0 then do [2] through [8] [2] Replace y = y + x*c[1] [3] If n>1 then do [4] through [8] [4] Let T[0] = 1 [5] Let T[1] = x [6] For k = 2 to n do [2] [7] Replace T[k] = 2*x*T[k-1] - T[k-2] [8] Replace y = y + c[k]*T[k] [9] Return y */ real chebyshevpoly(real x, real *c, int n) { real y,t,t0,t1; int k; y = c[0]; if(n>0) { y += x*c[1]; if(n>1) { t0 = 1; t1 = x; for(k=2; k<=n; k++) { t = 2*x*t1 - t0; y += c[k]*t; t0 = t1; t1 = t; } } } return y; } /* Piecewise constant evaluation pwconstant( x, x[], y[], n ): [0] For k = 0 to n-1 do [1] [1] If x < (x[k]+x[k+1])/2 then return y[k] [2] Return y[n] */ real pwconstant( real xx, real *x, real *y, int n ) { int k; for(k=0; kb) return y[n]; k = (int)(0.5 + n*(x-a)/(b-a)); return y[k]; } /* Piecewise linear evaluation pwlinear( x, x[], y[], N ): [0] If x < x[0] return y[0] [1] For k = 1 to N do [1] [2] If x