[Scipy-tickets] [SciPy] #840: Potential subtractive cancellation bug in scipy/interpolate/Rbf.py

SciPy scipy-tickets@scipy....
Fri Jan 9 00:14:24 CST 2009


#840: Potential subtractive cancellation bug in scipy/interpolate/Rbf.py
-------------------------------+--------------------------------------------
 Reporter:  fghorow            |       Owner:  somebody
     Type:  defect             |      Status:  new     
 Priority:  normal             |   Milestone:  0.7.1   
Component:  scipy.interpolate  |     Version:  devel   
 Severity:  normal             |    Keywords:  Rbf     
-------------------------------+--------------------------------------------
 In the function `_euclidean_norm()` of Rbf.py (line 100 in the current
 SVN) there is a one liner:

 `        return sqrt( ((x1 - x2)**2).sum(axis=0) )`

 A problem can arise when both `x1` and `x2` are single precision, since
 the difference will then be computed in single precision, and the usual
 "subtractive cancellation" problem of numerical analysis can occur. (Yes,
 I was just bitten by this bug by a slight local modification of this code,
 where I was trying to compute in-place to save some memory.)

 The fix is simple. Both `x1` and `x2` need to be cast into double
 precision before subtraction.

 In the old Numeric, I would simply have coded `x1.astype(Numeric.Float64)`
 etc. and moved on. Here in the numpy environment, I'm unsure of the side
 effects, and think that using `asarray()` ''might'' be preferable.

 I'll leave the coding of the actual fix to someone more familiar with the
 `numpy` intricacies than I am.

-- 
Ticket URL: <http://scipy.org/scipy/scipy/ticket/840>
SciPy <http://www.scipy.org/>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list