【问题标题】:How to calculate the log2 of integer in C as precisely as possible with bitwise operations如何使用按位运算尽可能精确地计算 C 中整数的 log2
【发布时间】:2018-12-16 02:12:44
【问题描述】:

我需要计算熵,并且由于我的系统的限制,我需要使用受限的 C 功能(无循环,无浮点支持),并且我需要尽可能高的精度。从here 我知道如何使用按位运算来估计整数的下限 log2。不过,我需要提高结果的精度。由于不允许浮点运算,有没有什么办法可以用x < y计算log2(x/y),这样结果就会像log2(x/y)*10000一样,旨在通过算术整数获得我需要的精度?

【问题讨论】:

标签: c bit-manipulation bitwise-operators entropy


【解决方案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,将其存储为整数的比率。

示例代码

这是一个玩具代码,但它适用于 xy 的小值,例如 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 = 12Ky = 1 &lt;&lt; 20 将大大提高准确性,但为此您需要更简洁的代码或超过 64 位。

节俭代码

想要花费更多精力编写的节俭代码可能代表质因数中的分子和分母。例如,它可能将 500 表示为 [2 0 3],意思是 (22)(30)(53)。

您的想象力可能会进一步改进。

另一种方法

对于另一种方法,虽然它可能无法完全满足您的要求,但如果您的程序是我的,@phuclv 给出了我倾向于遵循的建议:反过来解决问题,猜测 一个值 c/d 用于对数,然后计算 2^(c/d),大概是通过 Newton-Raphson 迭代。就个人而言,我更喜欢 Newton-Raphson 方法。见节。 4.8 here(我的原创)。

数学背景

包括我已经链接的多个来源解释了第一种方法的泰勒级数和第二种方法的牛顿-拉夫森迭代。不幸的是,数学是不平凡的,但你有它。祝你好运。

【讨论】:

  • 我一直在尝试遵循您的建议,但到目前为止没有成功,我想我无法完全理解您的推理(例如找到Kalog(x/y) 系列) .即使不理解它,我也尝试对其进行编码,但结果是错误的。您能否进一步解释一下,如果可能的话,提供一个编码示例?
  • 当然,为什么不呢?你是在方便的时候问的,碰巧。我已经很多年没有写过这种代码了,所以这个练习有助于更新旧技能。代码发布。给你。
  • @phuclv :关于信息,我的回答采纳了你的评论。谢谢。
  • 这正是我所需要的。非常感谢!一旦我弄清楚如何修改代码以满足更大数量的估计,我将发布更新版本。再次感谢!
猜你喜欢
  • 2012-07-07
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2020-02-25
  • 2020-10-24
相关资源
最近更新 更多