From 031e1ee091a1f62756e8e56a019564bf3385308f Mon Sep 17 00:00:00 2001 From: Rusty Russell Date: Sat, 4 Apr 2009 11:59:54 +1030 Subject: [PATCH] New polynomial module (not quite CCAN format). --- .../polynomial_adt.h | 430 ++++++++++++++++++ .../polynomial_adt_test.c | 132 ++++++ 2 files changed, 562 insertions(+) create mode 100644 junkcode/iasoule32@gmail.com-polynomial/polynomial_adt.h create mode 100644 junkcode/iasoule32@gmail.com-polynomial/polynomial_adt_test.c diff --git a/junkcode/iasoule32@gmail.com-polynomial/polynomial_adt.h b/junkcode/iasoule32@gmail.com-polynomial/polynomial_adt.h new file mode 100644 index 00000000..97c42830 --- /dev/null +++ b/junkcode/iasoule32@gmail.com-polynomial/polynomial_adt.h @@ -0,0 +1,430 @@ +/* Polynomial ADT +** A polynomial module with +** ability to add,sub,mul derivate/integrate, compose ... polynomials +** ..expansion in progress ... + * Copyright (c) 2009 I. Soule + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * +** iasoule32@gmail.com +*/ + +#ifndef __POLYNOMIAL_ADT +#define __POLYNOMIAL_ADT + +#include +#include +#include //C99 compliance +#include + +#define max(a, b) (a) > (b) ? (a) : (b) +#define sgn(a) (a) < 0 ? '+' : '-' //for quadratic factored form + +typedef struct node { + int exp; + float coeff; + struct node *next; +}Node; + +typedef struct polynomial_adt { + Node *head; + int terms, hp; //hp highest power +}PolyAdt; + +void display_poly(const PolyAdt *a); +/** +* create_adt - create a polynomial on the heap +* @hp: the highest power in the polynomial +*/ +PolyAdt *create_adt(int hp) +{ + PolyAdt *pAdt = malloc(sizeof(PolyAdt)); + assert(pAdt != NULL); + + pAdt->head = NULL; + pAdt->terms = 0; + pAdt->hp = hp; + + return pAdt; +} +/** +* create_node - creates a Node (exponent, constant and next pointer) on the heap +* @constant: the contant in the term +* @exp: the exponent on the term +* @next: the next pointer to another term in the polynomial +* +* This should not be called by client code (hence marked static) +* used to assist insert_term() +*/ +static Node *create_node(float constant, int exp, Node *next) +{ + Node *nNode = malloc(sizeof(Node)); + assert(nNode != NULL); + + nNode->exp = exp; + nNode->coeff = constant; + nNode->next = next; + return nNode; +} + +/** +* insert_term - inserts a term into the polynomial +* @pAdt: the polynomial +* @c: constant value on the term +* @e: the exponent on the term +*/ +void insert_term(PolyAdt *pAdt, float c, int e) +{ + assert(pAdt != NULL); //assume client code didnt call create_adt() + Node *n = malloc(sizeof(Node)); + + if(pAdt->head == NULL) + pAdt->head = create_node(c, e, pAdt->head); + else + for(n = pAdt->head; n->next != NULL; n = n->next); //go to the end of list + n->next = create_node(c, e, NULL); + + pAdt->terms++; +} + +/** +* polyImage - returns an image (direct) copy of the polynomial +* @orig: the polynomial to be duplicated +*/ +PolyAdt *polyImage(const PolyAdt *orig) +{ + PolyAdt *img = create_adt(orig->hp); + Node *origHead = orig->head; + + for(; origHead; origHead = origHead->next) + insert_term(img, origHead->coeff, origHead->exp); + return img; +} + + +/** +* add - adds two polynomials together, and returns their sum (as a polynomial) +* @a: the 1st polynomial +* @b: the 2nd polynomial +*/ +PolyAdt *add(const PolyAdt *a, const PolyAdt *b) +{ + PolyAdt *sum; + Node *n, *np; + _Bool state = true; + + assert(a != NULL && b != NULL); + + int hpow = max(a->hp, b->hp); + sum = create_adt(hpow); //create space for it + + /* using state machine to compare the poly with the most terms to + ** the poly with fewer, round robin type of effect comparison of + ** exponents => 3 Cases: Equal, Less, Greater + */ + n = a->head; np = b->head; + while(state) { + /* compare the exponents */ + if(n->exp == np->exp){ + insert_term(sum, n->coeff + np->coeff, n->exp); + n = n->next; + np = np->next; + } + + else if(n->exp < np->exp){ + insert_term(sum, np->coeff, np->exp); + np = np->next; //move to next term of b + } + + else { //greater than + insert_term(sum, n->coeff, n->exp); + n = n->next; + } + /* check whether at the end of one list or the other */ + if(np == NULL && state == true){ //copy rest of a to sum + for(; n != NULL; n = n->next) + insert_term(sum, n->coeff, n->exp); + state = false; + } + + if(n == NULL && state == true){ + for(; np != NULL; np = np->next) + insert_term(sum, np->coeff, np->exp); + state = false; + } + } + return sum; +} + +/** +* sub - subtracts two polynomials, and returns their difference (as a polynomial) +* @a: the 1st polynomial +* @b: the 2nd polynomial +* Aids in code reuse by negating the terms (b) and then calls the add() function +*/ +PolyAdt *subtract(const PolyAdt *a, const PolyAdt *b) +{ + assert(a != NULL && b != NULL); + + PolyAdt *tmp = create_adt(b->hp); + Node *bptr; + + for(bptr = b->head; bptr != NULL; bptr = bptr->next) + insert_term(tmp,-bptr->coeff,bptr->exp); //negating b's coeffs + return add(a,tmp); +} + +/** +* multiply - multiply two polynomials, and returns their product (as a polynomial) +* @a: the 1st polynomial +* @b: the 2nd polynomial +*/ +PolyAdt *multiply(const PolyAdt *a, const PolyAdt *b) +{ + assert(a != NULL && b != NULL); + + //the polys are inserted in order for now + PolyAdt *prod = create_adt(a->head->exp + b->head->exp); + Node *n = a->head, *np = b->head; + Node *t = b->head; + + if(a->terms < b->terms){ + n = b->head; + np = t = a->head; + } + + for(; n != NULL; n = n->next){ + np = t; //reset to the beginning + for(; np != NULL; np = np->next){ //always the least term in this loop + insert_term(prod, n->coeff * np->coeff, n->exp + np->exp); + } + } + + return prod; +} + +/** +* derivative - computes the derivative of a polynomial and returns the result +* @a: the polynomial to take the derivative upon +*/ +PolyAdt *derivative(const PolyAdt *a) +{ + assert(a != NULL); + + PolyAdt *deriv = create_adt(a->head->exp - 1); + Node *n = a->head; + + for(; n != NULL; n = n->next){ + if(n->exp == 0) break; + insert_term(deriv, n->coeff * n->exp, n->exp-1); + } + return deriv; +} +/** +* integrate - computes the integral of a polynomial and returns the result +* @a: the polynomial to take the integral of +* +* Will compute an indefinite integral over a +*/ +PolyAdt *integrate(const PolyAdt *a) +{ + assert(a != NULL); + + PolyAdt *integrand = create_adt(a->head->exp + 1); + Node *n; + + for(n = a->head; n != NULL; n = n->next) //very simple term by term + insert_term(integrand, (float)n->coeff/(n->exp+1.0F), n->exp + 1); + + return integrand; +} +/** +* quadratic_roots - finds the roots of the polynomial ax^2+bx+c, a != 0 && b != 0 +* @a: the polynomial +* @real: a pointer to float of the real(R) part of a +* @cplx: a pointer to float of the imaginary(I) part of a +* +* Usage: +* Two options can be done by the client +* 1. Either pass NULL to real and cplx +* this will display the roots by printf +* quadratic_roots(myPolynomial, NULL, NULL); +* +* 2. Pass in pointers** to type float of the real and complex +* if the discriminant is >0 cplx = -ve root of X +* quadratic_roots(myPolynomial, &realPart, &complexPart); +*/ +void quadratic_roots(const PolyAdt *a, float *real, float *cplx) +{ + assert(a != NULL); + + float dscrmnt, _a, b, c; + float u, v; + + Node *n = a->head; + _a = n->coeff; b = n->next->coeff; c = n->next->next->coeff; + + dscrmnt = (b*b) - 4*_a*c; + u = -b/(2*_a); v = sqrt((double)fabs(dscrmnt))/(2*_a); + + if(real && !cplx || !real && cplx) + assert(true); + + if(real == NULL && cplx == NULL){ + if(a->hp != 2 && a->terms < 3){ + printf("Invalid Quadratic*, A and B must be non-zero"); + return; + } + + if(dscrmnt != 0) + printf("X = %.2f +/- %.2f%c\n",u,v,dscrmnt < 0 ? 'I':' '); + else{ + printf("(X %c %.2f)(X %c %.2f)\n",sgn(u),fabs(u),sgn(u),fabs(u)); + printf("X1,2 = %.2f\n",u); + } + } + //save values in pointers + else { + if(dscrmnt < 0){ //x = u +/- vI Re(x) = u, Im(x) = +v + *real = u; + *cplx = v; //understand +/- is not representable + } + else if(dscrmnt == 0){ + *real = u; + *cplx = 0.00F; + } + else{ + *real = u + v; + *cplx = u - v; + } + } +} + +/** +* exponentiate - computes polynomial exponentiation (P(x))^n, n E Z* +* @a: the polynomial +* @n: the exponent +* Works fast for small n (n < 8) currently runs ~ O(n^2 lg n) +*/ +PolyAdt *exponentiate(const PolyAdt *a, int n) +{ + assert(a != NULL); + + PolyAdt *expn = create_adt(a->hp * n); + PolyAdt *aptr = polyImage(a); + int hl = n / 2; + + //check default cases before calculation + if(n == 0){ + insert_term(expn, 1, 0); + return expn; + } + else if(n == 1){ + return aptr; + } + + for(; hl ; hl--) + aptr = multiply(aptr, aptr); + + if(n % 2) //odd exponent do a^(n-1) * a = a^n + expn = multiply(aptr, a); + else + expn = aptr; + return expn; +} +/** +* compose - computes the composition of two polynomials P(Q(x)) and returns the composition +* @p: polynomial P(x) which will x will be equal to Q(x) +* @q: polynomial Q(x) which is the argument to P(x) +*/ +PolyAdt *compose(const PolyAdt *p, const PolyAdt *q) +{ + assert(p && q); + + PolyAdt *comp = create_adt(p->head->exp * q->head->exp); + PolyAdt *exp; + + Node *pp = p->head; + Node *qq = q->head; + + int swap = 0; + + if(p->terms < q->terms){ + pp = q->head; + qq = p->head; + swap = 1; + } + + /* going through, exponentiate each term with the exponent of p */ + for(; pp != NULL; pp = pp->next){ + exp = exponentiate(swap ? p: q, pp->exp); + insert_term(comp, pp->coeff * exp->head->coeff, exp->head->exp); + } + + return comp; +} +/** +* destroy_poly - completely frees the polynomial from the heap and resets all values +* @poly: the polynomial to release memory back to the heap +* Usage: +* destroy_poly(myPoly); //puts polynomial on free list +*/ +void destroy_poly(PolyAdt *poly) +{ + Node *ps = poly->head; + Node *tmp = NULL; + while(ps != NULL){ + tmp = ps; + free(tmp); + ps = ps->next; + } + poly->hp = poly->terms = 0; + poly->head = NULL; +} +/** +* display_poly - displays the polynomial to the console in nice format +* @a: the polynomial to display +*/ +void display_poly(const PolyAdt *a) +{ + assert(a != NULL); + Node *n; + + for(n = a->head; n != NULL; n = n->next){ + + n->coeff < 0 ? putchar('-') : putchar('+'); + if(n->exp == 0) + printf(" %.2f ",fabs(n->coeff)); + else if(n->coeff == 1) + printf(" X^%d ",n->exp); + else if(n->exp == 1) + printf(" %.2fX ",fabs(n->coeff)); + else if(n->coeff == 0) + continue; + else + printf(" %.2fX^%d ",fabs(n->coeff),n->exp); + } + printf("\n\n"); +} + +#endif diff --git a/junkcode/iasoule32@gmail.com-polynomial/polynomial_adt_test.c b/junkcode/iasoule32@gmail.com-polynomial/polynomial_adt_test.c new file mode 100644 index 00000000..f7c9c8a5 --- /dev/null +++ b/junkcode/iasoule32@gmail.com-polynomial/polynomial_adt_test.c @@ -0,0 +1,132 @@ +/* +** polynomial_adt_test.c +** Test (minimalistic) for the polynomial module + * More of a display of functionality + * Copyright (c) 2009 I. Soule + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * +** iasoule32@gmail.com +*/ + +#include +#include "polynomial_adt.h" + +int main() +{ + PolyAdt *p = create_adt(5), *q = create_adt(4); + PolyAdt *sum, *diff, *prod; + + insert_term(p,1,5); + insert_term(p,3,4); + insert_term(p,1,3); + insert_term(p,9,2); + insert_term(p,8,1); + + insert_term(q,2,4); + insert_term(q,8,3); + insert_term(q,7,2); + insert_term(q,6,1); + + + printf("Displaying Polynomials ...\n"); + display_poly(p); + display_poly(q); + + sum = add(p,q); + printf("P(x) + Q(x) = "); + display_poly(sum); + + diff = subtract(p,q); + printf("P(x) - Q(x) = "); + display_poly(diff); + + prod = multiply(p,q); + printf("P(x)*Q(x) = "); + display_poly(prod); + + PolyAdt *quad = create_adt(2); + insert_term(quad, 10, 2); + insert_term(quad, 30, 1); + insert_term(quad, 2, 0); + + quadratic_roots(quad, NULL, NULL); //print out the roots + + float real, cplx; + quadratic_roots(quad, &real, &cplx); + + printf("X1 = %f, X2 = %f\n\n", real, cplx); + + PolyAdt *deriv, *integral; + + deriv = derivative(p); + printf("The derivitive of p = "); + display_poly(deriv); + integral = integrate(q); + + printf("The integral of q = "); + display_poly(integral); + + printf("\n Computing P(x)^3\n"); + + PolyAdt *expo; + expo = exponentiate(p, 3); + display_poly(expo); + printf("\n"); + + printf("Computing Integral[Q(x)^2]\n"); + expo = exponentiate(q, 2); + integral = integrate(expo); + display_poly(integral); + + + printf(" Differentiating and Integrating P\n"); + display_poly(integrate(derivative(p))); + + PolyAdt *L, *M; + + L = create_adt(3), M = create_adt(2); + + insert_term(L, 4, 3); + insert_term(L, 10, 2); + insert_term(L, 15, 1); + + insert_term(M, 4, 2); + printf("L = "); + display_poly(L); + printf("M = "); + display_poly(M); + + + printf("Computing composition L(M(X))\n"); + display_poly(compose(L, M)); + + printf("Freed memory back to heap for allocated polys'"); + destroy_poly(sum); + destroy_poly(diff); + destroy_poly(prod); + destroy_poly(L); destroy_poly(M); + destroy_poly(q); destroy_poly(p); + + return 0; +} -- 2.39.2