[Numpy-discussion] double-precision sqrt?

[Charles R Harris]

On Sat, Oct 17, 2009 at 12:08 PM, Adam Ginsburg
<adam.ginsburg at colorado.edu>wrote:

> Hi folks,
>   I'm trying to write a ray-tracing code for which high precision is
> required.  I also "need" to use square roots.  However, math.sqrt and
> numpy.sqrt seem to only use single-precision floats.  Is there a
> simple way to make sqrt use higher precision?  Alternately, am I
> simply being obtuse?
>
> Thanks,
> Adam
>
> Example code:
> from scipy.optimize.minpack import fsolve
> from numpy import sqrt
>
> sqrt(float64(1.034324523462345))
> # 1.0170174646791199
> f=lambda x: x**2-float64(1.034324523462345)**2
>
> f(sqrt(float64(1.034324523462345)))
> # -0.03550269637326231
>
> fsolve(f,1.01)
> # 1.0343245234623459
>
> f(fsolve(f,1.01))
> # 1.7763568394002505e-15
>
> fsolve(f,1.01) - sqrt(float64(1.034324523462345))
> # 0.017307058783226026
> ____


The routines *are* in double precision, but why are you using fsolve?
optimize.zeros.brentq would probably be a better choice. Also, you are using
differences with squared terms that will lose you precision. The last time I
wrote a ray tracing package was 30 years ago, but I think much will depend
on how you represent the curved surfaces. IIRC, there were also a lot of
quadratic equations to solve, and using the correct formula for the root you
want (the usual formula is poor for at least one of the roots) will also
make a difference. In other words, you probably need to take a lot of care
with how you set up the problem in order to minimize roundoff error.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20091017/b2fcc3fa/attachment.html>