Listing 1 (s)

/*
** bc program to calculate Chebyshef economized polynomial
** for evaluation of sin(x) */
/* use bc -1 to get c() and s() functions */
define t(x) {                           /* sin(x)/x */
      if(x==0)return(1.); /* derivative of s function */
   return (s(x) / x); /* put function to be fit here */ }
define b(x) {
   if (x < 0) return (-x);
   return (x); }
define m(x, y) {
   if (x > y) return (x);
   return (y); }
n = 22; /* number of Chebyshef terms */
scale = 40;
p = a(1.) * 4; /* pi */
b = p * .5; /* upper end of curve fit interval */
a = -b; /* lower end of interval */
/* chebft adapted from Press Flannery et al */
/* "Numerical Recipes" FORTRAN version */
for (k = 1; k <= n; ++k) {
c[k] = 0;
f[k] = t(c((k - .5) * p / n) * (b - a) * .5 + (b + a) * .5);
}
/* because of symmetry, even c[] are 0 */
for (j = 1; j <= n; j += 2) {
   s = 0;
   q = (j - 1) * p / n;
   for (k = 1; k <= n; ++k) s += c(q * (k - .5)) * f[k];
   (c[j] = 2 / n * s); }
/* skip even terms, which are 0 */
for (n = 5; n <= 19; n += 2) {
 /* chebpc */
   for (j = 1; j <= n; ++j) d[j] = e[j] = 0;
   d[1] = c[n];
   for (j = n - 1; j >= 2; -j) {
      for (k = n - j + 1; k >= 2; -k) {
         s = d[k];
         d[k] = d[k - 1] * 2 - e[k];
         e[k] = s; }
      s = d[1];
      d[1] = c[j] - e[1];
      e[1] = s; }
   for (j = n; j >= 2; -j) d[j] = d[j - 1] - e[j];
   d[1] = c[1] * .5 - e[1];
 /* pcshft */
   g = 2 / (b - a);
   for (j = 2; j <= n; ++j) {
      d[j] *= g;
      g *= 2 / (b - a); }
   for (j = 1; j < n; ++j) {
      h = d[n];
      for (k = n - 1; k >= j; -k) {
         h = d[k] - (a + b) * .5 * h;
         d[k] = h; }
   }
   "Chebyshev Sin fit |x|<Pi/2 coefficients"
      " Maximum Rel Error:"
      m(b(c[n + 2]), b(c[2])) / t(b);
   for (i = 1; i <= n; i += 2) d[i];
}
/* End of File */