0

I have been trying to compute the square root from a fixed point data-type <24,8>.
Unfortunately nothing seems to work.
Does anyone know how to do this fast and efficient in C(++)?

Navnath Godse
  • 2,233
  • 2
  • 23
  • 32
Alex van Rijs
  • 803
  • 5
  • 17
  • 39

1 Answers1

3

Here is a prototype in python showing how to do a square root in fixed point using Newton's method.

import math

def sqrt(n, shift=8):
     """
     Return the square root of n as a fixed point number.  It uses a
     second order Newton-Raphson convergence.  This doubles the number
     of significant figures on each iteration.

     Shift is the number of bits in the fractional part of the fixed
     point number.
     """
     # Initial guess - could do better than this
     x = 1 << shift // 32 bit type
     n_one = n << shift // 64 bit type
     while 1:
         x_old = x
         x = (x + n_one // x) // 2
         if x == x_old:
             break
     return x

def main():
    a = 4.567
    print "Should be", math.sqrt(a)
    fp_a = int(a * 256)
    print "With fixed point", sqrt(fp_a)/256.

if __name__ == "__main__":
    main()

When converting this to C++ be really careful about the types - in particular n_one needs to be a 64 bit type or otherwise it will overflow on the <<8 bit step. Note also that // is an integer divide in python.

Nick Craig-Wood
  • 52,955
  • 12
  • 126
  • 132
  • weird. Newton rapson method as i understand works as f(x) = x0 - f(x)/f'(x). I don't see that function here. I see your x0 is 1. What are you doing with n_one? Can you post some explanation of this code? – newbie_old Jun 25 '15 at 04:43
  • @newbie_old if you work the differential of sqrt out and plug into the equation above and simplify, you'll get `x(i+1) = (x(i) + n/x(i))/2` where `n` is the number you are trying to square root. `n_one` is `n` shifted by the `shift` bits, which represents `n` in the fixed point world. – Nick Craig-Wood Jun 25 '15 at 15:02
  • but n_one is already in fixed point format no? I see you multiply that with 256 which is 1<<8 before feeding it to the function. So why again shifting by shift which is 8? – newbie_old Jun 25 '15 at 21:41
  • @newbie_old I see what you mean, my explanation is incorrect above. `n_one` is actually `n * one²`, ie `n_float << 16` so when you divide it by `x` which is `x_float * one` = `x_float << 8` you get `(n_float/x_float) * one` = `(n_float/x_float) << 8` which is correctly normalized, but with the full precision of the division. I hope that makes sense now! – Nick Craig-Wood Jun 26 '15 at 12:35