2010/11/12

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)
12
gosh> (sqrt 169/81)
13/9
gosh> (sqrt 340282366920938463463374607431768211456)
18446744073709551616
gosh> (sqrt 32317006071311007300714876688669951960444102669
71548403213034542752465513886789089319720141152291346368871
79609218980194941195591504909210950881523864482831206308773
67300996091750197750389652106796057638384067568276792218642
61975616183809433847617047058164585203630504288757589154106
58086075523991239303855219143333896683424206849747865645694
94856176035326322058077805659331026192708460314150258592864
17711672594360371846185735759835115230164590440369761323328
72312271256847108202097251571017269313234696785425806566979
35045997268352998638215525166389437335543602135433229604645
318478604952148193555853611059596230656)
17976931348623159077293051907890247336179769789423065727343
00811577326758055009631327084773224075360211201138798713933
57658789768814416622492847430639474124377767893424865485276
30221960124609411945308295208500576883815068234246288147391
31105408272371633505106845862982399472459384797163048353563
29624224137216

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