gcd.c

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

/*
FILE : GCD.C

Implementation of the Euclidean algorithm for Laurent polynomials.

(C) C. Valens

Created     : 16/09/1999
Last update : 24/09/1999
*/


//#define debug


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

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


/*
 * Tree stuff
 */
typedef struct __gcd_tree {
  zpoly_pair divmod;
  struct __gcd_tree *left;
  struct __gcd_tree *right;
} gcd_tree;


/*
 * Local prototypes.
 */
gcd_list *gcd_remove_multiples(gcd_list *l);
gcd_list *gcd_tree_to_list(gcd_tree *root);


gcd_tree *gcd_create_tree(void)
{
  gcd_tree *p=NULL;
  if ((p=malloc(sizeof(gcd_tree)))==NULL) {
    redir_printf("gcd_create_tree() : ");
    zpoly_error(1);
  }
  p->divmod.p1 = NULL;
  p->divmod.p2 = NULL;
  p->left = NULL;
  p->right = NULL;
  return p;
}


/*
 * Internal function.
 */
void __gcd_destroy_tree(gcd_tree *root)
{
  if (root) {
    __gcd_destroy_tree(root->left);
    free(root->left);
    __gcd_destroy_tree(root->right);
    free(root->right);
  }
}


void gcd_destroy_tree(gcd_tree *root)
{
  __gcd_destroy_tree(root);
  free(root);
}


/*
 * Builds a tree of all possible div-mod combinations of a division of
 * a polynomial in z.
 */
gcd_tree *gcd_divmod_tree(zpoly *poly1, zpoly *poly2)
{
  zpoly_pair temp;
  gcd_tree *root = NULL;
  int d1 = zpoly_degree(poly1);

  root = gcd_create_tree();

#ifdef debug
  redir_printf("{gcd([");
  zpoly_write(poly1);
  redir_printf("],[");
  zpoly_write(poly2);
  redir_printf("]}\n");
#endif

  if ((d1>0) && (d1>=zpoly_degree(poly2))) {

    temp = zpoly_divmod_fp_step(poly1,poly2);
    root->left = gcd_divmod_tree(temp.p2,poly2);
    root->left->divmod = temp;

    temp = zpoly_divmod_lp_step(poly1,poly2);
    root->right = gcd_divmod_tree(temp.p2,poly2);
    root->right->divmod = temp;
  }

  return root;
}


void gcd_write_tree(gcd_tree *root)
{
  if (root) {
    redir_printf("{div} "); zpoly_writeln(root->divmod.p1,0);
    redir_printf("{mod} "); zpoly_writeln(root->divmod.p2,0);
    gcd_write_tree(root->left);
    gcd_write_tree(root->right);
  }
}


gcd_list *gcd_create_list(void)
{
  gcd_list *p=NULL;
  if ((p=malloc(sizeof(gcd_list)))==NULL) {
    redir_printf("gcd_create_list() : ");
    zpoly_error(1);
  }
  p->divmod.p1 = NULL;
  p->divmod.p2 = NULL;
  p->next = NULL;
  return p;
}


void gcd_destroy_list(gcd_list *root)
{
  gcd_list *p1 = root;
  gcd_list *p2 = root;

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

}


/*
 * Internal function.
 */
static gcd_list *__gcd_tree_to_list(gcd_tree *root)
{
  gcd_tree *t = root;
  gcd_list *list = NULL;
  gcd_list *l1;
  gcd_list *l2;
  zpoly *p;
  int path_count;
  int depth = -1;
  unsigned int mask;

  if (root==NULL) return NULL;

  while (t) {
    depth++;
    t = t->left;
  }
  path_count = 1<<depth;

  while (path_count) {

    l1 = gcd_create_list();
    mask = 1<<(depth-1);
    t = root;

    while (t) {

      /*
       * By using the create/add_to method to build the quotient,
       * we assure that the resulting polynomial will be ordered correctly.
       */
      if (l1->divmod.p1==NULL) {
        l1->divmod.p1 = zpoly_copy(t->divmod.p1);
      }
      else {
        l1->divmod.p1 = zpoly_add_to(l1->divmod.p1,t->divmod.p1);
      }

      /*
       * Keep the last remainder.
       */
      p = t->divmod.p2;

      /*
       * Select next sub-tree.
       */
      if ((path_count-1)&mask) {
        t = t->left;
      }
      else {
        t = t->right;
      }
      mask >>= 1;

    }

    /*
     * Copy the remainder of the division.
     */
    l1->divmod.p2 = zpoly_copy(p);

    /*
     * Add to list.
     */
    if (list==NULL) {
      list = l1;
      l2 = list;
    }
    else {
      l2->next = l1;
      l2 = l2->next;
    }

    path_count--;
  }

  return list;
}


/*
 * Converts a zpoly tree to a list of quotients and remainders.
 */
gcd_list *gcd_tree_to_list(gcd_tree *root)
{
  gcd_list *list=NULL;
  gcd_list *l=NULL;

  if (root==NULL) return NULL;

  list = gcd_create_list();
  list->divmod.p1 = zpoly_copy(root->divmod.p1);
  list->divmod.p2 = zpoly_copy(root->divmod.p2);

  /*
   * Because the root of a tree contains the input polynomials
   * we use this construction to skip it further on.
   */
  l = list;

  /*
   * Add left sub-tree.
   */
  l->next = __gcd_tree_to_list(root->left);

  /*
   * Go to the end of the list.
   */
  while (l->next) {
    l = l->next;
  }

  /*
   * Add right sub-tree.
   */
  l->next = __gcd_tree_to_list(root->right);

  return list;
}


/*
 * Display a gcd_list.
 */
void gcd_write_list(gcd_list *list)
{
  gcd_list *l = list;

  while (l) {
    redir_printf("{div} ");
    zpoly_writeln(l->divmod.p1,0);
    redir_printf("{mod} ");
    zpoly_writeln(l->divmod.p2,0);
    l = l->next;
  }

}


/*
 * Returns a list of all div-mod pairs of poly1 & poly2.
 * Double entries are removed.
 */
gcd_list *gcd_euclides(zpoly *poly1, zpoly *poly2)
{
  gcd_tree *divmod_tree;
  gcd_list *divmod_list;
  gcd_list *temp;

  divmod_tree = gcd_divmod_tree(poly1,poly2);

#ifdef debug
  redir_printf("{GCD TREE BEGIN}\n");
  gcd_write_tree(divmod_tree);
  redir_printf("{GCD TREE END}\n");
#endif

  temp = gcd_tree_to_list(divmod_tree);
  gcd_destroy_tree(divmod_tree);
  divmod_list = temp->next;
  free(temp);

  divmod_list = gcd_remove_multiples(divmod_list);

  return divmod_list;
}


/*
 * Remove multiple entries from a divmod list.
 */
gcd_list *gcd_remove_multiples(gcd_list *l)
{
  gcd_list *t1, *t2, *t3;

  t1 = l;
  while (t1) {
    t3 = t1;
    t2 = t1->next;
    while (t2) {
      if ((zpoly_equal(t1->divmod.p1,t2->divmod.p1)) &&
          (zpoly_equal(t1->divmod.p2,t2->divmod.p2))) {
        t3->next = t2->next;
        free(t2);
        t2 = t3;
      }
      t3 = t2;
      t2 = t2->next;
    }
    t1 = t1->next;
  }

  return l;
}