Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : CHEBY.C

 
Output of file : CHEBY.C contained in archive : CEPHES22.ZIP
/* 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 */

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];
double cos(), log(), exp(), sqrt();

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 );
}
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : CHEBY.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/