From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/polmisc.c | 309 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 309 insertions(+) create mode 100644 libm/double/polmisc.c (limited to 'libm/double/polmisc.c') diff --git a/libm/double/polmisc.c b/libm/double/polmisc.c new file mode 100644 index 000000000..7d517ae69 --- /dev/null +++ b/libm/double/polmisc.c @@ -0,0 +1,309 @@ + +/* Square root, sine, cosine, and arctangent of polynomial. + * See polyn.c for data structures and discussion. + */ + +#include +#include +#ifdef ANSIPROT +extern double atan2 ( double, double ); +extern double sqrt ( double ); +extern double fabs ( double ); +extern double sin ( double ); +extern double cos ( double ); +extern void polclr ( double *a, int n ); +extern void polmov ( double *a, int na, double *b ); +extern void polmul ( double a[], int na, double b[], int nb, double c[] ); +extern void poladd ( double a[], int na, double b[], int nb, double c[] ); +extern void polsub ( double a[], int na, double b[], int nb, double c[] ); +extern int poldiv ( double a[], int na, double b[], int nb, double c[] ); +extern void polsbt ( double a[], int na, double b[], int nb, double c[] ); +extern void * malloc ( long ); +extern void free ( void * ); +#else +double atan2(), sqrt(), fabs(), sin(), cos(); +void polclr(), polmov(), polsbt(), poladd(), polsub(), polmul(); +int poldiv(); +void * malloc(); +void free (); +#endif + +/* Highest degree of polynomial to be handled + by the polyn.c subroutine package. */ +#define N 16 +/* Highest degree actually initialized at runtime. */ +extern int MAXPOL; + +/* Taylor series coefficients for various functions + */ +double patan[N+1] = { + 0.0, 1.0, 0.0, -1.0/3.0, 0.0, + 1.0/5.0, 0.0, -1.0/7.0, 0.0, 1.0/9.0, 0.0, -1.0/11.0, + 0.0, 1.0/13.0, 0.0, -1.0/15.0, 0.0 }; + +double psin[N+1] = { + 0.0, 1.0, 0.0, -1.0/6.0, 0.0, 1.0/120.0, 0.0, + -1.0/5040.0, 0.0, 1.0/362880.0, 0.0, -1.0/39916800.0, + 0.0, 1.0/6227020800.0, 0.0, -1.0/1.307674368e12, 0.0}; + +double pcos[N+1] = { + 1.0, 0.0, -1.0/2.0, 0.0, 1.0/24.0, 0.0, + -1.0/720.0, 0.0, 1.0/40320.0, 0.0, -1.0/3628800.0, 0.0, + 1.0/479001600.0, 0.0, -1.0/8.7179291e10, 0.0, 1.0/2.0922789888e13}; + +double pasin[N+1] = { + 0.0, 1.0, 0.0, 1.0/6.0, 0.0, + 3.0/40.0, 0.0, 15.0/336.0, 0.0, 105.0/3456.0, 0.0, 945.0/42240.0, + 0.0, 10395.0/599040.0 , 0.0, 135135.0/9676800.0 , 0.0 +}; + +/* Square root of 1 + x. */ +double psqrt[N+1] = { + 1.0, 1./2., -1./8., 1./16., -5./128., 7./256., -21./1024., 33./2048., + -429./32768., 715./65536., -2431./262144., 4199./524288., -29393./4194304., + 52003./8388608., -185725./33554432., 334305./67108864., + -9694845./2147483648.}; + +/* Arctangent of the ratio num/den of two polynomials. + */ +void +polatn( num, den, ans, nn ) + double num[], den[], ans[]; + int nn; +{ + double a, t; + double *polq, *polu, *polt; + int i; + + if (nn > N) + { + mtherr ("polatn", OVERFLOW); + return; + } + /* arctan( a + b ) = arctan(a) + arctan( b/(1 + ab + a**2) ) */ + t = num[0]; + a = den[0]; + if( (t == 0.0) && (a == 0.0 ) ) + { + t = num[1]; + a = den[1]; + } + t = atan2( t, a ); /* arctan(num/den), the ANSI argument order */ + polq = (double * )malloc( (MAXPOL+1) * sizeof (double) ); + polu = (double * )malloc( (MAXPOL+1) * sizeof (double) ); + polt = (double * )malloc( (MAXPOL+1) * sizeof (double) ); + polclr( polq, MAXPOL ); + i = poldiv( den, nn, num, nn, polq ); + a = polq[0]; /* a */ + polq[0] = 0.0; /* b */ + polmov( polq, nn, polu ); /* b */ + /* Form the polynomial + 1 + ab + a**2 + where a is a scalar. */ + for( i=0; i<=nn; i++ ) + polu[i] *= a; + polu[0] += 1.0 + a * a; + poldiv( polu, nn, polq, nn, polt ); /* divide into b */ + polsbt( polt, nn, patan, nn, polu ); /* arctan(b) */ + polu[0] += t; /* plus arctan(a) */ + polmov( polu, nn, ans ); + free( polt ); + free( polu ); + free( polq ); +} + + + +/* Square root of a polynomial. + * Assumes the lowest degree nonzero term is dominant + * and of even degree. An error message is given + * if the Newton iteration does not converge. + */ +void +polsqt( pol, ans, nn ) + double pol[], ans[]; + int nn; +{ + double t; + double *x, *y; + int i, n; +#if 0 + double z[N+1]; + double u; +#endif + + if (nn > N) + { + mtherr ("polatn", OVERFLOW); + return; + } + x = (double * )malloc( (MAXPOL+1) * sizeof (double) ); + y = (double * )malloc( (MAXPOL+1) * sizeof (double) ); + polmov( pol, nn, x ); + polclr( y, MAXPOL ); + + /* Find lowest degree nonzero term. */ + t = 0.0; + for( n=0; n 0 ) + { + if (n & 1) + { + printf("error, sqrt of odd polynomial\n"); + return; + } + /* Divide by x^n. */ + y[n] = x[n]; + poldiv (y, nn, pol, N, x); + } + + t = x[0]; + for( i=1; i<=nn; i++ ) + x[i] /= t; + x[0] = 0.0; + /* series development sqrt(1+x) = 1 + x / 2 - x**2 / 8 + x**3 / 16 + hopes that first (constant) term is greater than what follows */ + polsbt( x, nn, psqrt, nn, y); + t = sqrt( t ); + for( i=0; i<=nn; i++ ) + y[i] *= t; + + /* If first nonzero coefficient was at degree n > 0, multiply by + x^(n/2). */ + if (n > 0) + { + polclr (x, MAXPOL); + x[n/2] = 1.0; + polmul (x, nn, y, nn, y); + } +#if 0 +/* Newton iterations */ +for( n=0; n<10; n++ ) + { + poldiv( y, nn, pol, nn, z ); + poladd( y, nn, z, nn, y ); + for( i=0; i<=nn; i++ ) + y[i] *= 0.5; + for( i=0; i<=nn; i++ ) + { + u = fabs( y[i] - z[i] ); + if( u > 1.0e-15 ) + goto more; + } + goto done; +more: ; + } +printf( "square root did not converge\n" ); +done: +#endif /* 0 */ + +polmov( y, nn, ans ); +free( y ); +free( x ); +} + + + +/* Sine of a polynomial. + * The computation uses + * sin(a+b) = sin(a) cos(b) + cos(a) sin(b) + * where a is the constant term of the polynomial and + * b is the sum of the rest of the terms. + * Since sin(b) and cos(b) are computed by series expansions, + * the value of b should be small. + */ +void +polsin( x, y, nn ) + double x[], y[]; + int nn; +{ + double a, sc; + double *w, *c; + int i; + + if (nn > N) + { + mtherr ("polatn", OVERFLOW); + return; + } + w = (double * )malloc( (MAXPOL+1) * sizeof (double) ); + c = (double * )malloc( (MAXPOL+1) * sizeof (double) ); + polmov( x, nn, w ); + polclr( c, MAXPOL ); + polclr( y, nn ); + /* a, in the description, is x[0]. b is the polynomial x - x[0]. */ + a = w[0]; + /* c = cos (b) */ + w[0] = 0.0; + polsbt( w, nn, pcos, nn, c ); + sc = sin(a); + /* sin(a) cos (b) */ + for( i=0; i<=nn; i++ ) + c[i] *= sc; + /* y = sin (b) */ + polsbt( w, nn, psin, nn, y ); + sc = cos(a); + /* cos(a) sin(b) */ + for( i=0; i<=nn; i++ ) + y[i] *= sc; + poladd( c, nn, y, nn, y ); + free( c ); + free( w ); +} + + +/* Cosine of a polynomial. + * The computation uses + * cos(a+b) = cos(a) cos(b) - sin(a) sin(b) + * where a is the constant term of the polynomial and + * b is the sum of the rest of the terms. + * Since sin(b) and cos(b) are computed by series expansions, + * the value of b should be small. + */ +void +polcos( x, y, nn ) + double x[], y[]; + int nn; +{ + double a, sc; + double *w, *c; + int i; + double sin(), cos(); + + if (nn > N) + { + mtherr ("polatn", OVERFLOW); + return; + } + w = (double * )malloc( (MAXPOL+1) * sizeof (double) ); + c = (double * )malloc( (MAXPOL+1) * sizeof (double) ); + polmov( x, nn, w ); + polclr( c, MAXPOL ); + polclr( y, nn ); + a = w[0]; + w[0] = 0.0; + /* c = cos(b) */ + polsbt( w, nn, pcos, nn, c ); + sc = cos(a); + /* cos(a) cos(b) */ + for( i=0; i<=nn; i++ ) + c[i] *= sc; + /* y = sin(b) */ + polsbt( w, nn, psin, nn, y ); + sc = sin(a); + /* sin(a) sin(b) */ + for( i=0; i<=nn; i++ ) + y[i] *= sc; + polsub( y, nn, c, nn, y ); + free( c ); + free( w ); +} -- cgit v1.2.3