Существует ли нециклическая беззнаковая 32-битная целочисленная функция квадратного корня C

Я видел взлом битов с плавающей запятой для получения квадратного корня, как показано здесь быстрый квадратный корень с плавающей запятой, но этот метод работает для чисел с плавающей запятой.

Есть ли аналогичный метод для нахождения целого квадратного корня без циклов 32-битного целого числа без знака? Я искал в Интернете один, но не видел ни одного

(я думаю, что чистое двоичное представление не имеет достаточно информации для этого, но, поскольку оно ограничено 32-битным, я бы предположил иначе)


person yosmo78    schedule 01.02.2021    source источник
comment
Конечно, вы можете вычислить максимальное количество необходимых итераций и полностью развернуть цикл.   -  person Nate Eldredge    schedule 01.02.2021
comment
uint32_t y = sqrt(x); выполнит эту работу путем преобразования в double, извлечения квадратного корня и обратного преобразования.   -  person Eric Postpischil    schedule 01.02.2021


Ответы (2)


Этот ответ предполагает, что целевая платформа не поддерживает операции с плавающей запятой или поддерживает очень медленную работу с плавающей запятой (возможно, посредством эмуляции).

Как было указано в комментариях, инструкция подсчета ведущих нулей (CLZ) может использоваться для обеспечения быстрой функциональности log2, которая обеспечивается через экспоненциальную часть операндов с плавающей запятой. CLZ также можно с достаточной эффективностью эмулировать на платформах, которые не предоставляют функциональные возможности через встроенные функции, как показано ниже.

Начальное приближение, хорошее для нескольких битов, можно получить из таблицы поиска (LUT), которую можно дополнительно уточнить с помощью итераций Ньютона, как и в случае с плавающей запятой. Для извлечения квадратного корня из 32-битного целого числа обычно достаточно одной-двух итераций. Приведенный ниже код ISO-C99 показывает работающую примерную реализацию, включая исчерпывающий тест.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

uint8_t clz (uint32_t a); // count leading zeros
uint32_t umul_16_16 (uint16_t a, uint16_t b); // 16x16 bit multiply
uint16_t udiv_32_16 (uint32_t x, uint16_t y); // 32/16 bit division

/* LUT for initial square root approximation */
static const uint16_t sqrt_tab[32] = 
{ 
    0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
    0x85ff, 0x8cff, 0x94ff, 0x9aff, 0xa1ff, 0xa7ff, 0xadff, 0xb3ff,
    0xb9ff, 0xbeff, 0xc4ff, 0xc9ff, 0xceff, 0xd3ff, 0xd8ff, 0xdcff, 
    0xe1ff, 0xe6ff, 0xeaff, 0xeeff, 0xf3ff, 0xf7ff, 0xfbff, 0xffff
};

/* table lookup for initial guess followed by division-based Newton iteration */
uint16_t my_isqrt (uint32_t x)
{
    uint16_t q, lz, y, i, xh;

    if (x == 0) return x; // early out, code below can't handle zero

    // initial guess based on leading 5 bits of argument normalized to 2.30
    lz = clz (x);
    i = ((x << (lz & ~1)) >> 27);
    y = sqrt_tab[i] >> (lz >> 1);
    xh = x >> 16; // use for overflow check on divisions

    // first Newton iteration, guard against overflow in division
    q = 0xffff;
    if (xh < y) q = udiv_32_16 (x, y);
    y = (q + y) >> 1;

    if (lz < 10) {
        // second Newton iteration, guard against overflow in division
        q = 0xffff;
        if (xh < y) q = udiv_32_16 (x, y);
        y = (q + y) >> 1;
    }

    if (umul_16_16 (y, y) > x) y--; // adjust quotient if too large

    return y; // (uint16_t)sqrt((double)x)
}

static const uint8_t clz_tab[32] = 
{
    31, 22, 30, 21, 18, 10, 29,  2, 20, 17, 15, 13, 9,  6, 28, 1,
    23, 19, 11,  3, 16, 14,  7, 24, 12,  4,  8, 25, 5, 26, 27, 0
};

/* count leading zeros (for non-zero argument); a machine instruction on many architectures */
uint8_t clz (uint32_t a)
{
    a |= a >> 16;
    a |= a >> 8;
    a |= a >> 4;
    a |= a >> 2;
    a |= a >> 1;
    return clz_tab [0x07c4acdd * a >> 27];
}

/* 16x16->32 bit unsigned multiply; machine instruction on many architectures */
uint32_t umul_16_16 (uint16_t a, uint16_t b)
{
    return (uint32_t)a * b;
}

/* 32/16->16 bit division. Note: Will overflow if x[31:16] >= y */
uint16_t udiv_32_16 (uint32_t x, uint16_t y)
{
    uint16_t r = x / y;
    return r;
}

int main (void)
{
    uint32_t x;
    uint16_t res, ref;
    
    printf ("testing 32-bit integer square root\n");
    x = 0;
    do {
        ref = (uint16_t)sqrt((double)x);
        res = my_isqrt (x);
        if (res != ref) {
            printf ("error: x=%08x  res=%08x  ref=%08x\n", x, res, ref);
            printf ("exhaustive test FAILED\n");
            return EXIT_FAILURE;
        }
        x++;
    } while (x);
    printf ("exhaustive test PASSED\n");
    return EXIT_SUCCESS;
}
person njuffa    schedule 01.02.2021
comment
Что ж, впечатляет. - person Joshua; 01.02.2021

Нет. Вам нужно где-то ввести журнал; быстрый квадратный корень с плавающей запятой работает из-за журнала в битовом представлении.

Вероятно, самым быстрым методом является таблица поиска n -> floor(sqrt(n)). Вы не сохраняете все значения в таблице, а только значения, для которых изменяется квадратный корень. Используйте бинарный поиск, чтобы найти результат в таблице за время log(n).

person Joshua    schedule 01.02.2021
comment
С другой стороны, журналы целых чисел не сложны, если у вас есть внутренний счет с ведущими нулями. - person Nate Eldredge; 01.02.2021
comment
Бинарный поиск без зацикливания? - person Ajay Brahmakshatriya; 01.02.2021
comment
@AjayBrahmakshatriya: не более 16 поисковых запросов; вы можете развернуть его. - person Joshua; 01.02.2021