Mailman 3 double-precision sqrt? - NumPy-Discussion

double-precision sqrt?

Adam Ginsburg

17 Oct 2009 17 Oct '09

6:08 p.m.

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

Show replies by date

Adam Ginsburg

17 Oct 17 Oct

6:17 p.m.

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.

Anne Archibald

6:37 p.m.

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

Adam Ginsburg

10:03 p.m.

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 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 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.

Dag Sverre Seljebotn

6:25 p.m.

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)) -- Dag Sverre

Nadav Horesh

6:27 p.m.

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

Charles R Harris

6:31 p.m.

On Sat, Oct 17, 2009 at 12:08 PM, 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?

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

5298

Age (days ago)

5298

Last active (days ago)

List overview

Download

6 comments

5 participants

participants (5)

Adam Ginsburg
Anne Archibald
Charles R Harris
Dag Sverre Seljebotn
Nadav Horesh

double-precision sqrt?

Adam Ginsburg

Adam Ginsburg

Anne Archibald

Adam Ginsburg

Dag Sverre Seljebotn

Nadav Horesh

Charles R Harris

tags

participants (5)