Talk:Square root algorithms

Reciprocal of the square root

This piece of code is a composite of a quirky square root starting estimate and Newton's method iterations. Newton's method is covered elsewhere; there's a section on Rough estimate already - the estimation code should go there. Also, I first saw this trick in Burroughs B6700 system intrinsics about 1979, and it predated my tenure there, so it's been around a long time. That's well before the drafting of IEEE 754 in 1985. Since the trick is based on linear approximation of an arc seqment of x² which in the end, is how all estimates must be made, I'm certain that the method has been reinvented a number of times.

There're numerous issues with this code:

• the result of type punning via pointer dereferencing in C/C++ is undefined
• the result of bit-twiddling floating point numbers, bypassing the API, is undefined
• we don't know whether values like zero, infinity and unnormalized floating point numbers, and big/little endian formats are correctly handled.
• It definitely won't work on architectures like Unisys mainframes which are 48/96 bit, or 64-bit IEEE floats. Restructuring the expression to make it work in those cases is non-trivial.
• since I can readily find, by testing incremental offsets from the original, a constant which reduces the maximum error, the original constant isn't optimal, probably resulting from trial and error. How does one verify something that's basically a plausible random number? That it works for a range of typical values is cold comfort. (Because its only use is as an estimate, maybe we don't actually care that enumerable cases aren't handled(?)... they'll just converge slowly.)
• because it requires a multiply to get back a square root, on architectures without fast multiply, it won't be such a quick estimate relative to others (if multiply and divide are roughly comparable, it'll be no faster than a random seed plus a Newton iteration).

I think we should include at least an outline of the derivation of the estimate expression, thus: a normalized floating point number is basically some power of the base multiplied by 1+k, where 0 <= k < 1. The '1' is not represented in the mantissa, but is implicit. The square root of a number around 1, i.e. 1+k, is (as a linear approximation) 1+k/2. Shifting the fp number represented as an integer down by 1 effectively divides k by 2, since the '1' is not represented. Subtracting that from a constant fixes up the 'smeared' mantissa and exponent, and leaves the sign bit flipped, so the result is an estimate of the reciprocal square root, which requires a multiply to re-vert back to the square root. It's just a quick way of guessing that gives you 1+ digits of precision, not an algorithm.

That cryptic constant is actually a composite of three bitfields, and twiddling it requires some understanding of what those fields are. It would be clearer, but a few more operations, to do that line as a pair of bitfield extract/inserts. But we're saving divides in the subsequent iterations, so the extra 1-cycle operations are a wash.

Erroneous introduction

The first paragraph of the article is this:

"Methods of computing square roots are numerical analysis algorithms for finding the principal, or non-negative, square root (usually denoted √S, ²√S, or S^1/2) of a real number. Arithmetically, it means given S, a procedure for finding a number which when multiplied by itself, yields S; algebraically, it means a procedure for finding the non-negative root of the equation x² - S = 0; geometrically, it means given the area of a square, a procedure for constructing a side of the square."

I do not know much about numerical analysis. But I do know that this is totally misleading!

A "method for computing a square root" of a number does not mean a method for *finding* the square root of that number.

Instead, it means a method for approximating the square root of the number, usually to various degrees of accuracy.

The distinction is essential to understanding what a "method for computing a square root" means. So the article should not mislead readers with an erroneous first paragraph.66.37.241.35 (talk) 16:39, 7 October 2020 (UTC)

Is this not dealt with adequately in the second paragraph? JRSpriggs (talk) 17:05, 7 October 2020 (UTC)

I'd say no, it's not, because however well it is dealt with in the second paragraph, that's the SECOND paragraph, and it's the FIRST paragraph which should mention that it's an approximation, since that's what is being accomplished. Therefore I think the word 'finding' in the first paragraph should be changed to approximating. UnderEducatedGeezer (talk) 02:48, 8 December 2021 (UTC)

{=3 } =4

In procedure int32_t isqrt(int32_t n)

• 3 left curly brackets
• 4 right curly brackets Jumpow (talk) 20:37, 21 December 2021 (UTC)

I tried to fix it. Is it OK now? JRSpriggs (talk) 22:57, 22 December 2021 (UTC)

Yes, now OK Jumpow (talk) 21:33, 4 January 2022 (UTC)

Undefined behaviour

The examples using unions are invalid C, as they invoke undefined behaviour. An easy solution that is probably even clearer for the purpose of example code would be to use memcpy, e.g.

  float f;
  uint32_t u;
  memcpy (&u, &f, sizeof (u));

37.49.68.13 (talk) 13:08, 1 April 2022 (UTC)

binary method in c

I believe that the example that is given in C is what is referred to as the Chinese Abacus Method - confirmation from this article and code that appears to be the same algorithm.