zpoly.c

FIR2LIFT is a program to factor a wavelet FIR filter pair into lifting steps.

/*
FILE : ZPOLY.C

Routines to handle Laurent polynomials (in z). These differ from normal 
polynomials in that they can have negative powers.

(C) C. Valens

Created     : 10/09/1999
Last update : 25/09/1999
*/


#include "redir.h"
#include "zpoly.h"

#include <stdlib.h>
#include <stdio.h>
#include <math.h>


void zpoly_error(const int error)
{
  redir_printf("Exiting due to error %d.\n", error);
  exit(0);
}


/*
 * Print a polynomial.
 */
void zpoly_write(zpoly *poly)
{
  zpoly *p = poly;
  char first_sign=1;

  while (p) {

    /*
     * Do not show coefficients that are too small.
     */
    if (fabs(p->coeff)>=EPSILON) {

      /*
       * Show sign.
       */
      if (first_sign) {
        if (p->coeff<0.0) redir_printf("-");
        first_sign = 0;
      }
      else {
        if (p->coeff<0.0) redir_printf("- ");
        else redir_printf("+ ");
      }

      /*
       * Show coefficient, but don't show 1.0.
       */
      if ((fabs(p->coeff)<=(1.0-EPSILON)) ||
          (fabs(p->coeff)>=(1.0+EPSILON)) ||
          (p->power==0)) {
        redir_printf("%f",fabs(p->coeff));
        if (p->power!=0) {
          redir_printf(" ");
        }
      }

      /*
       * Show power.
       */
      if (p->power<0){
        redir_printf("z^(%d)",p->power);
      }
      else if (p->power==1) {
        redir_printf("z");
      }
      else if (p->power>1){
        redir_printf("z^%d",p->power);
      }
    }

    if (p->next) {
      redir_printf(" ",p->power);
    }

    p = p->next;
  }
}


/*
 * Print a polynomial with optional degree.
 */
void zpoly_writeln(zpoly *poly, const int show_degree)
{
  zpoly *p = poly;
  char first_sign=1;

  if (show_degree) {
    if (p) {
      redir_printf("(degree=%d) : ",zpoly_degree(p));
    }
    else {
      redir_printf("(degree=-1) : 0");
    }
  }

  zpoly_write(p);


  redir_printf("\n");
}


/*
 * Create a new polynomial.
 */
zpoly *zpoly_create(const int power, const double coeff)
{
  zpoly *p=NULL;
  if ((p=malloc(sizeof(zpoly)))==NULL) {
    redir_printf("zpoly_create() : ");
    zpoly_error(1);
  }
  p->power = power;
  p->coeff = coeff;
  p->next = NULL;
  return p;
}


/*
 * Destroy a polynomial.
 */
void zpoly_destroy(zpoly *poly)
{
  zpoly *p1 = poly;
  zpoly *p2;

  while (p1) {
    p2 = p1;
    p1 = p1->next;
    free(p2);
  }

}


/*
 * Create a copy of poly.
 */
zpoly *zpoly_copy(zpoly *poly)
{
  zpoly *p1 = poly;
  zpoly *p2 = NULL;

  if (p1) {
    p2 = zpoly_create(p1->power,p1->coeff);
    p1 = p1->next;
    p2 = zpoly_add_to(p2,p1);
  }

  return p2;
}


/*
 * Returns 1 if the two polynomials are equal, 0 if not.
 */
int zpoly_equal(zpoly *poly1, zpoly *poly2)
{
  zpoly *p1 = poly1;
  zpoly *p2 = poly2;
  int equal = 1;

  while ((p1) && (p2)) {
    if (p1->power!=p2->power) {
      equal = 0;
      break;
    }
    if (fabs(p1->coeff-p2->coeff)>=EPSILON) {
      equal = 0;
      break;
    }
    p1 = p1->next;
    p2 = p2->next;
  }

  if ((p1) || (p2)) return 0;
  else return equal;
}


/*
 * Returns 1 if the polynomial is equal to 0.
 */
int zpoly_equals_zero(zpoly *poly)
{
  zpoly *p = poly;
  int equal = 1;

  while (p) {
    if (fabs(p->coeff)>=EPSILON) {
      equal = 0;
      break;
    }
    p = p->next;
  }

  if (p) return 0;
  else return equal;
  return 0;
}


/*
 * All elements with a coefficient smaller than EPSILON are removed
 * from poly.
 */
zpoly *zpoly_remove_zeroes(zpoly *poly)
{
  zpoly *p1 = poly;
  zpoly *p2;
  zpoly *pnew = NULL;

  p2 = p1;
  while (p1) {
    if (fabs(p1->coeff)<EPSILON) {
      if (p1->power==p2->power) {
        pnew = p1->next;
      }
      p2 = p1->next;
      free(p1);
      p1 = p2;
    }
    else {
      p2 = p1;
      p1 = p1->next;
    }
  }
  return pnew;
}


/*
 * The degree of a polynomial is defined as the absolute value of the
 * difference between its highest power and its lowest power. A monomial
 * therefore has degree 0.
 */
int zpoly_degree(zpoly *poly)
{
  zpoly *p = poly;
  int min, max;

  if (p) {
    min = p->power;
    while (p) {
      max = p->power;
      p = p->next;
    }
    return abs(max-min);
  }
  return -1;
}


/*
 * Reverses a polynomial. If reverse_time==1 all powers are multiplied
 * by -1, which corresponds to the reversal of time for polynomials in z.
 */
zpoly *zpoly_reverse(zpoly *poly, const int reverse_time)
{
  zpoly *p1 = poly;
  zpoly *p2;
  zpoly *p3 = NULL;

  while (p1) {
    if (reverse_time==1) {
      p1->power *= -1;
    }
    p2 = p1;
    p1 = p1->next;
    p2->next = p3;
    p3 = p2;
  }

  return p3;
}


/*
 * Creates a reversed version of poly. If reverse_time==1 all powers
 * will be negated.
 */
zpoly *zpoly_copy_reverse(zpoly *poly, int reverse_time)
{
  zpoly *p1 = poly;
  zpoly *p2;
  zpoly *p3 = NULL;

  while (p1) {
    if (reverse_time) {
      p2 = zpoly_create(-p1->power,p1->coeff);
    }
    else {
      p2 = zpoly_create(p1->power,p1->coeff);
    }
    if (p2) {
      p2->next = p3;
      p3 = p2;
    }
    p1 = p1->next;
  }

  return p3;
}


/*
 * Adds two polynomials and puts the result in a new polynomial.
 *
 * Actually this is a list merge, except that elements with the same
 * power will have their coefficients added together.
 */
zpoly *zpoly_add(zpoly *poly1, zpoly *poly2)
{
  zpoly *p1 = poly1;
  zpoly *p2 = poly2;
  zpoly *p3;
  zpoly *temp;
  zpoly *result = NULL;

  while ((p1) || (p2)) {

    temp = NULL;

    if ((p1) && (p2)) {
      if (p1->power<p2->power) {
        if (fabs(p1->coeff)>=EPSILON) {
          temp = zpoly_create(p1->power,p1->coeff);
        }
        p1 = p1->next;
      }
      else if (p2->power<p1->power) {
        if (fabs(p2->coeff)>=EPSILON) {
          temp = zpoly_create(p2->power,p2->coeff);
        }
        p2 = p2->next;
      }
      else {
        if (fabs(p1->coeff+p2->coeff)>=EPSILON) {
          temp = zpoly_create(p1->power,p1->coeff+p2->coeff);
        }
        p1 = p1->next;
        p2 = p2->next;
      }
    }
    else if (p1) {
      if (fabs(p1->coeff)>=EPSILON) {
        temp = zpoly_create(p1->power,p1->coeff);
      }
      p1 = p1->next;
    }
    else if (p2) {
      if (fabs(p2->coeff)>=EPSILON) {
        temp = zpoly_create(p2->power,p2->coeff);
      }
      p2 = p2->next;
    }

    if (result==NULL) {
      result = temp;
      p3 = result;
    }
    else if (temp) {
      p3->next = temp;
      p3 = p3->next;
    }

  }

  return result;
}


/*
 * Adds polynomial poly2 to poly1.
 * usage : poly1 = polynomial_add_to(poly1,poly2);
 */
zpoly *zpoly_add_to(zpoly *poly1, zpoly *poly2)
{
  zpoly *p;

  p = zpoly_add(poly1,poly2);
  zpoly_destroy(poly1);
  return p;
}


/*
 * Adds an element to polynomial poly.
 * usage : poly1 = polynomial_add_to(poly,exponent,coefficient);
 */
zpoly *zpoly_add_element_to(zpoly *poly, int power, double coeff)
{
  zpoly *p1, *p2;

  p1 = zpoly_create(power,coeff);
  p2 = zpoly_add(poly,p1);
  zpoly_destroy(poly);
  zpoly_destroy(p1);
  return p2;
}


/*
 * Subtracts poly2 from poly1 and puts the result in a new polynomial.
 */
zpoly *zpoly_sub(zpoly *poly1, zpoly *poly2)
{
  zpoly *p1 = poly1;
  zpoly *p2 = poly2;
  zpoly *result = NULL;
  zpoly *temp;

  while ((p1) || (p2)) {

    temp = NULL;

    if ((p1) && (p2)) {
      if (p1->power<p2->power) {
        if (fabs(p1->coeff)>=EPSILON) {
          temp = zpoly_create(p1->power,p1->coeff);
        }
        p1 = p1->next;
      }
      else if (p1->power>p2->power) {
        if (fabs(p2->coeff)>=EPSILON) {
          temp = zpoly_create(p2->power,-p2->coeff);
        }
        p2 = p2->next;
      }
      else {
        if (fabs(p1->coeff-p2->coeff)>=EPSILON) {
          temp = zpoly_create(p1->power,p1->coeff-p2->coeff);
        }
        p1 = p1->next;
        p2 = p2->next;
      }
    }
    else if (p1) {
      if (fabs(p1->coeff)>=EPSILON) {
        temp = zpoly_create(p1->power,p1->coeff);
      }
      p1 = p1->next;
    }
    else if (p2) {
      if (fabs(p2->coeff)>=EPSILON) {
        temp = zpoly_create(p2->power,-p2->coeff);
      }
      p2 = p2->next;
    }

    if (result==NULL) {
      result = temp;
    }
    else {
      result = zpoly_add_to(result,temp);
      zpoly_destroy(temp);
    }

  }

  return result;
}


/*
 * Subtracts poly2 from poly1.
 * usage : poly1 = polynomial_sub_from(poly1,poly2);
 */
zpoly *zpoly_sub_from(zpoly *poly1, zpoly *poly2)
{
  zpoly *p;

  p = zpoly_sub(poly1,poly2);
  zpoly_destroy(poly1);
  return p;
}


/*
 * Multiplies a polynomial by a constant.
 */
zpoly *zpoly_scale(zpoly *poly1, const double s)
{
  zpoly *p = poly1;

  while (p) {
    p->coeff *= s;
    p = p->next;
  }

  return poly1;
}


/*
 * Multiplies two polynomials and puts the result in a new polynomial.
 */
zpoly *zpoly_mul(zpoly *poly1, zpoly *poly2)
{
  zpoly *p1 = poly1;
  zpoly *p2 = poly2;
  zpoly *temp;
  zpoly *result = NULL;
  double c;

  while (p1) {
    p2 = poly2;
    while (p2) {

      temp = NULL;

      c = p1->coeff*p2->coeff;
      if (fabs(c)>=EPSILON) {
        temp = zpoly_create(p1->power+p2->power,c);
      }

      if (result==NULL) {
        result = temp;
      }
      else {
        result = zpoly_add_to(result,temp);
        zpoly_destroy(temp);
      }

      p2 = p2->next;
    }

    p1 = p1->next;
  }

  return result;
}


/*
 * Multiplies polynomial poly1 by poly2.
 * usage : poly1 = polynomial_mul_with(poly1,poly2);
 */
zpoly *zpoly_mul_by(zpoly *poly1, zpoly *poly2)
{
  zpoly *p;

  p = zpoly_mul(poly1,poly2);
  zpoly_destroy(poly1);
  return p;
}


/*
 * Does one step of a polynomial div/mod.
 * Matches the first powers of the two polynomials.
 */
zpoly_pair zpoly_divmod_fp_step(zpoly *poly1, zpoly *poly2)
{
  zpoly *p1 = poly1;
  zpoly *p2 = poly2;
  zpoly *temp = NULL;
  zpoly_pair divmod = {NULL,NULL};
  double c;

  if ((p1==NULL) || (p2==NULL)) {
    return divmod;
  }

  /*
   * poly2 has to be of lower degree than poly1.
   */
  if (zpoly_degree(p2) > zpoly_degree(p1)) {
    return divmod;
  }

  c = p2->coeff;
  if (c) {
    c = p1->coeff/c;
    if (fabs(c)>=EPSILON) {
      divmod.p1 = zpoly_create(p1->power-p2->power,c);
      temp = zpoly_mul(p2,divmod.p1);
      divmod.p2 = zpoly_sub(p1,temp);
    }
  }

  zpoly_destroy(temp);

  return divmod;
}


/*
 * Does one step of a polynomial div/mod.
 * Matches the first powers of the two polynomials.
 */
zpoly_pair zpoly_divmod_lp_step(zpoly *poly1, zpoly *poly2)
{
  zpoly *p1 = zpoly_copy_reverse(poly1,1);
  zpoly *p2 = zpoly_copy_reverse(poly2,1);
  zpoly *temp = NULL;
  zpoly_pair divmod = {NULL,NULL};
  double c;

  if ((p1==NULL) || (p2==NULL)) {
    return divmod;
  }

  /*
   * poly2 has to be of lower degree than poly1.
   */
  if (zpoly_degree(p2) > zpoly_degree(p1)) {
    return divmod;
  }

  c = p2->coeff;
  if (c) {
    c = p1->coeff/c;
    if (fabs(c)>=EPSILON) {
      divmod.p1 = zpoly_create(p1->power-p2->power,c);
      temp = zpoly_mul(p2,divmod.p1);
      divmod.p2 = zpoly_sub(p1,temp);
    }
  }

  zpoly_destroy(p1);
  zpoly_destroy(p2);
  zpoly_destroy(temp);

  divmod.p1 = zpoly_reverse(divmod.p1,1);
  divmod.p2 = zpoly_reverse(divmod.p2,1);

  return divmod;
}


/*
 * Divides two polynomials and returns the result as a polynomial pair:
 *
 *   pair.p1 = poly1 div poly2
 *   pair.p2 = poly1 mod poly2
 *
 * If match==0 the first power of poly1 is matched,
 * else the last power of poly1 is matched.
 *
 * Make sure that the polynomials are correctly ordered and that the
 * degree of poly2 is smaller than the degree of poly1.
 */
zpoly_pair zpoly_divmod(zpoly *poly1, zpoly *poly2, const int match)
{
  zpoly_pair divmod = {NULL,NULL};
  zpoly_pair temp;

  /*
   * poly2 has to be of lower degree than poly1.
   */
  if (zpoly_degree(poly2) > zpoly_degree(poly1)) {
    return divmod;
  }

  divmod.p1 = NULL;
  divmod.p2 = zpoly_copy(poly1);

  while ((divmod.p2) && (zpoly_degree(divmod.p2))) {

    if (match==0) {
      temp = zpoly_divmod_fp_step(divmod.p2,poly2);
    }
    else {
      temp = zpoly_divmod_lp_step(divmod.p2,poly2);
    }

    if (divmod.p1==NULL) {
      divmod.p1 = zpoly_create(temp.p1->power,temp.p1->coeff);
    }
    else {
      divmod.p1 = zpoly_add_to(divmod.p1,temp.p1);
      zpoly_destroy(temp.p1);
    }

    zpoly_destroy(divmod.p2);
    divmod.p2 = temp.p2;
  }

  return divmod;
}