Talk:Square root algorithms

The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

Archives

Archive 1

This page has archives. Sections older than 360 days may be auto-archived by if there are more than 5.

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.

Merge "Approximations that depend on the floating point representation" into "Initial estimate"

I believe the section "Approximations that depend on the floating point representation" should be merged into "Initial estimate", since it is a special case of "Binary estimates". Merging would clear up the fact that the floating point trick gives an initial rough approximation, which is then typically iteratively improved.

I also believe the "Initial estimate" section should appear after the section on Heron's method, as the reader is likely more interested in the general idea of iterative refinement than in the details of how to obtain a good initial estimate in all possible ways.

Additionally, in my opinion the entirety of the article could benefit from some trimming/rewriting, as many sections contain redundant information, unnecessary details, and awkward formulations. BlueRavel (talk) 14:54, 4 December 2023 (UTC)[reply]

Your proposition makes sense to me, and I dont necessarily disagree. That said though, as a pure mathematician, I am uninclined to blur the lines between programmatical issues and mathematical problems. I think maintaining a distinction is appropriate. An analysis of the pure mathematical problem of initial estimation in these abstract reiterative processes is a decidedly distinct discussion from considerations in this programming language, or that programming language, or this architecture, or that architecture. The former is future-proofed, the latter is not. CogitoErgoCogitoSum (talk) 21:09, 11 February 2024 (UTC)[reply]

Useful addition??

Not sure if its useful, but I have found that, in general, ${\sqrt {x+2}}\approx {\frac {x+1}{\sqrt {x}}}$ , and if $x = n 2$ we get ${\sqrt {n^{2}+2}}\approx n+{\frac {1}{n}}$ .

Similarly ${\sqrt {x+4}}\approx {\frac {x+2}{\sqrt {x}}}$ .

I sometimes use this for quick pencil and paper calculations, if Im close enough to a convenient value.

Not sure if this is a known or established property, proven, bounded, or if its already in the article in some alternative capacity, or if its even appropriate for this article. I do know the taylor series approximation with two terms connects these expressions. CogitoErgoCogitoSum (talk) 21:05, 11 February 2024 (UTC)[reply]

There is nothing special about 2 and 4:

{\sqrt {x+2c}}\approx {\frac {x+c}{\sqrt {x}}}

provided that c is small compared to x. This is, in fact, just the first two terms of the series given in the article under the section heading "Taylor series". JBW (talk) 01:45, 13 February 2024 (UTC)[reply]

I don't think they are useful. In the first, you have replaced a square root and an addition with a square root, an addition, and a division to get an approximate answer. Bubba73 ^{You talkin' to me?} 08:02, 13 February 2024 (UTC)[reply]

Python tweak needed

I noticed the interesting Python example and thought I had better look at it before a gnome removes it as OR or whatever. There is a syntax error in the following line: the inner "..." should use apostrophes not quotes because of the quotes on the overall string.

print(f"2) {sqrt_Heron(Decimal("3.1415926535897932384626433832795028841971693993"))}")

Johnuniq (talk) 03:42, 21 October 2025 (UTC)[reply]

Fixed. Changed to

print(f"2) {sqrt_Heron(Decimal('3.1415926535897932384626433832795028841971693993'))}")

Marc Schroeder (talk) 13:25, 22 October 2025 (UTC)[reply]

Halley's method versus Heron's method

Halley's method should not be mixed into Heron's method, but put into a separate section.

Halley's method $x_{n+1}=x_{n}\cdot {\frac {x_{n}^{2}+3S}{3x_{n}^{2}+S}}$ requires one division and three multiplications per iteration. I do not count the factor of 3 in the denominator because that can be done by a shift and an addition which are negligible compared to multiplication. The 3S in the numerator need only be done once at the outset, so it is negligible also. After the initial division, the quotient is known to sufficient precision that we only need to do one step of the division algorithm which amounts to two multiplications. Thus the cost of one iteration of Halley's method amounts to five multiplications. Two iterations would be ten multiplications for a nine fold increase in precision.

By contrast, Heron's method only requires one division. (Multiplication by one half is just a shift and the addition is also negligible.) So it costs two multiplications per iteration. Four iterations of Heron's method would thus cost eight multiplications for a sixteen fold increase in precision. So Heron's method costs less while giving you more precision. JRSpriggs (talk) 15:49, 7 November 2025 (UTC)[reply]