===================================================================
RCS file: /home/cvs/OpenXM_contrib/gnuplot/Attic/interpol.c,v
retrieving revision 1.1.1.1
retrieving revision 1.1.1.3
diff -u -p -r1.1.1.1 -r1.1.1.3
--- OpenXM_contrib/gnuplot/Attic/interpol.c	2000/01/09 17:00:52	1.1.1.1
+++ OpenXM_contrib/gnuplot/Attic/interpol.c	2003/09/15 07:09:25	1.1.1.3
@@ -1,5 +1,5 @@
 #ifndef lint
-static char *RCSid = "$Id: interpol.c,v 1.1.1.1 2000/01/09 17:00:52 maekawa Exp $";
+static char *RCSid = "$Id: interpol.c,v 1.1.1.3 2003/09/15 07:09:25 ohara Exp $";
 #endif
 
 /* GNUPLOT - interpol.c */
@@ -188,7 +188,7 @@ 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));
 
 
 /*
@@ -251,6 +251,15 @@ 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.
@@ -268,10 +277,12 @@ 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--)
@@ -322,6 +333,7 @@ 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) */
@@ -344,18 +356,20 @@ 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 dsr_to_the_di = s_pow(1 - sr, n);
-	double reciprocal_dsr = 1.0 / (1 - sr);
+	double log_dsr_to_the_n = n * log(1 -sr);
+	double log_sr_over_dsr = log(sr) - log(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;
@@ -909,10 +923,13 @@ struct curve_points *plot;
  * (MGR 1992)
  */
 
-static int compare_points(p1, p2)
-struct coordinate *p1;
-struct coordinate *p2;
+int compare_points(argp1, argp2)
+    SORTFUNC_ARGS argp1;
+    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)
@@ -921,7 +938,7 @@ struct coordinate *p2;
 }
 
 void sort_points(plot)
-struct curve_points *plot;
+    struct curve_points *plot;
 {
     int first_point, num_points;
 
@@ -929,7 +946,7 @@ 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;