diff options
Diffstat (limited to 'libm/double/revers.c')
-rw-r--r-- | libm/double/revers.c | 156 |
1 files changed, 0 insertions, 156 deletions
diff --git a/libm/double/revers.c b/libm/double/revers.c deleted file mode 100644 index 370bdb5d6..000000000 --- a/libm/double/revers.c +++ /dev/null @@ -1,156 +0,0 @@ -/* revers.c - * - * Reversion of power series - * - * - * - * SYNOPSIS: - * - * extern int MAXPOL; - * int n; - * double x[n+1], y[n+1]; - * - * polini(n); - * revers( y, x, n ); - * - * Note, polini() initializes the polynomial arithmetic subroutines; - * see polyn.c. - * - * - * DESCRIPTION: - * - * If - * - * inf - * - i - * y(x) = > a x - * - i - * i=1 - * - * then - * - * inf - * - j - * x(y) = > A y , - * - j - * j=1 - * - * where - * 1 - * A = --- - * 1 a - * 1 - * - * etc. The coefficients of x(y) are found by expanding - * - * inf inf - * - - i - * x(y) = > A > a x - * - j - i - * j=1 i=1 - * - * and setting each coefficient of x , higher than the first, - * to zero. - * - * - * - * RESTRICTIONS: - * - * y[0] must be zero, and y[1] must be nonzero. - * - */ - -/* -Cephes Math Library Release 2.8: June, 2000 -Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier -*/ - -#include <math.h> - -extern int MAXPOL; /* initialized by polini() */ - -#ifdef ANSIPROT -/* See polyn.c. */ -void polmov ( double *, int, double * ); -void polclr ( double *, int ); -void poladd ( double *, int, double *, int, double * ); -void polmul ( double *, int, double *, int, double * ); -void * malloc ( long ); -void free ( void * ); -#else -void polmov(), polclr(), poladd(), polmul(); -void * malloc(); -void free (); -#endif - -void revers( y, x, n) -double y[], x[]; -int n; -{ -double *yn, *yp, *ysum; -int j; - -if( y[1] == 0.0 ) - mtherr( "revers", DOMAIN ); -/* printf( "revers: y[1] = 0\n" );*/ -j = (MAXPOL + 1) * sizeof(double); -yn = (double *)malloc(j); -yp = (double *)malloc(j); -ysum = (double *)malloc(j); - -polmov( y, n, yn ); -polclr( ysum, n ); -x[0] = 0.0; -x[1] = 1.0/y[1]; -for( j=2; j<=n; j++ ) - { -/* A_(j-1) times the expansion of y^(j-1) */ - polmul( &x[j-1], 0, yn, n, yp ); -/* The expansion of the sum of A_k y^k up to k=j-1 */ - poladd( yp, n, ysum, n, ysum ); -/* The expansion of y^j */ - polmul( yn, n, y, n, yn ); -/* The coefficient A_j to make the sum up to k=j equal to zero */ - x[j] = -ysum[j]/yn[j]; - } -free(yn); -free(yp); -free(ysum); -} - - -#if 0 -/* Demonstration program - */ -#define N 10 -double y[N], x[N]; -double fac(); - -main() -{ -double a, odd; -int i; - -polini( N-1 ); -a = 1.0; -y[0] = 0.0; -odd = 1.0; -for( i=1; i<N; i++ ) - { -/* sin(x) */ -/* - if( i & 1 ) - { - y[i] = odd/fac(i); - odd = -odd; - } - else - y[i] = 0.0; -*/ - y[i] = 1.0/fac(i); - } -revers( y, x, N-1 ); -for( i=0; i<N; i++ ) - printf( "%2d %.10e %.10e\n", i, x[i], y[i] ); -} -#endif |