[BACK]Return to interpol.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gnuplot

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

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.
Line 268  int points;
Line 277  int points;
   
     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) */
Line 344  double *c;
Line 356  double *c;
         *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;
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.1.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>