If one digs long enough, Kahan pops up. Inverse square root.

[pier4r]

We all know that W. Kahan (the author of many numerical solutions for the HP calculators) was great, but I was surprised to learn that he was behind the inverse square root solution.

https://www.netlib.org/fdlibm/e_sqrt.c (point B) , this is related to the FISR problem.

What I find funny (and sad) is that somehow people think that back then people were plain stupid or without equipment. A comment on youtube notice this when another user pointed out Kahan "It's still quite a feat to acquire that knowledge back in the day." (as if today would be easy to discover the solution)

I'd argue it is harder today to find such solutions, because nowadays people just throw hardware at the problem before optimizing (I am guilty of this). In the past they needed to optimize more given the limited HW.

[EdS2]

Nice find! From the comment block in the linked code:

Code:

(This is a copy of a drafted paper by Prof W. Kahan 
and K.C. Ng, written in May, 1986)
...
   (1)    Initial approximation

    Let x0 and x1 be the leading and the trailing 32-bit words of
    a floating point number x (in IEEE double format) respectively 

        1    11             52                  ...widths
       ------------------------------------------------------
    x: |s|      e     |          f                |
       ------------------------------------------------------
          msb    lsb  msb                      lsb ...order

 
         ------------------------           ------------------------
    x0:  |s|   e    |    f1     |     x1: |          f2           |
         ------------------------           ------------------------

    By performing shifts and subtracts on x0 and x1 (both regarded
    as integers), we obtain an 8-bit approximation of sqrt(x) as
    follows.

        k  := (x0>>1) + 0x1ff80000;
        y0 := k - T1[31&(k>>15)].    ... y ~ sqrt(x) to 8 bits
    Here k is a 32-bit integer and T1[] is an integer array containing
    correction terms. Now magically the floating value of y (y's
    leading 32-bit word is y0, the value of its trailing word is 0)
    approximates sqrt(x) to almost 8-bit.

    Value of T1:
    static int T1[32]= {
    0,    1024,    3062,    5746,    9193,    13348,    18162,    23592,
    29598,    36145,    43202,    50740,    58733,    67158,    75992,    85215,
    83599,    71378,    60428,    50647,    41945,    34246,    27478,    21581,
    16499,    12183,    8588,    5674,    3403,    1742,    661,    130,};

[Joe Horn]

(10-03-2023 07:37 AM)pier4r Wrote:  We all know that W. Kahan (the author of many numerical solutions for the HP calculators) was great ...

... and still is: https://en.wikipedia.org/wiki/William_Kahan

[pier4r]

(10-03-2023 02:18 PM)Joe Horn Wrote:  ... and still is: https://en.wikipedia.org/wiki/William_Kahan

Yes sorry, I wanted to mean for the work he has done. Not that he was dead (I always mixed up this part of the English language)

[Wes Loewer]

(10-03-2023 07:37 AM)pier4r Wrote:  What I find funny (and sad) is that somehow people think that back then people were plain stupid or without equipment.

I had the opposite reaction from a Geometry student once. When I was explaining how the Greeks figured it all out without calculators or computers, one student replied, "Wow, they were really smart back then."

Indeed they were.