On Sat, Oct 17, 2009 at 12:08 PM, Adam Ginsburg <adam.ginsb...@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
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
