1 files changed, 149 insertions, 0 deletions
diff --git a/libm/double/cheby.c b/libm/double/cheby.c
new file mode 100644
index 000000000..8da9b350e
--- /dev/null
+++ b/libm/double/cheby.c
@@ -0,0 +1,149 @@
+/*	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 );
+	}
+}