/* This is -*- C -*- */
/* vim: set sw=2: */
/* $Id: guppi-polynomial.c,v 1.7 2002/01/14 05:01:23 trow Exp $ */

/*
 * guppi-polynomial.c
 *
 * Copyright (C) 2001 The Free Software Foundation

 * Contains code
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * Otherwise developed by Jon Trowbridge <trow@gnu.org>
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA
 */

#include <config.h>
#include "guppi-polynomial.h"

#include <gtk/gtk.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <guppi-memory.h>

static GtkObjectClass *parent_class = NULL;

enum {
  CHANGED,
  LAST_SIGNAL
};

static guint guppi_polynomial_signals[LAST_SIGNAL] = { 0 };

typedef struct _GuppiPolynomialPrivate GuppiPolynomialPrivate;
struct _GuppiPolynomialPrivate {
  gint d;     /* degree */
  gint N;     /* allocated degree */
  double *c;  /* coefficients */

  gint num_roots;
  double *roots;
  
  gint num_minmax;
  double *minmax;

  gboolean process_lock;
  gint freeze_count;
  gboolean touched;
};

#define priv(x) ((GuppiPolynomialPrivate*)(GUPPI_POLYNOMIAL((x))->opaque_internals))

#define POLYEPSILON 1e-12

static void
guppi_polynomial_finalize (GtkObject *obj)
{
  GuppiPolynomial *x = GUPPI_POLYNOMIAL(obj);
  GuppiPolynomialPrivate *p = priv (x);

  guppi_free0 (p->c);
  guppi_free0 (p->roots);
  guppi_free0 (p->minmax);

  g_free (x->opaque_internals);
  x->opaque_internals = NULL;

  guppi_finalized (obj);

  if (parent_class->finalize)
    parent_class->finalize (obj);
}

static void
guppi_polynomial_class_init (GuppiPolynomialClass *klass)
{
  GtkObjectClass *object_class = (GtkObjectClass *)klass;

  parent_class = gtk_type_class (GTK_TYPE_OBJECT);

  object_class->finalize = guppi_polynomial_finalize;

  /* Signal definition template */
  guppi_polynomial_signals[CHANGED] =
    gtk_signal_new ("changed",
                    GTK_RUN_FIRST,
                    object_class->type,
                    GTK_SIGNAL_OFFSET (GuppiPolynomialClass, changed),
                    gtk_marshal_NONE__NONE, GTK_TYPE_NONE, 0);

  gtk_object_class_add_signals (object_class, guppi_polynomial_signals,
                                LAST_SIGNAL);
}

static void
guppi_polynomial_init (GuppiPolynomial *obj)
{
  GuppiPolynomialPrivate *p = g_new0 (GuppiPolynomialPrivate, 1);
  obj->opaque_internals = p;

  p->num_roots = -1;
  p->num_minmax = -1;

  p->d = 0;
  p->N = 2;
  p->c = guppi_new0 (double, p->N+1);
  p->c[0] = 1.0;
}

GtkType
guppi_polynomial_get_type (void)
{
  static GtkType guppi_polynomial_type = 0;
  if (!guppi_polynomial_type) {
    static const GtkTypeInfo guppi_polynomial_info = {
      "GuppiPolynomial",
      sizeof (GuppiPolynomial),
      sizeof (GuppiPolynomialClass),
      (GtkClassInitFunc)guppi_polynomial_class_init,
      (GtkObjectInitFunc)guppi_polynomial_init,
      NULL, NULL, (GtkClassInitFunc)NULL
    };
    guppi_polynomial_type = gtk_type_unique (GTK_TYPE_OBJECT, &guppi_polynomial_info);
  }
  return guppi_polynomial_type;
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

static inline void
guppi_polynomial_reset_cache (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p = priv (poly);

  p->num_roots = -1;
  guppi_free0 (p->roots);

  p->num_minmax = -1;
  guppi_free0 (p->minmax);
}

static inline void
guppi_polynomial_changed (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p = priv (poly);

  if (p->freeze_count > 0) {
    p->touched = TRUE;
  } else {
    gtk_signal_emit (GTK_OBJECT (poly), guppi_polynomial_signals[CHANGED]);
    p->touched = FALSE;
  }
}

static void
guppi_polynomial_sanitize (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p = priv (poly);

  while (p->d > 0 && fabs (p->c[p->d]) < POLYEPSILON) {
    --p->d;
    p->touched = TRUE;
  }
    
}

static void
guppi_polynomial_grow (GuppiPolynomial *poly, gint N)
{
  GuppiPolynomialPrivate *p = priv (poly);

  if (N <= p->N)
    return;

  p->c = guppi_realloc (p->c, sizeof (double) * (N+1));
  p->N = N;
}

static void
guppi_polynomial_polish_cached_roots (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p = priv (poly);
  gint i;

  if (p->num_roots <= 0 || p->roots == NULL)
    return;

  for (i=0; i < p->num_roots; ++i) {
    p->roots[i] = guppi_polynomial_newton_polish (poly, p->roots[i], 5, POLYEPSILON);
  }
}

static gint
double_cmp (const void *a, const void *b)
{
  return (*(double *)a < *(double *)b) - (*(double *)a > *(double *)b);
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

GuppiPolynomial *
guppi_polynomial_new (gint degree, ...)
{
  GuppiPolynomial *poly;
  GuppiPolynomialPrivate *p;
  va_list args;
  gint i;

  g_return_val_if_fail (degree >= 0, NULL);

  poly = gtk_type_new (GUPPI_TYPE_POLYNOMIAL);
  p = priv (poly);

  guppi_polynomial_grow (poly, degree);
  p->d = degree;

  va_start (args, degree);
  for (i=0; i<degree+1; ++i) {
    p->c[i] = va_arg (args, double);
  }
  va_end (args);
  
  guppi_polynomial_sanitize (poly);

  return poly;
}

GuppiPolynomial *
guppi_polynomial_newv (gint degree, const double *vec)
{
  GuppiPolynomial *poly;
  GuppiPolynomialPrivate *p;
  gint i;

  g_return_val_if_fail (degree >= 0, NULL);
  g_return_val_if_fail (vec, NULL);

  poly = gtk_type_new (GUPPI_TYPE_POLYNOMIAL);
  p = priv (poly);

  guppi_polynomial_grow (poly, degree);
  p->d = degree;

  for (i=0; i<degree+1; ++i) {
    p->c[i] = vec[i];
  }

  guppi_polynomial_sanitize (poly);

  return poly;
}

GuppiPolynomial *
guppi_polynomial_new_by_roots (gint degree, ...)
{
  GuppiPolynomial *poly;
  va_list args;
  gint i;

  g_return_val_if_fail (degree >= 0, NULL);

  poly = guppi_polynomial_new_constant (1);
  guppi_polynomial_freeze (poly);
  guppi_polynomial_grow (poly, degree);

  va_start (args, degree);
  for (i=0; i<degree; ++i) {
    guppi_polynomial_inflate (poly, va_arg (args, double));
  }
  va_end (args);

  guppi_polynomial_thaw (poly);

  return poly;
}

GuppiPolynomial *
guppi_polynomial_new_by_rootsv (gint degree, const double *vec) 
{
  GuppiPolynomial *poly;
  gint i;

  g_return_val_if_fail (degree >= 0, NULL);

  poly = guppi_polynomial_new_constant (1);
  guppi_polynomial_freeze (poly);
  guppi_polynomial_grow (poly, degree);

  for (i=0; i<degree; ++i) {
    guppi_polynomial_inflate (poly, vec[i]);
  }

  guppi_polynomial_thaw (poly);

  return poly;
}


GuppiPolynomial *
guppi_polynomial_new_constant (double c0)
{
  return guppi_polynomial_new (0, c0);
}

GuppiPolynomial *
guppi_polynomial_new_linear (double c0, double c1)
{
  return guppi_polynomial_new (1, c0, c1);
}

GuppiPolynomial *
guppi_polynomial_new_quadratic (double c0, double c1, double c2)
{
  return guppi_polynomial_new (2, c0, c1, c2);
}

GuppiPolynomial *
guppi_polynomial_copy (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p;
  GuppiPolynomial *cpy;
  GuppiPolynomialPrivate *cpy_p;

  if (poly == NULL)
    return NULL;

  p = priv (poly);

  cpy = guppi_polynomial_newv (p->d, p->c);
  cpy_p = priv (cpy);

  cpy_p->num_roots = p->num_roots;
  if (p->num_roots > 0) {
    cpy_p->roots = guppi_new0 (double, p->num_roots);
    memcpy (cpy_p->roots, p->roots, sizeof (double) * p->num_roots);
  }

  cpy_p->num_minmax = p->num_minmax;
  if (p->num_minmax > 0) {
    cpy_p->minmax = guppi_new0 (double, p->num_minmax);
    memcpy (cpy_p->minmax, p->minmax, sizeof (double) * p->num_minmax);
  }

  return cpy;
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

void
guppi_polynomial_freeze (GuppiPolynomial *poly)
{
  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  ++ priv (poly)->freeze_count;
}

void
guppi_polynomial_thaw (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p;
  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  p = priv (poly);
  g_return_if_fail (p->freeze_count > 0);
  --p->freeze_count;
  if (p->freeze_count == 0 && p->touched)
    guppi_polynomial_changed (poly);
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

gint
guppi_polynomial_degree (GuppiPolynomial *poly)
{
  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), -1);

  return priv (poly)->d;
}

double
guppi_polynomial_coefficient (GuppiPolynomial *poly, gint i)
{
  GuppiPolynomialPrivate *p;
  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), 0);
  p = priv (poly);

  return (0 <= i && i <= p->d) ? p->c[i] : 0;
}

void
guppi_polynomial_set_coefficient (GuppiPolynomial *poly, gint i, double c)
{
  GuppiPolynomialPrivate *p;
  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  g_return_if_fail (i >= 0);
  p = priv (poly);

  if (fabs (c - (i <= p->d ? p->c[i] : 0)) < POLYEPSILON)
    return;

  if (i > p->d && fabs (c) > POLYEPSILON) {
    guppi_polynomial_grow (poly, i);
    p->d = i;
  }

  if (i <= p->d)
    p->c[i] = c;

  guppi_polynomial_reset_cache (poly);
  if (i == p->d)
    guppi_polynomial_sanitize (poly);
  guppi_polynomial_changed (poly);
}

double
guppi_polynomial_eval (GuppiPolynomial *poly, double x)
{
  GuppiPolynomialPrivate *p;
  double run;
  gint i;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), 0);
  p = priv (poly);

  i = p->d;
  run = p->c[i];
  --i;
  while (i >= 0) {
    run = run * x + p->c[i];
    --i;
  }
  return run;
}

double
guppi_polynomial_eval_D (GuppiPolynomial *poly, double x)
{
  GuppiPolynomialPrivate *p;
  double run;
  gint i;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), 0);
  p = priv (poly);

  if (p->d == 0)
    return 0;
  
  i = p->d;
  run = i * p->c[i];
  --i;
  while (i >= 1) {
    run = run * x + i * p->c[i];
    --i;
  }
  return run;
}

double
guppi_polynomial_eval_DD (GuppiPolynomial *poly, double x)
{
  GuppiPolynomialPrivate *p;
  double run;
  gint i;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), 0);
  p = priv (poly);

  if (p->d < 2)
    return 0;
  
  i = p->d;
  run = i * (i-1) * p->c[i];
  --i;
  while (i >= 2) {
    run = run * x + i * (i-1) * p->c[i];
    --i;
  }
  return run;
}

void
guppi_polynomial_sample (GuppiPolynomial *poly, gint N,
			 const double *x, gint x_stride,
			 double *y, gint y_stride)
{
  GuppiPolynomialPrivate *p;
  gint i, j;
  double run, x0;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  if (N == 0)
    return;
  g_return_if_fail (N > 0);
  g_return_if_fail (x != NULL);
  g_return_if_fail (y != NULL);

  p = priv (poly);

  for (j=0; j < N; ++j) {
    x0 = *x;

    i = p->d;
    run = p->c[i];
    --i;

    while (i >= 0) {
      run = run * x0 + p->c[i];
      --i;
    }

    *y = run;

    x  = (const double *)( ((gchar *)x) + x_stride );
    y  =       (double *)( ((gchar *)y) + y_stride );
  }
}

void
guppi_polynomial_sample_uniformly (GuppiPolynomial *poly,
				   double a, double b, gsize N,
				   double *x, gint x_stride,
				   double *y, gint y_stride)
{
  GuppiPolynomialPrivate *p;
  gint i, j;
  double run, x0;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  if (N == 0)
    return;
  g_return_if_fail (N > 0);

  if (x == NULL && y == NULL)
    return;

  p = priv (poly);
  
  for (j=0; j < N; ++j) {
    x0 = a + j*(b-a)/(N-1);

    i = p->d;
    run = p->c[i];
    --i;

    while (i >= 0) {
      run = run * x0 + p->c[i];
      --i;
    }

    if (x) {
      *x = x0;
      x  = (double *)( ((gchar *)x) + x_stride );
    }
    if (y) {
      *y = run;
      y  = (double *)( ((gchar *)y) + y_stride );
    }

  }
}

#define APPROX_PATH_MAX 1000
#define ROOT_BUCKET_MAX 200
ArtVpath *
guppi_polynomial_approximate_path (GuppiPolynomial *poly,
				   double a, double b,
				   double min_y, double max_y,
				   double x_error, double y_error,
				   double min_angle,
				   double scale_x, double scale_y)
{
  GuppiPolynomial *Dpoly;
  ArtVpath path_seg[APPROX_PATH_MAX], *path;
  double root_bucket[ROOT_BUCKET_MAX];
  gint j, p=0, pos, N, DN;
  GList *subdiv = NULL;
  gboolean done = FALSE;
  GList *i, *next_i;
  double midx, midy, x0, y0, x1, y1, dx, dy, tx, ty, normt;
  double rxa, rya, rxb, ryb, s, s0, s1, s2, s3, ms, Ms, m, M;
  double trx0, try0, trx1, try1; /* transformed coordinates */

  Dpoly = guppi_polynomial_copy (poly);
  guppi_polynomial_D (Dpoly);
  
  path_seg[p].x = a;
  path_seg[p].y = guppi_polynomial_eval (poly, a);
  subdiv = g_list_append (subdiv, &path_seg[p]);
  ++p;

  /*
    We always want to sample at the maxima and minima.  This keeps the plotted curve
    from "wobbling" when we translate it.
  */
  if (guppi_polynomial_degree (Dpoly) < ROOT_BUCKET_MAX) {
    DN = guppi_polynomial_find_real_roots (Dpoly, root_bucket);
    for (j=0; j<DN; ++j) {
      if (j > 0 && fabs (root_bucket[j] - root_bucket[j-1]) < POLYEPSILON)
	continue;

      if (a <= root_bucket[j] && root_bucket[j] <= b) {
	path_seg[p].x = root_bucket[j];
	path_seg[p].y = guppi_polynomial_eval (poly, root_bucket[j]);
	if (min_y <= path_seg[p].y && path_seg[p].y <= max_y) {
	  subdiv = g_list_append (subdiv, &path_seg[p]);
	  ++p;
	}
      }
    }
  }

  path_seg[p].x = b;
  path_seg[p].y = guppi_polynomial_eval (poly, b);
  subdiv = g_list_append (subdiv, &path_seg[p]);
  ++p;

  while (!done) {
    done = TRUE;

    pos = 0;
    i = subdiv;
    while ((next_i = g_list_next (i)) && pos < APPROX_PATH_MAX) {

      x0 = ((ArtVpath *) i->data)->x;
      y0 = ((ArtVpath *) i->data)->y;

      x1 = ((ArtVpath *) next_i->data)->x;
      y1 = ((ArtVpath *) next_i->data)->y;

      trx0 = scale_x * x0;
      try0 = scale_y * y0;

      trx1 = scale_x * x1;
      try1 = scale_y * y1;

      if (fabs (trx1 - trx0) > x_error) {
	
	guppi_polynomial_minmax_on_range (poly, x0, x1, &m, &M);

	if (min_y <= M && m <= max_y) {
	
	  midx = (trx0 + trx1)/2;
	  midy = (try0 + try1)/2;

	  dx = trx1 - trx0;
	  dy = try1 - try0;

	  tx = -dy;
	  ty =  dx;
	  normt = sqrt(tx*tx + ty*ty);
	  tx *= x_error / normt;
	  ty *= y_error / normt;

	  rxa = midx + tx;
	  rya = midy + ty;

	  rxb = midx - tx;
	  ryb = midy - ty;
      
	  s  = (try1 - try0) / (trx1 - trx0);
	  s0 = (rya - try0)  / (rxa - trx0);
	  s1 = (rya - try1)  / (rxa - trx1);
	  s2 = (ryb - try0)  / (rxb - trx0);
	  s3 = (ryb - try1)  / (rxb - trx1);
	  
	  ms = Ms = s;
	  if (s0 < ms) ms = s0;
	  if (Ms < s0) Ms = s0;
	  if (s1 < ms) ms = s1;
	  if (Ms < s1) Ms = s1;
	  if (s2 < ms) ms = s2;
	  if (Ms < s2) Ms = s2;
	  if (s3 < ms) ms = s3;
	  if (Ms < s3) Ms = s3;
	  
	  guppi_polynomial_minmax_on_range (Dpoly, x0, x1, &m, &M);

	  m *= scale_y / scale_x;
	  M *= scale_y / scale_x;
	  
	  if (m < ms || Ms < M) {
	    path_seg[p].x = (x0+x1)/2;
	    path_seg[p].y = guppi_polynomial_eval (poly, (x0+x1)/2);
	    ++pos;
	    subdiv = g_list_insert (subdiv, &path_seg[p], pos);
	    ++p;
	    done = FALSE;
	  }
	}
      }      
      i = next_i;
      ++pos;
    }
  }

  N = g_list_length (subdiv);
  path = g_new0 (ArtVpath, N+1);
  i = subdiv;
  j = 0;
  while (i != NULL && j < N) {
    memcpy (&path[j], i->data, sizeof (ArtVpath));
    path[j].code = j ? ART_LINETO : ART_MOVETO;
    i = g_list_next (i);
    ++j;
  }
  path[j].code = ART_END;
  g_print ("samples: %d\n", N);

  g_list_free (subdiv);
  guppi_unref0 (Dpoly);
  
  return path;
}

void
guppi_polynomial_scale (GuppiPolynomial *poly, double c)
{
  GuppiPolynomialPrivate *p;
  gint i;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  p = priv (poly);

  if (fabs (c-1) < POLYEPSILON)
    return;

  for (i = 0; i <= p->d; ++i)
    p->c[i] *= c;

  /* roots and minima/maxima are not affected by scaling. */
  guppi_polynomial_sanitize (poly);
  guppi_polynomial_changed (poly);
}

void
guppi_polynomial_D (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p;
  gint i;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  p = priv (poly);

  if (p->d == 0) {
    gboolean changed = fabs (p->c[0]) > POLYEPSILON;
    p->c[0] = 0;
    if (changed)
      goto done;
    return;
  }

  for (i = 1; i <= p->d; ++i) {
    p->c[i-1] = i * p->c[i];
  }
   
  --p->d;

 done:
  guppi_polynomial_reset_cache (poly);
  guppi_polynomial_changed (poly);
}

void
guppi_polynomial_normalize (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p;
  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  p = priv (poly);

  guppi_polynomial_scale (poly, 1 / p->c[p->d]);
}

void
guppi_polynomial_modulo (GuppiPolynomial *poly, GuppiPolynomial *mod)
{
  GuppiPolynomialPrivate *p;
  GuppiPolynomialPrivate *mod_p;
  double factor;
  gint i;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  g_return_if_fail (mod && GUPPI_IS_POLYNOMIAL (mod));

  p = priv (poly);
  mod_p = priv (mod);

  if (mod_p->d == 0) {
    p->c[0] = 0;
    p->d = 0;
  }

  if (p->d < mod_p->d)
    return;
  
  while (p->d >= mod_p->d) {
    factor = p->c[p->d] / mod_p->c[mod_p->d];

    if (fabs (factor) > POLYEPSILON) {
      for (i=0; i<=mod_p->d; ++i)
	p->c[p->d - mod_p->d + i] -= factor * mod_p->c[i];
    }

    --p->d;
    guppi_polynomial_sanitize (poly);
  }

  guppi_polynomial_reset_cache (poly);
  guppi_polynomial_changed (poly);
}

void
guppi_polynomial_deflate (GuppiPolynomial *poly, double x0)
{
  GuppiPolynomialPrivate *p;
  gint i;
  double old, t;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  p = priv (poly);

  if (p->d == 0)
    return;

  old = p->c[p->d-1];
  p->c[p->d-1] = p->c[p->d];
  for (i = p->d-2; i >= 0; --i) {
    t = p->c[i];
    p->c[i] = p->c[i+1] * x0 + old;
    old = t;
  }
  --p->d;

  guppi_polynomial_sanitize (poly);
  guppi_polynomial_reset_cache (poly);
  guppi_polynomial_changed (poly);
}

void
guppi_polynomial_deflate_complex (GuppiPolynomial *poly, double re, double im)
{
  GuppiPolynomialPrivate *p;
  gint i;
  double a, b, older, old, t;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  p = priv (poly);
  if (p->d < 2)
    return;

  /* We don't handle this case yet. */
  g_assert (p->d != 2);

  a = 2*re;
  b = -(re*re + im*im);

  older = p->c[p->d-2];
  old   = p->c[p->d-3];

  p->c[p->d-2] = p->c[p->d];
  p->c[p->d-3] = p->c[p->d-1] + a * p->c[p->d-2];

  for (i = p->d-4; i >= 0; --i) {
    t = p->c[i];

    p->c[i] = older + a * p->c[i+1] + b * p->c[i+2];
    
    older = old;
    old = t;
  }

  p->d -= 2;

  guppi_polynomial_sanitize (poly);
  guppi_polynomial_reset_cache (poly);
  guppi_polynomial_changed (poly);
}

void
guppi_polynomial_inflate (GuppiPolynomial *poly, double x0)
{
  GuppiPolynomialPrivate *p;
  gint i;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  p = priv (poly);

  guppi_polynomial_grow (poly, p->d+1);

  p->c[p->d+1] = p->c[p->d];
  for (i = p->d; i > 0; --i) {
    p->c[i] = p->c[i-1] - x0 * p->c[i];
  }
  p->c[0] *= -x0;

  ++p->d;

  guppi_polynomial_reset_cache (poly);
  guppi_polynomial_changed (poly);
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

/* There two functions are taken from GSL.  Thanks to Brian Gough. */

/* Solve a x^2 + b x + c == 0 */

static gint
solve_quadratic (double a, double b, double c,
		 double *x)
{
  double disc = b * b - 4 * a * c;

  if (disc > 0) {
      if (b == 0) {
          double r = fabs (0.5 * sqrt (disc) / a);
          x[0] = -r;
          x[1] =  r;
      } else {
          double sgnb = (b > 0 ? 1 : -1);
          double temp = -0.5 * (b + sgnb * sqrt (disc));
          double r1 = temp / a;
          double r2 = c / temp;

          if (r1 < r2) {
	    x[0] = r1;
	    x[1] = r2;
	  } else {
	    x[0] = r2;
	    x[1] = r1;
	  }
      }
      return 2;
  } else if (disc == 0) {
    x[0] = -0.5 * b / a;
    return 1;
  } else {
    return 0;
  }
}

#define GSL_SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0)

/* Solve x^3 + a x^2 + b x + c == 0. */

static gint
solve_cubic (double a, double b, double c,double *x)
{
  double q = (a * a - 3 * b);
  double r = (2 * a * a * a - 9 * a * b + 27 * c);

  double Q = q / 9;
  double R = r / 54;

  double Q3 = Q * Q * Q;
  double R2 = R * R;

  double CR2 = 729 * r * r;
  double CQ3 = 2916 * q * q * q;

  if (R == 0 && Q == 0) {

    x[0] = - a / 3;
    return 1;

  } else if (CR2 == CQ3) {
    /* this test is actually R2 == Q3, written in a form suitable
       for exact computation with integers */

    /* Due to finite precision some double roots may be missed, and
       considered to be a pair of complex roots z = x +/- epsilon i
       close to the real axis. */

    double sqrtQ = sqrt (Q);

    if (R > 0) {
      x[0] = -2 * sqrtQ  - a / 3;
      x[1] = sqrtQ - a / 3;
    } else {
      x[0] = - sqrtQ  - a / 3;
      x[1] = 2 * sqrtQ - a / 3;
    }
    return 2;
  
  } else if (CR2 < CQ3) { /* equivalent to R2 < Q3 */

    double sqrtQ = sqrt (Q);
    double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
    double theta = acos (R / sqrtQ3);
    double norm = -2 * sqrtQ;
    x[0] = norm * cos (theta / 3) - a / 3;
    x[1] = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3;
    x[2] = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3;
    
    /* Sort *x0, *x1, *x2 into increasing order */

    if (x[0] > x[1])
      GSL_SWAP(x[0], x[1]);
    
    if (x[1] > x[2]) {
      GSL_SWAP(x[1], x[2]);
	
      if (x[0] > x[1])
	GSL_SWAP(x[0], x[1]);
    }
    
    return 3;

  } else {

    double sgnR = (R >= 0 ? 1 : -1);
    double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0/3.0);
    double B = Q / A;
    x[0] = A + B - a / 3;
    return 1;

  }

}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

/*** Special handling for low-degree polynomials ***/

static void
process_deg0 (GuppiPolynomialPrivate *p)
{
  p->num_roots = 0;
  p->num_minmax = 0;
}

static void
process_deg1 (GuppiPolynomialPrivate *p)
{
  p->num_roots = 1;
  p->roots = guppi_new0(double, 1);
  p->roots[0] = -p->c[0] / p->c[1];

  p->num_minmax = 0;
}

static void
process_deg2 (GuppiPolynomialPrivate *p)
{
  p->roots = guppi_new0 (double, 2);
  p->num_roots = solve_quadratic (p->c[2], p->c[1], p->c[0], 
				       p->roots);

  p->num_minmax = 1;
  p->minmax = guppi_new0 (double, 1);
  p->minmax[0] = -p->c[0] / (2*p->c[1]);
}

static void
process_deg3 (GuppiPolynomialPrivate *p)
{
  double C = p->c[3];
  p->roots = guppi_new0 (double, 3);

  p->num_roots = solve_cubic (p->c[2] / C,
				   p->c[1] / C,
				   p->c[0] / C,
				   p->roots);

  p->minmax = guppi_new0 (double, 2);
  p->num_minmax = solve_quadratic (3 * p->c[3],
				   2 * p->c[1],
				   p->c[0],
				   p->minmax);
}

static void
process_deg4 (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p;
  GuppiPolynomial *q;
  GuppiPolynomialPrivate *q_p;
  gboolean have_root;
  double x0, C;

  p = priv (poly);

  /* Find the first root numerically, then deflate and solve the cubic. */
  
  have_root = guppi_polynomial_find_one_real_root (poly, &x0);
  if (have_root) {

    p->roots = guppi_new0 (double, 4);
    p->roots[0] = x0;
    p->num_roots = 1;

    q = guppi_polynomial_copy (poly);
    q_p = priv (q);
    guppi_polynomial_deflate (q, x0);

    C = priv (q)->c[3];
    p->num_roots += solve_cubic (q_p->c[2] / C,
				 q_p->c[1] / C,
				 q_p->c[0] / C,
				 p->roots+1);

    guppi_unref (q);

    guppi_polynomial_polish_cached_roots (poly);
    
  } else {

    p->num_roots = 0;

  }

  p->minmax = guppi_new0 (double, 3);
  C = 4 * p->c[4];
  p->num_minmax = solve_cubic (3 * p->c[3] /  C,
			       2 * p->c[2] /  C,
			       p->c[1] / C,
			       p->minmax);
}

static gboolean
process_switch (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p;
  gboolean rv = FALSE;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), FALSE);
  p = priv (poly);

  if (p->process_lock)
    return FALSE;

  p->process_lock = TRUE;

  switch (p->d) {

  case 0:
    process_deg0 (p);
    rv = TRUE;
    break;

  case 1:
    process_deg1 (p);
    rv = TRUE;
    break;

  case 2:
    process_deg2 (p);
    rv = TRUE;
    break;

  case 3:
    process_deg3 (p);
    rv = TRUE;
    break;

  case 4:
    process_deg4 (poly);
    rv = TRUE;
    break;
  }

  p->process_lock = FALSE;

  return rv;
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

double
guppi_polynomial_newton_polish (GuppiPolynomial *poly, double x, gint iter, double epsilon)
{
  GuppiPolynomialPrivate *p;
  double best_x, best_abs_y;
  double fx, Dfx, abs_fx;
  gint i;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), x);
  p = priv (poly);

  if (epsilon <= 0)
    epsilon = POLYEPSILON;

  best_x = x;
  best_abs_y = fabs (guppi_polynomial_eval (poly, best_x));

  while (iter > 0 && best_abs_y > epsilon) {
    
    /* We simultaneously evaluate the polynomial
       and its first derivative. */
    i = p->d;
    fx = p->c[i];
    Dfx = i * p->c[i];
    --i;
    while (i >= 1) {
      fx = fx * x + p->c[i];
      Dfx = Dfx * x + i * p->c[i];
      --i;
    }
    fx = fx * x + p->c[0];

    abs_fx = fabs (fx);

    if (abs_fx < best_abs_y) {
      best_abs_y = abs_fx;
      best_x = x;
    }

    if (fabs (Dfx) < epsilon)
      break;

    x -= fx / Dfx;

    --iter;
  }

  return best_x;
}

double
guppi_polynomial_gershgorin_radius (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p;
  double max_c, top_c, q;
  gint i;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), 0);
  p = priv (poly);
  
  if (p->d == 0)
    return 0;

  max_c = 0;
  top_c = p->c[p->d];
  for (i = 0; i < p->d; ++i) {
    q = fabs (p->c[i] / top_c);
    if (q > max_c)
      max_c = q;
  }
  
  return max_c + 1;
}

gint
guppi_polynomial_real_roots (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p;
  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), -1);
  p = priv (poly);
  
  if (p->num_roots < 0) {

    if (! process_switch (poly)) {
      double R = guppi_polynomial_gershgorin_radius (poly);
      p->num_roots = guppi_polynomial_real_roots_in_interval (poly, -R, R);
    }

  }

  return p->num_roots;
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

static gint
sign_changes (GuppiPolynomial **sturm_seq, double a)
{
  gint i, count=0;
  double x, prev_x=0;
  gboolean started = FALSE;

  for (i = 0; sturm_seq[i]; ++i) {
    x = guppi_polynomial_eval (sturm_seq[i], a);
    if (fabs (x) > POLYEPSILON) {
      if (started && x * prev_x < 0)
	++count;
      prev_x = x;
      started = TRUE;
    }
  }

  return count;
}

static GuppiPolynomial **
build_sturm_sequence (GuppiPolynomial *poly)
{
  GuppiPolynomial **sturm_seq;
  gint i;

  sturm_seq = guppi_new0 (GuppiPolynomial *, guppi_polynomial_degree (poly)+2);

  sturm_seq[0] = guppi_polynomial_copy (poly);

  if (guppi_polynomial_degree (poly) <= 1)
    return sturm_seq;

  sturm_seq[1] = guppi_polynomial_copy (poly);
  guppi_polynomial_D (sturm_seq[1]);

  for (i=2; guppi_polynomial_degree (sturm_seq[i-1]) > 0; ++i) {
    sturm_seq[i] = guppi_polynomial_copy (sturm_seq[i-2]);
    guppi_polynomial_modulo (sturm_seq[i], sturm_seq[i-1]);
    guppi_polynomial_scale (sturm_seq[i], -1);
  }

  return sturm_seq;
}

static void
free_sturm_sequence (GuppiPolynomial **seq)
{
  gint i;
  
  if (seq) {
    for (i = 0; seq[i]; ++i) 
      guppi_unref (seq[i]);
  }

  guppi_free (seq);
}

gint
guppi_polynomial_real_roots_in_interval (GuppiPolynomial *poly, double a, double b)
{
  GuppiPolynomialPrivate *p;
  gint Fa, Fb;
  GuppiPolynomial **sturm_seq;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), -1);
  p = priv (poly);

  if (p->roots == NULL)
    process_switch (poly);

  if (p->roots != NULL) {
    gint i, count=0;
    for (i = 0; i < p->num_roots; ++i) {
      if (a <= p->roots[i] && p->roots[i] <= b)
	++count;
    }
    return count;
  }

  sturm_seq = build_sturm_sequence (poly);

  Fa = sign_changes (sturm_seq, a);
  Fb = sign_changes (sturm_seq, b);

  free_sturm_sequence (sturm_seq);

  return abs (Fa - Fb);
}

gboolean
guppi_polynomial_find_one_real_root (GuppiPolynomial *poly, double *root)
{
  GuppiPolynomialPrivate *p;
  GuppiPolynomial **sturm_seq = NULL;
  gboolean found = FALSE, bracketed = FALSE;
  double epsilon_width = 1e-3;
  double epsilon_error;
  double a, b, x=0, y, ya=0, yb=0;
  gint Fa, Fb, Fx;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), FALSE);
  p = priv (poly);

  /* Handle special cases. */

  if (p->num_roots < 0)
    process_switch (poly);

  if (p->num_roots == 0) {
    return FALSE;
  }

  if (p->roots) {
    if (root)
      *root = p->roots[0];
    return TRUE;
  }

  /* Do things the hard way. */

  sturm_seq = build_sturm_sequence (poly);
  b = guppi_polynomial_gershgorin_radius (poly);
  a = -b;

  epsilon_error = MIN (b/100, POLYEPSILON);

  Fa = sign_changes (sturm_seq, a);
  Fb = sign_changes (sturm_seq, b);

  if (Fa == Fb)
    goto done;

  /* Yes, there is a root to be found... */  
  found = TRUE;

  if (root == NULL)
    return TRUE;

  while (b-a > epsilon_width) {

    x = (a+b)/2;
    y = guppi_polynomial_eval (poly, x);
    if (fabs (y) < epsilon_error) {
      *root = x;
      goto done;
    }

    Fx = sign_changes (sturm_seq, x);

    if (Fa != Fx) {
      b = x;
      Fb = Fx;
    } else if (Fb != Fx) {
      a = x;
      Fa = Fx;
    } else {
      g_assert_not_reached ();
    }

    /* Check to see if we've managed to bracket a root with different
       signs on either side. */
    if (abs (Fa - Fb) == 1) {
      ya = guppi_polynomial_eval (poly, a);
      yb = guppi_polynomial_eval (poly, b);
      if (ya * yb < 0) {
	bracketed = TRUE;
	break;
      }
    }

  }

  /* Finish up by subdividing to find our bracketed root. */
  if (bracketed) {
    while (b - a > epsilon_width) {

      x = (a+b)/2;
      y = guppi_polynomial_eval (poly, x);
      if (fabs (y) < epsilon_error) {
	*root = y;
	goto done;
      }
      
      if (ya * y > 0) {
	a = x;
	ya = y;
      } else {
	b = x;
	yb = y;
      }
    }
  }

  /* Polish up our root */
  *root = guppi_polynomial_newton_polish (poly, x, 10, epsilon_error);

 done:
  free_sturm_sequence (sturm_seq);
  return found;
}

gint
guppi_polynomial_find_real_roots (GuppiPolynomial *poly, double *roots)
{
  GuppiPolynomialPrivate *p;
  GuppiPolynomial *q;
  GuppiPolynomialPrivate *q_p;
  double x0;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), -1);
  p = priv (poly);

  if (roots == NULL)
    return guppi_polynomial_real_roots (poly);

  if (p->num_roots < 0 || (p->num_roots > 0 && p->roots == NULL))
    process_switch (poly);

  if (p->num_roots == 0) {
    return 0;
  } else if (p->num_roots > 0 && p->roots) {
    memcpy (roots, p->roots, sizeof (double) * p->num_roots);
    return p->num_roots;
  }

  q = guppi_polynomial_copy (poly);
  q_p = priv (q);

  p->num_roots = 0;
  while (guppi_polynomial_find_one_real_root (q, &x0)) {

    /* If we've got all of the roots, grab them all rather than limiting ourselves to
       just one. */
    if (q_p->num_roots == 0 || (q_p->num_roots > 0 && q_p->roots)) {
      memcpy (roots+p->num_roots, q_p->roots, sizeof (double) * q_p->num_roots);
      p->num_roots += q_p->num_roots;
      break;
    }

    roots[p->num_roots] = x0;
    ++p->num_roots;
    guppi_polynomial_deflate (q, x0);
  }
  guppi_unref (q);

  /* Cache root information. */
  if (p->num_roots > 0 && p->roots == NULL) {
    p->roots = guppi_new (double, p->num_roots);
    memcpy (p->roots, roots, sizeof (double) * p->num_roots);
  }

  /* Make sure our roots are all nice and shiny, then copy the polished roots back out
     to the user's array. */
  if (p->num_roots > 0) {
    guppi_polynomial_polish_cached_roots (poly);
    memcpy (roots, p->roots, sizeof (double) * p->num_roots);
  }

  /* sort our roots */
  if (p->roots)
    qsort (p->roots, p->num_roots, sizeof (double), double_cmp);

  return p->num_roots;
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

static void
guppi_polynomial_cache_minmax (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p = priv (poly);
  GuppiPolynomial *Dp;
  
  if (p->num_minmax == 0 || (p->num_minmax > 0 && p->minmax != NULL))
    return;
  
  Dp = guppi_polynomial_copy (poly);
  guppi_polynomial_D (Dp);

  p->minmax = guppi_new (double, p->d - 1);
  p->num_minmax = guppi_polynomial_find_real_roots (Dp, p->minmax);

  guppi_unref (Dp);
}

void
guppi_polynomial_minmax_on_range (GuppiPolynomial *poly, double a, double b,
				  double *min, double *max)
{
  GuppiPolynomialPrivate *p;
  double va, vb, m, M;
  gint i;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  p = priv (poly);

  /* Optimizations for low degees. */

  if (p->d == 0) {
    if (min) *min = p->c[0];
    if (max) *max = p->c[0];
    return;
  }

  if (p->d == 1) {
    double v1 = p->c[0] + a * p->c[1];
    double v2 = p->c[0] + b * p->c[1];
    if (min) *min = MIN (v1, v2);
    if (max) *max = MAX (v1, v2);
    return;
  }

  /* The general case. */

  guppi_polynomial_cache_minmax (poly);

  va = guppi_polynomial_eval (poly, a);
  vb = guppi_polynomial_eval (poly, b);

  m = MIN (va, vb);
  M = MAX (va, vb);

  for (i=0; i<p->num_minmax; ++i) {
    double x = p->minmax[i], y;
    if (a < x && x < b) {
      y = guppi_polynomial_eval (poly, x);
      if (y < m)
	m = y;
      if (M < y)
	M = y;
    }
  }

  if (min)
    *min = m;
  if (max)
    *max = M;
}

gboolean
guppi_polynomial_find_bounded_range (GuppiPolynomial *poly, double x, double min, double max,
				     double *a, double *b)
{
  GuppiPolynomialPrivate *p;
  double y;
  gint i, mm_i;
  const gint max_iter = 10;

  double a0, a1, b0, b1;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), FALSE);
  p = priv (poly);

  if (guppi_polynomial_degree (poly) == 0)
    return FALSE;

  y = guppi_polynomial_eval (poly, x);
  if (y < min || y > max)
    return FALSE;

  guppi_polynomial_cache_minmax (poly);

  mm_i = 0;
  while (mm_i < p->num_minmax && x <= p->minmax[mm_i]) {
      ++mm_i;
  }

  a0 = a1 = b0 = b1 = x;

  if (a) {
    for (i = mm_i; i > 0; --i) {
      a1 = a0;
      a0 = p->minmax[i];
      y = guppi_polynomial_eval (poly, a0);
      if (y < min || max < y)
	break;
    }

    if (i <= 0) {
      double q = 1;
      while (1) {
	a1 = a0;
	a0 = a0 - q;
	q *= 2;
	y = guppi_polynomial_eval (poly, a0);
	if (y < min || max < y)
	  break;
      }
    }

    for (i=0; i<max_iter && a1 - a0 > POLYEPSILON; ++i) {
      double mid_x = (a0+a1)/2;
      y = guppi_polynomial_eval (poly, mid_x);

      if (min <= y && y <= max) {
	a1 = mid_x;
      } else {
	a0 = mid_x;
      }
    }

    *a = a1;
  }

  if (b) {

    for (i = mm_i+1; i < p->num_minmax; ++i) {
      b0 = b1;
      b1 = p->minmax[i];
      y = guppi_polynomial_eval (poly, b1);
      if (y < min || max < y)
	break;
    }

    if (i == p->num_minmax) {
      double q = 1;
      while (1) {
	b0 = b1;
	b1 = b1 + q;
	q *= 2;
	y = guppi_polynomial_eval (poly, b1);
	if (y < min || max < y)
	  break;
      }
    }

    for (i=0; i<max_iter && b1 - b0 > POLYEPSILON; ++i) {
      double mid_x = (b0+b1)/2;
      y = guppi_polynomial_eval (poly, mid_x);
      
      if (min <= y && y <= max) {
	b0 = mid_x;
      } else {
	b1 = mid_x;
      }
    }

    *b = b0;

  }

  return TRUE;
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

void
guppi_polynomial_spew (GuppiPolynomial *poly)
{
  GuppiPolynomialPrivate *p;
  gint i;

  g_return_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly));
  p = priv (poly);

  for (i = p->d; i >= 0; --i) {

    if (fabs (p->c[i]) > POLYEPSILON) {
      
      if (fabs (p->c[i] - 1) > POLYEPSILON || i == 0)
	g_print ("%g ", p->c[i]);

      if (i) {
	if (i == 1)
	  g_print ("x + ");
	else
	  g_print ("x^%d + ", i);
      }
    }
  }

  g_print ("\n");
}

xmlNodePtr
guppi_polynomial_export_xml (GuppiPolynomial *poly, GuppiXMLDocument *doc)
{
  GuppiPolynomialPrivate *p;
  xmlNodePtr poly_node;
  xmlNodePtr term_node;
  gint i;

  g_return_val_if_fail (poly && GUPPI_IS_POLYNOMIAL (poly), NULL);
  g_return_val_if_fail (doc != NULL, NULL);

  p = priv (poly);

  poly_node = guppi_xml_new_node (doc, "Polynomial");
  guppi_xml_set_propertyf (poly_node, "degree", "%d", p->d);
  if (p->num_roots >= 0)
    guppi_xml_set_propertyf (poly_node, "roots", "%d", p->num_roots);
  if (p->num_minmax >= 0)
    guppi_xml_set_propertyf (poly_node, "minmax", "%d", p->num_minmax);

  for (i = 0; i <= p->d; ++i) {
    if (fabs (p->c[i]) >= POLYEPSILON) {
      term_node = guppi_xml_new_text_nodef (doc, "term", "%g", p->c[i]);
      guppi_xml_set_propertyf (term_node, "degree", "%d", i);
      xmlAddChild (poly_node, term_node);
    }
  }

  for (i = 0; i < p->num_roots; ++i) {
    term_node = guppi_xml_new_text_nodef (doc, "root", "%g", p->roots[i]);
    xmlAddChild (poly_node, term_node);
  }

  for (i = 0; i < p->num_minmax; ++i) {
    term_node = guppi_xml_new_text_nodef (doc, "minmax", "%g", p->minmax[i]);
    xmlAddChild (poly_node, term_node);
  }

  return poly_node;
}

GuppiPolynomial *
guppi_polynomial_import_xml (GuppiXMLDocument *doc, xmlNodePtr node)
{
  GuppiPolynomial *poly;
  GuppiPolynomialPrivate *p;
  gint i, deg;
  gchar *buf;
  
  g_return_val_if_fail (doc != NULL, NULL);
  g_return_val_if_fail (node != NULL, NULL);

  if (strcmp (node->name, "Polynomial"))
    return NULL;

  buf = xmlGetProp (node, "degree");
  deg = buf ? atoi (buf) : 0;
  free (buf);

  poly = gtk_type_new (GUPPI_TYPE_POLYNOMIAL);
  p = priv (poly);

  guppi_polynomial_freeze (poly);

  guppi_polynomial_grow (poly, deg);

  node = node->xmlChildrenNode;

  while (node) {

    if (!strcmp (node->name, "term")) {

      buf = xmlGetProp (node, "degree");
      i = buf ? atoi (buf) : 0;
      free (buf);

      buf = xmlNodeListGetString (doc->doc, node->xmlChildrenNode, 1);
      if (0 <= i && i <= deg)
	p->c[i] = atof (buf);
      free (buf);
    }
    
    node = node->next;
  }

  guppi_polynomial_thaw (poly);

  return poly;
}