如何纠正 Floor(t / dt) * dt != t 中的数值不稳定？

Question

我有两个

float

t

和

dt

。典型值为

t = 1

和

dt = .01f

。我想从

foo

到

0

的每个“时间”执行一个函数

t

，步长为

dt

。如果碰巧

t

不能被

dt

整除，我也需要为最后一个时间点

foo

调用

t

。

这就是我正在做的：

std::size_t const k = static_cast<std::size_t>(std::floor(t / dt));
for (std::size_t i = 0; i < k; ++i)
    foo(i * dt, dt);

float const s = k * dt,
    h = t - s;
if (h > 0)
    foo(s, h);

如您所见，我还需要将增量

dt

（分别为

h

）传递给

foo

。

现在我观察到，由于数值不准确，即使

h > 0

可以被

t

整除，

dt

，在这种情况下，

foo

也会在第二个参数中用一个很小的

h

来执行。这是我非常想避免的事情。是否有可能以数字更稳定的方式编写代码，以便当

h > 0

能被

t

整除时，它不会执行

dt

分支？

Answer 1

问题在于

0.01

无法用浮点数完美表示，因此当您将其乘以

100

时，它将等于

0.9999999776482582

。为了解决这个问题，您可以与机器 epsilon 除以

2

进行比较（因为指数是

-1

）。

这是一些 Python 示例代码：

from math import floor

from numpy import finfo, float32

t = float32(1)
dt = float32(.01)
k = floor(t / dt)
for i in range(k):
    print("Loop:", i * dt, dt)

s = k * dt
h = t - s
if h > finfo(float32).eps / 2:
    print("Leftover:", s, h)

其中

finfo(float32).eps / 2

=

5.960464477539063e-08

>

h

=

2.2351741790771484e-08