оптимизируя с фиксированной точкой Sqrt

Question

Я сделал то, что я думаю, что это хорошая 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); 
} 

Однако, после того, как профилирование (в режиме выпуска) я обнаружил, что деление здесь занимает 83% времени этот алгоритм проводит. Я могу ускорить его 6 раз, заменив деление на умножение, но это просто неправильно. Я узнал, что целочисленное деление гораздо медленнее, чем умножение, к сожалению. Есть ли способ оптимизировать это?

В случае, если эта таблица необходима.

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 просто SQRT_32_EVEN_LOOKUP делится на sqrt(2).

Answer 1

Переосмыслить колесо, действительно, и не в хорошем смысле. Правильное решение состоит в том, чтобы вычислить 1/sqrt(x) с использованием NR, а затем умножить один раз, чтобы получить x/sqrt(x) - просто проверьте на x==0 спереди.

Причина, по которой это намного лучше, заключается в том, что шаг NR для y=1/sqrt(x) равен y = (3y-x*y*y)*y/2. Это простое умножение.