Delete moved file.
[ccan] / 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 didnt 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 }