Henry Baker <hbaker1@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! */