[math-fun] Origin of SQRT hack?
hgb>When making an initial guess of sqrt(x) for a Newton iteration, where x is a positive floating point number, a good initial guess is the bits of the entire number (both exponent & mantissa) shifted right by 1. This is the sort of thing that would have been done on the 7090 or the PDP-10. Does anyone here recall seeing this trick in the 1950's or 1960's ? I checked, and it isn't in HAKMEM. I've used this trick myself (at least on machines where moving numbers between the integer & floating point units isn't prohibitive). ------ rwg>It's really old, and can probably be found in the 704 and pdp-6 library sqrts. I think Knuth told me he found it as an undergrad, but that basketball youtube shows him using a biquinary machine (650?). --rwg
Gosper is correct, it's surely quite old--how old I can't say. "First you forget to zip up. Then you forget to zip down. Finally you forget your assembly hacks." I first encountered it in the early 1970's in the code for PDP-10 Space War, and surely it long predates that (again, as RWG suggests, a sufficiently motivated historian might look at the PDP-6 and earlier libraries--perhaps the Computer History Museum has them?). By the time I saw it it clearly bespoke evolving sophistication. Not only did it whack the approximation twice with unrolled copies of the iteration (no loop & test--twice being provably sufficient for the desired accuracy) but it had a strange twiddle before that (an OR or something, I forget), right after the initial shift, involving that alien bit from across the border between the exponent and mantissa. Famously, the only comment in that entire page of highly mysterious source code was "Watch out for the bum who wrote this!" I think I adapted it for the S-1, as unremarkable best practice du jour. Since Knuth was hanging around the Lab around then too he also might have run across it, whether or not it crept into TAOCP et seq. Indeed, it's plausible he independently invented it (heh, clever lad, that Don). At Apple (in pre-68000 days) and elsewhere I've often been annoyed hearing it (erroneously, I trust) ascribed to Kahan. "Those who cannot learn from history are doomed to reimplement it." --MLB
="Bill Gosper" <billgosper@gmail.com>
hgb>When making an initial guess of sqrt(x) for a Newton iteration, where x is a
positive floating point number, a good initial guess is the bits of the entire number (both exponent & mantissa) shifted right by 1.
This is the sort of thing that would have been done on the 7090 or the PDP-10. Does anyone here recall seeing this trick in the 1950's or 1960's ?
I checked, and it isn't in HAKMEM.
I've used this trick myself (at least on machines where moving numbers between the integer & floating point units isn't prohibitive). ------ rwg>It's really old, and can probably be found in the 704 and pdp-6 library sqrts. I think Knuth told me he found it as an undergrad, but that basketball youtube shows him using a biquinary machine (650?). --rwg _______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
Suppose you go t'other way: given flonum x approximate x^2 by left shifting it (assume no overflow). How close? Can you get it even closer using only a fixed number of "linear" operations (add, shift, bitwise logical ops etc)--short, of course, of "morally" multiplying via shift-and-add etc?
I'll bite. Let's come up with a Newton approximation for x^2, which is the inverse function for sqrt(x). Find x, s.t. f(x) = y-sqrt(x) = 0. x_i+1 = x_i - f(x_i)/f'(x_i) = x_i - (y-sqrt(x_i))/(-1/(2*sqrt(x_i))) = x_i + 2*sqrt(x_i)*(y-sqrt(x_i)) = sqrt(x_i)*(2*y-sqrt(x_i)) So (leaving aside the problem of how to compute sqrt(x)!), we would need to figure out how many iterations are required to get a decent approximation to x^2. BTW, I have an analogous problem with my proposed "asinh" representation for numbers. Traditionally simple operations like x^2 aren't that simple, so a lot of "simple" operations are going to involve Newton and/or other approximations. If you look carefully at asinh(x), the "high order" bits form the "exponent" of x, while the "low order" bits form the "mantissa" of x. Asinh(x) is thus seen to be a perfectly smooth floating point system without the "bumps" introduced by the arbitrary split between the exponent & mantissa of usual floating point systems. I.e., asinh(x) is already a representation with "gradual underflow" a la Kahan. The asinh "number system" deserves an entry in a 21st Century HAKMEM ("HakMem" ??). BTW, asinh isn't just my proposal; it is used by astronomers as an "improvement" for measuring magnitudes: http://www.sdss.org/dr7/algorithms/photometry.html#asinh http://adsabs.harvard.edu/cgi-bin/bib_query?1999AJ....118.1406L At 05:46 PM 4/14/2011, Marc LeBrun wrote:
Suppose you go t'other way: given flonum x approximate x^2 by left shifting it (assume no overflow). How close? Can you get it even closer using only a fixed number of "linear" operations (add, shift, bitwise logical ops etc)--short, of course, of "morally" multiplying via shift-and-add etc?
Richard Fateman contributes: sure you may quote me; i would not be surprised if it is described in one of his files at www.cs.berkeley.edu/~wkahan<http://www.cs.berkeley.edu/%7Ewkahan> but I can't spot it. W. Kahan. Implementation of algorithms (lecture notes by W. S. Haugeland and D. Hough). Technical Report 20, Department of Computer Science, University of California, Berkeley, CA, 1973. might contain it, but I can't spot it online.but see this, Bentley column.. http://delivery.acm.org/10.1145/320000/315707/p1155-bentley.pdf?key1=315707&... On Thu, Apr 14, 2011 at 8:32 AM, R Fateman <fateman@eecs.berkeley.edu>wrote:
w kahan taught this trick; I don't know if he invented it. He "proved" his 7090 code was optimal by essentially enumerating all 1-argument subroutines of fewer (plausible) instructions. Something like 12 instructions.
RJF
------------- Tom Duff>It worked pretty well on old enough hardware, but it mostly can't keep up
on modern machines with built-in sqrt (last I checked, which wasn't recently.) rwg>Eh? The last hardware sqrt I recall was on the Packard Bell 250 (1963?). (Integer sqrt, with remainder.) I have been griping for years at its absence from modern order codes, since it is far simpler than divide. So they've finally come around to reviving it? --rwg mrob> The next ("ultimate") editions of Knuth's books, if he manages to finish them will use MMIX[1] which uses IEEE floating-point, but I think it's unlikely Knuth would have held on to anecdotes about this right-shift hack just to include them 40 years later... rwg>Did you check Vol 4a? It *is* (in part) about bit bumming. On Wed, Apr 13, 2011 at 11:35 PM, Bill Gosper <billgosper@gmail.com> wrote:
hgb>When making an initial guess of sqrt(x) for a Newton iteration, where x is a
positive floating point number, a good initial guess is the bits of the entire number (both exponent & mantissa) shifted right by 1.
This is the sort of thing that would have been done on the 7090 or the PDP-10. Does anyone here recall seeing this trick in the 1950's or 1960's ?
I checked, and it isn't in HAKMEM.
I've used this trick myself (at least on machines where moving numbers between the
integer & floating point units isn't prohibitive). ------ rwg>It's really old, and can probably be found in the 704 and pdp-6 library sqrts. I think Knuth told me he found it as an undergrad, but that basketball youtube
shows him using a biquinary machine (650?). --rwg
This is not worth the $$ to me. Can anybody just tell us if it helps answer the original question, or even give the relevant quote?
rwg>Did you check Vol 4a? It *is* (in part) about bit bumming
Oh yeah, tAoCPv4. A lot of great stuff in the exercises of fascicle 1A. I just lost a good amount of time looking at that (-: Anyway, doesn't seem to be there either. The broadword stuff is really cool though. On Thu, Apr 14, 2011 at 19:29, Bill Gosper <billgosper@gmail.com> wrote:
Richard Fateman contributes: [...] but see this, Bentley column..
http://delivery.acm.org/10.1145/320000/315707/p1155-bentley.pdf?key1=315707&...
-- Robert Munafo -- mrob.com Follow me at: fb.com/mrob27 - twitter.com/mrob_27 - mrob27.wordpress.com- youtube.com/user/mrob143 - rilybot.blogspot.com
participants (4)
-
Bill Gosper -
Henry Baker -
Marc LeBrun -
Robert Munafo