Самый быстрый Integer Square Root в наименьшем количестве инструкций

Ответ 1

Посмотрите здесь.

Например, в 3 (a) существует этот метод, который тривиально адаптируется для 64- > 32-битного квадратного корня, а также тривиально транскрибируется ассемблеру:

/* by Jim Ulery */
static unsigned julery_isqrt(unsigned long val) {
    unsigned long temp, g=0, b = 0x8000, bshft = 15;
    do {
        if (val >= (temp = (((g << 1) + b)<<bshft--))) {
           g += b;
           val -= temp;
        }
    } while (b >>= 1);
    return g;
}

Нет делений, нет умножений, бит сдвигов. Тем не менее, время будет несколько непредсказуемым, особенно если вы используете ветку (в условных инструкциях ARM RISC будет работать).

В общем, эта страница перечисляет способы вычисления квадратных корней. Если вам захочется создать быстрый обратный квадратный корень (т.е. x**(-0.5)) или просто заинтересованы в удивительных способах оптимизации кода, посмотрите этот, и .

Ответ 2

Это то же самое, что и у вас, но с меньшим количеством опций. (Я рассчитываю 9 операций в цикле в вашем коде, включая тест и приращение i в цикле for и 3 назначениях, но, возможно, некоторые из них исчезают при кодировании в ASM? В коде ниже 6 операций, если вы считаете g*g>n как два (без присвоения)).

int isqrt(int n) {
  int g = 0x8000;
  int c = 0x8000;
  for (;;) {
    if (g*g > n) {
      g ^= c;
    }
    c >>= 1;
    if (c == 0) {
      return g;
    }
    g |= c;
  }
}

Я получил здесь. Вы можете исключить сравнение, если вы разворачиваете цикл и переходите в соответствующее место на основе самого высокого ненулевого бита на входе.

Обновление

Я больше думал об использовании метода Ньютона. Теоретически число бит точности должно удваиваться для каждой итерации. Это означает, что метод Ньютона намного хуже любого другого предложения, когда в ответе есть несколько правильных бит; однако ситуация меняется, когда в ответе есть много правильных бит. Учитывая, что большинство предложений, по-видимому, занимают 4 цикла на бит, это означает, что одна итерация метода Ньютона (16 циклов для деления + 1 для добавления + 1 для сдвига = 18 циклов) не стоит, если она не дает более 4 бит.

Итак, мое предложение состоит в том, чтобы построить 8 бит ответа одним из предложенных методов (8 * 4 = 32 цикла), затем выполнить одну итерацию метода Ньютона (18 циклов), чтобы удвоить количество бит до 16. Это всего 50 циклов (плюс, возможно, дополнительные 4 цикла, чтобы получить 9 бит до применения метода Ньютона, плюс, возможно, 2 цикла для преодоления перерегулирования, иногда испытываемого методом Ньютона). Это максимум 56 циклов, которые, насколько я вижу, соперничают с любыми другими предложениями.

Второе обновление

Я закодировал идею гибридного алгоритма. Метод Ньютона не имеет накладных расходов; вы просто применяете и удваиваете число значимых цифр. Проблема заключается в том, чтобы иметь предсказуемое количество значимых цифр, прежде чем применять метод Ньютона. Для этого нам нужно выяснить, где появится самый старший бит ответа. Используя модификацию быстрого метода последовательности DeBruijn, заданного другим плакатом, мы можем выполнить этот расчет примерно за 12 циклов в моей оценке. С другой стороны, знание позиции msb ответа ускоряет все методы (средний, не худший случай), поэтому кажется целесообразным независимо от того, что.

После вычисления msb ответа я запускаю несколько раундов предложенного выше алгоритма, а затем завершаю его одним или двумя раундами метода Ньютона. Мы решаем, когда запускать метод Ньютона следующим вычислением: один бит ответа занимает около 8 циклов в соответствии с расчетом в комментариях; один раунд метода Ньютона занимает около 18 циклов (деление, сложение и сдвиг и, возможно, назначение), поэтому мы должны запускать только метод Ньютона, если мы собираемся извлечь из него как минимум три бита. Итак, для 6-битных ответов мы можем запустить линейный метод 3 раза, чтобы получить 3 значащих бита, а затем запустить метод Ньютона 1 раз, чтобы получить еще 3. Для 15-битных ответов мы запускаем линейный метод 4 раза, чтобы получить 4 бита, затем Ньютон метод дважды, чтобы получить еще 4, затем еще 7. И так далее.

Эти вычисления зависят от того, как точно знать, сколько циклов требуется, чтобы получить бит по линейному методу, и сколько из них требуется по методу Ньютона. Если "экономика" изменяется, например, путем обнаружения более быстрого способа создания битов линейным образом, решение о том, когда вызывать метод Ньютона, изменится.

Я развернул петли и реализовал решения как переключатели, которые, я надеюсь, перейдут в быстрый поиск таблиц в сборке. Я не совсем уверен, что у меня есть минимальное количество циклов в каждом случае, так что возможно дальнейшая настройка возможна. Например, для s = 10 вы можете попытаться получить 5 бит, а затем применить метод Ньютона один раз вместо двух.

Я тщательно протестировал алгоритм для правильности. Некоторые дополнительные незначительные ускорения возможны, если вы согласитесь принять несколько неверные ответы в некоторых случаях. По крайней мере два цикла используются после применения метода Ньютона для исправления ошибки, которая возникает с номерами формы m^2-1. И цикл используется для тестирования ввода 0 в начале, так как алгоритм не может обрабатывать этот ввод. Если вы знаете, что никогда не возьмете квадратный корень из нуля, вы можете устранить этот тест. Наконец, если вам нужно только 8 значащих бит в ответе, вы можете отказаться от одного из расчетов метода Ньютона.

#include <inttypes.h>
#include <stdint.h>
#include <stdbool.h>
#include <stdio.h>

uint32_t isqrt1(uint32_t n);

int main() {
  uint32_t n;
  bool it_works = true;
  for (n = 0; n < UINT32_MAX; ++n) {
    uint32_t sr = isqrt1(n);
    if ( sr*sr > n || ( sr < 65535 && (sr+1)*(sr+1) <= n )) {
      it_works = false;
      printf("isqrt(%" PRIu32 ") = %" PRIu32 "\n", n, sr);
    }
  }
  if (it_works) {
    printf("it works\n");
  }
  return 0;
}

/* table modified to return shift s to move 1 to msb of square root of x */
/*
static const uint8_t debruijn32[32] = {
    0, 31, 9, 30, 3,  8, 13, 29,  2,  5,  7, 21, 12, 24, 28, 19,
    1, 10, 4, 14, 6, 22, 25, 20, 11, 15, 23, 26, 16, 27, 17, 18
};
*/

static const uint8_t debruijn32[32] = {
  15,  0, 11, 0, 14, 11, 9, 1, 14, 13, 12, 5, 9, 3, 1, 6,
  15, 10, 13, 8, 12,  4, 3, 5, 10,  8,  4, 2, 7, 2, 7, 6
};

/* based on CLZ emulation for non-zero arguments, from
 * http://stackoverflow.com/questions/23856596/counting-leading-zeros-in-a-32-bit-unsigned-integer-with-best-algorithm-in-c-pro
 */
uint8_t shift_for_msb_of_sqrt(uint32_t x) {
  x |= x >>  1;
  x |= x >>  2;
  x |= x >>  4;
  x |= x >>  8;
  x |= x >> 16;
  x++;
  return debruijn32 [x * 0x076be629 >> 27];
}

uint32_t isqrt1(uint32_t n) {
  if (n==0) return 0;

  uint32_t s = shift_for_msb_of_sqrt(n);
  uint32_t c = 1 << s;
  uint32_t g = c;

  switch (s) {
  case 9:
  case 5:
    if (g*g > n) {
      g ^= c;
    }
    c >>= 1;
    g |= c;
  case 15:
  case 14:
  case 13:
  case 8:
  case 7:
  case 4:
    if (g*g > n) {
      g ^= c;
    }
    c >>= 1;
    g |= c;
  case 12:
  case 11:
  case 10:
  case 6:
  case 3:
    if (g*g > n) {
      g ^= c;
    }
    c >>= 1;
    g |= c;
  case 2:
    if (g*g > n) {
      g ^= c;
    }
    c >>= 1;
    g |= c;
  case 1:
    if (g*g > n) {
      g ^= c;
    }
    c >>= 1;
    g |= c;
  case 0:
    if (g*g > n) {
      g ^= c;
    }
  }

  /* now apply one or two rounds of Newton method */
  switch (s) {
  case 15:
  case 14:
  case 13:
  case 12:
  case 11:
  case 10:
    g = (g + n/g) >> 1;
  case 9:
  case 8:
  case 7:
  case 6:
    g = (g + n/g) >> 1;
  }

  /* correct potential error at m^2-1 for Newton method */
  return (g==65536 || g*g>n) ? g-1 : g;
}

При тестировании на моей машине (который, по общему признанию, не похож на ваш), новая процедура isqrt1 работает в среднем на 40% быстрее, чем предыдущая процедура isqrt, которую я дал.

Ответ 3

Если умножение является той же скоростью (или быстрее, чем!) добавлением и сдвигом, или если вам не хватает быстрой инструкции по сокращению по сумме, содержащейся в регистре, то следующее не будет полезно. В противном случае:

Вы вычисляете temp*temp заново в каждом цикле цикла, но temp = res | add, что совпадает с res + add, так как их биты не перекрываются, и (a) вы уже вычислили res*res на предыдущий цикл цикла и (b) add имеет специальную структуру (это всегда один бит). Таким образом, используя тот факт, что (a+b)^2 = a^2 + 2ab + b^2, и у вас уже есть a^2, а b^2 - это всего лишь один бит, сдвинутый в два раза влево, как одноразрядный b, а 2ab - это просто a сдвинуто влево на 1 позицию, чем местоположение одного бита в b, вы можете избавиться от умножения:

unsigned short int int_sqrt32(unsigned int x)
{
    unsigned short int res = 0;
    unsigned int res2 = 0;
    unsigned short int add = 0x8000;   
    unsigned int add2 = 0x80000000;   
    int i;
    for(i = 0; i < 16; i++)
    {
        unsigned int g2 = res2 + (res << i) + add2;
        if (x >= g2)
        {
            res |= add;
            res2 = g2;
        }
        add >>= 1;
        add2 >>= 2;
    }
    return res;
}

Также я бы предположил, что лучше использовать один и тот же тип (unsigned int) для всех переменных, поскольку в соответствии со стандартом C для всей арифметики требуется продвижение (преобразование) более узких целочисленных типов к самому широкому типу, включенному до выполняется арифметическая операция, а затем последующее обратное преобразование, если необходимо. (Разумеется, это может быть оптимизировано достаточно интеллектуальным компилятором, но зачем рисковать?)

Ответ 4

Из комментариев след, кажется, что RISC-процессор обеспечивает только 32x32- > 32-битное умножение и 16x16- > 32-битное умножение. Расширенное 32x-32- > 64-битное умножение или инструкция MULHI, возвращающая верхние 32 бита 64-разрядного продукта, не предоставляется.

Это, по-видимому, исключает подходы, основанные на итерации Ньютона-Рафсона, которые, вероятно, будут неэффективными, поскольку они обычно требуют либо инструкции MULHI, либо расширения умножения для промежуточной арифметики с фиксированной точкой.

В приведенном ниже коде C99 используется итеративный подход, требующий только умножения 16x16- > 32 бит, но сходится несколько линейно, требуя до шести итераций. Для этого подхода требуется CLZ функциональность, чтобы быстро определить исходное предположение для итераций. Аскер заявил в комментариях, что используемый RISC-процессор не обеспечивает функциональность CLZ. Таким образом, требуется эмуляция CLZ, и поскольку эмуляция добавляет к количеству хранения и подсчета команд, это может сделать этот подход неконкурентоспособным. Я выполнил поиск грубой силы, чтобы определить таблицу поиска deBruijn с наименьшим множителем.

Этот итерационный алгоритм обеспечивает исходные результаты, близкие к желаемым результатам, т.е. (int)sqrt(x), но всегда несколько с высокой стороны из-за усеченного характера целочисленной арифметики. Чтобы достигнуть конечного результата, результат условно уменьшается, пока квадрат результата не будет меньше или равен исходному аргументу.

Использование квалификатора volatile в коде только для того, чтобы установить, что все именованные переменные фактически могут быть распределены как 16-битные данные, не влияя на функциональность. Я не знаю, дает ли это какое-либо преимущество, но заметил, что OP специально использовал 16-битные переменные в своем коде.

Обратите внимание, что для большинства процессоров шаги коррекции в конце не должны включать какого-либо разветвления. Продукт y*y можно вычесть из x с выполнением (или заимствованием), затем y корректируется вычитанием с переносом (или заимствованием). Таким образом, каждый шаг должен быть последовательностью MUL, SUBcc, SUBC.

Поскольку реализация итерации по циклу несет существенные накладные расходы, я решил полностью развернуть цикл, но предоставить две ранние проверки. Выполняя подсчет операций вручную, я рассчитываю 46 операций для самого быстрого случая, 54 операции для среднего случая и 60 операций для наихудшего случая.

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

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};

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];
}

/* 16 x 16 -> 32 bit unsigned multiplication; should be single instruction */
uint32_t umul16w (uint16_t a, uint16_t b)
{
    return (uint32_t)a * b;
}

/* Reza Hashemian, "Square Rooting Algorithms for Integer and Floating-Point
   Numbers", IEEE Transactions on Computers, Vol. 39, No. 8, Aug. 1990, p. 1025
*/
uint16_t isqrt (uint32_t x)
{
    volatile uint16_t y, z, lsb, mpo, mmo, lz, t;

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

    lz = clz (x);         // #leading zeros, 32-lz = #bits of argument
    lsb = lz & 1;
    mpo = 17 - (lz >> 1); // m+1, result has roughly half the #bits of argument
    mmo = mpo - 2;        // m-1
    t = 1 << mmo;         // power of two for two complement of initial guess
    y = t - (x >> (mpo - lsb)); // initial guess for sqrt
    t = t + t;            // power of two for two complement of result
    z = y;

    y = (umul16w (y, y) >> mpo) + z;
    y = (umul16w (y, y) >> mpo) + z;
    if (x >= 0x40400) {
        y = (umul16w (y, y) >> mpo) + z;
        y = (umul16w (y, y) >> mpo) + z;
        if (x >= 0x1002000) {
            y = (umul16w (y, y) >> mpo) + z;
            y = (umul16w (y, y) >> mpo) + z;
        }
    }

    y = t - y; // raw result is 2 complement of iterated solution
    y = y - umul16w (lsb, (umul16w (y, 19195) >> 16)); // mult. by sqrt(0.5) 

    if ((int32_t)(x - umul16w (y, y)) < 0) y--; // iteration may overestimate 
    if ((int32_t)(x - umul16w (y, y)) < 0) y--; // result, adjust downward if 
    if ((int32_t)(x - umul16w (y, y)) < 0) y--; // necessary 

    return y; // (int)sqrt(x)
}

int main (void)
{
    uint32_t x = 0;
    uint16_t res, ref;

    do {
        ref = (uint16_t)sqrt((double)x);
        res = isqrt (x);
        if (res != ref) {
            printf ("!!!! x=%08x  res=%08x  ref=%08x\n", x, res, ref);
            return EXIT_FAILURE;
        }
        x++;
    } while (x);
    return EXIT_SUCCESS;
}

Другая возможность - использовать итерацию Ньютона для квадратного корня, несмотря на высокую стоимость деления. Для небольших входов потребуется только одна итерация. Несмотря на то, что об этом не спрашивал об этом, на основании времени выполнения 16 циклов для операции DIV я сильно подозреваю, что это фактически бит-бит 32/16->16, который требует дополнительного защитного кода, чтобы избежать переполнения, определенного как частное, которое делает не вписывается в 16 бит. Я добавил соответствующие гарантии для моего кода на основе этого предположения.

Так как итерация Newton удваивает количество хороших бит каждый раз, когда она применяется, нам нужна только начальная догадка с низкой точностью, которую можно легко извлечь из таблицы на основе пяти ведущих бит аргумента. Чтобы перенести их, мы сначала нормализуем аргумент в 2.30 формат с фиксированной запятой с дополнительным неявным масштабным коэффициентом 2 ^{32- (lz и ~ 1)} где lz - число начальных нулей в аргументе. Как и в предыдущем подходе, итерация не всегда дает точный результат, поэтому необходимо приложить корректировку, если предварительный результат будет слишком большим. Я считаю 49 циклов для быстрого пути, 70 циклов для медленного пути (в среднем 60 циклов).

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
};

/* 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;
}

/* table lookup for initial guess followed by division-based Newton iteration*/ 
uint16_t isqrt (uint32_t x)
{
    volatile 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); // needed for overflow check on division

    // 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 (umul16w (y, y) > x) y--; // adjust quotient if too large

    return y; // (int)sqrt(x)
}

Ответ 5

Вот менее инкрементная версия описанной техники @j_random_hacker. По крайней мере, на одном процессоре это было немного быстрее, когда я возился с этим пару лет назад. Понятия не имею почему.

// assumes unsigned is 32 bits
unsigned isqrt1(unsigned x) {
  unsigned r = 0, r2 = 0; 
  for (int p = 15; p >= 0; --p) {
    unsigned tr2 = r2 + (r << (p + 1)) + (1u << (p + p));
    if (tr2 <= x) {
      r2 = tr2;
      r |= (1u << p);
    }
  }
  return r;
}

/*
gcc 6.3 -O2
isqrt(unsigned int):
        mov     esi, 15
        xor     r9d, r9d
        xor     eax, eax
        mov     r8d, 1
.L3:
        lea     ecx, [rsi+1]
        mov     edx, eax
        mov     r10d, r8d
        sal     edx, cl
        lea     ecx, [rsi+rsi]
        sal     r10d, cl
        add     edx, r10d
        add     edx, r9d
        cmp     edx, edi
        ja      .L2
        mov     r11d, r8d
        mov     ecx, esi
        mov     r9d, edx
        sal     r11d, cl
        or      eax, r11d
.L2:
        sub     esi, 1
        cmp     esi, -1
        jne     .L3
        rep ret
*/

Если вы включите оптимизацию gcc 9 x86, он полностью развернет цикл и свернет константы. В результате осталось всего около 100 инструкций.

Ответ 6

Я не знаю, как превратить его в эффективный алгоритм, но когда я исследовал это в 80-х годах, появился интересный шаблон. При округлении квадратных корней есть еще два целых числа с этим квадратным корнем, чем предыдущий (после нуля).

Итак, одно число (ноль) имеет квадратный корень из нуля, два имеют квадратный корень из 1 (1 и 2), 4 имеют квадратный корень из двух (3, 4, 5 и 6) и т.д. Вероятно, это не полезный ответ, но тем не менее интересный.