summaryrefslogtreecommitdiff
path: root/libm/double/revers.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/revers.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/revers.c')
-rw-r--r--libm/double/revers.c156
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