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

Question

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

Answer 1

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

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 + ...

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

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

ОБРАЗЕЦ КОДА

Это игрушечный код, но он работает для небольших значений 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, но для этого вам нужен либо более экономичный код, либо более 64 бит.

ТРИФТЕРНЫЙ КОД

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

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

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

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

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

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