Уловка деления константы (степени двойки) на целое число

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

• Ничто не может сравниться с предварительно вычисленной таблицей поиска по скорости.
• Существуют быстрые алгоритмы приближения квадратного корня (вероятно, уже в вашем DSP), и вы можете улучшить приближение одной или двумя итерациями Ньютона-Рафсона. Если вычисления с числами с плавающей запятой достаточно точны для вас, вы, вероятно, сможете превзойти 64-битное целочисленное деление с точки зрения скорости (но не с точки зрения ясности кода).

Вы упомянули, что результат будет преобразован в числа с плавающей запятой позже, может быть полезно вообще не вычислять целочисленное деление, а использовать ваше оборудование с плавающей запятой.

Я разработал версию Matlab, используя арифметику с фиксированной точкой.

Этот метод предполагает, что целочисленная версия log2(x) может быть эффективно вычислена, что верно для dsPIC30 / 33F и TI C6000, у которых есть инструкция для определения наиболее значимой единицы целого числа.

По этой причине этот код сильно зависит от ISA и не может быть написан на переносимом / стандартном C и может быть улучшен с помощью таких инструкций, как умножение и сложение, умножение и сдвиг, поэтому я не буду пытаться переводить его на C.

nrdiv.m

function [ y ] = nrdiv( q, x, lut) 
                          % assume q>31, lut = 2^31/[1,1,2,...255]
    p2 = ceil(log2(x));   % available in TI C6000, instruction LMBD
                          % available in Microchip dsPIC30F/33F, instruction FF1L 
    if p2<8
        pre_shift=0;
    else
        pre_shift=p2-8;
    end                                  % shr = (p2-8)>0?(p2-8):0;

    xn = shr(x, pre_shift);              % xn  = x>>pre_shift;
    y  = shr(lut(xn), pre_shift);        % y   = lut[xn]>pre_shift; 
    y  = shr(y * (2^32 - y*x), 30);      % basic iteration
                                         % step up from q31 to q32
    y  = shr(y * (2^33 - y*x), (64-q));  % step up from q32 to desired q
    if q>39
        y = shr(y * (2^(1+q) - y*x), (q));  % when q>40, additional 
                                            % iteration is required, 
    end                                     % no step up is performed
end
function y = shr(x, r)
    y=floor(x./2^r);             % simulate operator >>
end

test.m

test_number = (2^22-12345);
test_q      = 48;

lut_q31 = round(2^31 ./ [1,[1:1:255]]);
display(sprintf('tested 2^%d/%d, diff=%f\n',test_q, test_number,...
                 nrdiv( 39, (2^22-5), lut_q31) - 2^39/(2^22-5)));

образец вывода

tested 2^48/4181959, diff=-0.156250

ссылка:

деление Ньютона – Рафсона

Немного поздно, но вот мое решение.

Сначала несколько предположений:

Проблема:

X = N / D, где N - постоянная и степень двойки.

Все 32-битные беззнаковые целые числа.

X неизвестен, но у нас есть хорошая оценка (предыдущее, но уже не точное решение).

Точного решения не требуется.

Примечание: из-за целочисленного усечения это неточный алгоритм!

Итеративное решение - это нормально (улучшается с каждым циклом).

Деление намного дороже умножения:

Для 32-битного целого числа без знака для Arduino UNO:

'+/-' ~0.75us

'*' ~3.5us

'/' ~ 36us 4 Мы стремимся заменить в основном метод Ньютона:

Xnew=Xold-f(x)/(f`(x)

где f (x) = 0 для искомого решения.

Решив это, я получаю:

Xnew=XNew*(C-X*D)/N

где C = 2 * N

Первый трюк:

Теперь, когда числитель (константа) теперь является делителем (константой), одно решение здесь (которое не требует, чтобы N было степенью 2):

Xnew=XNew*(C-X*D)*A>>M

где C = 2 * N, A и M - константы (ищите деление на постоянные трюки).

или (оставаясь с методом Ньютона):

Xnew=XNew*(C-X*D)>>M

где C = 2 >> M, где M - мощность.

Итак, у меня есть 2 '*' (7.0us), '-' (0.75us) и '>>' (0.75us?) Или всего 8,5us (а не 36us), исключая другие накладные расходы.

Ограничения:

Поскольку тип данных - 32-битный беззнаковый, 'M' не должен превышать 15, иначе возникнут проблемы с переполнением (вы, вероятно, можете обойти это, используя 64-битный промежуточный тип данных).

N> D (иначе алгоритм взорвется! Хотя бы с целым числом без знака)

Очевидно, алгоритм будет работать со знаковыми и плавающими типами данных)

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
int main(void)
{
  unsigned long c,d,m,x;
  // x=n/d where n=1<<m
  m=15;
  c=2<<m;
  d=10;
  x=10;
  while (true)
  {
    x=x*(c-d*x)>>m;
    printf("%ld",x);
    getchar();
  }
  return(0);
}

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

/*
 * Converts the reference frequency count for a specific signal frequency
 * to a ratio.
 *   Xs = Ns * 2^32 / Nr
 *   Where:
 *   2^32 is a constant scaling so that the maximum accuracy can be achieved.
 *   Ns is the number of signal counts (fixed at 0x4000 by hardware).
 *   Nr is the number of reference counts, passed in W1:W0.
 * @param  W1:W0    The number of reference frequency pulses.
 * @return W1:W0    The scaled ratio.
 */
    .align  2
    .global _signalToReferenceRatio
    .type   _signalToReferenceRatio, @function

    ; This is the position of the most significant bit of the fixed Ns (0x4000).
    .equ    LOG2_DIVIDEND,  14
    .equ    DIVISOR_LIMIT,  LOG2_DIVIDEND+1
    .equ    WORD_SIZE,      16

_signalToReferenceRatio:
    ; Create a dividend, MSB-aligned with the divisor, in W2:W3 and place the
    ; number of iterations required for the MSW in [W14] and the LSW in [W14+2].
    LNK     #4
    MUL.UU  W2, #0, W2
    FF1L    W1, W4
    ; If MSW is zero the argument is out of range.
    BRA     C, .returnZero
    SUBR    W4, #WORD_SIZE, W4
    ; Find the number of quotient MSW loops.
    ; This is effectively 1 + log2(dividend) - log2(divisor).
    SUBR    W4, #DIVISOR_LIMIT, [W14]
    BRA     NC, .returnZero
    ; Since the SUBR above is always non-negative and the C flag set, use this
    ; to set bit W3<W5> and the dividend in W2:W3 = 2^(16+W5) = 2^log2(divisor).
    BSW.C   W3, W4
    ; Use 16 quotient LSW loops.
    MOV     #WORD_SIZE, W4
    MOV     W4, [W14+2]

    ; Set up W4:W5 to hold the divisor and W0:W1 to hold the result.
    MOV.D   W0, W4
    MUL.UU  W0, #0, W0

.checkLoopCount:
    ; While the bit count is non-negative ...
    DEC     [W14], [W14]
    BRA     NC,  .nextWord

.alignQuotient:
    ; Shift the current quotient word up by one bit.
    SL      W0, W0
    ; Subtract divisor from the current dividend part.
    SUB     W2, W4, W6
    SUBB    W3, W5, W7
    ; Check if the dividend part was less than the divisor.
    BRA     NC, .didNotDivide
    ; It did divide, so set the LSB of the quotient.
    BSET    W0, #0
    ; Shift the remainder up by one bit, with the next zero in the LSB.
    SL      W7, W3
    BTSC    W6, #15
    BSET    W3, #0
    SL      W6, W2
    BRA     .checkLoopCount
.didNotDivide:
    ; Shift the next (zero) bit of the dividend into the LSB of the remainder.
    SL      W3, W3
    BTSC    W2, #15
    BSET    W3, #0
    SL      W2, W2
    BRA     .checkLoopCount

.nextWord:
    ; Test if there are any LSW bits left to calculate.
    MOV     [++W14], W6
    SUB     W6, #WORD_SIZE, [W14--]
    BRA     NC, .returnQ
    ; Decrement the remaining bit counter before writing it back.
    DEC     W6, [W14]
    ; Move the working part of the quotient up into the MSW of the result.
    MOV     W0, W1
    BRA     .alignQuotient

.returnQ:
    ; Return the quotient in W0:W1.
    ULNK
    RETURN

.returnZero:
    MUL.UU  W0, #0, W0
    ULNK
    RETURN
.size   _signalToReferenceRatio, .-_signalToReferenceRatio