/* 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 * * 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 #include "guppi-polynomial.h" #include #include #include #include #include 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; ic[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; ic[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= 0, NULL); poly = guppi_polynomial_new_constant (1); guppi_polynomial_freeze (poly); guppi_polynomial_grow (poly, degree); for (i=0; id, 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 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; inum_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 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 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; }