аппроксимация деления, почему не используется Опции

[alexkarnaukhov]

У меня стоит задача нормировки потока чисел (поток пикселей изображения) на некоторое переменное значение (например разницу максимальной и минимальной яркости изображения). Соответственно без деления не обойтись. Те алгоритмы, которые используются чаще всего, делятся на два типа: либо тупо LUT с забитой таблицей 1/х, либо итеративные алгоритмы, имитирующие деление в столбик, с различной степенью оптимизации. Еще можно складывать делитель, пока результат не превзойдет делимое и посчитать количество таких операций, но это совсем плохой алгоритм.

Почему никто не использует тупо полиномиальную аппроксимацию 1/х? Вот здесь всего навсего полином 3-й степени, 100 слайсов на 3-ем Спартане и при этом точность 11 бит. Мне так кажется, что я здесь что-то не понимаю, это же так просто, почему так не делают?
Еще и конкретно для моей задачи - даже при точности 11 бит, я мог бы получить приемлемый результат, забив на линейность нормировки.

Спасибо.

[vladec]

Для вычисления 1/x есть еще итерационные алгоритмы, для них характерна очень быстрая сходимость, в том числе и на больших разрядностях

[Maverick]

Для вычисления 1/x есть еще итерационные алгоритмы, для них характерна очень быстрая сходимость, в том числе и на больших разрядностях

например? (хорошие по Вашему мнению - 32/64 бит)

[blackfin]

например? (хорошие по Вашему мнению - 32/64 бит)

Ну можно, наверное, методом Ньютона попробовать..

Главное - выбрать правильное начальное приближение..

Например, мы хотим найти частное от деления: A/B = D, где A = 555555, B = 11 и, соответственно, D = 50505.

Тогда:

A/B = [(A*2³²/B)>>32] = D.

Введем обозначение:

x = 2³²/B.

Это равносильно:

2³²/x = B.

Или:

f(x) == 2³²/x - B = 0.

Находим производную:

f'(x) = -2³²/x².

Теперь подставляем всё это в формулу Ньютона и находим:

x_n+1 = x_n - f(x_n)/f'(x_n), или

x_n+1 = x_n + [2³²/x_n - B]/[2³²/x_n²], или

x_n+1 = 2*x_n - B*x_n²/2³², или

x_n+1 = 2*x_n - [(B*x_n*x_n)>>32].

Начальное приближение для x₀ выбираем равным:

x₀ = 2^32-m, где m - максимальная степень двойки в двоичном представлении числа B..

То есть, для B = 11 это дает:

B = 11_d = 1011_b = 1*2³ + 0*2² + 1*2¹ + 1*2⁰.

Откуда находим, что m = 3, и:

x₀ = 2²⁹.

Тогда:

x₁ = 2³⁰ - 11*2²⁶ = 335544320;

x₂ = 2*335544320 - (11*335544320*335544320)>>32 = 382730240;

x₃ = 2*382730240 - (11*382730240*382730240)>>32 = 390298880;

x₄ = 2*390298880 - (11*390298880*390298880)>>32 = 390451513;

x₅ = 2*390451513 - (11*390451513*390451513)>>32 = 390451572;

x₆ = 2*390451572 - (11*390451572*390451572)>>32 = 390451572;

После чего находим:

A/B = [(A*x)>>32] = [(555555*390451572)>>32] = 50505 = D.

Но в целом, все эти "быстро сходящиеся итерации" не слишком то оптимальны для FPGA, КМК..

Как-то так..