← Back to branch summary

~ubuntu-branches/ubuntu/karmic/grace/karmic

« back to all changes in this revision

Viewing changes to cephes/revers.c

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-19 14:19:58 UTC
  • Revision ID: james.westby@ubuntu.com-20020319141958-5gxna6vo1ek3zjml
Tags: upstream-5.1.7
Import upstream version 5.1.7

Show diffs side-by-side

added added

removed removed

Lines of Context:
cephes/revers.c
 
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
 
 
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