home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / polyn / polmisc.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  3.8 KB  |  196 lines
  1. /* square root, sine, and arctangent of polynomial
  2.  * See polyn.c for data structures and discussion.
  3.  */
  4.  
  5.  
  6. /* Highest degree of polynomial to be handled
  7.  * by the polyn.c subroutine package */
  8. #define N 16
  9.  
  10. /* Taylor series coefficients for various functions
  11.  */
  12. double patan[N+1] = {
  13. 0.0,     1.0,      0.0, -1.0/3.0,     0.0,
  14. 1.0/5.0, 0.0, -1.0/7.0,      0.0, 1.0/9.0, 0.0, -1.0/11.0,
  15. 0.0, 1.0/13.0, 0.0, -1.0/15.0, 0.0 };
  16.  
  17. double psin[N+1] = {
  18. 0.0, 1.0, 0.0,   -1.0/6.0,  0.0, 1.0/120.0,  0.0,
  19. -1.0/5040.0, 0.0, 1.0/362880.0, 0.0, -1.0/39916800.0,
  20. 0.0, 1.0/6227020800.0, 0.0, -1.0/1.307674368e12, 0.0};
  21.  
  22. double pcos[N+1] = {
  23. 1.0, 0.0,   -1.0/2.0,  0.0, 1.0/24.0,  0.0,
  24. -1.0/720.0, 0.0, 1.0/40320.0, 0.0, -1.0/3628800.0, 0.0,
  25. 1.0/479001600.0, 0.0, -1.0/8.7179291e10, 0.0, 1.0/2.0922789888e13};
  26.  
  27. double pasin[N+1] = {
  28. 0.0,     1.0,  0.0, 1.0/6.0,  0.0,
  29. 3.0/40.0, 0.0, 15.0/336.0, 0.0, 105.0/3456.0, 0.0, 945.0/42240.0,
  30. 0.0, 10395.0/599040.0 , 0.0, 135135.0/9676800.0 , 0.0
  31. };
  32.  
  33.  
  34. static double polt[N+1];
  35. static double polq[N+1];
  36. static double polu[N+1];
  37.  
  38.  
  39.  
  40. /* Arctangent of the ratio num/den of two polynomials.
  41.  */
  42. polatn( den, num, ans )
  43. double num[], den[], ans[];
  44. {
  45. double a, t;
  46. int i;
  47. double atan2();
  48.  
  49. /*
  50.  * arctan( a + b ) = arctan(a) + arctan( b/(1 + ab + a**2) )
  51.  */
  52. t = num[0];
  53. a = den[0];
  54. if( (t == 0.0) && (a == 0.0 ) )
  55.     {
  56.     t = num[1];
  57.     a = den[1];
  58.     }
  59. t = atan2( a, t ); /* arctan(a) */
  60. polclr( polq, N );
  61. i = poldiv( den, N, num, N, polq );
  62. a = polq[0]; /* a */
  63. polq[0] = 0.0; /* b */
  64. polmov( polq, N, polu ); /* b */
  65. /*
  66.  * Form the polynomial
  67.  *     1 + ab + a**2
  68.  * where a is a scalar.
  69.  */
  70. for( i=0; i<=N; i++ )
  71.     polu[i] *= a;
  72. polu[0] += 1.0 + a * a;
  73. poldiv( polu, N, polq, N, polt ); /* divide into b */
  74. polsbt( polt, N, patan, N, polu ); /* arctan(b)  */
  75. polu[0] += t; /* plus arctan(a) */
  76. }
  77.  
  78.  
  79.  
  80. /* Square root of a polynomial.
  81.  * Assumes the lowest degree nonzero term is dominant
  82.  * and of even degree.  An error message is given
  83.  * if the Newton iteration does not converge.
  84.  */
  85. polsqt( pol, ans )
  86. double pol[], ans[];
  87. {
  88. double x[N+1], y[N+1], z[N+1];
  89. double t, u;
  90. int i, j, n0;
  91. double sqrt(), fabs();
  92.  
  93. polmov( pol, N, x );
  94. polclr( y, N );
  95. /* Initial guess is square root of lowest order
  96.  * nonzero term
  97.  */
  98. t = 0.0;
  99. for( i=0; i<N; i++ )
  100.     {
  101.     if( x[i] != 0.0 )
  102.         goto nzero;
  103.     }
  104. polmov( y, N, ans );
  105. return;
  106.  
  107. nzero:
  108.  
  109. if( i > 0 )
  110.     {
  111.     t = x[i];
  112.     i /= 2;
  113.     y[i] = sqrt( t );
  114.     }
  115. else
  116.     {
  117.     t = x[0];
  118.     n0 = 0;
  119.     y[0] = 1.0;
  120.     for( i=1; i<=N; i++ )
  121.         x[i] /= t;
  122.     x[0] = 0.0;
  123. /* series development sqrt(1+x) = 1  +  x / 2  -  x**2 / 8  +  x**3 / 16
  124.  * assumes first (constant) term is greater than what follows
  125.  */
  126.     polmov( x, N, z );
  127.     for( i=0; i<=N; i++ )
  128.         z[i] *= 0.5;
  129.     poladd( z, N, y, N, y );
  130.     polmul( x, N, x, N, z );
  131.     for( i=0; i<=N; i++ )
  132.         z[i] *= 0.25;
  133.     polsub( z, N, y, N, y );
  134.     polmul( x, N, z, N, z );
  135.     for( i=0; i<=N; i++ )
  136.         z[i] *= 0.5;
  137.     poladd( z, N, y, N, y );
  138.     t = sqrt( t );
  139.     for( i=0; i<=N; i++ )
  140.         y[i] *= t;
  141.     }
  142. /* Newton iterations */
  143. for( j=0; j<10; j++ )
  144.     {
  145.     poldiv( y, N, pol, N, z );
  146.     poladd( y, N, z, N, y );
  147.     for( i=0; i<=N; i++ )
  148.         y[i] *= 0.5;
  149.     for( i=0; i<=N; i++ )
  150.         {
  151.         u = fabs( y[i] - z[i] );
  152.         if( u > 1.0e-15 )
  153.             goto more;
  154.         }
  155.     goto done;
  156. more:    ;
  157.     }
  158. printf( "square root did not converge\n" );
  159. done:
  160. polmov( y, N, ans );
  161. }
  162.  
  163.  
  164.  
  165. /* Sine of a polynomial.
  166.  * The computation uses
  167.  *     sin(a+b) = sin(a) cos(b) + cos(a) sin(b)
  168.  * where a is the constant term of the polynomial and
  169.  * b is the sum of the rest of the terms.
  170.  * Since sin(b) and cos(b) are computed by series expansions,
  171.  * the value of b should be small.
  172.  */
  173. polsin( x, y )
  174. double x[], y[];
  175. {
  176. double w[N+1], c[N+1];
  177. double a, sc;
  178. int i;
  179. double sin(), cos();
  180.  
  181. polmov( x, N, w );
  182. polclr( c, N );
  183. polclr( y, N );
  184. a = w[0];
  185. w[0] = 0.0;
  186. polsbt( w, N, pcos, N, c );
  187. sc = sin(a);
  188. for( i=0; i<=N; i++ )
  189.     c[i] *= sc;
  190. polsbt( w, N, psin, N, y );
  191. sc = cos(a);
  192. for( i=0; i<=N; i++ )
  193.     y[i] *= sc;
  194. poladd( c, N, y, N, y );
  195. }
  196.