Gauche Devlog

< More control on redirection in run-process | Partial continuations and fun of stack fiddling >


Exact sqrt

Several months ago one of our users asked me why (sqrt 4) yields 2.0; if the argument is an exact square number, should the result be exact as well?

Yes, it should, although coercing the result to inexact numbers is allowed in R5RS, I believe. It was just as it was because I was lazy: "Fix it when it becomes a problem." It seemed that the time had come.

It turned out to be an interesting little exercise to implement a 'proper' sqrt that returns an exact result when the input is a square of an exact number. That includes rational numbers. Maybe it's a nice quiz for CS freshmen. If you haven't implemented one, just think how you'll do it.

Here are some results in the current Gauche (the long digit numbers are wrapped around):

gosh> (sqrt 144)
gosh> (sqrt 169/81)
gosh> (sqrt 340282366920938463463374607431768211456)
gosh> (sqrt 32317006071311007300714876688669951960444102669

Let's forget about non-integer rationals and think about exact integers. I ended up handling three cases depending on the range of the input.

  • If the input is less than 2^52, the integer can be represented exactly in IEEE double floating point number, and if FP sqrt yields an integer we know it's the result. FP sqrt which is fast on the modern hardware. (Edit: This part is corrected on 2013/03/20 02:14:15 UTC, thanks to Mark H Weaver. See this entry.)
  • Otherwise we run Newton-Rhapson with exact arithmetic.
    • If the input is representable by IEEE double (e.g. up to approx. 2.0^1024), we can get an approximated result by using FP sqrt, and we use an integer close to it as the initial approximation.
    • Otherwise we use 2^floor((log2(input)+1)/2) as the initial approximation. That's because log2(input) may be calculated quickly.

Other than utilizing FP sqrt, this is pretty straightforward. If you know clever implementation of exact integer sqrt, let me know.

For exact rationals, I naively apply exact integer sqrt to both its numerator and its denominator. I guess that's faster than running Newton-Rhapson directly with exact rationals, for Gauche's rational artihmetic isn't optimized at all.

Oh, btw, R6RS's exact-integer-sqrt also returns the remainder if the input is not an exact square of an integer. That adds another little fun to the exercise.

Tags: sqrt, exact-integer-sqrt

Past comment(s)

okagawa (2010/11/14 02:38:23):

Maybe not for your use, but there is an interesting inverse sqrt approx. implementation in Quake 3. Hope this makes you enjoy.

shiro (2010/11/14 03:54:08):

Yeah, whenever people talk about cool hacks I bring it as the coolest one I've ever seen.

Mark H Weaver (2013/03/19 16:13:20):

There is a mistake in this post. Unfortunately, it is not true that if x < 2^53 and FP sqrt(x) is an integer (using IEEE 64-bit doubles), then we know the result is exact. sqrt(9007199136250226) rounds to an integer, but that result is not exact.

shiro (2013/03/19 23:24:31):

Mark, thanks for letting me know! Indeed, (sqrt 9007199136250226.) => 94906265.0. I think I did back-of-the-envelope error estimation when I wrote this, but that was wrong.

I'll write up a new entry for error estimation and come back to this entry to fix.

Post a comment