4

I am trying to implement Fast Inverse Square Root for a fixed point number, but I'm not getting anywhere.

I am trying to follow exactly the same principle as the article, except instead of writing the number in the floating point format x = (-1) ^ s * (1 + M) * 2 ^ (E-127), I am using the format x = M * 2 ^ -16, which is a 32-bit fixed point number with 16 decimal bits and 16 fractional bits.

The problem is that I cannot find the value of the "magic constant". According to my calculations, it doesn’t exist, but I’m not a mathematician and I think I’m doing everything wrong.

To solve Y = 1 / sqrt (x), I used the following reasoning (I don't know if it is correct).

In the original code we have that Y0 for approximation of newton is given by:

i = 0x5f3759df - (i >> 1);

Which means that we will have as a result a floating point number given by:

y0 = (1 + R2 - M / 2) * 2 ^ (R1 - E / 2);

This is because the operation >> divides exponent and mantissa by 2, and then we perform a subtraction of the numbers as integers.

Following the steps shown in the article, I set the format of x to:

x = M * 2 ^ -16

In an attempt to perform the same logic, I try to define Y0 for:

Y0 = (R2 - M / 2) * 2 ^ (R1 - (-16/2));

I'm trying to find a number, which can minimize the error given by:

error = (Y - Y0) / Y

Regardless of the value of R1, I can do shift operations to correct the exponent value of my final result, having the correct result at a fixed point.

Where am I wrong?

John Kugelman
  • 307,513
  • 65
  • 473
  • 519
vicente cesar
  • 41
  • 1
  • 3
  • Perhaps this question belongs on the math exchange if it's math that is the problem. – Christian Gibbons Apr 08 '21 at 16:20
  • 4
    @ChristianGibbons: Certainly not; the “Fast Inverse Square Root” uses specifics of the IEEE-754 binary32 floating-point representation as well as characteristics of floating-point arithmetic, and the task of adapting the algorithm to fixed-point arithmetic is a software engineering issue. – Eric Postpischil Apr 08 '21 at 16:48

1 Answers1

3

It can't be done.

The fast inverse sqrt is due to the floating point representation, that has already split the number into powers of two (exponent) and the significant.

It can be done.

With the same tricks as done for floating points, it's possible to convert your fixed point into 2^exp * x. Given uint32_t a, uint8_t exp = bias- builtin_count_leading_zeros(a); uint32_t b = a << exp, with the constants (and domain of a) so carefully chosen, that there will be no underflows or overflows.

Thus, you will actually have a custom floating point representation, which is tailored for this specific purpose, omitting the sign bit at least and having the best possible number of bits for the exponent, which might as well be 8.

Aki Suihkonen
  • 15,929
  • 1
  • 30
  • 50
  • I saw a technique that used this same principle, however I will implement this algorithm in an FPGA in the future, and "builtin_count_leading_zeros" would not be a good thing. This would involve a "Priority Encoder", and as far as I know it will limit the clock speed that I could use on the circuit. Is this algorithm the wrong way to go? – vicente cesar Apr 08 '21 at 17:55
  • 1
    @vincente cesar Would have been better to mention that in the question. We don't know the specifics of your FPGA. You could try a bit-by-bit algorithm, either restoring or non-restoring variant, and maybe pipeline it. If you have a count-leading-zeros circuit, multipliers, and a bit of ROM space, you could try a Newton-Raphson based algorithm as I showed in [this answer](https://stackoverflow.com/a/32337283/780717) – njuffa Apr 08 '21 at 18:17
  • The method from [John Carmack's, Fast Inverse Square Root](https://en.wikipedia.org/wiki/Fast_inverse_square_root) utilizing carefully chosen constants is discussed there. The beauty is the solution generally converges within 3 iterations making it quite efficient within the chosen accuracy. – David C. Rankin Apr 08 '21 at 19:04
  • 1
    In FPGA, I would probably normalise the value by shifting left by 16,8,4,2,1 bits at a time, getting one bit for the exponent at each iteration. It possibly doesn't limit the clock speed but adds to latency. – Aki Suihkonen Apr 09 '21 at 03:15
  • @AkiSuihkonen I am trying to study the algorithm on the link. Is it possible to adapt it to calculate the reciprocal of the number? I am still trying to understand. I'm reading as many solutions as possible and its algorithm looks the best both for calculating the reciprocal and the reciprocal of the square root. I am using an FPGA Spartan 6 xc6slx16. – vicente cesar Apr 09 '21 at 16:48
  • Yes, a reciprocal of a number of 2^n *(1+m) is approximately 2^-n * (1-m). – Aki Suihkonen Apr 09 '21 at 17:41