[Gnuplot-bugs] gnuplot-4.2.3 (and earlier): serious floating-point bug, and a simple fix

[Nelson H. F. B. <be...@ma...>]

I finally tracked down a problem that has bothered me for some time.
It can be exhibited by a one-line test in ANY current version of
gnuplot (4.3.2 or earlier):

	gnuplot> print log10(1.0e-162)
		       ^
		 undefined value

This slightly larger argument works:

	gnuplot> print log10(1.0e-161)
	-161.002592673726

On impact of this bug is that it is impossible to graph a function
result over a logarithmic representation of the entire floating-point
range with, e.g.,

	plot "foo.dat" using (log10(abs($1))):($2)

The fix is, fortunately, easy:

diff -c -r gnuplot-4.2.3/src/eval.c gnuplot-4.2.3.1/src/eval.c
*** gnuplot-4.2.3/src/eval.c	Wed Mar 28 16:23:10 2007
--- gnuplot-4.2.3.1/src/eval.c	Thu Aug 28 09:27:44 2008
***************
*** 323,332 ****
--- 323,337 ----
      case INTGR:
  	return ((double) abs(val->v.int_val));
      case CMPLX:
+ #if 0
  	return (sqrt(val->v.cmplx_val.real *
  		     val->v.cmplx_val.real +
  		     val->v.cmplx_val.imag *
  		     val->v.cmplx_val.imag));
+ #else
+ 	return (hypot(val->v.cmplx_val.real,
+ 		     val->v.cmplx_val.imag));
+ #endif
      default:
  	int_error(NO_CARET, "unknown type in magnitude()");
      }

The correct way to compute a Euclidean norm is NOT sqrt(x*x + y*y),
but hypot(x,y).  The first slashes the floating-point range in half
because of premature overflow or underflow, but the second is valid
for all possible argument pairs.

As the diff line suggests, I made further changes in the configure.in
file's AC_INIT() line to change the version locally from 4.2.3 to
4.2.3.1.  However, versioning is the developers' job, so 4.2.3.1 is
unique to my systems.

I scanned for other instances of incorrect use of sqrt() instead of
correct hypot(), but did not fix them.  Many of them SHOULD be changed
too, but I leave that job to the developers for the next official
release of gnuplot:

% find . -name '*.[ch]' | xargs grep -n 'sqrt.*[*]'
./src/fit.c:345:            *lambda = sqrt(*lambda / num_data / num_params);
./src/os2/gclient.c:3833:        ds = sqrt(dx*dx + dy*dy);
./src/os2/gclient.c:3874:        ds = sqrt(dx*dx + dy*dy);
./src/datafile.c:2778:    C1 = sqrt(x*x + y*y + z*z);
./src/datafile.c:2779:    C2 = sqrt(x*x + y*y);
./src/eval.c:327:       return (sqrt(val->v.cmplx_val.real *
./src/term.c:1165:    double len_arrow = sqrt(dx * dx + dy * dy);
./src/specfun.c:347:    /* log( sqrt( 2*pi ) ) */
./src/specfun.c:377:    /* log( sqrt( 2*pi ) ) */
./src/specfun.c:408:    /* log( sqrt( 2*pi ) ) */
./src/specfun.c:440:    /* log( sqrt( 2*pi ) ) */
./src/specfun.c:1039: * z = sqrt( -2.0 * log(y) );  then the approximation is
./src/specfun.c:1043: * w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
./src/specfun.c:1070:/* sqrt(2pi) */
./src/specfun.c:1401:    x = sqrt(-2.0 * log(y));
./src/specfun.c:1859:    x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
./src/specfun.c:1860:    x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
./src/specfun.c:1861:    x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
./src/specfun.c:1890:   p = sqrt(2.0 * (exp(1.0) * x + 1.0));
./src/matrix.c:156:                 w = fsign(C[j][j]) * sqrt(C[j][j] * C[j][j] + C[i][j] * C[i][j]);
./src/matrix.c:210:             gamma = sqrt(1 - sigma * sigma);
./src/matrix.c:213:             sigma = fsign(rho) * sqrt(1 - gamma * gamma);
./src/standard.h:73:void f_sqrt __PROTO((union argument *x));
./src/hidden3d.c:631:    s = sqrt(plane[0] * plane[0] + plane[1] * plane[1] + plane[2] * plane[2]);
./src/hidden3d.c:656:   s = sqrt(plane[0] * plane[0] + plane[1] * plane[1] + plane[2] * plane[2]);
./src/pm3d.c:100:    /* pow(x, 0.25) is slightly faster than sqrt(sqrt(x)) */
./src/standard.c:585:   push(Gcomplex(&a, 0.0, -log(-y + sqrt(y * y + 1)) / ang2rad));
./src/standard.c:587:   beta = sqrt((x + 1) * (x + 1) + y * y) / 2 - sqrt((x - 1) * (x - 1) + y * y) / 2;
./src/standard.c:590:   alpha = sqrt((x + 1) * (x + 1) + y * y) / 2 + sqrt((x - 1) * (x - 1) + y * y) / 2;
./src/standard.c:591:   push(Gcomplex(&a, asin(beta) / ang2rad, -log(alpha + sqrt(alpha * alpha - 1)) / ang2rad));
./src/standard.c:609:   double alpha = sqrt((x + 1) * (x + 1) + y * y) / 2
./src/standard.c:610:                  + sqrt((x - 1) * (x - 1) + y * y) / 2;
./src/standard.c:611:   double beta = sqrt((x + 1) * (x + 1) + y * y) / 2
./src/standard.c:612:                 - sqrt((x - 1) * (x - 1) + y * y) / 2;
./src/standard.c:618:                     log(alpha + sqrt(alpha * alpha - 1)) / ang2rad));
./src/standard.c:749:   push(Gcomplex(&a, log(y + sqrt(y * y + 1)) / ang2rad, 0.0));
./src/standard.c:751:   beta = sqrt((x + 1) * (x + 1) + y * y) / 2 - sqrt((x - 1) * (x - 1) + y * y) / 2;
./src/standard.c:752:   alpha = sqrt((x + 1) * (x + 1) + y * y) / 2 + sqrt((x - 1) * (x - 1) + y * y) / 2;
./src/standard.c:753:   push(Gcomplex(&a, log(alpha + sqrt(alpha * alpha - 1)) / ang2rad, -asin(beta) / ang2rad));
./src/standard.c:770:   push(Gcomplex(&a, log(x + sqrt(x * x - 1)) / ang2rad, 0.0));
./src/standard.c:772:   alpha = sqrt((x + 1) * (x + 1) + y * y) / 2
./src/standard.c:773:           + sqrt((x - 1) * (x - 1) + y * y) / 2;
./src/standard.c:774:   beta = sqrt((x + 1) * (x + 1) + y * y) / 2
./src/standard.c:775:          - sqrt((x - 1) * (x - 1) + y * y) / 2;
./src/standard.c:776:   push(Gcomplex(&a, log(alpha + sqrt(alpha * alpha - 1)) / ang2rad,
./src/standard.c:868:f_sqrt(union argument *arg)
./src/standard.c:882:   /* -pi < ang < pi, so real(sqrt(z)) >= 0 */
./src/standard.c:1072:  return (sqrt(TWO_ON_PI / x) *
./src/standard.c:1085:  return (sqrt(TWO_ON_PI / x) *
./src/standard.c:1168:  w = sqrt(TWO_ON_PI / x) *
./src/standard.c:1185:  return (sqrt(TWO_ON_PI / x) *
./src/mouse.c:586:      rho = sqrt((x - rx) * (x - rx) + (y - ry) * (y - ry)); /* distance */
./src/mouse.c:1522:         int dist = sqrt((double)(dist_x * dist_x + dist_y * dist_y));
./src/graphics.c:3239:                  x1 = xlowM + bar_size * tic / sqrt(1.0 + slope * slope);
./src/graphics.c:3240:                  x2 = xlowM - bar_size * tic / sqrt(1.0 + slope * slope);
./src/graphics.c:3248:                  x1 = xhighM + bar_size * tic / sqrt(1.0 + slope * slope);
./src/graphics.c:3249:                  x2 = xhighM - bar_size * tic / sqrt(1.0 + slope * slope);
./src/plot3d.c:403:         K[i][j] = K[j][i] = -splines_kernel(sqrt(dx * dx + dy * dy));
./src/plot3d.c:453:             z = z - b[k] * splines_kernel(sqrt(dx * dx + dy * dy));

-------------------------------------------------------------------------------
- Nelson H. F. Beebe                    Tel: +1 801 581 5254                  -
- University of Utah                    FAX: +1 801 581 4148                  -
- Department of Mathematics, 110 LCB    Internet e-mail: be...@ma...  -
- 155 S 1400 E RM 233                       be...@ac...  be...@co... -
- Salt Lake City, UT 84112-0090, USA    URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------