OpenXM_contrib/gnuplot/interpol.c - diff

Return to interpol.c CVS log

Up to [local] / OpenXM_contrib / gnuplot

Diff for /OpenXM_contrib/gnuplot/Attic/interpol.c between version 1.1.1.1 and 1.1.1.3

version 1.1.1.1, 2000/01/09 17:00:52

version 1.1.1.3, 2003/09/15 07:09:25

Line 188 static int solve_five_diag __PROTO((five_diag m[], dou

static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points));

static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points));

static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest));

static int compare_points __PROTO((struct coordinate * p1, struct coordinate * p2));

int compare_points __PROTO((SORTFUNC_ARGS p1, SORTFUNC_ARGS p2));

Line 251 struct curve_points *plot;

/* HBB 990205: rewrote the 'bezier' interpolation routine,

* to prevent numerical overflow and other undesirable things happening

* for large data files (num_data about 1000 or so), where binomial

* coefficients would explode, and powers of 'sr' (0 < sr < 1) become

* extremely small. Method used: compute logarithms of these

* extremely large and small numbers, and only go back to the

* real numbers once they've cancelled out each other, leaving

* a reasonable-sized one. */

* cp_binomial() computes the binomial coefficients needed for BEZIER stuff

* and stores them into an array which is hooked to sdat.

Line 268 int points;

Line 277 int points;

n = points - 1;

e = n / 2;

coeff[0] = 1.0;

/* HBB 990205: calculate these in 'logarithmic space',

* as they become _very_ large, with growing n (4^n) */

coeff[0] = 0.0;

for (k = 0; k < e; k++) {

coeff[k + 1] = coeff[k] * ((double) (n - k)) / ((double) (k + 1));

coeff[k + 1] = coeff[k] + log(((double) (n-k)) / ((double) (k+1)));

}

for (k = n; k >= e; k--)

Line 322 unsigned int exponent;

Line 333 unsigned int exponent;

return (y);

}

static void eval_bezier(cp, first_point, num_points, sr, px, py, c)

struct curve_points *cp;

int first_point; /* where to start in plot->points (to find x-range) */

Line 344 double *c;

Line 356 double *c;

*px = this_points[n].x;

*py = this_points[n].y;

} else {

/* HBB 990205: do calculation in 'logarithmic space',

* to avoid over/underflow errors, which would exactly cancel

* out each other, anyway, in an exact calculation

unsigned int i;

double lx = 0.0, ly = 0.0;

double sr_to_the_i = 1.0;

double log_dsr_to_the_n = n * log(1 -sr);

double dsr_to_the_di = s_pow(1 - sr, n);

double log_sr_over_dsr = log(sr) - log(1 - sr);

double reciprocal_dsr = 1.0 / (1 - sr);

for (i = 0; i <= n; i++) {

double u = c[i] * sr_to_the_i * s_pow(1 - sr, n - i);

double u = exp(c[i] + log_dsr_to_the_n + i * log_sr_over_dsr);

lx += this_points[i].x * u;

ly += this_points[i].y * u;

sr_to_the_i *= sr;

dsr_to_the_di *= reciprocal_dsr;

}

*px = lx;

Line 909 struct curve_points *plot;

Line 923 struct curve_points *plot;

* (MGR 1992)

static int compare_points(p1, p2)

int compare_points(argp1, argp2)

struct coordinate *p1;

SORTFUNC_ARGS argp1;

struct coordinate *p2;

SORTFUNC_ARGS argp2;

{

const struct coordinate *p1 = argp1;

const struct coordinate *p2 = argp2;

if (p1->x > p2->x)

return (1);

if (p1->x < p2->x)

Line 921 struct coordinate *p2;

Line 938 struct coordinate *p2;

}

void sort_points(plot)

struct curve_points *plot;

{

int first_point, num_points;

Line 929 struct curve_points *plot;

Line 946 struct curve_points *plot;

while ((num_points = next_curve(plot, &first_point)) > 0) {

/* Sort this set of points, does qsort handle 1 point correctly? */

qsort((char *) (plot->points + first_point), num_points,

sizeof(struct coordinate), (sortfunc) compare_points);

sizeof(struct coordinate), compare_points);

first_point += num_points;

}

return;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>