Как максимально точно вычислить log2 целого числа в C с помощью побитовых операций - PullRequest
Новая
       36

Как максимально точно вычислить log2 целого числа в C с помощью побитовых операций

0 голосов
/ 16 декабря 2018

Мне нужно рассчитать энтропию, и из-за ограничений моей системы мне нужно использовать ограниченные функции C (без циклов, без поддержки с плавающей запятой) и мне нужна как можно большая точность.Начиная с здесь я выясняю, как оценить пол log2 целого числа, используя побитовые операции.Тем не менее, мне нужно повысить точность результатов.Поскольку никакие операции с плавающей запятой не допускаются, есть ли способ вычислить log2(x/y) с помощью x < y, чтобы результат был примерно таким же, как log2(x/y)*10000, с целью получения необходимой мне точности с помощью арифметического целого числа?

1 Ответ

0 голосов
/ 16 декабря 2018

Вы будете основывать алгоритм на формуле 

log2(x/y) = K*(-log(x/y));

, где 

 K        = -1.0/log(2.0); // you can precompute this constant before run-time
 a        = (y-x)/y;
-log(x/y) = a + a^2/2 + a^3/3 + a^4/4 + a^5/5 + ...

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

(y^N*(1*2*3*4*5*...*N)) * (-log(x/y))
  = y^(N-1)*(2*3*4*5*...*N)*(y-x) + y^(N-2)*(1*3*4*5*...*N)*(y-x)^2 + ...

Конечно, ^, оператор power, связывающий более сильно, чем *, не является оператором C, но вы можетеэффективно реализовать это в контексте вашего (возможно, развернутого) цикла как работающего продукта.

* N - это целое число, достаточно большое, чтобы обеспечить желаемую точность, но не настолько большое, чтобы оно превышало количество битов, которые у вас естьимеется в наличии.Если не уверены, попробуйте, например, N = 6.Что касается K, вы можете возразить, что это число с плавающей запятой, но это не проблема для вас, потому что вы собираетесь предварительно вычислить K, сохранив его как отношение целых чисел.

SAMPLE CODE

Это игрушечный код, но он работает для небольших значений x и y, таких как 5 и 7, таким образом, достаточно для подтверждения концепции.В игрушечном коде большие значения могут молча переполнять стандартные 64-битные регистры.Чтобы сделать код более надежным, потребовалось бы больше работы. 

#include <stddef.h>
#include <stdlib.h>
// Your program will not need the below headers, which are here
// included only for comparison and demonstration.
#include <math.h>
#include <stdio.h>

const size_t     N = 6;
const long long Ky = 1 << 10; // denominator of K
// Your code should define a precomputed value for Kx here.

int main(const int argc, const char *const *const argv)
{
    // Your program won't include the following library calls but this
    // does not matter.  You can instead precompute the value of Kx and
    // hard-code its value above with Ky.
    const long long Kx = lrintl((-1.0/log(2.0))*Ky); // numerator of K
    printf("K == %lld/%lld\n", Kx, Ky);

    if (argc != 3) exit(1);

    // Read x and y from the command line.
    const long long x0 = atoll(argv[1]);
    const long long y  = atoll(argv[2]);
    printf("x/y == %lld/%lld\n", x0, y);
    if (x0 <= 0 || y <= 0 || x0 > y) exit(1);

    // If 2*x <= y, then, to improve accuracy, double x repeatedly
    // until 2*x > y. Each doubling offsets the log2 by 1. The offset
    // is to be recovered later.
    long long               x = x0;
    int integral_part_of_log2 = 0;
    while (1) {
        const long long trial_x = x << 1;
        if (trial_x > y) break;
        x = trial_x;
        --integral_part_of_log2;
    }
    printf("integral_part_of_log2 == %d\n", integral_part_of_log2);

    // Calculate the denominator of -log(x/y).
    long long yy = 1;
    for (size_t j = N; j; --j) yy *= j*y;

    // Calculate the numerator of -log(x/y).
    long long xx = 0;
    {
        const long long y_minus_x = y - x;
        for (size_t i = N; i; --i) {
            long long term = 1;
            size_t j       = N;
            for (; j > i; --j) {
                term *= j*y;
            }
            term *= y_minus_x;
            --j;
            for (; j; --j) {
                term *= j*y_minus_x;
            }
            xx += term;
        }
    }

    // Convert log to log2.
    xx *= Kx;
    yy *= Ky;

    // Restore the aforementioned offset.
    for (; integral_part_of_log2; ++integral_part_of_log2) xx -= yy;

    printf("log2(%lld/%lld) == %lld/%lld\n", x0, y, xx, yy);
    printf("in floating point, this ratio of integers works out to %g\n",
      (1.0*xx)/(1.0*yy));
    printf("the CPU's floating-point unit computes the log2 to be  %g\n",
      log2((1.0*x0)/(1.0*y)));

    return 0;
}

Запустив это на моем компьютере с аргументами командной строки 5 7, вы получите: 

K == -1477/1024
x/y == 5/7
integral_part_of_log2 == 0
log2(5/7) == -42093223872/86740254720
in floating point, this ratio of integers works out to -0.485279
the CPU's floating-point unit computes the log2 to be  -0.485427

Точность будетЗначительно улучшены на N = 12 и Ky = 1 << 20, но для этого вам нужен либо код Thriftier, либо более 64 бит.

КОД THRIFTIER

Код Thriftier, требующий большегоусилия, чтобы написать, может представлять числитель и знаменатель в простых факторах.Например, он может представлять 500 как [2 0 3], что означает (2 2 ) (3 0 ) (5 3 ).

Еще больше улучшений может произойти в вашем воображении.

АЛЬТЕРНАТИВНЫЙ ПОДХОД

Для альтернативного подхода, хотя он может не соответствовать вашим требованиям точно так, как вы их сформулировали@phuclv дал предложение, которому я был бы склонен следовать, если бы ваша программа была моей: обработайте задачу в обратном порядке, угадывая значение c/d для логарифма и затем вычисляя 2^(c/d), предположительно черезИтерация Ньютона-Рафсона.Лично мне больше нравится подход Ньютона-Рафсона.См. Раздел4.8 здесь (мой оригинал).

МАТЕМАТИЧЕСКИЙ ФОН

Несколько источников, включая мой, уже связаны, объясняют ряд Тейлора, лежащий в основе первого подхода, и Ньютон-Рапсоновская итерация второго подхода.Математика, к сожалению, нетривиальна, но она у вас есть.Удачи.

...