Вычисление логарифмического выражения без плавающей точки арифметики или войти

Question

мне нужно вычислить математическое выражение floor(ln(u)/ln(1-p)) для 0 < u < 1 и 0 < p < 1 в C на встроенный процессор с не с плавающей точкой арифметики и не ln функцию. Результат - положительное целое число. Я знаю о предельных случаях (р = 0), я буду иметь дело с ними позже ...

Я полагаю, что решение подразумевает наличие u и p диапазона над 0..UINT16_MAX и обратиться к справочной таблице для логарифмов, но я не могу понять, как именно: на что выглядит таблица поиска?

Результат не должен быть на 100% точным, приближения в порядке.

Спасибо!

Answer 1

Поскольку логарифм используется как в дивиденде, так и в делителе, нет необходимости использовать log(); мы можем использовать log2(). Из-за ограничений на входы u и p логарифмы, как известно, оба отрицательные, поэтому мы можем ограничить себя вычислением положительной величины -log2().

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

Этот процесс известен в литературе как мультипликативная нормализация или псевдоделение, а некоторые ранние публикации, описывающие его, - это работы Де Лугиша и Меггита. Последнее указывает, что происхождение в основном является методом Генри Бриггса для вычисления общих логарифмов.

B. de Lugish. «Класс алгоритмов автоматической оценки функций и вычислений на цифровом компьютере». Кандидатская диссертация, кафедра компьютерных наук, Университет Иллинойса, Урбана, 1970.

J. E. Meggitt. «Процессы псевдоделения и псевдопроцесса». IBM Journal of Research and Development, Vol. 6, № 2, апрель 1962, стр. 210-226

В качестве выбранного набора факторов включает в себя 2 ^я и необходимые умножения могут быть выполнены без необходимости (1 +-^-i) инструкция умножения: продукты могут быть вычислены сдвигом или сдвигом плюс добавление.

Поскольку входы u и p являются чисто дробными числами с 16 битами, мы можем выбрать результат с фиксированной точкой 5.16 для логарифма. Просто разделив два значения логарифма, мы удалим масштабный коэффициент с фиксированной точкой и применим операцию floor() в то же время, потому что для положительных чисел floor(x) идентичен trunc(x)10, а целочисленное деление усекается.

Обратите внимание, что вычисление логарифма с фиксированной запятой приводит к большой относительной ошибке для входов около 1.Это, в свою очередь, означает, что вся функция, рассчитанная с использованием арифметики с фиксированной точкой, может выдавать результаты, значительно отличающиеся от ссылки, если p невелик. Примером этого является следующий тестовый пример: u=55af p=0052 res=848 ref=874.

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

/* input x is a 0.16 fixed-point number in [0,1) 
    function returns -log2(x) as a 5.16 fixed-point number in (0, 16] 
*/ 
uint32_t nlog2_16 (uint16_t x) 
{ 
    uint32_t r = 0; 
    uint32_t t, a = x; 

    /* try factors 2**i with i = 8, 4, 2, 1 */ 
    if ((t = a << 8  ) < 0x10000) { a = t; r += 0x80000; } 
    if ((t = a << 4  ) < 0x10000) { a = t; r += 0x40000; } 
    if ((t = a << 2  ) < 0x10000) { a = t; r += 0x20000; } 
    if ((t = a << 1  ) < 0x10000) { a = t; r += 0x10000; } 
    /* try factors (1+2**(-i)) with i = 1, .., 16 */ 
    if ((t = a + (a >> 1)) < 0x10000) { a = t; r += 0x095c0; } 
    if ((t = a + (a >> 2)) < 0x10000) { a = t; r += 0x0526a; } 
    if ((t = a + (a >> 3)) < 0x10000) { a = t; r += 0x02b80; } 
    if ((t = a + (a >> 4)) < 0x10000) { a = t; r += 0x01664; } 
    if ((t = a + (a >> 5)) < 0x10000) { a = t; r += 0x00b5d; } 
    if ((t = a + (a >> 6)) < 0x10000) { a = t; r += 0x005ba; } 
    if ((t = a + (a >> 7)) < 0x10000) { a = t; r += 0x002e0; } 
    if ((t = a + (a >> 8)) < 0x10000) { a = t; r += 0x00171; } 
    if ((t = a + (a >> 9)) < 0x10000) { a = t; r += 0x000b8; } 
    if ((t = a + (a >> 10)) < 0x10000) { a = t; r += 0x0005c; } 
    if ((t = a + (a >> 11)) < 0x10000) { a = t; r += 0x0002e; } 
    if ((t = a + (a >> 12)) < 0x10000) { a = t; r += 0x00017; } 
    if ((t = a + (a >> 13)) < 0x10000) { a = t; r += 0x0000c; } 
    if ((t = a + (a >> 14)) < 0x10000) { a = t; r += 0x00006; } 
    if ((t = a + (a >> 15)) < 0x10000) { a = t; r += 0x00003; } 
    if ((t = a + (a >> 16)) < 0x10000) { a = t; r += 0x00001; } 
    return r; 
} 

/* Compute floor(log(u)/log(1-p)) for 0 < u < 1 and 0 < p < 1, 
    where 'u' and 'p' are represented as 0.16 fixed-point numbers 
    Result is an integer in range [0, 1048676] 
*/ 
uint32_t func (uint16_t u, uint16_t p) 
{ 
    uint16_t one_minus_p = 0x10000 - p; // 1.0 - p 
    uint32_t log_u = nlog2_16 (u); 
    uint32_t log_p = nlog2_16 (one_minus_p); 
    uint32_t res = log_u/log_p; // divide and floor in one go 
    return res; 
} 

Answer 2

Максимальное значение этой функции в основном зависит от предел точности; то есть как сколь угодно близко к пределам (u -> 0) или (1 - p -> 1) значения фиксированной точки могут быть.

Если мы предположим (k) дробные биты, например, с пределами: u = (2^-k) и 1 - p = 1 - (2^-k),
, то максимальное значение: k/(k - log2(2^k - 1))

(как отношение натуральных логарифмов, мы свободны использовать любые основания, например , lb(x) или log2)

---

в отличие от njuffa's ответа, я пошел с табличного подхода, оседая на k = 10 дробных битов для представления 0 < frac(u) < 1024 и 0 < frac(p) < 1024. Для этого требуется таблица журналов с записями 2^k. Используя 32-битные значения таблиц, мы смотрим только на таблицу 4KiB.

Не более того, и вы используете достаточно памяти, чтобы серьезно рассмотреть возможность использования соответствующих частей библиотеки «soft-float». например, k = 16 даст 256KiB LUT.

Мы вычисляем значения - log2(i/1024.0) для 0 < i < 1024. Поскольку эти значения находятся в открытом интервале (0, k), нам нужно только 4 двоичных разряда для хранения интегральной части. Таким образом, мы сохраняем LUT предварительно вычисленное в 32-битном формате [4.28] с фиксированной точкой:

uint32_t lut[1024]; /* never use lut[0] */ 

for (uint32_t i = 1; i < 1024; i++) 
    lut[i] = (uint32_t) (- (log2(i/1024.0) * (268435456.0)); 

Дано: u, p представлены [0.10] значений с фиксированной точкой в ​​[1, 1023]:

uint32_t func (uint16_t u, uint16_t p) 
{ 
    /* assert: 0 < u, p < 1024 */ 

    return lut[u]/lut[1024 - p]; 
} 

Мы можем легко проверить все действительные (u, p) пары против оценки «с наименьшими числами»:

floor(log(u/1024.0)/log(1.0 - p/1024.0))

и получить только несоответствие (+1 слишком высока) на следующих случаях:

u = 193, p = 1 : 1708 vs 1707 (1.7079978488147417e+03) 
u = 250, p = 384 :  3 vs  2 (2.9999999999999996e+00) 
u = 413, p = 4 : 232 vs 231 (2.3199989016957960e+02) 
u = 603, p = 1 : 542 vs 541 (5.4199909906444600e+02) 
u = 680, p = 1 : 419 vs 418 (4.1899938077226307e+02) 

---

Наконец, оказывается, что используя натуральный логарифм в формате [3.29] с фиксированной точкой дает нам даже более высокая точность, где:

lut[i] = (uint32_t) (- (log(i/1024.0) * (536870912.0)); 

только дает один «несоответствие», хотя 'bignum' точности предполагает, что это правильно:

u = 250, p = 384 :  3 vs  2 (2.9999999999999996e+00)