Listing 1 bc program to calculate coefficients for approximations up to 28 digits precision

/* bc program to calculate Chebyshef economized polynomial
** for evaluation of sin(x) */
/* use bc -l < listing1.bc to get s() function
** replace t() with c() to fit cos(x)
** a(x)/x for atan() */
define t(x) {                           /* sin(x)/x */
       if(x==0)return(1.); /* l'Hopital's rule */
   return (s(x) / x); /* put function to be fit here */ }
define b(x) { /* abs(x) */
   if (x < 0) return (-x);
   return (x); }
define m(x, y) { /* max(x,y) */
   if (x > y) return (x);
   return (y); }
n = 22; /* number of Chebyshef terms - need at least 3
           more than ** are used below */
scale = 40;
/* set precision for dc -must be more than desired result
*/
p = a(1.) * 4; /* pi */
b = p * .25;
/* upper end of curve fit interval; b = 7/16 for atan()
*/
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];
    /* old style bc requires =+ */
    (c[j] = 2 / n * s); } /* display results */
/* 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; 1 == 1; --k) {
          s = e[k];
          e[k] = d[k];
          if(k==1)break;
          d[k] = d[k - 1] * 2 - s; }
       d[1] = c[j] - 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;
 /*  bc doesn't support multiple assignment */
       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/4 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 */