给定范围内数字的浮点仿真或定点答案

【问题标题】：Floating point emulation or Fixed Point for numbers in a given range给定范围内数字的浮点仿真或定点
【发布时间】：2013-04-18 11:10:17
【问题描述】：

我有一个不支持浮点的协处理器。我尝试使用 32 位固定点，但它无法处理非常小的数字。我的数字范围从 1 到 1e-18。一种方法是使用浮点仿真，但它太慢了。在我们知道数字不会大于 1 且小于 1e-18 的情况下，我们能否让它更快。或者有没有办法让修复点在非常小的数字上起作用。

【问题讨论】：

像往常一样，定义您要解决的问题。给我们您要计算的公式。
用数字计算（乘法、减法和加法），其中初始值、中间值和最终值始终在 1e-18 到 1 的范围内，即永远不会大于 1 且永远不会小于 1e -18.
你想要什么精度？
您可以在定点算术中将点放置在您想要的任何位置。唯一的窍门是，如果这个点是出路的，就像你的情况一样，提供标准的乘法或除法是没有意义的——这会自动产生超出范围的结果。您可以使用 constmul(a,b) 和 constdiv(a,b) 来代替，它们返回 a * b * const 和 a / (b * const) 的定点表示，用于精心选择的 const。
@Pascal Cuoq - 是的，我已经为我的固定点计算移动了点，以便它最适合小于 1 的数字，但 1e-18 仍然太小而无法以 32 位表示整数。而且像你提到的缩放问题会使数字大于 1，我们将再次不得不改变这一点，所以我认为它不会起作用。

标签： c floating-point fixed-point

【解决方案1】：

32 位定点编码不可能表示从 10^–18 到 1 的数字。从 10^{-18 的跨度这一事实可以立即看出这一点} 是 10¹⁸ 的比率，但是 32 位整数的非零编码跨度小于 2³² 的比率，这要小得多超过 10¹⁸。因此，定点编码的尺度选择无法提供所需的跨度。

因此 32 位定点编码将不起作用，您必须使用其他技术。

在某些应用中，可能适合使用多个定点编码。也就是说，各种输入值将使用定点编码进行编码，但每个值都有适合它的比例，中间值和输出也将具有自定义比例。显然，只有在设计时可以确定合适的比例时，这才有可能。否则，您应该放弃 32 位定点编码并考虑替代方案。

【讨论】：

他可以使用定点对数系统，尽管在这种系统中加法很烦人。
64 位定点算法应该可以工作。以我的经验，软件中的 LNS 仅适用于计算加法相对较少且精度要求较低的情况，否则加减法所需的表会变得很大，这对于没有 FPU 的低端处理器（即也往往有小缓存）。研究 LNS 的一个合理起点是“FOCUS 微型计算机数字系统”，CACM，1979。
@njuffa：尚不清楚 64 位定点算术是否可行，因为 MetallicPriest 尚未回答有关所需精度的问题。 64 位的跨度可以表示 1e-18 的幅度，但它几乎没有为额外的精度留下空间。
很公平，请注意我使用了黄鼠狼这个词“应该”:-)。另一方面，我没有看到明确的“定点不起作用”适用。根据具体情况，在操作数较小的情况下进行一些动态重新缩放的定点可能会起作用。当然，重新缩放的实例应该相对较少，否则性能会超出窗口。
@njuffa：我更新了答案以将其限制为 32 位定点。我是在问题“我尝试使用 32 位固定点”的上下文中写的。

【解决方案2】：

简化的 24 位浮点是否足够快和足够准确？：

#include <stdio.h>
#include <limits.h>

#if UINT_MAX >= 0xFFFFFFFF
typedef unsigned myfloat;
#else
typedef unsigned long myfloat;
#endif

#define MF_EXP_BIAS 0x80

myfloat mfadd(myfloat a, myfloat b)
{
  unsigned ea = a >> 16, eb = b >> 16;
  if (ea > eb)
  {
    a &= 0xFFFF;
    b = (b & 0xFFFF) >> (ea - eb);
    if ((a += b) > 0xFFFF)
      a >>= 1, ++ea;
    return a | ((myfloat)ea << 16);
  }
  else if (eb > ea)
  {
    b &= 0xFFFF;
    a = (a & 0xFFFF) >> (eb - ea);
    if ((b += a) > 0xFFFF)
      b >>= 1, ++eb;
    return b | ((myfloat)eb << 16);
  }
  else
  {
    return (((a & 0xFFFF) + (b & 0xFFFF)) >> 1) | ((myfloat)++ea << 16);
  }
}

myfloat mfmul(myfloat a, myfloat b)
{
  unsigned ea = a >> 16, eb = b >> 16, e = ea + eb - MF_EXP_BIAS;
  myfloat p = ((a & 0xFFFF) * (b & 0xFFFF)) >> 16;
  return p | ((myfloat)e << 16);
}

myfloat double2mf(double x)
{
  myfloat f;
  unsigned e = MF_EXP_BIAS + 16;
  if (x <= 0)
    return 0;
  while (x < 0x8000)
    x *= 2, --e;
  while (x >= 0x10000)
    x /= 2, ++e;
  f = x;
  return f | ((myfloat)e << 16);
}

double mf2double(myfloat f)
{
  double x;
  unsigned e = (f >> 16) - 16;
  if ((f & 0xFFFF) == 0)
    return 0;
  x = f & 0xFFFF;
  while (e > MF_EXP_BIAS)
    x *= 2, --e;
  while (e < MF_EXP_BIAS)
    x /= 2, ++e;
  return x;
}

int main(void)
{
  double testConvData[] = { 1e-18, .25, 0.3333333, .5, 1, 2, 3.141593, 1e18 };
  unsigned i;
  for (i = 0; i < sizeof(testConvData) / sizeof(testConvData[0]); i++)
    printf("%e -> 0x%06lX -> %e\n",
           testConvData[i],
           (unsigned long)double2mf(testConvData[i]),
           mf2double(double2mf(testConvData[i])));

  printf("300 * 5 = %e\n", mf2double(mfmul(double2mf(300),double2mf(5))));
  printf("500 + 3 = %e\n", mf2double(mfadd(double2mf(500),double2mf(3))));
  printf("1e18 * 1e-18 = %e\n", mf2double(mfmul(double2mf(1e18),double2mf(1e-18))));
  printf("1e-18 + 2e-18 = %e\n", mf2double(mfadd(double2mf(1e-18),double2mf(2e-18))));
  printf("1e-16 + 1e-18 = %e\n", mf2double(mfadd(double2mf(1e-16),double2mf(1e-18))));

  return 0;
}

输出（ideone）：

1.000000e-18 -> 0x459392 -> 9.999753e-19
2.500000e-01 -> 0x7F8000 -> 2.500000e-01
3.333333e-01 -> 0x7FAAAA -> 3.333282e-01
5.000000e-01 -> 0x808000 -> 5.000000e-01
1.000000e+00 -> 0x818000 -> 1.000000e+00
2.000000e+00 -> 0x828000 -> 2.000000e+00
3.141593e+00 -> 0x82C90F -> 3.141541e+00
1.000000e+18 -> 0xBCDE0B -> 9.999926e+17
300 * 5 = 1.500000e+03
500 + 3 = 5.030000e+02
1e18 * 1e-18 = 9.999390e-01
1e-18 + 2e-18 = 2.999926e-18
1e-16 + 1e-18 = 1.009985e-16

减法留作练习。同上以获得更好的转换例程。

【讨论】：

【解决方案3】：

使用 64 位定点并完成它。

与 32 位定点相比，它的乘法速度会慢四倍，但仍然比浮点仿真效率高得多。

【讨论】：

【解决方案4】：

在嵌入式系统中，我建议使用 16+32、16+16、8+16 或 8+24 位冗余浮点表示，其中每个数字只是 M * 2^exp。

在这种情况下，您可以选择用 M=0 和 exp=0 来表示零；每个 2 的幂有 16-32 种表示形式——这主要使比较比通常更难。也可以推迟标准化，例如减法后。

【讨论】：