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

   1 #include "polynomial_adt.h"
   2
   3 PolyAdt *create_adt(int hp)
   4 {
   5     PolyAdt *pAdt = malloc(sizeof(PolyAdt));
   6     assert(pAdt != NULL);
   7
   8     pAdt->head = NULL;
   9     pAdt->terms = 0;
  10     pAdt->hp = hp;
  11
  12     return pAdt;
  13 }
  14
  15 void insert_term(PolyAdt *pAdt, float c, int e)
  16 {
  17     assert(pAdt != NULL); //assume client code didn't call create_adt()
  18     Node *n = malloc(sizeof(Node));
  19
  20     if(pAdt->head == NULL)
  21         pAdt->head = create_node(c, e, pAdt->head);
  22     else
  23         for(n = pAdt->head; n->next != NULL; n = n->next); //go to the end of list
  24             n->next = create_node(c, e, NULL);
  25
  26     pAdt->terms++;
  27 }
  28
  29 PolyAdt *polyImage(const PolyAdt *orig)
  30 {
  31     PolyAdt *img = create_adt(orig->hp);
  32     Node *origHead = orig->head;
  33
  34     for(; origHead; origHead = origHead->next)
  35              insert_term(img, origHead->coeff, origHead->exp);
  36     return img;
  37 }
  38
  39 PolyAdt *add(const PolyAdt *a, const PolyAdt *b)
  40 {
  41     PolyAdt *sum;
  42     Node *n, *np;
  43     int state = 1;
  44
  45     assert(a != NULL && b != NULL);
  46
  47     int hpow = max(a->hp, b->hp);
  48     sum = create_adt(hpow); //create space for it
  49
  50     /* using state machine to compare the poly with the most terms to
  51     ** the poly with fewer, round robin type of effect comparison of
  52     ** exponents => 3 Cases: Equal, Less, Greater
  53     */
  54         n = a->head; np = b->head;
  55         while(state) {
  56             /* compare the exponents */
  57             if(n->exp == np->exp){
  58                 insert_term(sum, n->coeff + np->coeff, n->exp);
  59                 n = n->next;
  60                 np = np->next;
  61             }
  62
  63             else if(n->exp < np->exp){
  64                 insert_term(sum, np->coeff, np->exp);
  65                 np = np->next; //move to next term of b
  66             }
  67
  68             else { //greater than
  69                 insert_term(sum, n->coeff, n->exp);
  70                 n = n->next;
  71             }
  72             /* check whether at the end of one list or the other */
  73             if(np == NULL && state){ //copy rest of a to sum
  74                 for(; n != NULL; n = n->next)
  75                     insert_term(sum, n->coeff, n->exp);
  76                 state = 0;
  77             }
  78
  79            if(n == NULL && state){
  80                 for(; np != NULL; np = np->next)
  81                     insert_term(sum, np->coeff, np->exp);
  82                 state = 0;
  83             }
  84      }
  85     return sum;
  86 }
  87
  88 PolyAdt *subtract(const PolyAdt *a, const PolyAdt *b)
  89 {
  90         assert(a != NULL && b != NULL);
  91
  92     PolyAdt *tmp = create_adt(b->hp);
  93     Node *bptr;
  94
  95     for(bptr = b->head; bptr != NULL; bptr = bptr->next)
  96         insert_term(tmp,-bptr->coeff,bptr->exp);  //negating b's coeffs
  97     return add(a,tmp);
  98 }
  99
 100 PolyAdt *multiply(const PolyAdt *a, const PolyAdt *b)
 101 {
 102         assert(a != NULL && b != NULL);
 103
 104     //the polys are inserted in order for now
 105     PolyAdt *prod = create_adt(a->head->exp + b->head->exp);
 106     Node *n = a->head, *np = b->head;
 107     Node *t = b->head;
 108
 109     if(a->terms < b->terms){
 110         n = b->head;
 111         np = t = a->head;
 112     }
 113
 114     for(; n != NULL; n = n->next){
 115         np = t; //reset to the beginning
 116         for(; np != NULL; np = np->next){ //always the least term in this loop
 117                 insert_term(prod, n->coeff * np->coeff, n->exp + np->exp);
 118         }
 119     }
 120
 121     return prod;
 122 }
 123
 124 PolyAdt *derivative(const PolyAdt *a)
 125 {
 126         assert(a != NULL);
 127
 128         PolyAdt *deriv = create_adt(a->head->exp - 1);
 129         Node *n = a->head;
 130
 131         for(; n != NULL; n = n->next){
 132                 if(n->exp == 0) break;
 133                 insert_term(deriv, n->coeff * n->exp, n->exp-1);
 134         }
 135         return deriv;
 136 }
 137
 138 PolyAdt *integrate(const PolyAdt *a)
 139 {
 140         assert(a != NULL);
 141
 142         PolyAdt *integrand = create_adt(a->head->exp + 1);
 143         Node *n;
 144
 145         for(n = a->head; n != NULL; n = n->next) //very simple term by term
 146         insert_term(integrand, (float)n->coeff/(n->exp+1.0F), n->exp + 1);
 147
 148         return integrand;
 149 }
 150
 151 void quadratic_roots(const PolyAdt *a, float *real, float *cplx)
 152 {
 153         assert(a != NULL);
 154
 155         float dscrmnt, _a, b, c;
 156         float u, v;
 157
 158     Node *n = a->head;
 159     _a = n->coeff; b = n->next->coeff; c = n->next->next->coeff;
 160
 161         dscrmnt = (b*b) - 4*_a*c;
 162     u = -b/(2*_a); v = sqrt((double)fabs(dscrmnt))/(2*_a);
 163
 164         if((real && !cplx) || (!real && cplx))
 165                 assert(1);
 166
 167         if(real == NULL && cplx == NULL){
 168                 if(a->hp != 2 && a->terms < 3){
 169                   printf("Invalid Quadratic*, A and B must be non-zero");
 170                         return;
 171         }
 172
 173                 if(dscrmnt != 0)
 174                         printf("X = %.2f +/- %.2f%c\n",u,v,dscrmnt < 0 ? 'I':' ');
 175                 else{
 176                         printf("(X %c %.2f)(X %c %.2f)\n",sgn(u),fabs(u),sgn(u),fabs(u));
 177                         printf("X1,2 = %.2f\n",u);
 178                 }
 179         }
 180         //save values in pointers
 181         else {
 182                 if(dscrmnt < 0){ //x = u +/- vI Re(x) = u, Im(x) = +v
 183                         *real = u;
 184                         *cplx = v; //understand +/- is not representable
 185                 }
 186                 else if(dscrmnt == 0){
 187                         *real = u;
 188                         *cplx = 0.00F;
 189                 }
 190                 else{
 191                         *real = u + v;
 192                         *cplx = u - v;
 193                 }
 194         }
 195 }
 196
 197 PolyAdt *exponentiate(const PolyAdt *a, int n)
 198 {
 199         assert(a != NULL);
 200
 201         PolyAdt *expn = create_adt(a->hp *  n);
 202         PolyAdt *aptr = polyImage(a);
 203     int hl = n / 2;
 204
 205     //check default cases before calculation
 206     if(n == 0){
 207         insert_term(expn, 1, 0);
 208         return expn;
 209     }
 210     else if(n == 1){
 211         return aptr;
 212     }
 213
 214         for(; hl ; hl--)
 215         aptr = multiply(aptr, aptr);
 216
 217     if(n % 2) //odd exponent do a^(n-1) * a = a^n
 218         expn = multiply(aptr, a);
 219     else
 220         expn = aptr;
 221     return expn;
 222 }
 223
 224 PolyAdt *compose(const PolyAdt *p, const PolyAdt *q)
 225 {
 226     assert(p && q);
 227
 228         PolyAdt *comp = create_adt(p->head->exp * q->head->exp);
 229         PolyAdt *exp;
 230
 231     Node *pp = p->head;
 232     Node *qq = q->head;
 233
 234     int swap = 0;
 235
 236     if(p->terms < q->terms){
 237         pp = q->head;
 238         qq = p->head;
 239         swap = 1;
 240     }
 241
 242     /* going through, exponentiate each term with the exponent of p */
 243         for(; pp != NULL; pp = pp->next){
 244                 exp = exponentiate(swap ? p: q, pp->exp);
 245                 insert_term(comp, pp->coeff * exp->head->coeff, exp->head->exp);
 246         }
 247
 248     return comp;
 249 }
 250
 251 void destroy_poly(PolyAdt *poly)
 252 {
 253     Node *ps = poly->head;
 254     Node *tmp = NULL;
 255     while(ps != NULL){
 256         tmp = ps;
 257         free(tmp);
 258         ps = ps->next;
 259     }
 260     poly->hp = poly->terms = 0;
 261     poly->head = NULL;
 262 }
 263
 264 void display_poly(const PolyAdt *a)
 265 {
 266     assert(a != NULL);
 267     Node *n;
 268
 269     for(n = a->head; n != NULL; n = n->next){
 270
 271        n->coeff < 0 ? putchar('-') : putchar('+');
 272         if(n->exp == 0)
 273             printf(" %.2f ",fabs(n->coeff));
 274         else if(n->coeff == 1)
 275             printf(" X^%d ",n->exp);
 276         else if(n->exp == 1)
 277             printf(" %.2fX ",fabs(n->coeff));
 278         else if(n->coeff == 0)
 279             continue;
 280         else
 281             printf(" %.2fX^%d ",fabs(n->coeff),n->exp);
 282         }
 283     printf("\n\n");
 284 }