/* 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