Быстрое евклидово деление на С

Question

Я заинтересован в получении остатка от евклидова деления , то есть для пары целых чисел (i, n) найдите r, например:

i = k * n + r, 0 <= r < |k|

простое решение:

int euc(int i, int n)
{
    int r;

    r = i % n;
    if ( r < 0) {
        r += n;
    }
    return r;
}

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

• Разветвление, но также желательно быстрее.
• Приемлемо решение, которое работает только для положительного n (но оно должно работать для отрицательного i).
• n заранее неизвестен и может принимать любое значение> 0 и

Редактировать

На самом деле довольно легко ошибиться в результате, поэтому вот пример ожидаемых результатов:

• euc (0, 3) = 0
• euc (1, 3) = 1
• euc (2, 3) = 2
• euc (3, 3) = 0
• euc (-1, 3) = 2
• euc (-2, 3) = 1
• euc (-3,3) = 0

Некоторые люди также беспокоятся, что не имеет смысла оптимизировать это. Мне это нужно для многомерного итератора, в котором элементы за пределами границ заменяются элементами в «виртуальном массиве», который повторяет исходный массив. Поэтому, если мой массив x равен [1, 2, 3, 4], виртуальный массив равен [...., 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4] и, например, x [-2] равен x 1 и т. Д. *

Для второго массива измерения d мне нужно евклидово деление d для каждой точки. Если мне нужно сделать корреляцию между массивом n ^ d и ядром m ^ d, мне нужно евклидово деление n ^ d * m ^ d * d. Для трехмерного изображения 100x100x100 точек и ядра 5 * 5 * 5 точек это уже ~ 400 миллионов евклидовых делений.

Answer 1

Редактировать: Нет умножения или веток.

int euc(int i, int n)
{
    int r;

    r = i % n;
    r += n & (-(r < 0));

    return r;
}

Вот сгенерированный код. Согласно профилировщику инструментов MSVC ++ (мое тестирование) и тестированию OP, они работают почти одинаково.

; Original post
00401000  cdq              
00401001  idiv        eax,ecx 
00401003  mov         eax,edx 
00401005  test        eax,eax 
00401007  jge         euc+0Bh (40100Bh) 
00401009  add         eax,ecx 
0040100B  ret              

; Mine
00401020  cdq              
00401021  idiv        eax,ecx 
00401023  xor         eax,eax 
00401025  test        edx,edx 
00401027  setl        al   
0040102A  neg         eax  
0040102C  and         eax,ecx 
0040102E  add         eax,edx 
00401030  ret              

Answer 2

Я думаю, что у 280Z28 и Christopher ассемблерный гольф покрыт лучше, чем я, и это касается случайного доступа.

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

В этом случае, сначала применив подходящую проверку границ, вы можете выполнить внутренний цикл в том, что я буду называть «тире». Убедитесь, что следующие k приращений не приводят к переполнению наименьшего измерения в любом массиве, а затем «разбейте» k шагов, используя новый, еще более внутренний цикл, который просто увеличивает «физический» индекс на 1 каждый раз вместо делать другой идив. Вы или компилятор можете развернуть этот цикл, использовать устройство Даффа и т. Д.

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

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

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

В качестве альтернативы, если выполняемая операция представляет собой числовое равенство или какую-либо другую простую операцию (я не знаю, что такое «корреляция»), вы можете посмотреть SIMD-инструкции или что-то еще, в этом случае длина тире должна быть кратные из самых широких сравнений с одной инструкцией (или соответствующей операции SIMD) в вашей архитектуре. Это не то, с чем у меня есть опыт.

Answer 3

Без ветки, но немного суетливо:

int euc2(int i, int n)
{
    int r;
    r = i % n;
    r += (((unsigned int)r) >> 31) * n;
    return r;
}

Без умножения:

int euc2(int i, int n)
{
    int r;
    r = i % n;
    r += (r >> 31) & n;
    return r;
}

Это дает:

; _i$ = eax
; _n$ = ecx

cdq
idiv   ecx
mov eax, edx
sar eax, 31
and eax, ecx
add eax, edx

Answer 4

Целочисленное умножение намного быстрее деления. Для большого количества вызовов с известным N вы можете заменить деление на N умножением на псевдообратное значение N.

Я проиллюстрирую это на примере. Возьми N = 29. Затем вычислите один раз псевдообратное 2 ^ 16 / N: K = 2259 (усечено из 2259.86 ...). Я предполагаю, что я положительный, и я * K подходит для 32 бит.

Quo = (I*K)>>16;   // replaces the division, Quo <= I/N
Mod = I - Quo*N;   // Mod >= I%N
while (Mod >= N) Mod -= N;  // compensate for the approximation

В моем примере возьмем I = 753, получим Quo = 25 и Mod = 28. (компенсация не требуется)

EDIT.

В вашем примере трехмерной свертки большинство вызовов i% n будут с i в 0..n-1, поэтому в большинстве случаев это первая строка, такая как

if (i>=0 && i<n) return i;

обойдёт дорогостоящий и здесь бесполезный идив.

Кроме того, если у вас достаточно ОЗУ, просто выровняйте все измерения по степеням 2 и используйте битовые манипуляции (сдвиг и) вместо делений.

РЕДАКТИРОВАТЬ 2.

Я действительно попробовал это на 10 ^ 9 звонках. я% n: 2.93s, мой код: 1.38s. Просто имейте в виду, что это подразумевает ограничение на I (I * K должно соответствовать 32 битам).

Еще одна мысль: если ваши значения x + dx, с x в 0..n-1 и dx small, то следующее будет охватывать все случаи:

if (i<0) return i+n; else if (i>=n) return i-n;
return i;

Answer 5

Если у вас достаточно низкий диапазон, создайте таблицу соответствия - два dim-массива.Также вы можете сделать функцию Inline и убедиться в этом, посмотрев полученный код.

Answer 6

Мне очень нравится выражение:

r = ((i%n)+n)%n; 

Разборка очень короткая:

r = ((i% n) + n)% n;

004135AC  mov         eax,dword ptr [i] 
004135AF  cdq              
004135B0  idiv        eax,dword ptr [n] 
004135B3  add         edx,dword ptr [n] 
004135B6  mov         eax,edx 
004135B8  cdq              
004135B9  idiv        eax,dword ptr [n] 
004135BC  mov         dword ptr [r],edx 

У него нет переходов (2 idivs, которые могут быть дорогостоящими), и он может быть полностью встроенным, избегая затрат на вызов функции.

Что ты думаешь?

Answer 7

Я рассчитал все предложения в gcc -O3, используя TSC (кроме одного для постоянной N), и все они заняли одинаковое количество времени (в пределах 1%).

Я думал, что либо ((i% n) + n)% n (без разветвления), либо (i + (n << 16))% n (очевидно, что сбой при большом n или крайне отрицательном значении i) будет быстрее , но все они заняли одно и то же время. </p>

Answer 8

int euc(int i, int n)
{
    return (i % n) + (((i % n) < 0) * n);
}

Answer 9

Вот версия Кристофера с отступом до Джейсона , если смещение вправо не арифметическое.

#include <limits.h>
static inline int euc(int i, int n)
{
    // check for arithmetic shift
    #if (-1 >> 1) == -1
        #define OFFSET ((i % n >> (sizeof(int) * CHAR_BIT - 1)) & n)
    #else
        #define OFFSET ((i % n < 0) * n)
    #endif

    return i % n + OFFSET;
}

Резервная версия должна быть медленнее, поскольку она использует imul вместо and.

Answer 10

Если вы можете гарантировать, что размеры вашего массива всегда степени двух, то вы можете сделать это:

r = (i & (n - 1));

Если вы можете дополнительно гарантировать, что ваши размеры будут из заданного подмножества, вы можете сделать:

template<int n>
int euc(int i) {
    return (i & (n - 1));
}

int euc(int i, int n) {
    switch (n) {
        case 2: return euc<2>(i);
        case 4: return euc<4>(i);
    }
}