version 1.1.1.1, 2000/01/09 17:00:52 |
version 1.1.1.2, 2000/01/22 14:15:59 |
Line 251 struct curve_points *plot; |
|
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 |
* cp_binomial() computes the binomial coefficients needed for BEZIER stuff |
* and stores them into an array which is hooked to sdat. |
* and stores them into an array which is hooked to sdat. |
|
|
|
|
n = points - 1; |
n = points - 1; |
e = n / 2; |
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++) { |
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--) |
for (k = n; k >= e; k--) |
Line 322 unsigned int exponent; |
|
Line 333 unsigned int exponent; |
|
return (y); |
return (y); |
} |
} |
|
|
|
|
static void eval_bezier(cp, first_point, num_points, sr, px, py, c) |
static void eval_bezier(cp, first_point, num_points, sr, px, py, c) |
struct curve_points *cp; |
struct curve_points *cp; |
int first_point; /* where to start in plot->points (to find x-range) */ |
int first_point; /* where to start in plot->points (to find x-range) */ |
|
|
*px = this_points[n].x; |
*px = this_points[n].x; |
*py = this_points[n].y; |
*py = this_points[n].y; |
} else { |
} 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; |
unsigned int i; |
double lx = 0.0, ly = 0.0; |
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++) { |
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; |
lx += this_points[i].x * u; |
ly += this_points[i].y * u; |
ly += this_points[i].y * u; |
sr_to_the_i *= sr; |
|
dsr_to_the_di *= reciprocal_dsr; |
|
} |
} |
|
|
*px = lx; |
*px = lx; |