C 乘法或加法浮点结果 NaN答案

【问题标题】：C multiplication or addition floats results NaNC 乘法或加法浮点结果 NaN
【发布时间】：2012-12-10 21:36:39
【问题描述】：

我在程序中使用 libresample。一段时间后（大约 50 分钟），它在一个工作站的 lib 函数 lrsFilterUD() 中崩溃。

float lrsFilterUD(float Imp[],  /* impulse response */
              float ImpD[], /* impulse response deltas */
              UWORD Nwing,  /* len of one wing of filter */
              BOOL Interp,  /* Interpolate coefs using deltas? */
              float *Xp,    /* Current sample */
              double Ph,    /* Phase */
              int Inc,    /* increment (1 for right wing or -1 for left) */
              double dhb)
{
   float a;
   float *Hp, *Hdp, *End;
   float v, t;
   double Ho;

   v = 0.0; /* The output value */
   Ho = Ph*dhb;
   End = &Imp[Nwing];
   if (Inc == 1)        /* If doing right wing...              */
   {                      /* ...drop extra coeff, so when Ph is  */
      End--;            /*    0.5, we don't do too many mult's */
      if (Ph == 0)      /* If the phase is zero...           */
         Ho += dhb;     /* ...then we've already skipped the */
   }                         /*    first sample, so we must also  */
                        /*    skip ahead in Imp[] and ImpD[] */

   if (Interp)
      while ((Hp = &Imp[(int)Ho]) < End) {
         t = *Hp;       /* Get IR sample */
         Hdp = &ImpD[(int)Ho];  /* get interp bits from diff table*/
         a = Ho - floor(Ho);      /* a is logically between 0 and 1 */
         t += (*Hdp)*a; /* t is now interp'd filter coeff */
         t *= *Xp;      /* Mult coeff by input sample */
         v += t;            /* The filter output */
         Ho += dhb;     /* IR step */
         Xp += Inc;     /* Input signal step. NO CHECK ON BOUNDS */
      }
   else 
      while ((Hp = &Imp[(int)Ho]) < End) {
         dprintf("while begin: Hp = %p, *Hp = %a, (int)Ho = %d, Imp[(int)Ho] = %a, &Imp[(int)Ho] = %p", Hp, *Hp, (int)Ho, Imp[(int)Ho], &Imp[(int)Ho]);
         t = *Hp;       /* Get IR sample */
         dprintf("before t = %a, *Xp = %a, Xp = %p", t, *Xp, Xp);
         t *= *Xp;      /* Mult coeff by input sample */
         dprintf("after2 t = %a, v = %a", t, v);
         v += t;            /* The filter output */
         dprintf("v = %a", v);
         Ho += dhb;     /* IR step */
         Xp += Inc;     /* Input signal step. NO CHECK ON BOUNDS */
      }

   return v;
}

我在乘法前后记录了 t、*Xp、Xp 的值：

while begin: Hp = 0xaf5daa8, *Hp = -0.009034, (int)Ho = 16384, Imp[(int)Ho] = -0.009034, &Imp[(int)Ho] = 0xaf5daa8
before multiplication t = -0.009034, *Xp = 0.000000, Xp = 0xaebe9b8
after multiplication t = nan

这段代码运行多次，崩溃前的t和Xp值相同：

before multiplication t = -0.009034, *Xp = 0.000000, Xp = 0xaebe9c8
after multiplication t = -0.000000, v = 282.423676

或另一种情况加法：

before addition t = -460.799988, v = 0.000000
after addition v = nan

什么可能导致nan？这是在 Linux 上使用 gcc 4.1.2 编译的。

更新：将变量打印为 %a。结果：

//t = 0x1.2806bap+2
//Hp = 0xb3bb870
t = *Hp;
//t = nan

更新 2： 如果代码由 icpc 编译，则不会出现此类问题。那么是否存在编译器特定的问题？

【问题讨论】：

这种代码格式很烦人。
请发布您从中获取日志条目的实际代码。另外，请确保Ho 没有超出范围。（Ho是什么类型？）
添加了带有日志记录的整个函数的代码。

标签： c gcc floating-point nan

【解决方案1】：

显然，-0.009034•0.000000 不应产生 NaN。因此，要么问题中呈现的代码和数据不能准确表示实际计算，要么计算实现存在缺陷。

如果我们假设硬件和基本计算实现没有缺陷，那么需要调查的一些可能性包括：

t 和 *Xp 的记录未能在乘法之前立即记录 t 和 *Xp 的正确值或在乘法之后立即记录 t 的正确值。
t 或*Xp 的值显示不正确。例如，用于显示 *Xp 的格式显示“0.000000”，即使 *Xp 具有其他值，例如 NaN。
Xp 指向不合适的地方，导致 *Xp 不可靠（例如，被外部操作更改）。
显示的代码不准确。例如，它来自已更改的旧源，或者它是新源但正在执行之前编译的代码。

注意：使用浮点对象进行调试时，您应该不使用诸如“%f”之类的格式进行打印，尤其是不要使用数字位数的默认值。您应该使用“%a”打印，它使用十六进制表示打印浮点值的精确值。您也可以在许多情况下使用“%.99g”，前提是您的 C 实现可以很好地将浮点值转换为十进制。

【讨论】：

将变量值添加为“%a”以发布。
我今天遇到的另一个原因是这个警告（不是在海报代码中，而是我自己的）：“函数'fabs'的隐式声明”，因为我忘记包含 math.h

【解决方案2】：

Eric Postpischil 的出色回答没有提到第五种可能性：

乘法正在 x87 寄存器中执行，并且由于（可能不相关的）程序执行中的较早操作而发生了浮点堆栈溢出。当处理器处于这种故障状态时，所有在 x87 寄存器上执行的计算都会产生 NaN 结果。

这两个最常见的原因是调用返回浮点结果的函数在范围内没有原型（对于许多调用约定，这将导致调用者无法将结果从 FP 堆栈中弹出），以及不正确的手写（可能是内联）程序集。

故障仅在经过一段时间后才发生的事实为这种可能性提供了一些证据；如果有一个很少使用的代码路径泄漏了浮点堆栈的一个元素，则需要在故障出现之前使用它一定次数，这可能让它直到现在才被注意到。

要诊断或排除这种可能性，您需要查看浮点状态寄存器 (FPSR) 的第 6 位 (SF)。根据您使用的编译器，检查 FPSR 的确切方法可能会有所不同。

【讨论】：

“范围内没有原型”是什么意思？在涉及复杂浮点数的应用程序中尝试将 Python 与 C 接口时，我遇到了类似的行为。我发现，当使用分配有 malloc NaN 值的内存的 C 代码时，会在计算过程中的某个时刻出现。但是当我使用带有 extern "C" 声明的 g++ 编译器并从 malloc/free 切换到 new/delete [] 时，程序运行完美！

【解决方案3】：

从Wiki，可以返回NaN的操作有如下三种：

1. Operations with a NaN as at least one operand.
2. Indeterminate forms
      The divisions 0/0 and ±∞/±∞
      The multiplications 0×±∞ and ±∞×0
      The additions ∞ + (−∞), (−∞) + ∞ and equivalent subtractions
      The standard has alternative functions for powers:
      The standard pow function and the integer exponent pown function define 0pow(0), 1pow(∞),
      and ∞pow(0) as 1.
      The powr function defines all three indeterminate forms as invalid operations and 
      so returns NaN.
3. Real operations with complex results, for example:
      The square root of a negative number.
      The logarithm of a negative number
      The inverse sine or cosine of a number that is less than −1 or greater than +1.

现在这应该可以帮助您自己解决问题。

【讨论】：

我看到了，但没有适合我的情况。
确保您对值 t、*Xp、Xp 执行的算术运算正好在 0 和 FLT_MAX 之间

【解决方案4】：

您必须打印每个计算的子结果 - 或使用 isnan() 函数在常规位置进行检查，并追踪其来源。这要么是一些“糟糕”的数学，要么是你一开始就在喂垃圾（未初始化的变量很可能是 NaN）

【讨论】：