渲染引擎杂记（1）——浮点数误差分析及补偿

　　在各种渲染引擎中使用浮点数几乎可以肯定一定会产生浮点数误差，而渲染引擎的大量计算量也不允许使用其他特别高精度的浮点数，因而需要一定程度的精度补偿。事实上，浮点数并不适合用于精确计算，而适合进行科学计算。

float：2^23 = 8388608，一共七位，这意味着最多能有7位有效数字，但绝对能保证的为6位，也即float的精度为6~7位有效数字；
double：2^52 = 4503599627370496，一共16位，同理，double的精度为15~16位。

　　度量浮点数误差有两种形式：绝对误差（absolute）和相对误差（relative）。如果我们使用了浮点数计算得到了结果 a~ ，我们可以使用实际值和计算差的绝对值来计算绝对误差：

δa=|a~−a|

　　相对误差是指绝对误差和实际结果的比值，当 a≠0 时：δr=∣a~−aa∣=∣δaa∣　　定义了两种误差后，可以这样子表达：

a=a±δa=a(1±δr)　　下面尝试计算 a,b,c,d 四个数的和，如果我们这样计算：r=(((a+b)+c)+d)，加上误差可以表示为：(((a⊕b)⊕c)⊕d)⊂((((a+b)(1±ϵm))+c)(1±ϵm)+d)(1±ϵm)=(a+b)(1±ϵm)3+c(1±ϵm)2+d(1±ϵm)　　因为ϵm非常小，高阶的 ϵm 可以表示为 ϵm 的加法：(1±ϵm)n≤(1±(n+1)ϵm)　　（在实践中，(1±nϵm) 一般就可以表示误差的上界，因为高阶的 ϵm 下降得非常快）

　　所以，误差的方程可以改写成如下所示：(a+b)(1±4ϵm)+c(1±3ϵm)+d(1±2ϵm)=a+b+c+d+[±4ϵm(a+b)±3ϵmc±2ϵmd]　　就可以得到绝对误差的上界：4ϵm∣a+b∣+3ϵm∣c∣+2ϵm|d∣　　因此，很容易得到 a,b,c,d 之和的误差上界。

　　得到结果中，a+b 的误差贡献比 d的误差贡献更大，这也说明了如果将大量浮点数加在一起，从小到大排通常比任意排序的误差更小。

　　上述分析假定了编译器会根据表达式的顺序来计算和，因此必须要使编译器遵循运算顺序来最小化误差范围。

　　如果改变一下相加顺序得到 floatr=(a+b)+(c+d)的误差，同样适用上述过程，可以得到误差上界为3ϵm∣a+b∣+3ϵm∣c+d∣　　如果 |a+b| 相对比较大，此误差会比较小，但如果 |d| 相对较大，误差也可能更大。

　　这种分析误差的方法称为前向误差分析（forward error analysis）：给定输入，可以用相似的方法计算结果的误差上界。得到的误差上界一般会比实际误差更大，但是在实际运算中，误差项的符号往往是混合出现的，所以当相加时会产生误差抵消。

　　所以需要一种新的误差分析方法——后向误差分析，将计算结果暂时视为精确，然后给出相同输入产生结果波动的上界。

　　(1±ϵm)n的保守误差上界(1±(n+1)ϵm)在误差阶数较高时不太满足。Higham 得出了一种能够更加紧逼 (1±ϵm)乘积误差项。

　　如果有(1±ϵm)n，那么可以值的上界为 1+θn，当nϵm<1时有∣θn∣≤nϵm1−nϵm　　在上式中分母比1小n个ϵm误差，所以几乎没有增大nϵm 就得到了一个保守上界。如果用γn 来表示上界：γn=nϵm1−nϵm　　使用 γ 标记就可以得到四个值和的误差上界（这里转换更紧逼的(1±nϵm)），有：∣a+b∣γ3+∣c∣γ2+∣d∣γ1　　这种方法的一个优点是可以计算商的误差，如果有 (1±ϵm)m(1±ϵm)n，可以得到误差1±γm+n。

　　如果计算两个误差的乘积，可以得到误差a(1±γi)⊗b(1±γj)⊂ab(1±γi+j+1)　　故相对误差为∣abγi+j+1ab∣=γi+j+1　　这说明乘法的误差积累非常小。

　　可是对于加减法来说，误差会被积累得很大。对于a(1±γi)⊕b(1±γj)，有误差为∣a∣γi+1+∣b∣γj+1　　所以对于精度很敏感的部分，对加法一定要进行补偿。

　　在具体的实现中可以使用 constexpr 关键字来指定它为关键字常量：

inline constexpr Float gamma(int n) {
return (n * MachineEpsilon) / (1 - n * MachineEpsilon);
}

　　以下是 pbrt 中对光线变换坐标时做精度补偿的例子：

inline Ray Transform::operator()(const Ray &r) const {
Vector3f oError;
Point3f o = (*this)(r.o, &oError);
Vector3f d = (*this)(r.d);
// Offset ray origin to edge of error bounds and compute _tMax_
Float lengthSquared = d.LengthSquared();
Float tMax = r.tMax;
if (lengthSquared > 0) {
Float dt = Dot(Abs(d), oError) / lengthSquared;
o += d * dt;
tMax -= dt;
}
return Ray(o, d, tMax, r.time, r.medium);
}

　　使用了重载的小括号来计算变换后的原点 origin 的位置和浮点数绝对误差，计算误差的方式很简单。

template <typename T>
inline Point3<T> Transform::operator()(const Point3<T> &p,
Vector3<T> *pError) const {
T x = p.x, y = p.y, z = p.z;
// Compute transformed coordinates from point _pt_
T xp = m.m[0][0] * x + m.m[0][1] * y + m.m[0][2] * z + m.m[0][3];
T yp = m.m[1][0] * x + m.m[1][1] * y + m.m[1][2] * z + m.m[1][3];
T zp = m.m[2][0] * x + m.m[2][1] * y + m.m[2][2] * z + m.m[2][3];
T wp = m.m[3][0] * x + m.m[3][1] * y + m.m[3][2] * z + m.m[3][3];

// Compute absolute error for transformed point
T xAbsSum = (std::abs(m.m[0][0] * x) + std::abs(m.m[0][1] * y) +
std::abs(m.m[0][2] * z) + std::abs(m.m[0][3]));
T yAbsSum = (std::abs(m.m[1][0] * x) + std::abs(m.m[1][1] * y) +
std::abs(m.m[1][2] * z) + std::abs(m.m[1][3]));
T zAbsSum = (std::abs(m.m[2][0] * x) + std::abs(m.m[2][1] * y) +
std::abs(m.m[2][2] * z) + std::abs(m.m[2][3]));
*pError = gamma(3) * Vector3<T>(xAbsSum, yAbsSum, zAbsSum);
CHECK_NE(wp, 0);
if (wp == 1)
return Point3<T>(xp, yp, zp);
else
return Point3<T>(xp, yp, zp) / wp;
}

　　上述代码值得注意的是gamma函数的参数为3，表达的意思也就是3ϵm。

参考资料：Physically Based Rendering