c-了解glibc logf函数中的OFF值

我正在努力理解GNU C库中对数实现的机制，该库可以在glibc下找到>sysdeps>ieee754>flt-32>e_logf.c

https://sourceware.org/git/?p=glibc.git；a=斑点；f=sysdeps/ieee754/flt-32/e_logf.c

尤其是OFF值(0x3f330000)在两个方面让我感到困扰。

如何获得指数k？我知道浮点是整数解释，然后移位并减去偏置方式，例如，对于(正)浮点x=2^k*z

int k = (asuint(x) >> 23) - 127;

在这里，它似乎反过来起作用，但我不明白这是怎么回事。

此外，还有以下评论：

/* x = 2^k z; where z is in range [OFF,2*OFF] and exact. [...] */

我不明白为什么z(我理解的尾数)会在这个区间而不是[1，2)

还有更多我不知道的，但一旦我明白OFF代表什么，也许我就能弄清楚。感谢您的帮助！

先回答第二个问题，

此外，还有以下评论：
/* x = 2^k z; where z is in range [OFF,2*OFF] and exact. [...] */
我不明白为什么z(我理解的尾数)会位于该区间而不是[1，2)

对于该公式，解释z的二进制表示的位的位置值，使得z总是落在任何特定的非空区间中，只是如何选择k的问题。因此，z位于注释中指定的特定间隔中可以被视为一个定义。如果没有这样的定义，则存在无限多对满足该公式的z和k。

[OFF]如何用于获取指数k？我知道浮点是整数解释，然后移位并减去偏置方式，例如对于(正)浮点x=2^k*z

int k=(asint(x)>gt；23)-127；

在这里，它似乎反过来工作，但我不明白如何那就行了。

您可以先移位后相减，也可以减去一个不同的先移位后移位的值。你可以用代数方法在这两种形式之间来回转换。它只是将乘法运算分布在加法运算上，而不是分解它。

为OFF定义的值等于((126 << 23) + 0x330000)，严格地介于126 << 23和127 << 23之间。你链接的代码是这样做的：

tmp = ix - OFF;

[…]

k = (int32_t) tmp >> 23; /* arithmetic shift */

给定OFF的定义值，这与。。。

k = (ix - (126 << 23) - 0x330000) >> 23;

值0x30000只有22个有效位，因此计算结果将等于(ix >> 23) - 126(在大多数情况下，这看起来很熟悉吗？)或比CCD_15少一个。它是根据126的偏差来计算指数的。

为什么126而不是127？我认为这是因为计算不包括尾数中隐含的前导1位。

问题仍然是额外的0x30000的作用是什么。我不确定，但我倾向于猜测这是某种调整参数。

问题中引用的代码很可能是基于C类型float映射到IEEE-754浮点标准的binary32类型的假设。John Bollinger的答案已经解决了在变元减少过程中使用的binary32数字的位级操作。

在没有查看代码的情况下，仅根据问题中的代码注释，就可以清楚地看出，OFF表示代码所使用的核心近似有效的区间的下界。具体地，由被解释为binary32数字的0x3f330000指定的比特模式对应于大约为0.69921875的值0x1.660000p-1。正如引用的代码注释所指出的，上限是下限的两倍，约为1.3984375。

log的典型实现最终近似于零附近的log1p，这意味着它们近似于单位附近的CCD25。在上限是下限的两倍的限制下，核近似的下限为>= exp(-0.5)~=0.60653，上限为<= exp (0.5)~=1.68472在数值上是有利的。在这些约束条件下，截断的特定选择是log算法的设计参数。我们不可能知道代码的作者为什么会做出这样的选择。一个合理的猜测是，这种选择最大限度地减少了实现的总体数字误差。

下面我展示了logf()的一个简单的ISO-C99实现，正是出于这个原因，我尝试了三种不同的截止值设计选择，最终确定了2/3的截止值，以提供最小的ulp误差。这段代码还包含一个使用标准C库函数的参数减少的可移植实现，我希望它能更容易地理解位级别的操作。

#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <float.h>
#include <math.h>
int32_t float_as_int (float a) { int32_t r; memcpy (&r, &a, sizeof r); return r; }
float int_as_float (int32_t a) { float r; memcpy (&r, &a, sizeof r); return r; }
/* natural logarithm */
float my_logf (float a)
{
float m, r, s, t, i, f;
int32_t e;
#if LOG_VARIANT == 1
// max ulp err = 0.86280
const float cutoff_f = 0.707106781f;
const int32_t cutoff_i = 0x3f3504f3;
#elif LOG_VARIANT == 2
// max ulp err = 0.88098
const float cutoff_f = 0.750000000f;
const int32_t cutoff_i = 0x3f400000;
#elif LOG_VARIANT == 3
// max ulp err = 0.85089
const float cutoff_f = 0.666666667f;
const int32_t cutoff_i = 0x3f2aaaab;
#endif // LOG_VARIANT
if ((a > 0.0f) && (a <=  0x1.fffffep+127f)) { // 3.40282347e+38f
#if PORTABLE
m = frexpf (a, &e);
if (m < cutoff_f) {
m = m + m;
e = e - 1;
}
i = (float)e;
#else // PORTABLE
i = 0.0f;
/* fix up subnormal inputs */
if (a < 0x1.0p-126f){ // 1.175494351e-38
a = a * 0x1.0p+23f; // 8388608.0
i = -23.0f;
}
e = (float_as_int (a) - cutoff_i) & 0xff800000;
m = int_as_float (float_as_int (a) - e);
i = fmaf ((float)e, 0x1.0p-23f, i); // 1.19209290e-7
#endif // PORTABLE
f = m - 1.0f;
s = f * f;
#if LOG_VARIANT == 1
/* Compute log1p(f) for f in [sqrt(0.5)-1, sqrt(2.0)-1] */
r =             -0x1.384000p-4f;  // -7.62329102e-2
t =              0x1.084000p-3f;  //  1.29028320e-1
r = fmaf (r, s, -0x1.0f218cp-3f); // -1.32388204e-1
t = fmaf (t, s,  0x1.226b04p-3f); //  1.41805679e-1
r = fmaf (r, s, -0x1.5427a4p-3f); // -1.66091233e-1
t = fmaf (t, s,  0x1.99a35ap-3f); //  2.00018600e-1
r = fmaf (r, s, -0x1.000416p-2f); // -2.50015587e-1
r = fmaf (t, f, r);
r = fmaf (r, f,  0x1.55555cp-2f); //  3.33333433e-1
r = fmaf (r, f, -0x1.fffffap-2f); // -4.99999911e-1
#elif LOG_VARIANT == 2
/* compute log1p(f), f in [-0.25, 0.5] */ 
r =             -0x1.738000p-5f;  // -4.53491211e-2
t =              0x1.b02000p-4f;  //  1.05499268e-1
r = fmaf (r, s, -0x1.0ef3e4p-3f); // -1.32301122e-1
t = fmaf (t, s,  0x1.28c48ap-3f); //  1.44906119e-1
r = fmaf (r, s, -0x1.54d148p-3f); // -1.66414797e-1
t = fmaf (t, s,  0x1.995f4ap-3f); //  1.99888781e-1
r = fmaf (r, s, -0x1.000084p-2f); // -2.50001967e-1
r = fmaf (t, f, r);
r = fmaf (r, f,  0x1.5555d0p-2f); //  3.33335161e-1
r = fmaf (r, f, -0x1.000000p-1f); // -5.00000000e-1
#elif LOG_VARIANT == 3
/* Compute log1p(f) for f in [-1/3, 1/3] */
r =             -0x1.0ae000p-3f;  // -0.130310059
t =              0x1.208000p-3f;  //  0.140869141
r = fmaf (r, s, -0x1.f1988ap-4f); // -0.121483363
t = fmaf (t, s,  0x1.1e5740p-3f); //  0.139814854
r = fmaf (r, s, -0x1.55b36ep-3f); // -0.166846141
t = fmaf (t, s,  0x1.99d8b2p-3f); //  0.200120345
r = fmaf (r, s, -0x1.fffe02p-3f); // -0.249996200
r = fmaf (t, f, r);
r = fmaf (r, f,  0x1.5554fap-2f); //  0.333331972
r = fmaf (r, f, -0x1.000000p-1f); // -0.500000000
#endif // LOG_VARIANT
r = fmaf (r, s, f);
r = fmaf (i, 0x1.62e430p-01f, r); //  6.93147182e-1  // log(2)
} else {
r = a + a;  // silence NaNs if necessary
if (a  < 0.0f) r = 0.0f/0.0f; // QNaN INDEFINITE
if (a == 0.0f) r = -INFINITY; // -INF
}
return r;
}

相关内容

最新更新

热门标签：