home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / polyn / revers.c < prev   
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.2 KB  |  142 lines
  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. extern int MAXPOL; /* initialized by polini() */
  70.  
  71. revers( y, x, n)
  72. double y[], x[];
  73. int n;
  74. {
  75. double s, a;
  76. double *yn, *yp, *ysum;
  77. int i, j;
  78.  
  79. if( y[1] == 0.0 )
  80.     printf( "revers: y[1] = 0\n" );
  81. j = (MAXPOL + 1) * sizeof(double);
  82. yn = (double *)malloc(j);
  83. yp = (double *)malloc(j);
  84. ysum = (double *)malloc(j);
  85.  
  86. polmov( y, n, yn );
  87. polclr( ysum, n );
  88. x[0] = 0.0;
  89. x[1] = 1.0/y[1];
  90. for( j=2; j<=n; j++ )
  91.     {
  92. /* A_(j-1) times the expansion of y^(j-1)  */
  93.     polmul( &x[j-1], 0, yn, n, yp );
  94. /* The expansion of the sum of A_k y^k up to k=j-1 */
  95.     poladd( yp, n, ysum, n, ysum );
  96. /* The expansion of y^j */
  97.     polmul( yn, n, y, n, yn );
  98. /* The coefficient A_j to make the sum up to k=j equal to zero */
  99.     x[j] = -ysum[j]/yn[j];
  100.     }
  101. free(yn);
  102. free(yp);
  103. free(ysum);
  104. }
  105.  
  106.  
  107. #if 0
  108. /* Demonstration program
  109.  */
  110. #define N 10
  111. double y[N], x[N];
  112. double fac();
  113.  
  114. main()
  115. {
  116. double a, odd;
  117. int i;
  118.  
  119. polini( N-1 );
  120. a = 1.0;
  121. y[0] = 0.0;
  122. odd = 1.0;
  123. for( i=1; i<N; i++ )
  124.     {
  125. /* sin(x) */
  126. /*
  127.     if( i & 1 )
  128.         {
  129.         y[i] = odd/fac(i);
  130.         odd = -odd;
  131.         }
  132.     else
  133.         y[i] = 0.0;
  134. */
  135.     y[i] = 1.0/fac(i);
  136.     }
  137. revers( y, x, N-1 );
  138. for( i=0; i<N; i++ )
  139.     printf( "%2d %.10e %.10e\n", i, x[i], y[i] );
  140. }
  141. #endif
  142.