2

我正在尝试将浮点双精度值转换x为具有 12 个（正确舍入）有效数字的十进制值。我假设它x在 10^110 和 10^111 之间，因此它的十进制表示形式为x.xxxxxxxxxxxE110. 而且，只是为了好玩，我只尝试使用浮点运算。

我到达了下面的伪代码，其中所有操作都是双精度操作，符号1e98是最接近数学 10^98 的双精度数，1e98_2是最接近数学减法 10^98- 的结果的双精度数1e98。该符号适用于带有操作数, ,fmadd(X * Y + Z)的融合乘加运算。XYZ

  y = x * 2^-1074;    // exact
  q = y / 1e98;       // q is denormal and the significand of q interpreted
                      // as an integer is our candidate for the 12 decimal
                      // digits of x

  r = fmadd(q * 1e98 - y);  // close to 1e98 * (error made during the division)

  // If 1e98_2 >= 0, we divided by a number that was smaller than we wished
  // The correct answer may be q or q+1.

  if (r and 1e98_2 have opposite signs)
  {
    return the significand of q;
  }

  s = copysign(2^-1074, r);
  r1 = abs(r);
  r2 = abs(1e98_2);

  h = 1e98 * 0.5 * 2^-1074;

  Set rounding mode to downwards

  r3 = fmadd(r2 * q + r1);

  if (r3 < h)
  {
    return the significand of q;
  }
  else
  {
    return significand of (q + s)
  }

对于上述伪代码所造成的混乱，我深表歉意，但对我来说还不是很清楚，因此有以下问题：

  1. 第一个 fmadd 是否按预期工作（计算 1e98 *（除法期间出错））？

  2. 标志。我无法说服自己他们是对的。但我也无法说服自己他们是错的。

  3. 关于这个算法可能产生错误结果的频率的任何想法，也许是争论？

  4. 如果它确实有效，如果将“q = y / 1e98”更改为“q = y * 1e-98”（保持所有其他指令相同），算法是否有可能继续工作？

我没有测试过这个算法。我没有任何带有 fmadd 指令的计算机，尽管我希望能找到一台这样我就可以执行上述操作。

4

1 回答 1

2

让y/d是精确的操作，并将q=rnd(y/d)结果四舍五入到最接近的浮点数。
那么真正的误差乘以 d 是rt=(rnd(y/d)-y/d)*d=q*d-y，我们用 fmadd 执行的操作是r=rnd(q*d-y)
为什么q*d-y是精确的（fmadd 不进行最终舍入）不太清楚解释，但假设q*d位数有限（<nbits(q)+nbits(d)），的指数y是q*d(+/- 1) 并且由于错误是|rt|<0.5*ulp(q)*d，这意味着第一个nbits(q)正在消失...这就是问题 1 的答案。

所以q*1e98 - y = r，哪里|r|*2^1074 <= 0.5e98 < 5*10^98（第二个不等式是幸运的）

q*(10^98) - y = r + (10^98-1e98)*q其中|10^98-1e98|*q*2^1074 <= 0.5e95（假设至少 15 位精度，log(2^53)/log(10) > 15）

所以你问是否|q*(10^98)-y|*2^1074>5*10^97

你有一个近似值|q*(10^98)-y|是r+1e98_2*q

因为|r| < 5*10^98, 并且|r+(10^98-1e98)*q|<|r|如果符号相反，我认为这对问题 2 的回答是肯定的。但我不太确定 1e98_2 是否 < 0。

如果r和1e98_2具有相同的符号，它可能会超过5*10^97，因此您将进一步处理r3 = 1e98_2*q + rvs的讨论h=0.5e98*2^-1074

对于问题 3，乍一看，我想说两件事可能会使算法失败：

  • 1e98_2不准确（10^98-1e98-1e98_2 = -3.6e63大约）

  • 并且h比ht=0.5*10^98*2^-1074我们上面看到的要小一点。

真正的错误r3t大约是(1e98_2-3e63)*q + r < r3（并且只有 >0 的情况我们才感兴趣，因为 1e98_2>0）。

因此，当真实误差 r3t 低于真实关系 ht 时，误差 r3 的近似值落在近似关系 h 之上可能会导致不正确的舍入。是否有可能，如果是，您的问题 3 的频率如何？

为了减轻上述不平等风险，您尝试截断 r3 的大小，因此r3 <= 1e98_2*q + r. 对误差范围进行真正的分析我觉得有点累......

So I scanned for an error, and the first failing example I found was 1.0000000001835e110 (I assume correctly rounded to nearest double, but it is in fact 1000000000183.49999984153799821120915424942630528225695526491963291846957919215885146546696544423465444842668032e98).

在这种情况下，r和1e98_2具有相同的符号，并且

  • (x/1e98) > 1000000000183.50000215

  • q有效数字因此四舍五入为1000000000184

  • r3>h（r3*2^1074大约是 5.000001584620017e97）并且我们错误地增加q+s了，当它应该是q-s，绝对是一个错误。

我的回答是：

  1. 是的，r=fmadd(q * 1e98 - y)正好是 1e98*（除法时出错），但我们不关心除法，它只是提供一个猜测，重要的是减法是准确的。

  2. 是的，符号是正确的，因为|r| < 5*10^98,|r+(10^98-1e98)*q|<|r|如果符号相反。但我不太确定 1e98_2 是否 < 0。

  3. 以第一个失败的例子为例(1.0000000001835e110 - 1.0e110)/1.0e110 ulp -> 1.099632e6，一个非常非常天真的猜想是说百万分之一的情况下，r3 下降到 h... 所以一旦 q+s 修正为 qs，r3>hwhiler3t<ht的出现远小于 1/ 1,000,000 无论如何......在感兴趣的范围内有超过 10^15 双打，所以认为这不是一个严肃的答案......

  4. 是的，上面的讨论只是关于猜测 q，与它的产生方式无关，并且 1. 中的减法仍然是准确的......

于 2014-02-14T04:43:15.883 回答