通常nextafter通过以下方式实作：
double nextafter(double x, double y)
{
// handle corner cases
int delta = ((x > 0) == (x < y)) ? 1 : -1;
unsigned long long mant = __mant(x); // get mantissa as int
mant = delta;
...
}
这里，使用 获得二进制表示__mant(x)。
出于好奇：是否可以在nextafter不获得二进制表示的情况下实作？例如，使用一系列算术浮点运算。
uj5u.com热心网友回复：
下面的代码nextafter在上升方向上实作了 IEEE-754 算法的有限值，并使用舍入到最近关系到偶数。NaN、无穷大和下降方向的处理是显而易见的。
在不假设 IEEE-754 或舍入到最近的情况下，浮点属性在 C 2018 5.2.4.2.2 中有足够好的特征，我们可以通过这种方式实作 nextafter（再次在上升方向）：
- 如果输入是NaN，则回传，如果是信令NaN则报错。
- 如果输入是?∞，则回传
-DBL_MAX。 - 如果输入为
-DBL_TRUE_MIN，则回传零。 - 如果输入为零，则回传
DBL_TRUE_MIN。 - 如果输入是
DBL_MAX，则回传 ∞。 - 如果输入是 ∞，则回传 ∞。（请注意，这在完整
nextafter(x, y)实作中永远不会发生，因为它会将第一个自变量移动到第二个自变量的方向，所以我们永远不会从 ∞ 上升，因为我们永远不会收到大于 ∞ 的第二个自变量。） - 否则，如果它是正数，则使用
logb来查找指数e。如果e小于DBL_MIN，则回传输入加号DBL_TRUE_MIN（次法线和最低法线的 ULP）。如果e不小于DBL_MIN，则回传输入加号scalb(1, e 1 - DBL_MANT_DIG)（输入的特定 ULP）。舍入方法无关紧要，因为这些添加是准确的。 - 否则，输入为负。除非输入恰好是
FLT_RADIX（输入等于scalb(1, e)）的幂，否则使用上面的方法，将 的第二个自变量减scalb一（因为此nextafter步骤从较大的指数过渡到较低的指数）。
请注意，FLT_RADIX上面是正确的；没有DBL_RADIX; 所有浮点格式都使用相同的基数。
如果您想将logb和scalb视为操作浮点表示的函式，则可以将它们替换为普通算术。log可以找到可以快速改进为真实指数的快速近似值，并且scalb可以通过多种方式实作，可能只是查表。如果log仍然令人反感，那么试验比较就足够了。
上面处理有或没有次法线的格式，因为如果支持次法线，它会随着递减而进入它们，如果不支持次法线，最小法线幅度是DBL_TRUE_MIN，所以它在上面被认为是我们步进的点接下来归零。
有一个警告；C 标准允许“不确定”一个实作是否支持次正规，“如果浮点运算不能一致地将次正规表示解释为零，也不能解释为非零。” 在那种情况下，我没有看到标准指定了标准的nextafter作用，所以我们没有什么可以在我们的实作中匹配它。假设有时支持次正规，DBL_TRUE_MIN必须是次正规值，并且以上将尝试作业，好像次正规支持当前处于打开状态（例如，重绘 为零已关闭），如果关闭，您将获得任何您想要的得到。
#include <float.h>
#include <math.h>
/* Return the next floating-point value after the finite value q.
This was inspired by Algorithm 3.5 in Siegfried M. Rump, Takeshi Ogita, and
Shin'ichi Oishi, "Accurate Floating-Point Summation", _Technical Report
05.12_, Faculty for Information and Communication Sciences, Hamburg
University of Technology, November 13, 2005.
IEEE-754 and the default rounding mode,
round-to-nearest-ties-to-even, may be required.
*/
double NextAfter(double q)
{
/* Scale is .625 ULP, so multiplying it by any significand in [1, 2)
yields something in [.625 ULP, 1.25 ULP].
*/
static const double Scale = 0.625 * DBL_EPSILON;
/* Either of the following may be used, according to preference and
performance characteristics. In either case, use a fused multiply-add
(fma) to add to q a number that is in [.625 ULP, 1.25 ULP]. When this
is rounded to the floating-point format, it must produce the next
number after q.
*/
#if 0
// SmallestPositive is the smallest positive floating-point number.
static const double SmallestPositive = DBL_EPSILON * DBL_MIN;
if (fabs(q) < 2*DBL_MIN)
return q SmallestPositive;
return fma(fabs(q), Scale, q);
#else
return fma(fmax(fabs(q), DBL_MIN), Scale, q);
#endif
}
#if defined CompileMain
#include <stdio.h>
#include <stdlib.h>
#define NumberOf(a) (sizeof (a) / sizeof *(a))
int main(void)
{
int status = EXIT_SUCCESS;
static const struct { double in, out; } cases[] =
{
{ INFINITY, INFINITY },
{ 0x1.fffffffffffffp1023, INFINITY },
{ 0x1.ffffffffffffep1023, 0x1.fffffffffffffp1023 },
{ 0x1.ffffffffffffdp1023, 0x1.ffffffffffffep1023 },
{ 0x1.ffffffffffffcp1023, 0x1.ffffffffffffdp1023 },
{ 0x1.0000000000003p1023, 0x1.0000000000004p1023 },
{ 0x1.0000000000002p1023, 0x1.0000000000003p1023 },
{ 0x1.0000000000001p1023, 0x1.0000000000002p1023 },
{ 0x1.0000000000000p1023, 0x1.0000000000001p1023 },
{ 0x1.fffffffffffffp1022, 0x1.0000000000000p1023 },
{ 0x1.fffffffffffffp1, 0x1.0000000000000p2 },
{ 0x1.ffffffffffffep1, 0x1.fffffffffffffp1 },
{ 0x1.ffffffffffffdp1, 0x1.ffffffffffffep1 },
{ 0x1.ffffffffffffcp1, 0x1.ffffffffffffdp1 },
{ 0x1.0000000000003p1, 0x1.0000000000004p1 },
{ 0x1.0000000000002p1, 0x1.0000000000003p1 },
{ 0x1.0000000000001p1, 0x1.0000000000002p1 },
{ 0x1.0000000000000p1, 0x1.0000000000001p1 },
{ 0x1.fffffffffffffp-1022, 0x1.0000000000000p-1021 },
{ 0x1.ffffffffffffep-1022, 0x1.fffffffffffffp-1022 },
{ 0x1.ffffffffffffdp-1022, 0x1.ffffffffffffep-1022 },
{ 0x1.ffffffffffffcp-1022, 0x1.ffffffffffffdp-1022 },
{ 0x1.0000000000003p-1022, 0x1.0000000000004p-1022 },
{ 0x1.0000000000002p-1022, 0x1.0000000000003p-1022 },
{ 0x1.0000000000001p-1022, 0x1.0000000000002p-1022 },
{ 0x1.0000000000000p-1022, 0x1.0000000000001p-1022 },
{ 0x0.fffffffffffffp-1022, 0x1.0000000000000p-1022 },
{ 0x0.ffffffffffffep-1022, 0x0.fffffffffffffp-1022 },
{ 0x0.ffffffffffffdp-1022, 0x0.ffffffffffffep-1022 },
{ 0x0.ffffffffffffcp-1022, 0x0.ffffffffffffdp-1022 },
{ 0x0.0000000000003p-1022, 0x0.0000000000004p-1022 },
{ 0x0.0000000000002p-1022, 0x0.0000000000003p-1022 },
{ 0x0.0000000000001p-1022, 0x0.0000000000002p-1022 },
{ 0x0.0000000000000p-1022, 0x0.0000000000001p-1022 },
{ -0x1.fffffffffffffp1023, -0x1.ffffffffffffep1023 },
{ -0x1.ffffffffffffep1023, -0x1.ffffffffffffdp1023 },
{ -0x1.ffffffffffffdp1023, -0x1.ffffffffffffcp1023 },
{ -0x1.0000000000004p1023, -0x1.0000000000003p1023 },
{ -0x1.0000000000003p1023, -0x1.0000000000002p1023 },
{ -0x1.0000000000002p1023, -0x1.0000000000001p1023 },
{ -0x1.0000000000001p1023, -0x1.0000000000000p1023 },
{ -0x1.0000000000000p1023, -0x1.fffffffffffffp1022 },
{ -0x1.0000000000000p2, -0x1.fffffffffffffp1 },
{ -0x1.fffffffffffffp1, -0x1.ffffffffffffep1 },
{ -0x1.ffffffffffffep1, -0x1.ffffffffffffdp1 },
{ -0x1.ffffffffffffdp1, -0x1.ffffffffffffcp1 },
{ -0x1.0000000000004p1, -0x1.0000000000003p1 },
{ -0x1.0000000000003p1, -0x1.0000000000002p1 },
{ -0x1.0000000000002p1, -0x1.0000000000001p1 },
{ -0x1.0000000000001p1, -0x1.0000000000000p1 },
{ -0x1.0000000000000p-1021, -0x1.fffffffffffffp-1022 },
{ -0x1.fffffffffffffp-1022, -0x1.ffffffffffffep-1022 },
{ -0x1.ffffffffffffep-1022, -0x1.ffffffffffffdp-1022 },
{ -0x1.ffffffffffffdp-1022, -0x1.ffffffffffffcp-1022 },
{ -0x1.0000000000004p-1022, -0x1.0000000000003p-1022 },
{ -0x1.0000000000003p-1022, -0x1.0000000000002p-1022 },
{ -0x1.0000000000002p-1022, -0x1.0000000000001p-1022 },
{ -0x1.0000000000001p-1022, -0x1.0000000000000p-1022 },
{ -0x1.0000000000000p-1022, -0x0.fffffffffffffp-1022 },
{ -0x0.fffffffffffffp-1022, -0x0.ffffffffffffep-1022 },
{ -0x0.ffffffffffffep-1022, -0x0.ffffffffffffdp-1022 },
{ -0x0.ffffffffffffdp-1022, -0x0.ffffffffffffcp-1022 },
{ -0x0.0000000000004p-1022, -0x0.0000000000003p-1022 },
{ -0x0.0000000000003p-1022, -0x0.0000000000002p-1022 },
{ -0x0.0000000000002p-1022, -0x0.0000000000001p-1022 },
{ -0x0.0000000000001p-1022, -0x0.0000000000000p-1022 },
};
for (int i = 0; i < NumberOf(cases); i)
{
double in = cases[i].in, expected = cases[i].out;
double observed = NextAfter(in);
printf("NextAfter(%a) = %a.\n", in, observed);
if (! (observed == expected))
{
printf("\tError, expected %a.\n", expected);
status = EXIT_FAILURE;
}
}
return status;
}
#endif // defined CompileMain
uj5u.com热心网友回复：
考虑到 FP 可能/可能不支持次法线、无穷大、大多数值具有唯一值（例如使用 2float表示double）、支持 /= 0、没有 FP 编码的内在知识，使用假设mant = delta;是下一个值导致可移植性失败 - 即使使用二进制操作。仅使用 FP 操作需要对 FP 编码进行许多假设。
我认为更有用的方法是发布使用“一系列算术浮点运算”的候选代码，然后询问 1）它失败的条件 2）如何改进？
uj5u.com热心网友回复：
我只能说，这只有在融合乘加 (FMA) 可用时才有可能。对于我下面的示例 ISO-C99 实作，我使用float映射到 IEEE-754，binary32因为我可以通过这种方式获得更好的测验覆写率。假设 C 的 IEEE-754 系结有效（因此 C 浮点型别系结到 IEEE-754 二进制型别，支持次正规等），有效的舍入模式是舍入到最近或偶数, 和 ISO-C 标准 ( FE_INEXACT, FE_OVERFLOW, FE_UNDERFLOW)指定的例外信号要求被放弃（传递掩码回应）。
代码可能比它需要的更复杂；我只是简单地将各种操作数类分开并一一处理。我使用具有最严格浮点设定的 Intel 编译器进行编译。实施nextafterf()从英特尔数学库作为黄金参考。我认为我的实作存在与英特尔库实作中的错误相匹配的错误极不可能，但显然并非不可能。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <float.h>
#include <string.h>
#include <fenv.h>
#include <math.h>
float my_nextafterf (float a, float b)
{
const float FP32_MIN_NORMAL = 0x1.000000p-126f;
const float FP32_MAX_NORMAL = 0x1.fffffep 127f;
const float FP32_EPSILON = 0x1.0p-23f;
const float FP32_ONE = 1.0f;
const float FP32_HALF = 0.5f;
const float FP32_SUBNORM_SCALE = FP32_ONE / FP32_MIN_NORMAL;
const float FP32_INC = (FP32_ONE FP32_EPSILON) * FP32_EPSILON * FP32_HALF;
const float FP32_INT_SCALE = FP32_ONE / FP32_EPSILON;
float r;
if ((!(fabsf(a) <= INFINITY)) || (!(fabsf(b) <= INFINITY))) { // unordered
r = a b;
}
else if (a == b) { // equal
r = b;
}
else if (fabsf (a) == INFINITY) { // infinity
r = (a >= 0) ? FP32_MAX_NORMAL : (-FP32_MAX_NORMAL);
}
else if (fabsf (a) >= FP32_MIN_NORMAL) { // normal
float factor = ((a >= 0) != (a > b)) ? FP32_INC : (-FP32_INC);
r = fmaf (factor, a, a);
} else { // subnormal or zero
float scal = (a >= 0) ? FP32_INT_SCALE : (-FP32_INT_SCALE);
float incr = ((a >= 0) != (a > b)) ? FP32_ONE : (-FP32_ONE);
r = (a * scal * FP32_SUBNORM_SCALE incr) / scal / FP32_SUBNORM_SCALE;
}
return r;
}
float uint32_as_float (uint32_t a)
{
float r;
memcpy (&r, &a, sizeof r);
return r;
}
uint32_t float_as_uint32 (float a)
{
uint32_t r;
memcpy (&r, &a, sizeof r);
return r;
}
// Fixes via: Greg Rose, KISS: A Bit Too Simple. http://eprint.iacr.org/2011/007
static uint32_t kiss_z = 362436069;
static uint32_t kiss_w = 521288629;
static uint32_t kiss_jsr = 123456789;
static uint32_t kiss_jcong = 380116160;
#define znew (kiss_z=36969*(kiss_z&0xffff) (kiss_z>>16))
#define wnew (kiss_w=18000*(kiss_w&0xffff) (kiss_w>>16))
#define MWC ((znew<<16) wnew)
#define SHR3 (kiss_jsr^=(kiss_jsr<<13),kiss_jsr^=(kiss_jsr>>17),kiss_jsr^=(kiss_jsr<<5))
#define CONG (kiss_jcong=69069*kiss_jcong 13579)
#define KISS ((MWC^CONG) SHR3)
int main (void)
{
float a, b, res, ref;
uint32_t ia, ib, ires, iref;
unsigned long long int count = 0;
const uint32_t FP32_QNAN_BIT = 0x00400000;
const uint32_t FP32_SIGN_BIT = 0x80000000;
const uint32_t FP32_INFINITY = 0x7f800000;
printf ("Testing nextafterf()\n");
while (1) {
ia = KISS;
ib = KISS;
/* increase likelihood of zeros, infinities, equality */
if (!(KISS & 0xfff)) ia = ia & FP32_SIGN_BIT;
if (!(KISS & 0xfff)) ib = ib & FP32_SIGN_BIT;
if (!(KISS & 0xfff)) ia = ia | FP32_INFINITY;
if (!(KISS & 0xfff)) ib = ib | FP32_INFINITY;
if (!(KISS & 0xffffff)) ib = ia;
a = uint32_as_float (ia);
b = uint32_as_float (ib);
res = my_nextafterf (a, b);
ref = nextafterf (a, b);
ires = float_as_uint32 (res);
iref = float_as_uint32 (ref);
if (ires != iref) {
/* if both 'from' and 'to' are NaN, result may be either NaN, quietened */
if (!(isnan (a) && isnan (b) &&
((ires == (ia | FP32_QNAN_BIT)) || (ires == (ib | FP32_QNAN_BIT))))) {
printf ("error: a=x b=x res=x ref=x\n", ia, ib, ires, iref);
}
}
count ;
if (!(count & 0xffffff)) printf ("\rcount = 0x%llx", count);
}
return EXIT_SUCCESS;
}
0 评论