April 2011 - mathfun@mailman.xmission.com

Re: [math-fun] Origin of SQRT hack?
by Henry Baker 20 Apr '11

20 Apr '11

Very cool! Thx, Tom! I searched for some of the strings in this code to see what I could find on Google, but could only find _later_ versions of the same code, so your reference is the best so far. At 04:42 PM 4/18/2011, Tom Duff wrote: >Henry Baker <hbaker1(a)pipeline.com> writes: > >>After some additional thinking, it would appear that the "SQRT hack", under Robert Munafo's definition & my analysis, would appear to work _only on FP representations which have a "hidden bit"_. Therefore, the IBM & early DEC series machines need not apply. >> >>If we're going to find examples of this type of code, we're going to have to look at PDP11 & VAX routines, as they were apparently the first wide-spread use of "hidden bit" floating point units. >> >>This means that the early Unix hackers -- including those at Bell Labs & Berkeley -- may be responsible. > >As it happens, I have a copy of the 6th Edition UNIX source code right here on my disk. Here's the PDP-11 assembly code for the C library sqrt function, with the instructions that compute the initial guess at the square root marked *****. Sure enough, it pushes fr0 (the argument) on the stack, shifts it right one bit, adds a magic constant and pops the result back into fr0. This code was released by Bell Labs in May 1975, which puts an upper bound on the date of discovery. Most likely this routine was written by Bob Morris, who was responsible for much of the math library code in early UNIX editions. This code is probably also in earlier UNIX editions that you might be able to find on line. > >ldfps = 170100^tst >stfps = 170200^tst >/ >/ sqrt replaces the f.p. number in fr0 by its >/ square root. newton's method >/ >.globl sqrt, _sqrt >/ >/ >_sqrt: > mov r5,-(sp) > mov sp,r5 > movf 4(r5),fr0 > jsr pc,sqrt > mov (sp)+,r5 > rts pc > >sqrt: > tstf fr0 > cfcc > bne 1f > clc > rts pc /sqrt(0) >1: > bgt 1f > clrf fr0 > sec > rts pc / sqrt(-a) >1: > mov r0,-(sp) > stfps -(sp) > mov (sp),r0 > bic $!200,r0 / retain mode > ldfps r0 > movf fr1,-(sp) > movf fr2,-(sp) >/ > movf fr0,fr1 > movf fr0,-(sp) / ***** > asr (sp) / ***** > add $20100,(sp) / ***** > movf (sp)+,fr0 /initial guess / ***** > mov $4,r0 >1: > movf fr1,fr2 > divf fr0,fr2 > addf fr2,fr0 > mulf $half,fr0 / x = (x+a/x)/2 > sob r0,1b >2: > movf (sp)+,fr2 > movf (sp)+,fr1 > ldfps (sp)+ > mov (sp)+,r0 > clc > rts pc >/ >half = 40000 > >-- >Tom Duff. chan[2]='m'; /* Ooh, careful! */

3 3

[math-fun] Towards Newton's method a la mode.
by Marc LeBrun 20 Apr '11

20 Apr '11

Sorry if you've heard this before, but recent threads here inescapably rekindled this rant. It is exciting and empowering to learn to compute famous functions for yourself. Stranded on a desert island it's arguably more useful to know how your lost calculator could have computed square roots than knowing how to fondle the UI of the latest iPhone app. As denizens of the burgeoning technosphere, being bereft of our cultural computing heritage renders us much like dogs sniffing at an ATM. Newton's method provides a wonderful excuse to convey invigorating insights. Alas Sir Isaac's method is usually presented almost as a kind of side-show hat trick, perhaps rigorously dissected, clearly clever, but missing some vital spark of motivation and burdened with distracting allusions irrelevant to enabling a broader contemporary application and appreciation. I suggest that root-finding is a historical red herring that ought to be suppressed from the main narrative in modern presentations. Of course zeros are inextricably entwined with structural topics like the Fundamental Theorem of Algebra, but they are confusing and non-essential in the context of Newton's method, which is, at heart, about computational procedures. We might instead present Newton's method motivated top-down like this: "Given x, we wish to compute the value of some function f(x). However we may not have a way to conveniently compute f(x) directly. Therefore we try to construct another function g(y), composed from pieces we can more readily compute, with the property that iterates of g(y) will converge on a fixpoint equal to f(x)." Note zeros and root-finding don't arise; the focus is on computation, composition (both in the construction of g and in its iteration) and convergence on fixpoints (which connects it to au courant dynamics etc). Interesting developments proceed from there: We can derive g's familiar "classical Newton transform" Ng explicitly by stipulating that g() is "sufficiently well-approximated" by the first two terms of its Taylor series and solving the expansion of the fixpoint constraint y = g(y). Further observations and explorations include the following: * As N can be easily inverted ('cause z/dz = 1/d ln z) EVERY expression when iterated is the "Newton's method" for SOME function. * The reason iterating some expressions so magically converges whilest others suck goes back to that assumption of "sufficiently well-approximated" * We can compare iterations of g, gg, ggg... iterations of the transformed function Ng, NgNg, NgNgNg... and iterations of the transformation Ng, NNg, NNNg etc and seek limiting fixpoints for all these, along with interesting special identities (NgNg=NNg) etc. * We can consider higher-order expansions--solving the quadratic Taylor approximation requires square roots, so it's useless for SQRT, but could be interesting for CBRT etc, or for transcendentals (eg find log given exp). * We can consider solving other expansions, say Lambert or Fourier series, continued fractions, etc. All with nary a zero on the horizon...

2 2

Re: [math-fun] Origin of SQRT hack?
by Henry Baker 19 Apr '11

19 Apr '11

Ok, here's the _extremely short_ version, which I'll expand upon later: On a "hidden bit" floating point system, the positive number 2^e(1+m) (0<m<1) is represented by the pair [e,m]. If we have a number z=2^(2e)*(1+2x) (0<x<.5), represented by [2e,2x], then a good initial approximation to sqrt(z) is [e,x]. Why? Consider squaring s0=[e,x]=2^e*(1+x): [e,x]^2 = (2^e*(1+x))^2 = 2^(2e)*(1+2*x+x^2) ~~ 2^(2e)*(1+2*x) = [2e,2x]. So you get the idea. So, the initial guess isn't a lot of precision, but every little "bit" helps... (I don't have time to expand upon this right this second, but the oddp(e) case looks more like 2^(e-1)*(1-m), so the low bit shifted out of the exponent becomes the "sign" bit for the mantissa.) At 02:30 PM 4/17/2011, Robert Munafo wrote: >On Sun, Apr 17, 2011 at 11:59, Henry Baker <<mailto:hbaker1@pipeline.com>hbaker1(a)pipeline.com> wrote: >> I think we've established that shifting the exponent right by 1 for the initial approximation goes way far back. I suspect that it may go back much further -- perhaps to Von Neumann's original binary floating point proposal. >> >> The only question now is what are the contents of the mantissa. I may need to do some experiments with Maxima to determine how important various guesses at the mantissa are. >> >> [...] > >No, you asked "[...] the bits of the entire number (both exponent & mantissa) shifted right by 1." > >I still want the answer to that question, not any of the others. > >I wrote assembly code for the PDP-11, and apparently it used a floating point format in which an integer test "A>B" was suitable for comparing two positive floating point numbers, so the monolithic shift would have worked there. Maybe other PDP's were similar. But we're looking for actual evidence of it. > >- Robert > >-- > Robert Munafo -- <http://mrob.com>mrob.com > Follow me at: <http://fb.com/mrob27>fb.com/mrob27 - <http://twitter.com/mrob_27>twitter.com/mrob_27 - <http://mrob27.wordpress.com>mrob27.wordpress.com - <http://youtube.com/user/mrob143>youtube.com/user/mrob143 - <http://rilybot.blogspot.com>rilybot.blogspot.com

2 4

Re: [math-fun] Origin of SQRT hack?
by Henry Baker 19 Apr '11

19 Apr '11

After some additional thinking, it would appear that the "SQRT hack", under Robert Munafo's definition & my analysis, would appear to work _only on FP representations which have a "hidden bit"_. Therefore, the IBM & early DEC series machines need not apply. If we're going to find examples of this type of code, we're going to have to look at PDP11 & VAX routines, as they were apparently the first wide-spread use of "hidden bit" floating point units. This means that the early Unix hackers -- including those at Bell Labs & Berkeley -- may be responsible. I would look in the graphics & vision libraries, as they had the greatest need for speed at some cost in precision, and they certainly didn't care about denormalized numbers. At 10:02 AM 4/18/2011, Robert Munafo wrote: >Henry, thanks for all your work on this -- I hope my persistent questions weren't too frustrating. > >For now, I'm treating the 1997 source code by Blinn as the earliest-known code that meets the criteria (namely a single shift followed by Newton iteration and no conditions). > >I did note that it doesn't work with denormals (in fact that was in my article on 0x5f3759df from the start) and of course it doesn't work on infinity either. > >You can read it at mrob.com/pub/math/numbers-16.html#le009_16

2 3

Re: [math-fun] Origin of SQRT hack?
by Tom Duff 19 Apr '11

19 Apr '11

Henry Baker <hbaker1(a)pipeline.com> writes: > After some additional thinking, it would appear that the "SQRT hack", > under Robert Munafo's definition & my analysis, would appear to work > _only on FP representations which have a "hidden bit"_. Therefore, the > IBM & early DEC series machines need not apply. > > If we're going to find examples of this type of code, we're going to > have to look at PDP11 & VAX routines, as they were apparently the first > wide-spread use of "hidden bit" floating point units. > > This means that the early Unix hackers -- including those at Bell Labs & > Berkeley -- may be responsible. As it happens, I have a copy of the 6th Edition UNIX source code right here on my disk. Here's the PDP-11 assembly code for the C library sqrt function, with the instructions that compute the initial guess at the square root marked *****. Sure enough, it pushes fr0 (the argument) on the stack, shifts it right one bit, adds a magic constant and pops the result back into fr0. This code was released by Bell Labs in May 1975, which puts an upper bound on the date of discovery. Most likely this routine was written by Bob Morris, who was responsible for much of the math library code in early UNIX editions. This code is probably also in earlier UNIX editions that you might be able to find on line. ldfps = 170100^tst stfps = 170200^tst / / sqrt replaces the f.p. number in fr0 by its / square root. newton's method / .globl sqrt, _sqrt / / _sqrt: mov r5,-(sp) mov sp,r5 movf 4(r5),fr0 jsr pc,sqrt mov (sp)+,r5 rts pc sqrt: tstf fr0 cfcc bne 1f clc rts pc /sqrt(0) 1: bgt 1f clrf fr0 sec rts pc / sqrt(-a) 1: mov r0,-(sp) stfps -(sp) mov (sp),r0 bic $!200,r0 / retain mode ldfps r0 movf fr1,-(sp) movf fr2,-(sp) / movf fr0,fr1 movf fr0,-(sp) / ***** asr (sp) / ***** add $20100,(sp) / ***** movf (sp)+,fr0 /initial guess / ***** mov $4,r0 1: movf fr1,fr2 divf fr0,fr2 addf fr2,fr0 mulf $half,fr0 / x = (x+a/x)/2 sob r0,1b 2: movf (sp)+,fr2 movf (sp)+,fr1 ldfps (sp)+ mov (sp)+,r0 clc rts pc / half = 40000 -- Tom Duff. chan[2]='m'; /* Ooh, careful! */

2 1

[math-fun] "Rational Trigonometry" ??
by Henry Baker 18 Apr '11

18 Apr '11

Has anyone here dealt with Rational Trigonometry? I found it as a link through Chebyshev polynomials ?!? http://web.maths.unsw.edu.au/~norman/Rational1.htm http://en.wikipedia.org/wiki/Chebyshev_polynomial

2 3

Re: [math-fun] Origin of SQRT hack?
by Henry Baker 18 Apr '11

18 Apr '11

Thanks, Rich! It should have been written: "The function sinh(x) acts like exp(x)/2 for large x, and like x for small x". At 08:37 AM 4/18/2011, Schroeppel, Richard wrote: >Minor typo ... >sinh(x) ~ (e^x)/2 for large x; and asinh(x) ~ log(2x) for large x. ---Rich > >-----Original Message----- >From: math-fun-bounces(a)mailman.xmission.com [mailto:math-fun-bounces@mailman.xmission.com] On Behalf Of Henry Baker >Sent: Monday, April 18, 2011 9:27 AM >To: math-fun >Subject: Re: [math-fun] Origin of SQRT hack? > >I should point out the "denormalized" numbers completely blow out the "SQRT hack", regardless of the FP format. > >The function asinh(x) acts like exp(x)/2 for large x, and like x for small x -- i.e., like a "normalized" FP number for large x (albeit with a radix of e), and like a "denormalized" number for small x. So asinh(x) is the "smoothest" floating point format that includes denormalized numbers. > >So I tested Newton given a first approximation sinh(asinh(N)/2) to sqrt(N). > >In order to compute a full precision (in GNU Common Lisp) IEEE sqrt for N in the range 2^1 up to 2^512, exactly _five_ Newton iterations sufficed. > >However, when N started falling in the "denomalized" range -- i.e., N<1, the number of Newton iterations skyrocketed, basically requiring > >N Iters reqd >2^0 6 >2^-10 11 >2^-20 16 >2^-30 21 >2^-40 26 >2^-50 31 >etc. > >So, if you care about the square roots of denormalized numbers, you will have to treat them as a special case. > > >_______________________________________________ >math-fun mailing list >math-fun(a)mailman.xmission.com >http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun > > > >_______________________________________________ >math-fun mailing list >math-fun(a)mailman.xmission.com >http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun

1 0

Re: [math-fun] periodic nag
by Dan Asimov 18 Apr '11

18 Apr '11

Bill G. wrote: << Sum[k/n*Binomial[-n, n - k]*Fibonacci[k + 1], {k, 1, n}] == 1 ? >> I'm not entirely familiar with Mma syntax. May we assume there are parentheses around the k/n ? Thanks, Dan Sometimes the brain has a mind of its own.

3 2

[math-fun] 6205s
by Bill Gosper 17 Apr '11

17 Apr '11

quad>As a quick aside, Wikipedia says "The PDP-6 was infamous because of the 6205 board, a large (11 x 9 inches) board which contained 1 bit of AR, MB, and MQ (thus there were 36 of these). It had 88 transistors, a 2-sided PC etch, two 18-pin and two 22-pin connectors (two on each side of the module). Because of all these connectors, swapping this module was a major undertaking, and the mechanical coupling made it highly likely that fixing one fault would cause another." following with "There was also a great fear of powering off a PDP-6, since it would generally result in at least one 6205 board failing." Is this true? -------------- I think they were larger than a foot square. I might even have one somewhere. As I recall their commonest failure modes were cold solder joints and pulse transformers ("red turds"). And they were hugely expensive. Often, the failures would elude DEC's Field Service lame diagnostic routines, which they also insisted on running as "preventative maintenance" (An ounce of preventation is worth a pound of curatativeness). Once, Stew Nelson proved to them the machine was broken with a program which buzzed the teletype bell whenever you turned out the room lights. Even after being shown the program, devoid of input of any type, Field Service ran their diagnostics with lights both off and on, shrugged, and went home. Nelson etc eventually determined that Floating Subtract following Floating Divide got the wrong answer with the room lights off because their hum was needed to boost the Master Clear pulse enough to completely clear out the state left over by the FDVR(ounded) The only great fear of powering off that I remember was when an idiot Multics operator wandered in and powered off the apparently unattended PDP-6, whereupon the AMF hydraulic arm nearly removed Gerry Sussman's head in the next room. --rwg Field Service also refused to fix a problem that Nelson demonstrated by Scotch-taping a spring-return switch in the up position while resting the console telephone on another switch, and twiddling the speed control.

3 3

Re: [math-fun] Origin of SQRT hack?
by Henry Baker 17 Apr '11

17 Apr '11

http://pdp-10.trailing-edge.com/decuslib10-10/01/43,50516/basxct.mac.html ;VERSION 17E 2-OCT-74/NA ;VERSION 17D 4-MAY-73/KK ;VERSION 17C 2-JAN-73/KK ;VERSION 17B 25-JUL-72/KK ;VERSION 17A 10-FEB-1972/KK ;VERSION 17 15-OCT-1971/KK ;VERSION 16 5-APR-1971/KK ;VERSION 15 17-AUG-1970/KK ;VERSION 14 16-JUL-1970/AL/KK ;VERSION 13 15-SEP-1969 ... ;FLOATING POINT SINGLE PRECISION SQUARE ROOT FUNCTION ;THE SQUARE ROOT OF THE ABSOLUTE VALUE OF THE ARGUMENT IS ;CALCULATED. THE ARGUMENT IS WRITTEN IN THE FORM ; X= F*(2**2B) WHERE 0.LT.F.LT.1 ;SQRT(X) IS THEN CALCULATED AS (SQRT(X))*(2**B) ;SQRT(F) IS CALCULATED BY N LINEAR APPROXIMATION, THE NATURE ;OF WHICH DEPENDS ON WHETHER 1/4 .LT. F .LT. 1/2 OR 1/2 .LT. F .LT. 1, ;FOLLOWED BY TWO ITERATIONS OF NEWTON'S METHOD. SQRTB: MOVE T, N ;PICK UP THE ARGUMENT IN T JUMPL T,SQRMIN ;SQRT OF NEGATIVE NUMBER? JUMPE T,SQRT1 ;CHECK FOR ARGUMENT OF ZERO SQRTB0: ASHC T, -33 ;PUT EXPONENT IN T, FRACTION IN T1 SUBI T, 201 ;SUBTRACT 201 FROM EXPONENT ROT T, -1 ;CUT EXP IN HALF, SAVE ODD BIT HRRM T,EX1 ;SAVE FOR FUTURE SCALING OF ANS ;IN FSC N,. INSTRUCTION LSH T, -43 ;GET BIT SAVED BY PREVIOUS INST. ASH T1, -10 ;PUT FRACTION IN PROPER POSITION FSC T1, 177(T) ;PUT EXPONENT OF FRACT TO -1 OR 0 MOVEM T1, N ;SAVE IT. 1/4 .LT. F .LT. 1 FMP T1, SQCON1(T) ;LINEAR FIRST APPROX,DEPENDS ON FAD T1, SQCON2(T) ;WHETHER 1/4.LT.F.LT.1/2 OR 1/2.LT.F.LT.1. MOVE T, N ;START NEWTONS METHOD WITH FRAC FDV T, T1 ;CALCULATE X(0)/X(1) FAD T1, T ;X(1) + X(0)/X(1) FSC T1, -1 ;1/2(X(1) + X(0)/X(1)) FDV N, T1 ;X(0)/X(2) FADR N, T1 ;X(2) + X(0)/X(2) XCT EX1 SQRT1: POPJ P, ;EXIT SQCON1: 0.8125 ;CONSTANT, USED IF 1/4.LT.FRAC.LT.1/2 0.578125 ;CONSTANT, USED IF 1/2.LT.FRAC.LT.1 SQCON2: 0.302734 ;CONSTANT, USED IF 1/4.LT.FRAC.LT.1/2 0.421875 ;CONSTANT, USED IF 1/2.LT.FRAC.LT.1 SQRMIN: PUSH P,T ;SAVE ARG NFERR (54,0) PUSHJ P,INLMES ASCIZ / % SQRT of negative number/ PUSHJ P,GOSR3 ;PRINT LINE NUMBER POP P,T ;GET ARG MOVMS T JRST SQRTB0 ;USE ABSOLUTE VALUE At 04:52 PM 4/15/2011, rcs(a)xmission.com wrote: >I think I remember seeing this in the PDP6 library function for sqrt. >I believe the shift was followed by adding (integer add) a magic constant, >probably to adjust the range and repair the exponent offset. Then a >few rounds of Newton, unwound. > >(note to folks under 50) Before IEEE floating point came along, the >individual computer makers had their own formats for floating point. >Both the 7090 series and PDP6/10 used 36-bit words. >The PDP6 format had 8 bits of exponent, offset by 128, and 27 bits of >mantissa. A normalized number included the high-order 1 bit, which is >dropped in IEEE format. I don't recall if that bit had value .5 or 1.0. >0.0 was simply 0, and negatives were the (integer) two's complement of >the positives. This had the virtue that the ordering of (normalized) >floating numbers was the same as the ordering of integers, so the >compare instructions didn't need to come in two flavors. > >I think I'll register a disagreement with the comment that IEEE accuracy >requirements for sqrt (exactly correct rounding) effectively preclude >doing it in software: After a few rounds of Newton, the answer should be >within one bit, likely on the high side. You can check which of two >consecutive floating values is correct by comparing V * (V-epsilon) versus X. > >Rich > >---- >Quoting Bill Gosper <billgosper(a)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

2 2