1 /* revers.c
2 *
3 * Reversion of power series
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * extern int MAXPOL;
10 * int n;
11 * double x[n+1], y[n+1];
12 *
13 * polini(n);
14 * revers( y, x, n );
15 *
16 * Note, polini() initializes the polynomial arithmetic subroutines;
17 * see polyn.c.
18 *
19 *
20 * DESCRIPTION:
21 *
22 * If
23 *
24 * inf
25 * - i
26 * y(x) = > a x
27 * - i
28 * i=1
29 *
30 * then
31 *
32 * inf
33 * - j
34 * x(y) = > A y ,
35 * - j
36 * j=1
37 *
38 * where
39 * 1
40 * A = ---
41 * 1 a
42 * 1
43 *
44 * etc. The coefficients of x(y) are found by expanding
45 *
46 * inf inf
47 * - - i
48 * x(y) = > A > a x
49 * - j - i
50 * j=1 i=1
51 *
52 * and setting each coefficient of x , higher than the first,
53 * to zero.
54 *
55 *
56 *
57 * RESTRICTIONS:
58 *
59 * y[0] must be zero, and y[1] must be nonzero.
60 *
61 */
62
63 /*
64 Cephes Math Library Release 2.2: July, 1992
65 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
66 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
67 */
68
69 #include "mconf.h"
70 #include "cephes.h"
71
72 #include <stdlib.h>
73
74 extern int MAXPOL; /* initialized by polini() */
75
76 /* See polyn.c. */
77 void polmov(), polclr(), poladd(), polmul();
78
revers(y,x,n)79 void revers( y, x, n)
80 double y[], x[];
81 int n;
82 {
83 double *yn, *yp, *ysum;
84 int j;
85
86 if( y[1] == 0.0 )
87 mtherr( "revers", DOMAIN );
88 /* printf( "revers: y[1] = 0\n" );*/
89 j = (MAXPOL + 1) * sizeof(double);
90 yn = (double *)malloc(j);
91 yp = (double *)malloc(j);
92 ysum = (double *)malloc(j);
93
94 polmov( y, n, yn );
95 polclr( ysum, n );
96 x[0] = 0.0;
97 x[1] = 1.0/y[1];
98 for( j=2; j<=n; j++ )
99 {
100 /* A_(j-1) times the expansion of y^(j-1) */
101 polmul( &x[j-1], 0, yn, n, yp );
102 /* The expansion of the sum of A_k y^k up to k=j-1 */
103 poladd( yp, n, ysum, n, ysum );
104 /* The expansion of y^j */
105 polmul( yn, n, y, n, yn );
106 /* The coefficient A_j to make the sum up to k=j equal to zero */
107 x[j] = -ysum[j]/yn[j];
108 }
109 free(yn);
110 free(yp);
111 free(ysum);
112 }
113
114