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