1 files changed, 156 insertions, 0 deletions
diff --git a/libm/double/revers.c b/libm/double/revers.c
new file mode 100644
index 000000000..370bdb5d6
--- /dev/null
+++ b/libm/double/revers.c
@@ -0,0 +1,156 @@
+/*							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