xsqrt - Zeda/z80float GitHub Wiki

The algorithm used for xsqrt is as follows:

• If the input is negative, return NaN.
• If the input is 0, NaN, or inf, then return the input as the output.

Otherwise:

• The new exponent will be half the input.
  • The exponent has a bias of 0x4000. We basically need to do ((exp-0x4000)>>1)+0x4000 which is (exp+0x4000)>>1. Since we have a positive float, we don't have to worry about the sign bit messing stuff up.
  • If the shift does not set the carry flag , then we shift the mantissa right. We do this so that we don't need to multiply the final result by sqrt(1/2).

The actual algorithm is a convoluted implementation of Newton's Method. If you want a proof of the algorithm, just let me know.

• Take the top 16 bits and use sqrtHL to find the square root (x), with remainder (a). (ex. sqrt(19) returns 4 with remainder 3, since 19=4*4+3).
  • We now have the first 8 bits of the square root!
  • The bottom 48 bits of the input are also the bottom 48 bits of the "remainder".
• Next, we divide a (the remainder), by 2x, stopping at 8 bits and keeping the remainder. We will call the quotient, b.
  • b is now the next 8 bits of the square root! x is now 16 bit
  • Our remainder for the new division is our new a
• From our new a, subtract b*b. Remember, at this point, you can think of a as a 17.40 fixed-point number and b as an integer.
  • Keep track of the sign ! This cause a carry !
• Now we divide a by 2x again, but this time to 16 bits. The quotient is our new b
  • The quotient is now the next 16 bits of the square root. The remainder is our new a.
• subtract b*b from a
  • Note that this is a 16x16 multiply
• Finally, we divide our new a by 2x, getting 32 bits. No need to preserve a remainder. This quotient is the final 32 bits of the result.