# 'Optimizing Fixed-Point Sqrt

I made what I think is a good `fixed-point square root algorithm`

:

```
template<int64_t M, int64_t P>
typename enable_if<M + P == 32, FixedPoint<M, P>>::type sqrt(FixedPoint<M, P> f)
{
if (f.num == 0)
return 0;
//Reduce it to the 1/2 to 2 range (based around FixedPoint<2, 30> to avoid left/right shift branching)
int64_t num{ f.num }, faux_half{ 1 << 29 };
ptrdiff_t mag{ 0 };
while (num < (faux_half)) {
num <<= 2;
++mag;
}
int64_t res = (M % 2 == 0 ? SQRT_32_EVEN_LOOKUP : SQRT_32_ODD_LOOKUP)[(num >> (30 - 4)) - (1LL << 3)];
res >>= M / 2 + mag - 1; //Finish making an excellent guess
for (int i = 0; i < 2; ++i)
// \ | /
// \ | /
// _| V L
res = (res + (int64_t(f.num) << P) / res) >> 1; //Use Newton's method to improve greatly on guess
// 7 A r
// / | \
// / | \
// The Infamous Time Eater
return FixedPoint<M, P>(res, true);
}
```

However, after profiling (in release mode) I found out that the division here takes up 83% of the time this algorithm spends. I can speed it up 6x by replacing the division with multiplication, but that's just wrong. I found out that integer division is much slower than multiplication, unfortunately. Is there any way to optimize this?

In case this table is necessary.

```
const array<int32_t, 24> SQRT_32_EVEN_LOOKUP = {
0x2d413ccd, //magic numbers calculated by taking 0.5 + 0.5 * i / 8 from i = 0 to 23, multiplying by 2^30, and converting to hex
0x30000000,
0x3298b076,
0x3510e528,
0x376cf5d1,
0x39b05689,
0x3bddd423,
0x3df7bd63,
0x40000000,
0x41f83d9b,
0x43e1db33,
0x45be0cd2,
0x478dde6e,
0x49523ae4,
0x4b0bf165,
0x4cbbb9d6,
0x4e623850,
0x50000000,
0x5195957c,
0x532370b9,
0x54a9fea7,
0x5629a293,
0x57a2b749,
0x59159016
};
```

`SQRT_32_ODD_LOOKUP`

is just `SQRT_32_EVEN_LOOKUP`

divided by `sqrt(2)`

.

## Solution 1:^{[1]}

Reinventing the wheel, really, and not in a good way. The correct solution is to calculate `1/sqrt(x)`

using NR, and then multiply once to get `x/sqrt(x)`

- just check for `x==0`

up front.

The reason why this is so much better is that the NR step for `y=1/sqrt(x)`

is just `y = (3-x*y*y)*y/2`

. That's all straightforward multiplication.

## Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution | Source |
---|---|

Solution 1 | BlackBeans |