97c428305d8f288f62284211a680051569654e2d
[ccan] / 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