summaryrefslogtreecommitdiff
path: root/libm/double/cheby.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/cheby.c')
-rw-r--r--libm/double/cheby.c149
1 files changed, 0 insertions, 149 deletions
diff --git a/libm/double/cheby.c b/libm/double/cheby.c
deleted file mode 100644
index 8da9b350e..000000000
--- a/libm/double/cheby.c
+++ /dev/null
@@ -1,149 +0,0 @@
-/* cheby.c
- *
- * Program to calculate coefficients of the Chebyshev polynomial
- * expansion of a given input function. The algorithm computes
- * the discrete Fourier cosine transform of the function evaluated
- * at unevenly spaced points. Library routine chbevl.c uses the
- * coefficients to calculate an approximate value of the original
- * function.
- * -- S. L. Moshier
- */
-
-extern double PI; /* 3.14159... */
-extern double PIO2;
-double cosi[33] = {0.0,}; /* cosine array for Fourier transform */
-double func[65] = {0.0,}; /* values of the function */
-double cos(), log(), exp(), sqrt();
-
-main()
-{
-double c, r, s, t, x, y, z, temp;
-double low, high, dtemp;
-long n;
-int i, ii, j, n2, k, rr, invflg;
-short *p;
-char st[40];
-
-low = 0.0; /* low end of approximation interval */
-high = 1.0; /* high end */
-invflg = 0; /* set to 1 if inverted interval, else zero */
-/* Note: inverted interval goes from 1/high to 1/low */
-z = 0.0;
-n = 64; /* will find 64 coefficients */
- /* but use only those greater than roundoff error */
-n2 = n/2;
-t = n;
-t = PI/t;
-
-/* calculate array of cosines */
-puts("calculating cosines");
-s = 1.0;
-cosi[0] = 1.0;
-i = 1;
-while( i < 32 )
- {
- y = cos( s * t );
- cosi[i] = y;
- s += 1.0;
- ++i;
- }
-cosi[32] = 0.0;
-
-/* cheby.c 2 */
-
-/* calculate function at special values of the argument */
-puts("calculating function values");
-x = low;
-y = high;
-if( invflg && (low != 0.0) )
- { /* inverted interval */
- temp = 1.0/x;
- x = 1.0/y;
- y = temp;
- }
-r = (x + y)/2.0;
-printf( "center %.15E ", r);
-s = (y - x)/2.0;
-printf( "width %.15E\n", s);
-i = 0;
-while( i < 65 )
- {
- if( i < n2 )
- c = cosi[i];
- else
- c = -cosi[64-i];
- temp = r + s * c;
-/* if inverted interval, compute function(1/x) */
- if( invflg && (temp != 0.0) )
- temp = 1.0/temp;
-
- printf( "%.15E ", temp );
-
-/* insert call to function routine here: */
-/**********************************/
-
- if( temp == 0.0 )
- y = 1.0;
- else
- y = exp( temp * log(2.0) );
-
-/**********************************/
- func[i] = y;
- printf( "%.15E\n", y );
- ++i;
- }
-
-/* cheby.c 3 */
-
-puts( "calculating Chebyshev coefficients");
-rr = 0;
-while( rr < 65 )
- {
- z = func[0]/2.0;
- j = 1;
- while( j < 65 )
- {
- k = (rr * j)/n2;
- i = rr * j - n2 * k;
- k &= 3;
- if( k == 0 )
- c = cosi[i];
- if( k == 1 )
- {
- i = 32-i;
- c = -cosi[i];
- if( i == 32 )
- c = -c;
- }
- if( k == 2 )
- {
- c = -cosi[i];
- }
- if( k == 3 )
- {
- i = 32-i;
- c = cosi[i];
- }
- if( i != 32)
- {
- temp = func[j];
- temp = c * temp;
- z += temp;
- }
- ++j;
- }
-
- if( i != 32 )
- {
- temp /= 2.0;
- z = z - temp;
- }
- z *= 2.0;
- temp = n;
- z /= temp;
- dtemp = z;
- ++rr;
- sprintf( st, "/* %.16E */", dtemp );
- puts( st );
- }
-}