My code is actually wrong.... but I still have the problem I've identified that sqrt is leading to precision errors. Sorry about the earlier mistake. Adam On Sat, Oct 17, 2009 at 12:08 PM, Adam Ginsburg wrote:

sqrt(float64(1.034324523462345)) # 1.0170174646791199 f=lambda x: x**2-float64(1.034324523462345)**2

should be f=lambda x: x**2-float64(1.034324523462345) so the code I sent was not a legitimate test.

The default precision is double unless yue specify otherwise (float32 or long double (float128 or float96)) You can see this from: f(fsolve(f,1.01)) # 1.7763568394002505e-15 The last line should be:

...

...

fsolve(f,1.01) - float64(1.034324523462345) 8.8817841970012523e-16

Nadav -----הודעה מקורית----- מאת: numpy-discussion-bounces@scipy.org בשם Adam Ginsburg נשלח: ש 17-אוקטובר-09 20:08 אל: numpy-discussion@scipy.org נושא: [Numpy-discussion] double-precision sqrt? 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 _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

2009/10/17 Adam Ginsburg :

My code is actually wrong.... but I still have the problem I've identified that sqrt is leading to precision errors.  Sorry about the earlier mistake.

I think you'll find that numpy's sqrt is as good as it gets for double precision. You can try using numpy's float96 type, which at least on my machine, does give sa few more significant figures. If you really, really need accuracy, there are arbitrary-precision packages for python, which you could try. But I think you may find that your problem is not solved by higher precision. Something about ray-tracing just leads it to ferret out every little limitation of floating-point computation. For example, you can easily get "surface acne" when shooting shadow rays, where a ray shot from a surface to a light source accidentally intersects that same surface for some pixels but not for others. You can try to fix it by requiring some minimum intersection distance, but then you'll find lots of weird little quirks where your minimum distance causes problems. A better solution is one which takes into account the awkwardness of floating-point; for this particular case one trick is to mark the object you're shooting rays from as not a candidate for intersection. (This doesn't work, of course, if the object can cast shadows on itself...) I have even seen people advocate for using interval analysis inside ray tracers, to avoid this kind of problem. Anne