On Sat, Oct 17, 2009 at 12:08 PM, Adam Ginsburg <adam.ginsburg@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