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
8 * Redistribution and use in source and binary forms, with or without
\r
9 * modification, are permitted provided that the following conditions
\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
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
29 ** iasoule32@gmail.com
\r
32 #ifndef __POLYNOMIAL_ADT
\r
33 #define __POLYNOMIAL_ADT
\r
37 #include <stdbool.h> //C99 compliance
\r
40 #define max(a, b) (a) > (b) ? (a) : (b)
\r
41 #define sgn(a) (a) < 0 ? '+' : '-' //for quadratic factored form
\r
43 typedef struct node {
\r
49 typedef struct polynomial_adt {
\r
51 int terms, hp; //hp highest power
\r
54 void display_poly(const PolyAdt *a);
\r
56 * create_adt - create a polynomial on the heap
\r
57 * @hp: the highest power in the polynomial
\r
59 PolyAdt *create_adt(int hp)
\r
61 PolyAdt *pAdt = malloc(sizeof(PolyAdt));
\r
62 assert(pAdt != NULL);
\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
76 * This should not be called by client code (hence marked static)
\r
77 * used to assist insert_term()
\r
79 static Node *create_node(float constant, int exp, Node *next)
\r
81 Node *nNode = malloc(sizeof(Node));
\r
82 assert(nNode != NULL);
\r
85 nNode->coeff = constant;
\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
96 void insert_term(PolyAdt *pAdt, float c, int e)
\r
98 assert(pAdt != NULL); //assume client code didnt call create_adt()
\r
99 Node *n = malloc(sizeof(Node));
\r
101 if(pAdt->head == NULL)
\r
102 pAdt->head = create_node(c, e, pAdt->head);
\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
111 * polyImage - returns an image (direct) copy of the polynomial
\r
112 * @orig: the polynomial to be duplicated
\r
114 PolyAdt *polyImage(const PolyAdt *orig)
\r
116 PolyAdt *img = create_adt(orig->hp);
\r
117 Node *origHead = orig->head;
\r
119 for(; origHead; origHead = origHead->next)
\r
120 insert_term(img, origHead->coeff, origHead->exp);
\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
130 PolyAdt *add(const PolyAdt *a, const PolyAdt *b)
\r
134 _Bool state = true;
\r
136 assert(a != NULL && b != NULL);
\r
138 int hpow = max(a->hp, b->hp);
\r
139 sum = create_adt(hpow); //create space for it
\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
145 n = a->head; np = b->head;
\r
147 /* compare the exponents */
\r
148 if(n->exp == np->exp){
\r
149 insert_term(sum, n->coeff + np->coeff, n->exp);
\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
159 else { //greater than
\r
160 insert_term(sum, n->coeff, n->exp);
\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
170 if(n == NULL && state == true){
\r
171 for(; np != NULL; np = np->next)
\r
172 insert_term(sum, np->coeff, np->exp);
\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
185 PolyAdt *subtract(const PolyAdt *a, const PolyAdt *b)
\r
187 assert(a != NULL && b != NULL);
\r
189 PolyAdt *tmp = create_adt(b->hp);
\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
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
202 PolyAdt *multiply(const PolyAdt *a, const PolyAdt *b)
\r
204 assert(a != NULL && b != NULL);
\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
211 if(a->terms < b->terms){
\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
227 * derivative - computes the derivative of a polynomial and returns the result
\r
228 * @a: the polynomial to take the derivative upon
\r
230 PolyAdt *derivative(const PolyAdt *a)
\r
234 PolyAdt *deriv = create_adt(a->head->exp - 1);
\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
244 * integrate - computes the integral of a polynomial and returns the result
\r
245 * @a: the polynomial to take the integral of
\r
247 * Will compute an indefinite integral over a
\r
249 PolyAdt *integrate(const PolyAdt *a)
\r
253 PolyAdt *integrand = create_adt(a->head->exp + 1);
\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
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
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
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
277 void quadratic_roots(const PolyAdt *a, float *real, float *cplx)
\r
281 float dscrmnt, _a, b, c;
\r
285 _a = n->coeff; b = n->next->coeff; c = n->next->next->coeff;
\r
287 dscrmnt = (b*b) - 4*_a*c;
\r
288 u = -b/(2*_a); v = sqrt((double)fabs(dscrmnt))/(2*_a);
\r
290 if(real && !cplx || !real && cplx)
\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
300 printf("X = %.2f +/- %.2f%c\n",u,v,dscrmnt < 0 ? 'I':' ');
\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
306 //save values in pointers
\r
308 if(dscrmnt < 0){ //x = u +/- vI Re(x) = u, Im(x) = +v
\r
310 *cplx = v; //understand +/- is not representable
\r
312 else if(dscrmnt == 0){
\r
324 * exponentiate - computes polynomial exponentiation (P(x))^n, n E Z*
\r
325 * @a: the polynomial
\r
327 * Works fast for small n (n < 8) currently runs ~ O(n^2 lg n)
\r
329 PolyAdt *exponentiate(const PolyAdt *a, int n)
\r
333 PolyAdt *expn = create_adt(a->hp * n);
\r
334 PolyAdt *aptr = polyImage(a);
\r
337 //check default cases before calculation
\r
339 insert_term(expn, 1, 0);
\r
347 aptr = multiply(aptr, aptr);
\r
349 if(n % 2) //odd exponent do a^(n-1) * a = a^n
\r
350 expn = multiply(aptr, a);
\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
360 PolyAdt *compose(const PolyAdt *p, const PolyAdt *q)
\r
364 PolyAdt *comp = create_adt(p->head->exp * q->head->exp);
\r
367 Node *pp = p->head;
\r
368 Node *qq = q->head;
\r
372 if(p->terms < q->terms){
\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
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
390 * destroy_poly(myPoly); //puts polynomial on free list
\r
392 void destroy_poly(PolyAdt *poly)
\r
394 Node *ps = poly->head;
\r
401 poly->hp = poly->terms = 0;
\r
405 * display_poly - displays the polynomial to the console in nice format
\r
406 * @a: the polynomial to display
\r
408 void display_poly(const PolyAdt *a)
\r
413 for(n = a->head; n != NULL; n = n->next){
\r
415 n->coeff < 0 ? putchar('-') : putchar('+');
\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
425 printf(" %.2fX^%d ",fabs(n->coeff),n->exp);
\r