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
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 adam.ginsburg@colorado.edu 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.
2009/10/17 Adam Ginsburg adam.ginsburg@colorado.edu:
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
Hi again, I apologize, the mistake was entirely my own. Sqrt's do the right thing....
Adam
On Sat, Oct 17, 2009 at 12:17 PM, Adam Ginsburg adam.ginsburg@colorado.edu wrote:
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 adam.ginsburg@colorado.edu 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.
Adam Ginsburg 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?
How are you actually using the results of sqrt? When printing the results you may not get the full precision...try e.g.
print "%.50f" % np.sqrt(np.float64( 1.034324523462345))
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
On Sat, Oct 17, 2009 at 12:08 PM, Adam Ginsburg adam.ginsburg@colorado.eduwrote:
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
participants (5)
-
Adam Ginsburg -
Anne Archibald -
Charles R Harris -
Dag Sverre Seljebotn -
Nadav Horesh