lib/libm/polyr.c

   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  * `a', `b', `c' are arrays of fracts.
  29  * poleva( a, na, &x, &sum );   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 #include <stdio.h>
  63 #include "mconf.h"
  64 #ifndef NULL
  65 #define NULL 0
  66 #endif
  67 typedef struct{
  68         double n;
  69         double d;
  70         }fract;
  71
  72 #ifdef ANSIPROT
  73 extern void radd ( fract *, fract *, fract * );
  74 extern void rsub ( fract *, fract *, fract * );
  75 extern void rmul ( fract *, fract *, fract * );
  76 extern void rdiv ( fract *, fract *, fract * );
  77 void polmov ( fract *, int, fract * );
  78 void polmul ( fract *, int, fract *, int, fract * );
  79 int poldiv ( fract *, int, fract *, int, fract * );
  80 void * malloc ( long );
  81 void free ( void * );
  82 #else
  83 void radd(), rsub(), rmul(), rdiv();
  84 void polmov(), polmul();
  85 int poldiv();
  86 void * malloc();
  87 void free ();
  88 #endif
  89
  90 /* near pointer version of malloc() */
  91 /*
  92 #define malloc _nmalloc
  93 #define free _nfree
  94 */
  95 /* Pointers to internal arrays.  Note poldiv() allocates
  96  * and deallocates some temporary arrays every time it is called.
  97  */
  98 static fract *pt1 = 0;
  99 static fract *pt2 = 0;
 100 static fract *pt3 = 0;
 101
 102 /* Maximum degree of polynomial. */
 103 int MAXPOL = 0;
 104 extern int MAXPOL;
 105
 106 /* Number of bytes (chars) in maximum size polynomial. */
 107 static int psize = 0;
 108
 109
 110 /* Initialize max degree of polynomials
 111  * and allocate temporary storage.
 112  */
 113 void polini( maxdeg )
 114 int maxdeg;
 115 {
 116
 117 MAXPOL = maxdeg;
 118 psize = (maxdeg + 1) * sizeof(fract);
 119
 120 /* Release previously allocated memory, if any. */
 121 if( pt3 )
 122         free(pt3);
 123 if( pt2 )
 124         free(pt2);
 125 if( pt1 )
 126         free(pt1);
 127
 128 /* Allocate new arrays */
 129 pt1 = (fract * )malloc(psize); /* used by polsbt */
 130 pt2 = (fract * )malloc(psize); /* used by polsbt */
 131 pt3 = (fract * )malloc(psize); /* used by polmul */
 132
 133 /* Report if failure */
 134 if( (pt1 == NULL) || (pt2 == NULL) || (pt3 == NULL) )
 135         {
 136         mtherr( "polini", ERANGE );
 137         exit(1);
 138         }
 139 }
 140
 141
 142
 143 /* Print the coefficients of a, with d decimal precision.
 144  */
 145 static char *form = "abcdefghijk";
 146
 147 void polprt( a, na, d )
 148 fract a[];
 149 int na, d;
 150 {
 151 int i, j, d1;
 152 char *p;
 153
 154 /* Create format descriptor string for the printout.
 155  * Do this partly by hand, since sprintf() may be too
 156  * bug-ridden to accomplish this feat by itself.
 157  */
 158 p = form;
 159 *p++ = '%';
 160 d1 = d + 8;
 161 sprintf( p, "%d ", d1 );
 162 p += 1;
 163 if( d1 >= 10 )
 164         p += 1;
 165 *p++ = '.';
 166 sprintf( p, "%d ", d );
 167 p += 1;
 168 if( d >= 10 )
 169         p += 1;
 170 *p++ = 'e';
 171 *p++ = ' ';
 172 *p++ = '\0';
 173
 174
 175 /* Now do the printing.
 176  */
 177 d1 += 1;
 178 j = 0;
 179 for( i=0; i<=na; i++ )
 180         {
 181 /* Detect end of available line */
 182         j += d1;
 183         if( j >= 78 )
 184                 {
 185                 printf( "\n" );
 186                 j = d1;
 187                 }
 188         printf( form, a[i].n );
 189         j += d1;
 190         if( j >= 78 )
 191                 {
 192                 printf( "\n" );
 193                 j = d1;
 194                 }
 195         printf( form, a[i].d );
 196         }
 197 printf( "\n" );
 198 }
 199
 200
 201
 202 /* Set a = 0.
 203  */
 204 void polclr( a, n )
 205 fract a[];
 206 int n;
 207 {
 208 int i;
 209
 210 if( n > MAXPOL )
 211         n = MAXPOL;
 212 for( i=0; i<=n; i++ )
 213         {
 214         a[i].n = 0.0;
 215         a[i].d = 1.0;
 216         }
 217 }
 218
 219
 220
 221 /* Set b = a.
 222  */
 223 void polmov( a, na, b )
 224 fract a[], b[];
 225 int na;
 226 {
 227 int i;
 228
 229 if( na > MAXPOL )
 230         na = MAXPOL;
 231
 232 for( i=0; i<= na; i++ )
 233         {
 234         b[i].n = a[i].n;
 235         b[i].d = a[i].d;
 236         }
 237 }
 238
 239
 240 /* c = b * a.
 241  */
 242 void polmul( a, na, b, nb, c )
 243 fract a[], b[], c[];
 244 int na, nb;
 245 {
 246 int i, j, k, nc;
 247 fract temp;
 248 fract *p;
 249
 250 nc = na + nb;
 251 polclr( pt3, MAXPOL );
 252
 253 p = &a[0];
 254 for( i=0; i<=na; i++ )
 255         {
 256         for( j=0; j<=nb; j++ )
 257                 {
 258                 k = i + j;
 259                 if( k > MAXPOL )
 260                         break;
 261                 rmul( p, &b[j], &temp ); /*pt3[k] += a[i] * b[j];*/
 262                 radd( &temp, &pt3[k], &pt3[k] );
 263                 }
 264         ++p;
 265         }
 266
 267 if( nc > MAXPOL )
 268         nc = MAXPOL;
 269 for( i=0; i<=nc; i++ )
 270         {
 271         c[i].n = pt3[i].n;
 272         c[i].d = pt3[i].d;
 273         }
 274 }
 275
 276
 277
 278
 279 /* c = b + a.
 280  */
 281 void poladd( a, na, b, nb, c )
 282 fract a[], b[], c[];
 283 int na, nb;
 284 {
 285 int i, n;
 286
 287
 288 if( na > nb )
 289         n = na;
 290 else
 291         n = nb;
 292
 293 if( n > MAXPOL )
 294         n = MAXPOL;
 295
 296 for( i=0; i<=n; i++ )
 297         {
 298         if( i > na )
 299                 {
 300                 c[i].n = b[i].n;
 301                 c[i].d = b[i].d;
 302                 }
 303         else if( i > nb )
 304                 {
 305                 c[i].n = a[i].n;
 306                 c[i].d = a[i].d;
 307                 }
 308         else
 309                 {
 310                 radd( &a[i], &b[i], &c[i] ); /*c[i] = b[i] + a[i];*/
 311                 }
 312         }
 313 }
 314
 315 /* c = b - a.
 316  */
 317 void polsub( a, na, b, nb, c )
 318 fract a[], b[], c[];
 319 int na, nb;
 320 {
 321 int i, n;
 322
 323
 324 if( na > nb )
 325         n = na;
 326 else
 327         n = nb;
 328
 329 if( n > MAXPOL )
 330         n = MAXPOL;
 331
 332 for( i=0; i<=n; i++ )
 333         {
 334         if( i > na )
 335                 {
 336                 c[i].n = b[i].n;
 337                 c[i].d = b[i].d;
 338                 }
 339         else if( i > nb )
 340                 {
 341                 c[i].n = -a[i].n;
 342                 c[i].d = a[i].d;
 343                 }
 344         else
 345                 {
 346                 rsub( &a[i], &b[i], &c[i] ); /*c[i] = b[i] - a[i];*/
 347                 }
 348         }
 349 }
 350
 351
 352
 353 /* c = b/a
 354  */
 355 int poldiv( a, na, b, nb, c )
 356 fract a[], b[], c[];
 357 int na, nb;
 358 {
 359 fract *ta, *tb, *tq;
 360 fract quot;
 361 fract temp;
 362 int i, j, k, sing;
 363
 364 sing = 0;
 365
 366 /* Allocate temporary arrays.  This would be quicker
 367  * if done automatically on the stack, but stack space
 368  * may be hard to obtain on a small computer.
 369  */
 370 ta = (fract * )malloc( psize );
 371 polclr( ta, MAXPOL );
 372 polmov( a, na, ta );
 373
 374 tb = (fract * )malloc( psize );
 375 polclr( tb, MAXPOL );
 376 polmov( b, nb, tb );
 377
 378 tq = (fract * )malloc( psize );
 379 polclr( tq, MAXPOL );
 380
 381 /* What to do if leading (constant) coefficient
 382  * of denominator is zero.
 383  */
 384 if( a[0].n == 0.0 )
 385         {
 386         for( i=0; i<=na; i++ )
 387                 {
 388                 if( ta[i].n != 0.0 )
 389                         goto nzero;
 390                 }
 391         mtherr( "poldiv", SING );
 392         goto done;
 393
 394 nzero:
 395 /* Reduce the degree of the denominator. */
 396         for( i=0; i<na; i++ )
 397                 {
 398                 ta[i].n = ta[i+1].n;
 399                 ta[i].d = ta[i+1].d;
 400                 }
 401         ta[na].n = 0.0;
 402         ta[na].d = 1.0;
 403
 404         if( b[0].n != 0.0 )
 405                 {
 406 /* Optional message:
 407                 printf( "poldiv singularity, divide quotient by x\n" );
 408 */
 409                 sing += 1;
 410                 }
 411         else
 412                 {
 413 /* Reduce degree of numerator. */
 414                 for( i=0; i<nb; i++ )
 415                         {
 416                         tb[i].n = tb[i+1].n;
 417                         tb[i].d = tb[i+1].d;
 418                         }
 419                 tb[nb].n = 0.0;
 420                 tb[nb].d = 1.0;
 421                 }
 422 /* Call self, using reduced polynomials. */
 423         sing += poldiv( ta, na, tb, nb, c );
 424         goto done;
 425         }
 426
 427 /* Long division algorithm.  ta[0] is nonzero.
 428  */
 429 for( i=0; i<=MAXPOL; i++ )
 430         {
 431         rdiv( &ta[0], &tb[i], &quot ); /*quot = tb[i]/ta[0];*/
 432         for( j=0; j<=MAXPOL; j++ )
 433                 {
 434                 k = j + i;
 435                 if( k > MAXPOL )
 436                         break;
 437
 438                 rmul( &ta[j], &quot, &temp ); /*tb[k] -= quot * ta[j];*/
 439                 rsub( &temp, &tb[k], &tb[k] );
 440                 }
 441         tq[i].n = quot.n;
 442         tq[i].d = quot.d;
 443         }
 444 /* Send quotient to output array. */
 445 polmov( tq, MAXPOL, c );
 446
 447 done:
 448
 449 /* Restore allocated memory. */
 450 free(tq);
 451 free(tb);
 452 free(ta);
 453 return( sing );
 454 }
 455
 456
 457
 458
 459 /* Change of variables
 460  * Substitute a(y) for the variable x in b(x).
 461  * x = a(y)
 462  * c(x) = b(x) = b(a(y)).
 463  */
 464
 465 void polsbt( a, na, b, nb, c )
 466 fract a[], b[], c[];
 467 int na, nb;
 468 {
 469 int i, j, k, n2;
 470 fract temp;
 471 fract *p;
 472
 473 /* 0th degree term:
 474  */
 475 polclr( pt1, MAXPOL );
 476 pt1[0].n = b[0].n;
 477 pt1[0].d = b[0].d;
 478
 479 polclr( pt2, MAXPOL );
 480 pt2[0].n = 1.0;
 481 pt2[0].d = 1.0;
 482 n2 = 0;
 483 p = &b[1];
 484
 485 for( i=1; i<=nb; i++ )
 486         {
 487 /* Form ith power of a. */
 488         polmul( a, na, pt2, n2, pt2 );
 489         n2 += na;
 490 /* Add the ith coefficient of b times the ith power of a. */
 491         for( j=0; j<=n2; j++ )
 492                 {
 493                 if( j > MAXPOL )
 494                         break;
 495                 rmul( &pt2[j], p, &temp ); /*pt1[j] += b[i] * pt2[j];*/
 496                 radd( &temp, &pt1[j], &pt1[j] );
 497                 }
 498         ++p;
 499         }
 500
 501 k = n2 + nb;
 502 if( k > MAXPOL )
 503         k = MAXPOL;
 504 for( i=0; i<=k; i++ )
 505         {
 506         c[i].n = pt1[i].n;
 507         c[i].d = pt1[i].d;
 508         }
 509 }
 510
 511
 512
 513
 514 /* Evaluate polynomial a(t) at t = x.
 515  */
 516 void poleva( a, na, x, s )
 517 fract a[];
 518 int na;
 519 fract *x;
 520 fract *s;
 521 {
 522 int i;
 523 fract temp;
 524
 525 s->n = a[na].n;
 526 s->d = a[na].d;
 527 for( i=na-1; i>=0; i-- )
 528         {
 529         rmul( s, x, &temp ); /*s = s * x + a[i];*/
 530         radd( &a[i], &temp, s );
 531         }
 532 }
 533