数值分析作业 - 浮点运算和误差

Problem 1

说明

分析单精度计算 fl(9.4)−fl(9)−fl(0.4) 的结果，并进行计算机实践。

	Symbol	Exponent	mantissa
`9.4f`	+ `0`	3 + 127 = `10000010`	1`00101100110011001100110`
`9.0f`	+ `0`	3 + 127 = `10000010`	1`00100000000000000000000`
`9.4f-9.0f`	+ `0`	3 + 127 = `10000010`	0`00001100110011001100110`
Lsh 5	+ `0`	-2 + 127 = `01111101`	1`10011001100110011000000`
`0.4f`	+ `0`	-2 + 127 = `01111101`	1`10011001100110011001101`
`9.4f-9.0f-0.4f`	- `1`	-2 + 127 = `01111101`	0`00000000000000000001101`
Lsh 20	- `1`	-22 + 127 = `01101001`	1`10100000000000000000000`

得到计算结果为

编写辅助函数用于输出浮点数的比特

#define print_fpbits(x) print_fpbits_impl(x, #x)
template <typename T>
void print_fpbits_impl(T value, std::string expr) {
    auto y = reinterpret_cast<std::uint8_t*>(&value);
    std::cout << expr << ": " << std::endl;
    for (int i = sizeof(T) - 1; i >= 0; i--) {
        char byte = y[i];
        for (int j = 7; j >= 0; j--) {
            std::cout << ((byte >> j) & 1);
        }
    }
    std::cout << std::endl;
}

验证计算过程

print_fpbits(9.4f);
print_fpbits(9.0f);
print_fpbits(0.4f);
print_fpbits(9.4f - 9.0f);
print_fpbits(9.4f - 9.0f - 0.4f);
std::cout << std::format("{:.12f}\n", 9.4f - 9.0f - 0.4f);

得到输出

9.4f: 
01000001000101100110011001100110
9.0f:
01000001000100000000000000000000
0.4f:
00111110110011001100110011001101
9.4f - 9.0f:
00111110110011001100110011000000
9.4f - 9.0f - 0.4f:
10110100110100000000000000000000
-0.000000387430

Problem 2

说明

设计高效的多项式算法 , 计算机编程比较直接算法与优化算法的计算时间（采用双精度进行计算，分别循环次）

容易写出直接的多项式计算方法

const double coef[] = {1, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 4, 0, 0, 0, 5};

double pow_naive(double x, int n) {
    double res = 1.0;
    for (int i = 1; i <= n; i++) {
        res *= x;
    }
    return res;
}

double f_naive(double x) {
    double res = 0.;
    for (int i = 0; i <= 15; i++) {
        res += coef[i] * pow_naive(x, i);
    }
    return res;
}

和简单的优化

double f_fast(double x) {
    double res = 0.;
    double pow = 1.;
    for (int i = 0; i <= 15; i++) {
        res += coef[i] * pow;
        pow *= x;
    }
    return res;
}

需要注意的是, 针对 C++ 一类编译器具有强大静态优化能力的编程语言, 将简单的算法在同样的输入下单纯地重复运行并不是合理的性能测量方式, 因为编译器很容易将整个循环删除掉. 简单起见, 可以在循环中要求从 std::vector 中读取并存储数据, 以防止无用的循环被删除. 下面报告程序在 MSVC 19.33.31629.0 下, 简单循环次的用时

x = 2, naive: 2272.14ms
x = 2, fast: 418.473ms
x = 2.22... naive: 2245.62ms
x = 2.22..., fast: 426.92ms

启用 /O2 优化后, 性能有显著的提升

x = 2, naive: 686.996ms
x = 2, fast: 113.86ms
x = 2.22... naive: 670.233ms
x = 2.22..., fast: 108.462ms

Problem 3

说明

推导计算的递推公式,用正向和逆向递推计算，比较并分析两种迭代的误差和稳定性。

容易推导公式

初始值可以取

这个正向递推公式每次递推将误差放大五倍

1 2	forward = (NestList[{E - #[[2]] #[[1]], #[[2]] + 1} &, {1., 2}, 20] // Transpose)[[1]]; forward // InputForm

{1., 0.7182818284590451, 0.5634363430819098, 
 0.4645364561314058, 0.395599547802016, 0.34468454164694906, 
 0.30549003693040166, 0.27436153301583177, 0.24902803131655915, 
 0.22800151529345358, 0.21026516023105568, 0.19509990568637692, 
 0.18198305453614516, 0.17051906495301283, 0.1604958541638526, 
 0.15034816183740363, 0.16236307722318344, -0.2042535615582568, 
 6.599099498065924, -129.26370813285942, 2717.256152618507}

注意到在 16 项处由于误差的累积导致积分有界性已被破坏, 之后的迭代过程明显不收敛.

考虑反向迭代公式

取反向迭代

1
2
3

backward = (NestList[{1/#[[2]] (E - #[[1]]), #[[2]] - 1} &, {1., 30}, 
      29] // Transpose)[[1]] // Reverse;
backward // InputForm

{1., 0.7182818284590452, 0.5634363430819095, 
 0.4645364561314071, 0.3955995478020096, 0.3446845416469873, 
 0.30549003693013366, 0.27436153301797606, 0.24902803129726034, 
 0.2280015154864418, 0.21026515810818536, 0.19509993116082064, 
 0.18198272336837695, 0.17052370130176742, 0.16042630893253398, 
 0.15146088553850115, 0.14344677430452524, 0.13623989097759062, 
 0.12972389988482333, 0.12380383076257832, 0.11840138244490007, 
 0.1134514146712435, 0.10889929102044421, 0.10469884396838404, 
 0.10081072924944416, 0.09720286797349674, 0.09380439317463288, 
 0.09175881956932448, 0.057276060948634834, 1.}

与 NIntegrate 直接积分 (保证 20 位有效数字) 的结果比较并绘制图像

direct = 
 Table[NIntegrate[x^n E^x, {x, 0, 1}, AccuracyGoal -> 20, 
   WorkingPrecision -> 20], {n, 1, 20}];
ListLinePlot[{forward, backward, direct}]