Skip to content

C
ccan
Commit f5d28c40 authored by Rusty Russell's avatar Rusty Russell
Browse files
Options

Joseph Adams updates to polynomial.

parent 9c1974f6
Show whitespace changes
Inline Side-by-side
Showing with 623 additions and 430 deletions
junkcode/iasoule32@gmail.com-polynomial/_info.c 0 → 100644
View file @ f5d28c40
#include <string.h>
#include <stdio.h>
#include "config.h"
/**
* polynomial_adt - A polynomial module with ability to add,sub,mul derivate/integrate, compose ... polynomials
*
* ..expansion in progress ...
*
* Example:
* FULLY-COMPILABLE-INDENTED-TRIVIAL-BUT-USEFUL-EXAMPLE-HERE
*/
int main(int argc, char *argv[])
{
if (argc != 2)
return 1;
if (strcmp(argv[1], "depends") == 0) {
/* Nothing. */
return 0;
}
if (strcmp(argv[1], "libs") == 0) {
printf("m\n");
return 0;
}
return 1;
}
junkcode/iasoule32@gmail.com-polynomial/polynomial_adt.c 0 → 100644
View file @ f5d28c40
#include "polynomial_adt.h"
PolyAdt *create_adt(int hp)
{
PolyAdt *pAdt = malloc(sizeof(PolyAdt));
assert(pAdt != NULL);
pAdt->head = NULL;
pAdt->terms = 0;
pAdt->hp = hp;
return pAdt;
}
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++;
}
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;
}
PolyAdt *add(const PolyAdt *a, const PolyAdt *b)
{
PolyAdt *sum;
Node *n, *np;
int state = 1;
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){ //copy rest of a to sum
for(; n != NULL; n = n->next)
insert_term(sum, n->coeff, n->exp);
state = 0;
}
if(n == NULL && state){
for(; np != NULL; np = np->next)
insert_term(sum, np->coeff, np->exp);
state = 0;
}
}
return sum;
}
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);
}
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;
}
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;
}
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;
}
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(1);
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;
}
}
}
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;
}
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;
}
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;
}
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");
}
junkcode/iasoule32@gmail.com-polynomial/polynomial_adt.h
View file @ f5d28c40
......@@ -34,7 +34,7 @@
#include <assert.h>
#include <stdlib.h>
#include <stdbool.h> //C99 compliance
#include <stdio.h>
#include <math.h>
#define max(a, b) (a) > (b) ? (a) : (b)
......@@ -51,22 +51,12 @@ typedef struct polynomial_adt {
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;
PolyAdt *create_adt(int hp);
return pAdt;
}
/**
* create_node - creates a Node (exponent, constant and next pointer) on the heap
* @constant: the contant in the term
......@@ -76,8 +66,7 @@ PolyAdt *create_adt(int hp)
* 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)
{
static inline Node *create_node(float constant, int exp, Node *next) {
Node *nNode = malloc(sizeof(Node));
assert(nNode != NULL);
......@@ -93,33 +82,13 @@ static Node *create_node(float constant, int exp, Node *next)
* @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++;
}
void insert_term(PolyAdt *pAdt, float c, int e);
/**
* 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;
}
PolyAdt *polyImage(const PolyAdt *orig);
/**
......@@ -127,54 +96,7 @@ PolyAdt *polyImage(const PolyAdt *orig)
* @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;
}
PolyAdt *add(const PolyAdt *a, const PolyAdt *b);
/**
* sub - subtracts two polynomials, and returns their difference (as a polynomial)
......@@ -182,82 +104,29 @@ PolyAdt *add(const PolyAdt *a, const PolyAdt *b)
* @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);
}
PolyAdt *subtract(const PolyAdt *a, const PolyAdt *b);
/**
* 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;
}
PolyAdt *multiply(const PolyAdt *a, const PolyAdt *b);
/**
* 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 *derivative(const PolyAdt *a);
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;
PolyAdt *integrate(const PolyAdt *a);
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
......@@ -274,51 +143,7 @@ PolyAdt *integrate(const PolyAdt *a)
* 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;
}
}
}
void quadratic_roots(const PolyAdt *a, float *real, float *cplx);
/**
* exponentiate - computes polynomial exponentiation (P(x))^n, n E Z*
......@@ -326,105 +151,27 @@ void quadratic_roots(const PolyAdt *a, float *real, float *cplx)
* @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;
PolyAdt *exponentiate(const PolyAdt *a, int n);
//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 *compose(const PolyAdt *p, const PolyAdt *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;
}
void destroy_poly(PolyAdt *poly);
/**
* 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");
}
void display_poly(const PolyAdt *a);
#endif
junkcode/iasoule32@gmail.com-polynomial/test/run.c 0 → 100644
View file @ f5d28c40
/*
** 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 <stdio.h>
#include "../polynomial_adt.h"
#include "../polynomial_adt.c"
int main(void)
{
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'\n");
destroy_poly(sum);
destroy_poly(diff);
destroy_poly(prod);
destroy_poly(L); destroy_poly(M);
destroy_poly(q); destroy_poly(p);
return 0;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment
GitLab Nexedi Edition | About GitLab | About Nexedi | 沪ICP备2021021310号-2 | 沪ICP备2021021310号-7