summaryrefslogtreecommitdiff
path: root/libm/double/revers.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/revers.c')
-rw-r--r--libm/double/revers.c156
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