home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / polyn / polyn.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  7.4 KB  |  458 lines
  1. /*                            polyn.c
  2.  *                            polyr.c
  3.  * Arithmetic operations on polynomials
  4.  *
  5.  * In the following descriptions a, b, c are polynomials of degree
  6.  * na, nb, nc respectively.  The degree of a polynomial cannot
  7.  * exceed a run-time value MAXPOL.  An operation that attempts
  8.  * to use or generate a polynomial of higher degree may produce a
  9.  * result that suffers truncation at degree MAXPOL.  The value of
  10.  * MAXPOL is set by calling the function
  11.  *
  12.  *     polini( maxpol );
  13.  *
  14.  * where maxpol is the desired maximum degree.  This must be
  15.  * done prior to calling any of the other functions in this module.
  16.  * Memory for internal temporary polynomial storage is allocated
  17.  * by polini().
  18.  *
  19.  * Each polynomial is represented by an array containing its
  20.  * coefficients, together with a separately declared integer equal
  21.  * to the degree of the polynomial.  The coefficients appear in
  22.  * ascending order; that is,
  23.  *
  24.  *                                        2                      na
  25.  * a(x)  =  a[0]  +  a[1] * x  +  a[2] * x   +  ...  +  a[na] * x  .
  26.  *
  27.  *
  28.  *
  29.  * sum = poleva( a, na, x );    Evaluate polynomial a(t) at t = x.
  30.  * polprt( a, na, D );        Print the coefficients of a to D digits.
  31.  * polclr( a, na );        Set a identically equal to zero, up to a[na].
  32.  * polmov( a, na, b );        Set b = a.
  33.  * poladd( a, na, b, nb, c );    c = b + a, nc = max(na,nb)
  34.  * polsub( a, na, b, nb, c );    c = b - a, nc = max(na,nb)
  35.  * polmul( a, na, b, nb, c );    c = b * a, nc = na+nb
  36.  *
  37.  *
  38.  * Division:
  39.  *
  40.  * i = poldiv( a, na, b, nb, c );    c = b / a, nc = MAXPOL
  41.  *
  42.  * returns i = the degree of the first nonzero coefficient of a.
  43.  * The computed quotient c must be divided by x^i.  An error message
  44.  * is printed if a is identically zero.
  45.  *
  46.  *
  47.  * Change of variables:
  48.  * If a and b are polynomials, and t = a(x), then
  49.  *     c(t) = b(a(x))
  50.  * is a polynomial found by substituting a(x) for t.  The
  51.  * subroutine call for this is
  52.  *
  53.  * polsbt( a, na, b, nb, c );
  54.  *
  55.  *
  56.  * Notes:
  57.  * poldiv() is an integer routine; poleva() is double.
  58.  * Any of the arguments a, b, c may refer to the same array.
  59.  *
  60.  */
  61.  
  62. #ifndef NULL
  63. #define NULL 0
  64. #endif
  65. #include "mconf.h"
  66.  
  67.  
  68. /* near pointer version of malloc() */
  69. #define malloc _nmalloc
  70. #define free _nfree
  71.  
  72. /* Pointers to internal arrays.  Note poldiv() allocates
  73.  * and deallocates some temporary arrays every time it is called.
  74.  */
  75. static double *pt1 = 0;
  76. static double *pt2 = 0;
  77. static double *pt3 = 0;
  78.  
  79. /* Maximum degree of polynomial. */
  80. int MAXPOL = 0;
  81. extern int MAXPOL;
  82.  
  83. /* Number of bytes (chars) in maximum size polynomial. */
  84. static int psize = 0;
  85.  
  86.  
  87. /* Initialize max degree of polynomials
  88.  * and allocate temporary storage.
  89.  */
  90. polini( maxdeg )
  91. int maxdeg;
  92. {
  93. void *malloc();
  94. int m;
  95.  
  96. MAXPOL = maxdeg;
  97. psize = (maxdeg + 1) * sizeof(double);
  98.  
  99. /* Release previously allocated memory, if any. */
  100. if( pt3 )
  101.     free(pt3);
  102. if( pt2 )
  103.     free(pt2);
  104. if( pt1 )
  105.     free(pt1);
  106.  
  107. /* Allocate new arrays */
  108. pt1 = (double * )malloc(psize); /* used by polsbt */
  109. pt2 = (double * )malloc(psize); /* used by polsbt */
  110. pt3 = (double * )malloc(psize); /* used by polmul */
  111.  
  112. /* Report if failure */
  113. if( (pt1 == NULL) || (pt2 == NULL) || (pt3 == NULL) )
  114.     {
  115.     mtherr( "polini", ERANGE );
  116.     exit(1);
  117.     }
  118. }
  119.  
  120.  
  121.  
  122. /* Print the coefficients of a, with d decimal precision.
  123.  */
  124. static char *form = "abcdefghijk";
  125.  
  126. polprt( a, na, d )
  127. double a[];
  128. int na, d;
  129. {
  130. int i, j, d1;
  131. char *p;
  132.  
  133. /* Create format descriptor string for the printout.
  134.  * Do this partly by hand, since sprintf() may be too
  135.  * bug-ridden to accomplish this feat by itself.
  136.  */
  137. p = form;
  138. *p++ = '%';
  139. d1 = d + 8;
  140. sprintf( p, "%d ", d1 );
  141. p += 1;
  142. if( d1 >= 10 )
  143.     p += 1;
  144. *p++ = '.';
  145. sprintf( p, "%d ", d );
  146. p += 1;
  147. if( d >= 10 )
  148.     p += 1;
  149. *p++ = 'e';
  150. *p++ = ' ';
  151. *p++ = '\0';
  152.  
  153.  
  154. /* Now do the printing.
  155.  */
  156. d1 += 1;
  157. j = 0;
  158. for( i=0; i<=na; i++ )
  159.     {
  160. /* Detect end of available line */
  161.     j += d1;
  162.     if( j >= 78 )
  163.         {
  164.         printf( "\n" );
  165.         j = d1;
  166.         }
  167.     printf( form, a[i] );
  168.     }
  169. printf( "\n" );
  170. }
  171.  
  172.  
  173.  
  174. /* Set a = 0.
  175.  */
  176. polclr( a, n )
  177. register double *a;
  178. int n;
  179. {
  180. int i;
  181.  
  182. if( n > MAXPOL )
  183.     n = MAXPOL;
  184. for( i=0; i<=n; i++ )
  185.     *a++ = 0.0;
  186. }
  187.  
  188.  
  189.  
  190. /* Set b = a.
  191.  */
  192. polmov( a, na, b )
  193. register double *a, *b;
  194. int na;
  195. {
  196. int i;
  197.  
  198. if( na > MAXPOL )
  199.     na = MAXPOL;
  200.  
  201. for( i=0; i<= na; i++ )
  202.     {
  203.     *b++ = *a++;
  204.     }
  205. }
  206.  
  207.  
  208. /* c = b * a.
  209.  */
  210. polmul( a, na, b, nb, c )
  211. double a[], b[], c[];
  212. int na, nb;
  213. {
  214. int i, j, k, nc;
  215. double x;
  216.  
  217. nc = na + nb;
  218. polclr( pt3, MAXPOL );
  219.  
  220. for( i=0; i<=na; i++ )
  221.     {
  222.     x = a[i];
  223.     for( j=0; j<=nb; j++ )
  224.         {
  225.         k = i + j;
  226.         if( k > MAXPOL )
  227.             break;
  228.         pt3[k] += x * b[j];
  229.         }
  230.     }
  231.  
  232. if( nc > MAXPOL )
  233.     nc = MAXPOL;
  234. for( i=0; i<=nc; i++ )
  235.     c[i] = pt3[i];
  236. }
  237.  
  238.  
  239.  
  240.  
  241. /* c = b + a.
  242.  */
  243. poladd( a, na, b, nb, c )
  244. double a[], b[], c[];
  245. int na, nb;
  246. {
  247. int i, n;
  248.  
  249.  
  250. if( na > nb )
  251.     n = na;
  252. else
  253.     n = nb;
  254.  
  255. if( n > MAXPOL )
  256.     n = MAXPOL;
  257.  
  258. for( i=0; i<=n; i++ )
  259.     {
  260.     if( i > na )
  261.         c[i] = b[i];
  262.     else if( i > nb )
  263.         c[i] = a[i];
  264.     else
  265.         c[i] = b[i] + a[i];
  266.     }
  267. }
  268.  
  269. /* c = b - a.
  270.  */
  271. polsub( a, na, b, nb, c )
  272. double a[], b[], c[];
  273. int na, nb;
  274. {
  275. int i, n;
  276.  
  277.  
  278. if( na > nb )
  279.     n = na;
  280. else
  281.     n = nb;
  282.  
  283. if( n > MAXPOL )
  284.     n = MAXPOL;
  285.  
  286. for( i=0; i<=n; i++ )
  287.     {
  288.     if( i > na )
  289.         c[i] = b[i];
  290.     else if( i > nb )
  291.         c[i] = -a[i];
  292.     else
  293.         c[i] = b[i] - a[i];
  294.     }
  295. }
  296.  
  297.  
  298.  
  299. /* c = b/a
  300.  */
  301. poldiv( a, na, b, nb, c )
  302. double a[], b[], c[];
  303. int na, nb;
  304. {
  305. double quot;
  306. double *ta, *tb, *tq;
  307. int i, j, k, sing;
  308.  
  309. sing = 0;
  310.  
  311. /* Allocate temporary arrays.  This would be quicker
  312.  * if done automatically on the stack, but stack space
  313.  * may be hard to obtain on a small computer.
  314.  */
  315. ta = (double * )malloc( psize );
  316. polclr( ta, MAXPOL );
  317. polmov( a, na, ta );
  318.  
  319. tb = (double * )malloc( psize );
  320. polclr( tb, MAXPOL );
  321. polmov( b, nb, tb );
  322.  
  323. tq = (double * )malloc( psize );
  324. polclr( tq, MAXPOL );
  325.  
  326. /* What to do if leading (constant) coefficient
  327.  * of denominator is zero.
  328.  */
  329. if( a[0] == 0.0 )
  330.     {
  331.     for( i=0; i<=na; i++ )
  332.         {
  333.         if( ta[i] != 0.0 )
  334.             goto nzero;
  335.         }
  336.     mtherr( "poldiv", SING );
  337.     goto done;
  338.  
  339. nzero:
  340. /* Reduce the degree of the denominator. */
  341.     for( i=0; i<na; i++ )
  342.         ta[i] = ta[i+1];
  343.     ta[na] = 0.0;
  344.  
  345.     if( b[0] != 0.0 )
  346.         {
  347. /* Optional message:
  348.         printf( "poldiv singularity, divide quotient by x\n" );
  349. */
  350.         sing += 1;
  351.         }
  352.     else
  353.         {
  354. /* Reduce degree of numerator. */
  355.         for( i=0; i<nb; i++ )
  356.             tb[i] = tb[i+1];
  357.         tb[nb] = 0.0;
  358.         }
  359. /* Call self, using reduced polynomials. */
  360.     sing += poldiv( ta, na, tb, nb, c );
  361.     goto done;
  362.     }
  363.  
  364. /* Long division algorithm.  ta[0] is nonzero.
  365.  */
  366. for( i=0; i<=MAXPOL; i++ )
  367.     {
  368.     quot = tb[i]/ta[0];
  369.     for( j=0; j<=MAXPOL; j++ )
  370.         {
  371.         k = j + i;
  372.         if( k > MAXPOL )
  373.             break;
  374.         tb[k] -= quot * ta[j];
  375.         }
  376.     tq[i] = quot;
  377.     }
  378. /* Send quotient to output array. */
  379. polmov( tq, MAXPOL, c );
  380.  
  381. done:
  382.  
  383. /* Restore allocated memory. */
  384. free(tq);
  385. free(tb);
  386. free(ta);
  387. return( sing );
  388. }
  389.  
  390.  
  391.  
  392.  
  393. /* Change of variables
  394.  * Substitute a(y) for the variable x in b(x).
  395.  * x = a(y)
  396.  * c(x) = b(x) = b(a(y)).
  397.  */
  398.  
  399. polsbt( a, na, b, nb, c )
  400. double a[], b[], c[];
  401. int na, nb;
  402. {
  403. int i, j, k, n2;
  404. double x;
  405.  
  406. /* 0th degree term:
  407.  */
  408. polclr( pt1, MAXPOL );
  409. pt1[0] = b[0];
  410.  
  411. polclr( pt2, MAXPOL );
  412. pt2[0] = 1.0;
  413. n2 = 0;
  414.  
  415. for( i=1; i<=nb; i++ )
  416.     {
  417. /* Form ith power of a. */
  418.     polmul( a, na, pt2, n2, pt2 );
  419.     n2 += na;
  420.     x = b[i];
  421. /* Add the ith coefficient of b times the ith power of a. */
  422.     for( j=0; j<=n2; j++ )
  423.         {
  424.         if( j > MAXPOL )
  425.             break;
  426.         pt1[j] += x * pt2[j];
  427.         }
  428.     }
  429.  
  430. k = n2 + nb;
  431. if( k > MAXPOL )
  432.     k = MAXPOL;
  433. for( i=0; i<=k; i++ )
  434.     c[i] = pt1[i];
  435. }
  436.  
  437.  
  438.  
  439.  
  440. /* Evaluate polynomial a(t) at t = x.
  441.  */
  442. double poleva( a, na, x )
  443. double a[];
  444. int na;
  445. double x;
  446. {
  447. double s;
  448. int i;
  449.  
  450. s = a[na];
  451. for( i=na-1; i>=0; i-- )
  452.     {
  453.     s = s * x + a[i];
  454.     }
  455. return(s);
  456. }
  457.  
  458.