lib/libm/polynf.c

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