git.ozlabs.org Git - ccan/blob - junkcode/iasoule32@gmail.com-polynomial/polynomial_adt.h

   1 /* Polynomial ADT \r
   2 ** A polynomial module with\r
   3 ** ability to add,sub,mul derivate/integrate, compose ... polynomials \r
   4 ** ..expansion in progress ...\r
   5  * Copyright (c) 2009 I. Soule\r
   6  * All rights reserved.\r
   7  *\r
   8  * Redistribution and use in source and binary forms, with or without\r
   9  * modification, are permitted provided that the following conditions\r
  10  * are met:\r
  11  * 1. Redistributions of source code must retain the above copyright\r
  12  *    notice, this list of conditions and the following disclaimer.\r
  13  * 2. Redistributions in binary form must reproduce the above copyright\r
  14  *    notice, this list of conditions and the following disclaimer in the\r
  15  *    documentation and/or other materials provided with the distribution.\r
  16  *\r
  17  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND\r
  18  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE\r
  19  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE\r
  20  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE\r
  21  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL\r
  22  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS\r
  23  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)\r
  24  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT\r
  25  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY\r
  26  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF\r
  27  * SUCH DAMAGE.\r
  28  *\r
  29 **          iasoule32@gmail.com\r
  30 */\r
  31 \r
  32 #ifndef __POLYNOMIAL_ADT\r
  33 #define __POLYNOMIAL_ADT\r
  34 \r
  35 #include <assert.h>\r
  36 #include <stdlib.h>\r
  37 #include <stdbool.h> //C99 compliance\r
  38 #include <math.h>\r
  39 \r
  40 #define max(a, b) (a) > (b) ? (a) : (b)\r
  41 #define sgn(a)    (a) < 0 ? '+' : '-' //for quadratic factored form\r
  42 \r
  43 typedef struct node {\r
  44     int exp;\r
  45     float coeff;\r
  46     struct node *next;\r
  47 }Node;\r
  48 \r
  49 typedef struct polynomial_adt {\r
  50     Node *head;\r
  51     int terms, hp; //hp highest power\r
  52 }PolyAdt;\r
  53 \r
  54 void display_poly(const PolyAdt *a);\r
  55 /**\r
  56 * create_adt - create a polynomial on the heap\r
  57 * @hp: the highest power in the polynomial\r
  58 */\r
  59 PolyAdt *create_adt(int hp) \r
  60 {\r
  61     PolyAdt *pAdt = malloc(sizeof(PolyAdt));\r
  62     assert(pAdt != NULL);\r
  63     \r
  64     pAdt->head = NULL;\r
  65     pAdt->terms = 0;\r
  66     pAdt->hp = hp;\r
  67 \r
  68     return pAdt;\r
  69 }\r
  70 /**\r
  71 * create_node - creates a Node (exponent, constant and next pointer) on the heap \r
  72 * @constant: the contant in the term \r
  73 * @exp:      the exponent on the term\r
  74 * @next:     the next pointer to another term in the polynomial\r
  75 * \r
  76 * This should not be called by client code (hence marked static)\r
  77 * used to assist insert_term()\r
  78 */ \r
  79 static Node *create_node(float constant, int exp, Node *next)\r
  80 {\r
  81     Node *nNode = malloc(sizeof(Node));\r
  82     assert(nNode != NULL);\r
  83     \r
  84     nNode->exp = exp;\r
  85     nNode->coeff = constant;\r
  86     nNode->next = next;\r
  87     return nNode;\r
  88 }\r
  89 \r
  90 /**\r
  91 * insert_term - inserts a term into the polynomial\r
  92 * @pAdt: the polynomial \r
  93 * @c:    constant value on the term\r
  94 * @e:    the exponent on the term \r
  95 */\r
  96 void insert_term(PolyAdt *pAdt, float c, int e)\r
  97 {\r
  98     assert(pAdt != NULL); //assume client code didnt call create_adt()\r
  99     Node *n = malloc(sizeof(Node));\r
 100        \r
 101     if(pAdt->head == NULL)\r
 102         pAdt->head = create_node(c, e, pAdt->head);\r
 103     else\r
 104         for(n = pAdt->head; n->next != NULL; n = n->next); //go to the end of list\r
 105             n->next = create_node(c, e, NULL);\r
 106     \r
 107     pAdt->terms++;\r
 108 }\r
 109 \r
 110 /**\r
 111 * polyImage - returns an image (direct) copy of the polynomial\r
 112 * @orig: the polynomial to be duplicated\r
 113 */\r
 114 PolyAdt *polyImage(const PolyAdt *orig)\r
 115 {\r
 116     PolyAdt *img = create_adt(orig->hp);\r
 117     Node *origHead = orig->head;\r
 118     \r
 119     for(; origHead; origHead = origHead->next)\r
 120              insert_term(img, origHead->coeff, origHead->exp);\r
 121     return img;\r
 122 }\r
 123 \r
 124 \r
 125 /**\r
 126 * add - adds two polynomials together, and returns their sum (as a polynomial)\r
 127 * @a: the 1st polynomial \r
 128 * @b: the 2nd polynomial\r
 129 */\r
 130 PolyAdt *add(const PolyAdt *a, const PolyAdt *b)\r
 131 {\r
 132     PolyAdt *sum;\r
 133     Node *n, *np;\r
 134     _Bool state = true;\r
 135     \r
 136     assert(a != NULL && b != NULL);\r
 137     \r
 138     int hpow = max(a->hp, b->hp);\r
 139     sum = create_adt(hpow); //create space for it\r
 140     \r
 141     /* using state machine to compare the poly with the most terms to \r
 142     ** the poly with fewer, round robin type of effect comparison of\r
 143     ** exponents => 3 Cases: Equal, Less, Greater\r
 144     */\r
 145         n = a->head; np = b->head;\r
 146         while(state) {\r
 147             /* compare the exponents */\r
 148             if(n->exp == np->exp){\r
 149                 insert_term(sum, n->coeff + np->coeff, n->exp);\r
 150                 n = n->next;\r
 151                 np = np->next;\r
 152             }\r
 153             \r
 154             else if(n->exp < np->exp){\r
 155                 insert_term(sum, np->coeff, np->exp);\r
 156                 np = np->next; //move to next term of b\r
 157             }\r
 158             \r
 159             else { //greater than\r
 160                 insert_term(sum, n->coeff, n->exp);\r
 161                 n = n->next;\r
 162             }\r
 163             /* check whether at the end of one list or the other */\r
 164             if(np == NULL && state == true){ //copy rest of a to sum\r
 165                 for(; n != NULL; n = n->next)\r
 166                     insert_term(sum, n->coeff, n->exp);\r
 167                 state = false;\r
 168             }\r
 169             \r
 170            if(n == NULL && state == true){\r
 171                 for(; np != NULL; np = np->next)\r
 172                     insert_term(sum, np->coeff, np->exp);\r
 173                 state = false;\r
 174             }       \r
 175      }        \r
 176     return sum;               \r
 177 }            \r
 178 \r
 179 /**\r
 180 * sub - subtracts two polynomials, and returns their difference (as a polynomial)\r
 181 * @a: the 1st polynomial \r
 182 * @b: the 2nd polynomial\r
 183 * Aids in code reuse by negating the terms (b) and then calls the add() function\r
 184 */\r
 185 PolyAdt *subtract(const PolyAdt *a, const PolyAdt *b)\r
 186 {\r
 187         assert(a != NULL && b != NULL);\r
 188 \r
 189     PolyAdt *tmp = create_adt(b->hp);\r
 190     Node *bptr;\r
 191     \r
 192     for(bptr = b->head; bptr != NULL; bptr = bptr->next)\r
 193         insert_term(tmp,-bptr->coeff,bptr->exp);  //negating b's coeffs\r
 194     return add(a,tmp);\r
 195 }\r
 196 \r
 197 /**\r
 198 * multiply - multiply two polynomials, and returns their product (as a polynomial)\r
 199 * @a: the 1st polynomial \r
 200 * @b: the 2nd polynomial\r
 201 */\r
 202 PolyAdt *multiply(const PolyAdt *a, const PolyAdt *b)\r
 203 {\r
 204         assert(a != NULL && b != NULL);\r
 205 \r
 206     //the polys are inserted in order for now\r
 207     PolyAdt *prod = create_adt(a->head->exp + b->head->exp);\r
 208     Node *n = a->head, *np = b->head;\r
 209     Node *t = b->head; \r
 210     \r
 211     if(a->terms < b->terms){\r
 212         n = b->head;\r
 213         np = t = a->head;\r
 214     }\r
 215     \r
 216     for(; n != NULL; n = n->next){\r
 217         np = t; //reset to the beginning\r
 218         for(; np != NULL; np = np->next){ //always the least term in this loop\r
 219                 insert_term(prod, n->coeff * np->coeff, n->exp + np->exp);\r
 220         }\r
 221     }\r
 222 \r
 223     return prod;       \r
 224 }\r
 225 \r
 226 /**\r
 227 * derivative - computes the derivative of a polynomial and returns the result\r
 228 * @a: the polynomial to take the derivative upon\r
 229 */\r
 230 PolyAdt *derivative(const PolyAdt *a)\r
 231 {\r
 232         assert(a != NULL);\r
 233         \r
 234         PolyAdt *deriv = create_adt(a->head->exp - 1);\r
 235         Node *n = a->head;\r
 236 \r
 237         for(; n != NULL; n = n->next){\r
 238                 if(n->exp == 0) break;\r
 239                 insert_term(deriv, n->coeff * n->exp, n->exp-1);\r
 240         }\r
 241         return deriv;\r
 242 }\r
 243 /**\r
 244 * integrate - computes the integral of a polynomial and returns the result\r
 245 * @a: the polynomial to take the integral of\r
 246 * \r
 247 * Will compute an indefinite integral over a\r
 248 */\r
 249 PolyAdt *integrate(const PolyAdt *a)\r
 250 {\r
 251         assert(a != NULL);\r
 252         \r
 253         PolyAdt *integrand = create_adt(a->head->exp + 1);\r
 254         Node *n;\r
 255 \r
 256         for(n = a->head; n != NULL; n = n->next) //very simple term by term\r
 257         insert_term(integrand, (float)n->coeff/(n->exp+1.0F), n->exp + 1);\r
 258     \r
 259         return integrand;\r
 260 }\r
 261 /**\r
 262 * quadratic_roots - finds the roots of the polynomial ax^2+bx+c, a != 0 && b != 0\r
 263 * @a: the polynomial\r
 264 * @real: a pointer to float of the real(R) part of a\r
 265 * @cplx: a pointer to float of the imaginary(I) part of a\r
 266 *\r
 267 * Usage:\r
 268 * Two options can be done by the client \r
 269 * 1. Either pass NULL to real and cplx\r
 270 *    this will display the roots by printf\r
 271 *    quadratic_roots(myPolynomial, NULL, NULL);\r
 272 *\r
 273 * 2. Pass in pointers** to type float of the real and complex\r
 274 *    if the discriminant is >0 cplx = -ve root of X\r
 275 *    quadratic_roots(myPolynomial, &realPart, &complexPart);\r
 276 */\r
 277 void quadratic_roots(const PolyAdt *a, float *real, float *cplx)\r
 278 {\r
 279         assert(a != NULL);\r
 280         \r
 281         float dscrmnt, _a, b, c;\r
 282         float u, v;\r
 283     \r
 284     Node *n = a->head;\r
 285     _a = n->coeff; b = n->next->coeff; c = n->next->next->coeff;\r
 286     \r
 287         dscrmnt = (b*b) - 4*_a*c;\r
 288     u = -b/(2*_a); v = sqrt((double)fabs(dscrmnt))/(2*_a);\r
 289     \r
 290         if(real && !cplx || !real && cplx)\r
 291                 assert(true);\r
 292 \r
 293         if(real == NULL && cplx == NULL){\r
 294                 if(a->hp != 2 && a->terms < 3){\r
 295                   printf("Invalid Quadratic*, A and B must be non-zero");\r
 296                         return;\r
 297         }\r
 298         \r
 299                 if(dscrmnt != 0)\r
 300                         printf("X = %.2f +/- %.2f%c\n",u,v,dscrmnt < 0 ? 'I':' ');\r
 301                 else{\r
 302                         printf("(X %c %.2f)(X %c %.2f)\n",sgn(u),fabs(u),sgn(u),fabs(u));\r
 303                         printf("X1,2 = %.2f\n",u);\r
 304                 }\r
 305         }\r
 306         //save values in pointers\r
 307         else {\r
 308                 if(dscrmnt < 0){ //x = u +/- vI Re(x) = u, Im(x) = +v\r
 309                         *real = u; \r
 310                         *cplx = v; //understand +/- is not representable\r
 311                 }\r
 312                 else if(dscrmnt == 0){\r
 313                         *real = u; \r
 314                         *cplx = 0.00F;\r
 315                 }\r
 316                 else{\r
 317                         *real = u + v;\r
 318                         *cplx = u - v;\r
 319                 }\r
 320         }\r
 321 }\r
 322 \r
 323 /**\r
 324 * exponentiate - computes polynomial exponentiation (P(x))^n, n E Z*\r
 325 * @a: the polynomial\r
 326 * @n: the exponent\r
 327 * Works fast for small n (n < 8) currently runs ~ O(n^2 lg n)\r
 328 */\r
 329 PolyAdt *exponentiate(const PolyAdt *a, int n)\r
 330 {\r
 331         assert(a != NULL);\r
 332 \r
 333         PolyAdt *expn = create_adt(a->hp *  n);\r
 334         PolyAdt *aptr = polyImage(a);\r
 335     int hl = n / 2;\r
 336     \r
 337     //check default cases before calculation\r
 338     if(n == 0){\r
 339         insert_term(expn, 1, 0);\r
 340         return expn;\r
 341     }\r
 342     else if(n == 1){\r
 343         return aptr;\r
 344     }\r
 345         \r
 346         for(; hl ; hl--)\r
 347         aptr = multiply(aptr, aptr);\r
 348 \r
 349     if(n % 2) //odd exponent do a^(n-1) * a = a^n\r
 350         expn = multiply(aptr, a);\r
 351     else\r
 352         expn = aptr;\r
 353     return expn;\r
 354 }\r
 355 /**\r
 356 * compose - computes the composition of two polynomials P(Q(x)) and returns the composition\r
 357 * @p: polynomial P(x) which will x will be equal to Q(x)\r
 358 * @q: polynomial Q(x) which is the argument to P(x)\r
 359 */\r
 360 PolyAdt *compose(const PolyAdt *p, const PolyAdt *q)\r
 361 {\r
 362     assert(p && q);\r
 363     \r
 364         PolyAdt *comp = create_adt(p->head->exp * q->head->exp);\r
 365         PolyAdt *exp;\r
 366         \r
 367     Node *pp = p->head;\r
 368     Node *qq = q->head;\r
 369     \r
 370     int swap = 0;\r
 371     \r
 372     if(p->terms < q->terms){\r
 373         pp = q->head;\r
 374         qq = p->head;\r
 375         swap = 1;\r
 376     }\r
 377     \r
 378     /* going through, exponentiate each term with the exponent of p */\r
 379         for(; pp != NULL; pp = pp->next){\r
 380                 exp = exponentiate(swap ? p: q, pp->exp);\r
 381                 insert_term(comp, pp->coeff * exp->head->coeff, exp->head->exp);\r
 382         }\r
 383     \r
 384     return comp;\r
 385 }\r
 386 /** \r
 387 * destroy_poly - completely frees the polynomial from the heap and resets all values\r
 388 * @poly: the polynomial to release memory back to the heap\r
 389 * Usage:\r
 390 * destroy_poly(myPoly); //puts polynomial on free list\r
 391 */\r
 392 void destroy_poly(PolyAdt *poly)\r
 393 {\r
 394     Node *ps = poly->head;\r
 395     Node *tmp = NULL;\r
 396     while(ps != NULL){\r
 397         tmp = ps;\r
 398         free(tmp);\r
 399         ps = ps->next;\r
 400     }\r
 401     poly->hp = poly->terms = 0;\r
 402     poly->head = NULL;\r
 403 }\r
 404 /**\r
 405 * display_poly - displays the polynomial to the console in nice format\r
 406 * @a: the polynomial to display \r
 407 */\r
 408 void display_poly(const PolyAdt *a)\r
 409 {\r
 410     assert(a != NULL);\r
 411     Node *n;\r
 412     \r
 413     for(n = a->head; n != NULL; n = n->next){\r
 414         \r
 415        n->coeff < 0 ? putchar('-') : putchar('+'); \r
 416         if(n->exp == 0)\r
 417             printf(" %.2f ",fabs(n->coeff));\r
 418         else if(n->coeff == 1)\r
 419             printf(" X^%d ",n->exp);\r
 420         else if(n->exp == 1)\r
 421             printf(" %.2fX ",fabs(n->coeff));\r
 422         else if(n->coeff == 0)\r
 423             continue;\r
 424         else\r
 425             printf(" %.2fX^%d ",fabs(n->coeff),n->exp);\r
 426         }\r
 427     printf("\n\n");\r
 428 }\r
 429 \r
 430 #endif\r