New polynomial module (not quite CCAN format).
authorRusty Russell <rusty@rustcorp.com.au>
Sat, 4 Apr 2009 01:29:54 +0000 (11:59 +1030)
committerRusty Russell <rusty@rustcorp.com.au>
Sat, 4 Apr 2009 01:29:54 +0000 (11:59 +1030)
junkcode/iasoule32@gmail.com-polynomial/polynomial_adt.h [new file with mode: 0644]
junkcode/iasoule32@gmail.com-polynomial/polynomial_adt_test.c [new file with mode: 0644]

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