summaryrefslogtreecommitdiff
path: root/libm/double/polmisc.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/polmisc.c')
-rw-r--r--libm/double/polmisc.c309
1 files changed, 0 insertions, 309 deletions
diff --git a/libm/double/polmisc.c b/libm/double/polmisc.c
deleted file mode 100644
index 7d517ae69..000000000
--- a/libm/double/polmisc.c
+++ /dev/null
@@ -1,309 +0,0 @@
-
-/* Square root, sine, cosine, and arctangent of polynomial.
- * See polyn.c for data structures and discussion.
- */
-
-#include <stdio.h>
-#include <math.h>
-#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<nn; n++ )
- {
- if( x[n] != 0.0 )
- goto nzero;
- }
- polmov( y, nn, ans );
- return;
-
-nzero:
-
- if( 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 );
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-24 04:45:11 +0000