第6章数值考虑 - 浮点数的陷阱 (Numerical Considerations - Pitfalls of Floating-Point Arithmetic)

一、导读

6.1 本章学习目标

通过本章的学习，读者将能够：

深入理解IEEE 754浮点数标准的表示方法，包括单精度、双精度和半精度的格式

掌握浮点数运算的舍入误差来源，能够估算累积误差的范围

理解GPU上不同精度浮点运算的性能差异，合理选择数据类型

掌握数值稳定性的基本概念，能够识别和改进不稳定的算法

学会处理浮点数比较、除零、溢出和下溢等常见问题

理解归约操作中的数值误差累积问题，掌握 compensated summation 技术

了解混合精度计算的原理和应用，能够在保持精度的同时提升性能

掌握CUDA中特殊数学函数（intrinsics）的使用，在精度和性能之间权衡

6.2 为什么本章如此重要

数值计算是GPU应用的核心领域之一，从科学模拟到深度学习，从金融建模到图像处理，几乎所有GPU应用都涉及浮点数运算。然而，浮点数并非实数的完美表示，它存在精度限制、舍入误差、溢出和下溢等问题。忽视这些数值特性可能导致程序产生错误结果，甚至在极端情况下完全失效。

在GPU编程中，数值问题更加复杂。首先，GPU为了追求高性能，在某些情况下允许违反严格的IEEE 754合规性。例如，默认情况下CUDA可能使用快速近似数学函数，牺牲精度换取速度。其次，GPU的大规模并行性意味着大量的浮点运算同时进行，舍入误差的累积模式与串行程序不同。第三，混合精度计算在GPU上越来越普遍，如何在不同精度之间转换而不损失精度是一个重要技能。

对于科学计算应用，数值准确性直接关系到结果的可信度。一个看似微小的舍入误差，在经过数百万次迭代后可能被放大到不可接受的程度。对于深度学习应用，虽然神经网络对数值误差不太敏感，但不当的精度选择可能导致训练不稳定或推理精度下降。对于金融应用，数值错误可能导致直接的经济损失。

本章将系统地介绍浮点数表示、运算特性、误差分析和稳定性优化，帮助读者建立起扎实的数值计算基础。我们不仅会讲解理论知识，还会提供实用的编程技巧和代码示例，让读者能够在实际项目中避免常见的数值陷阱。

6.3 核心问题引导

IEEE 754浮点数标准是如何表示实数的？单精度、双精度和半精度的区别是什么？

什么是机器epsilon？它如何影响浮点运算的精度？

GPU上单精度和双精度运算的性能差异有多大？何时应该使用双精度？

什么是数值稳定性？如何判断一个算法是否数值稳定？

浮点数比较为什么不能直接使用==？正确的比较方法是什么？

什么是Kahan求和？它如何减少归约操作中的累积误差？

混合精度训练的原理是什么？如何在保持模型精度的同时使用低精度加速？

CUDA的快速数学函数（如__sinf、__expf）与标准函数有何区别？何时可以使用？

二、核心概念详解

6.1 IEEE 754浮点数标准

IEEE 754是浮点数算术的国际标准，定义了浮点数的表示格式、运算规则和特殊值处理。理解IEEE 754是进行正确数值计算的基础。

浮点数的表示：

一个浮点数由三部分组成：符号位（sign）、指数位（exponent）和尾数位（mantissa/fraction）。其值为：

$$value = (-1)^{sign} \times 1.mantissa \times 2^{exponent - bias}$$

其中bias是指数偏移量，用于表示负指数。

精度格式对比：

格式	总位数	符号位	指数位	尾数位	Bias	最小正正规数	最大有限数	机器epsilon
FP16 (半精度)	16	1	5	10	15	2^-14 ≈ 6.1×10^-5	65504	2^-10 ≈ 9.8×10^-4
FP32 (单精度)	32	1	8	23	127	2^-126 ≈ 1.2×10^-38	3.4×10^38	2^-23 ≈ 1.2×10^-7
FP64 (双精度)	64	1	11	52	1023	2^-1022 ≈ 2.2×10^-308	1.8×10^308	2^-52 ≈ 2.2×10^-16
TF32	32*	1	8	10+	127	同FP32	同FP32	2^-10 ≈ 9.8×10^-4

*TF32是NVIDIA Ampere引入的格式，使用FP32的存储但只有约10位有效数字。

特殊值：

IEEE 754定义了几个特殊值：

零：有正零（+0）和负零（-0），它们在大多数运算中行为相同，但在某些情况下（如除法）会产生不同的结果。

非规格化数（Subnormal/Denormal）：当指数为全0时，表示非常接近零的数。非规格化数的有效数字前导不是1而是0，这使得它们可以表示更小的数，但精度降低。处理非规格化数通常较慢，因为需要特殊的硬件路径。

无穷大（Infinity）：当指数为全1且尾数为全0时，表示无穷大。正无穷和负无穷分别由符号位决定。无穷大参与运算时有明确的规则：∞ + x = ∞，∞ × ∞ = ∞，x / ∞ = 0。

NaN（Not a Number）：当指数为全1且尾数非零时，表示NaN。NaN由无效运算产生，如0/0、√(-1)、∞ - ∞。NaN的传播规则是：任何涉及NaN的运算结果都是NaN。

GPU上的浮点表示：

在CUDA中，可以使用联合体来查看浮点数的二进制表示：

cpp

union FloatBits {
    float f;
    uint32_t bits;
};

// Print the binary representation of a float
void printFloatBits(float value) {
    FloatBits fb;
    fb.f = value;
    printf("Float %f: 0x%08X\n", value, fb.bits);
    
    uint32_t sign = (fb.bits >> 31) & 1;
    uint32_t exponent = (fb.bits >> 23) & 0xFF;
    uint32_t mantissa = fb.bits & 0x7FFFFF;
    
    printf("  Sign: %u, Exponent: %u (biased), Mantissa: 0x%06X\n", 
           sign, exponent, mantissa);
}

6.2 舍入误差与精度限制

舍入模式：

IEEE 754定义了四种舍入模式：

向最近舍入（Round to Nearest, Ties to Even）：默认模式。结果舍入到最接近的可表示值，如果恰好在中间，则舍入到尾数为偶数的那个值。这种模式统计上无偏。

向零舍入（Round toward Zero）：截断多余位，结果总是向零靠近。

向正无穷舍入（Round toward +∞）：结果总是向上舍入。

向负无穷舍入（Round toward -∞）：结果总是向下舍入。

在CUDA中，可以使用__fdivide_rn()、__fdivide_rz()等函数指定舍入模式。

机器Epsilon：

机器epsilon（ε）是1.0与大于1.0的最小可表示浮点数之间的差。它代表了浮点数的相对精度上限。

FP32: ε ≈ 1.2 × 10^-7（约7位十进制有效数字）
FP64: ε ≈ 2.2 × 10^-16（约16位十进制有效数字）

机器epsilon决定了浮点运算的精度极限。两个相差小于ε的数可能无法区分。

舍入误差的来源：

表示误差：许多十进制小数无法精确表示为二进制浮点数。例如，0.1在二进制中是无限循环小数，存储为FP32时会有舍入误差。

cpp

float x = 0.1f;
printf("%.20f\n", x);  // Output: 0.10000000149011611938

运算误差：每次浮点运算都可能引入舍入误差。加法、减法、乘法、除法都可能产生无法精确表示的结果。

抵消误差（Catastrophic Cancellation）：当两个相近的数相减时，有效数字大量丢失，相对误差急剧增大。这是数值计算中最危险的情况之一。

cpp

float a = 1.0000001f;
float b = 1.0000000f;
float c = a - b;  // Result: 0.0000001, but with significant relative error

在这个例子中，a和b都有约7位有效数字，但它们的差只有1位有效数字。相对误差从约10^-7放大到约10^-1。

误差传播：

舍入误差会在后续运算中传播和累积。对于加法，绝对误差相加；对于乘法，相对误差相加。

$$|fl(a + b) - (a + b)| \leq |a|\epsilon + |b|\epsilon$$

$$\frac{|fl(a \times b) - (a \times b)|}{|a \times b|} \leq 2\epsilon$$

其中fl(x)表示x的浮点表示。

6.3 GPU上的浮点性能

单精度 vs 双精度：

在消费级GPU上，双精度性能通常远低于单精度。这是因为NVIDIA刻意限制了消费级GPU的双精度吞吐量，以区分专业市场。

GPU型号	FP32 TFLOPS	FP64 TFLOPS	FP64/FP32比率
RTX 4090	82.6	1.3	1/64
A100	19.5	9.7	1/2
V100	15.7	7.8	1/2

对于科学计算应用，如果需要双精度，应选择Tesla/Ampere等专业GPU。

半精度的优势：

FP16的计算速度可以是FP32的两倍（使用Tensor Core时更高），内存占用减半。这在深度学习中特别有用，因为神经网络对数值精度要求不高。

Tensor Core与混合精度：

Volta及以后的GPU引入了Tensor Core，专门加速矩阵乘法。Tensor Core支持多种精度组合：

FP16输入 → FP32累加
BF16输入 → FP32累加
TF32输入 → FP32累加
INT8/INT4输入 → INT32累加

Tensor Core可以在一次操作中完成4×4矩阵的乘加，吞吐量远超普通CUDA Core。

6.4 数值稳定性

数值稳定性的定义：

一个算法是数值稳定的，如果它在存在舍入误差的情况下仍能产生接近精确解的结果。相反，数值不稳定的算法会将小的舍入误差放大为大的输出误差。

稳定性分析示例：

考虑计算方差的问题。有两种常见的方法：

方法1：两遍算法（稳定）

cpp

// First pass: compute mean
float sum = 0.0f;
for (int i = 0; i < N; i++) {
    sum += data[i];
}
float mean = sum / N;

// Second pass: compute variance
float varSum = 0.0f;
for (int i = 0; i < N; i++) {
    float diff = data[i] - mean;
    varSum += diff * diff;
}
float variance = varSum / N;

方法2：单遍公式（不稳定）

cpp

float sum = 0.0f;
float sumSq = 0.0f;
for (int i = 0; i < N; i++) {
    sum += data[i];
    sumSq += data[i] * data[i];
}
float variance = sumSq / N - (sum / N) * (sum / N);

方法2的问题是sumSq / N和(sum / N)^2可能非常接近，导致抵消误差。当数据值的范围很大时，这个问题尤其严重。

Welford在线算法（稳定）：

cpp

float mean = 0.0f;
float M2 = 0.0f;
for (int i = 0; i < N; i++) {
    float delta = data[i] - mean;
    mean += delta / (i + 1);
    float delta2 = data[i] - mean;
    M2 += delta * delta2;
}
float variance = M2 / N;

Welford算法在单次遍历中同时计算均值和方差，且数值稳定。

6.5 浮点数比较与特殊值处理

浮点数比较：

直接使用==比较浮点数是不可靠的，因为舍入误差可能导致理论上相等的数在实际表示中有微小差异。

正确的比较方法是检查差的绝对值是否小于某个容差：

cpp

// Absolute tolerance comparison
bool approximatelyEqual(float a, float b, float epsilon = 1e-6f) {
    return fabsf(a - b) < epsilon;
}

// Relative tolerance comparison (better for large numbers)
bool relativelyEqual(float a, float b, float relEpsilon = 1e-6f) {
    float maxVal = fmaxf(fabsf(a), fabsf(b));
    return fabsf(a - b) <= maxVal * relEpsilon;
}

// Combined comparison
bool isEqual(float a, float b, float absEps = 1e-8f, float relEps = 1e-6f) {
    float diff = fabsf(a - b);
    if (diff <= absEps) return true;  // Close to zero
    return diff <= fmaxf(fabsf(a), fabsf(b)) * relEps;
}

除零处理：

除以零会产生无穷大或NaN：

cpp

float safeDivide(float numerator, float denominator) {
    if (fabsf(denominator) < 1e-10f) {
        return 0.0f;  // Or handle as appropriate for your application
    }
    return numerator / denominator;
}

溢出和下溢：

cpp

// Check for overflow before multiplication
float safeMultiply(float a, float b) {
    float result = a * b;
    if (isinf(result)) {
        // Handle overflow
        return (a > 0 && b > 0) || (a < 0 && b < 0) ? FLT_MAX : -FLT_MAX;
    }
    if (result == 0.0f && a != 0.0f && b != 0.0f) {
        // Underflow detected
        return 0.0f;
    }
    return result;
}

检测特殊值：

cpp

#include <cmath>

if (isnan(value)) {
    // Handle NaN
}
if (isinf(value)) {
    // Handle infinity
}
if (isfinite(value)) {
    // Normal finite number
}

6.6 Kahan求和与补偿技术

问题：

在归约操作（如求和）中，简单的累加会导致舍入误差累积。当累加一个大数和一个小数时，小数的低位可能被舍去。

Kahan求和算法：

Kahan求和通过维护一个补偿变量来跟踪舍去的低位，显著减少累积误差。

cpp

// Serial Kahan summation
float kahanSum(const float* data, int N) {
    float sum = 0.0f;
    float compensation = 0.0f;  // Running compensation for lost low-order bits
    
    for (int i = 0; i < N; i++) {
        float y = data[i] - compensation;
        float t = sum + y;
        compensation = (t - sum) - y;  // Recover what was lost
        sum = t;
    }
    
    return sum;
}

补偿变量的工作原理：

y = data[i] - compensation：从当前数据中减去之前的补偿

t = sum + y：执行加法

compensation = (t - sum) - y：计算实际加到sum中的值与预期值的差，即舍去的低位

并行Kahan求和：

在GPU上实现并行Kahan求和较为复杂，因为补偿变量需要在归约树中传递。一种实用方法是：

每个线程使用Kahan求和计算局部和

使用标准归约合并局部和

cpp

__global__ void parallelKahanReduce(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x * 2 + threadIdx.x;
    
    // Each thread does Kahan sum of two elements
    float sum = 0.0f;
    float comp = 0.0f;
    
    for (int i = 0; i < 2 && idx + i * blockDim.x < N; i++) {
        float y = input[idx + i * blockDim.x] - comp;
        float t = sum + y;
        comp = (t - sum) - y;
        sum = t;
    }
    
    shared[tid] = sum;
    __syncthreads();
    
    // Standard tree reduction
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    if (tid == 0) {
        atomicAdd(output, shared[0]);
    }
}

虽然这不是完全的Kahan求和（归约阶段没有补偿），但比简单累加要好得多。

6.7 混合精度计算

混合精度训练：

在深度学习中，混合精度训练使用FP16进行前向和后向传播的大部分计算，但使用FP32保存权重和累加梯度。这种方法可以显著加速训练，同时保持模型精度。

关键技巧：

Loss Scaling：FP16的动态范围较小，梯度可能下溢为零。通过将loss乘以缩放因子（如2^16），将梯度值移入FP16的可表示范围。

FP32 Master Weights：权重更新使用FP32精度，避免FP16的精度不足导致权重停滞。

Mixed Precision Accumulation：矩阵乘法使用FP16输入，但累加使用FP32，保持数值稳定性。

CUDA中的混合精度：

cpp

// Using half precision in CUDA
#include <cuda_fp16.h>

__global__ void mixedPrecisionKernel(const __half* input, __half* output, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        // Convert to float for computation
        float val = __half2float(input[idx]);
        
        // Perform computation in float
        float result = val * val + 1.0f;
        
        // Convert back to half for storage
        output[idx] = __float2half(result);
    }
}

// Using half2 for vectorized operations
__global__ void half2Kernel(const __half2* input, __half2* output, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        // Process two half values at once
        __half2 val = input[idx];
        __half2 result = __hmul2(val, val);  // Element-wise multiply
        output[idx] = result;
    }
}

6.8 CUDA快速数学函数

CUDA提供了快速近似数学函数，牺牲精度换取速度：

标准函数	快速版本	精度	速度提升
`sinf()`	`__sinf()`	~2^-23	2-4x
`cosf()`	`__cosf()`	~2^-23	2-4x
`expf()`	`__expf()`	~2^-23	2-4x
`logf()`	`__logf()`	~2^-23	2-4x
`powf()`	`__powf()`	~2^-23	2-4x
`sqrtf()`	`__fsqrt_rn()`	IEEE compliant	1.5-2x

使用快速函数时需要评估精度要求。对于图形渲染和深度学习，快速函数通常足够精确；对于科学计算，可能需要标准函数。

启用--use_fast_math编译器选项可以将所有标准数学函数替换为快速版本，但需谨慎使用。

三、程序员的视角

3.1 在实际项目中的应用

科学模拟：在流体动力学模拟中，数值稳定性至关重要。使用双精度可以避免长时间积分后的误差累积。Kahan求和用于计算守恒量（如质量、能量），确保数值守恒。

深度学习：混合精度训练已成为标准做法，可以加速2-3倍而几乎不影响模型精度。Tensor Core的使用需要特定的数据布局和对齐要求。

金融计算：期权定价使用蒙特卡洛模拟，需要双精度确保定价准确性。舍入误差可能导致套利机会的错误估计。

图像处理：HDR图像处理使用FP32避免色调映射时的带状伪影。快速数学函数用于实时滤镜，精度损失不可察觉。

3.2 常见陷阱与误区

假设浮点数精确：认为0.1 + 0.2 == 0.3，实际上由于表示误差，这个等式为假。

忽视抵消误差：直接计算sumSq/N - mean*mean作为方差，当数据范围大时产生巨大误差。

过度使用双精度：在不需要高精度的场景使用双精度，浪费性能和内存。应先分析问题的条件数。

快速函数滥用：在需要高精度的地方使用__sinf等快速函数，导致结果偏差过大。

忽略非规格化数：非规格化数的处理速度慢，应在可能的情况下刷新为零（flush-to-zero）。

3.3 与其他层级的联系

与硬件的联系：不同GPU架构的浮点单元设计不同。Ampere的TF32格式在保持FP32易用性的同时提供接近FP16的性能。

与编译器的联系：编译器优化可能改变运算顺序，影响数值结果。使用-fp-model precise可以禁用可能影响精度的优化。

与算法设计的联系：数值稳定性应作为算法选择的标准之一。有时稍慢但更稳定的算法是更好的选择。

四、代码实战

4.1 浮点数精度演示

cpp

#include <cuda_runtime.h>
#include <stdio.h>
#include <cmath>
#include <cfloat>

// Demonstrate floating-point representation
void demonstrateRepresentation() {
    printf("=== Floating-Point Representation ===\n");
    
    // Show that 0.1 cannot be exactly represented
    float x = 0.1f;
    printf("0.1f stored as: %.20f\n", x);
    printf("Difference from 0.1: %.20e\n", x - 0.1);
    
    // Show machine epsilon
    printf("\nMachine Epsilon (FLT_EPSILON): %.20e\n", FLT_EPSILON);
    printf("1.0f + FLT_EPSILON: %.20f\n", 1.0f + FLT_EPSILON);
    printf("1.0f + FLT_EPSILON/2: %.20f (rounded to 1.0)\n", 1.0f + FLT_EPSILON/2);
    
    // Show subnormal numbers
    float subnormal = FLT_MIN / 2.0f;
    printf("\nSmallest normal float (FLT_MIN): %.20e\n", FLT_MIN);
    printf("Subnormal (FLT_MIN/2): %.20e\n", subnormal);
    
    // Show special values
    printf("\nPositive infinity: %f\n", INFINITY);
    printf("Negative infinity: %f\n", -INFINITY);
    printf("NaN: %f\n", NAN);
    printf("Is NaN == NaN? %s\n", (NAN == NAN) ? "true" : "false");
}

// Demonstrate catastrophic cancellation
void demonstrateCancellation() {
    printf("\n=== Catastrophic Cancellation ===\n");
    
    float a = 1.0000001f;
    float b = 1.0000000f;
    float c = a - b;
    
    printf("a = %.10f\n", a);
    printf("b = %.10f\n", b);
    printf("a - b = %.10e\n", c);
    printf("Expected: 1.0e-7\n");
    printf("Relative error: %.2f%%\n", 
           fabs(c - 1e-7) / 1e-7 * 100);
    
    // Better approach: use higher precision
    double da = 1.0000001;
    double db = 1.0000000;
    double dc = da - db;
    printf("\nUsing double:\n");
    printf("a - b = %.10e\n", dc);
}

// Demonstrate floating-point comparison
void demonstrateComparison() {
    printf("\n=== Floating-Point Comparison ===\n");
    
    float sum = 0.0f;
    for (int i = 0; i < 10; i++) {
        sum += 0.1f;
    }
    
    printf("Sum of ten 0.1s: %.20f\n", sum);
    printf("Expected: 1.00000000000000000000\n");
    printf("Direct equality (sum == 1.0f): %s\n", 
           (sum == 1.0f) ? "true" : "false");
    
    // Proper comparison
    bool approxEqual = fabsf(sum - 1.0f) < 1e-6f;
    printf("Approximate equality (|sum - 1| < 1e-6): %s\n",
           approxEqual ? "true" : "false");
}

int main() {
    demonstrateRepresentation();
    demonstrateCancellation();
    demonstrateComparison();
    return 0;
}

4.2 Kahan求和vs朴素求和

cpp

#include <cuda_runtime.h>
#include <stdio.h>
#include <cmath>

#define N 10000000

// Naive summation
__global__ void naiveSum(const float* data, float* result, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    float localSum = 0.0f;
    for (int i = idx; i < N; i += blockDim.x * gridDim.x) {
        localSum += data[i];
    }
    
    shared[tid] = localSum;
    __syncthreads();
    
    // Tree reduction
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    if (tid == 0) {
        atomicAdd(result, shared[0]);
    }
}

// Kahan summation per thread, then reduce
__global__ void kahanSum(const float* data, float* result, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Kahan summation for this thread's elements
    float sum = 0.0f;
    float comp = 0.0f;
    
    for (int i = idx; i < N; i += blockDim.x * gridDim.x) {
        float y = data[i] - comp;
        float t = sum + y;
        comp = (t - sum) - y;
        sum = t;
    }
    
    shared[tid] = sum;
    __syncthreads();
    
    // Tree reduction (without Kahan, for simplicity)
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    if (tid == 0) {
        atomicAdd(result, shared[0]);
    }
}

// CPU reference using double precision
double cpuReferenceSum(const float* data, int N) {
    double sum = 0.0;
    for (int i = 0; i < N; i++) {
        sum += data[i];
    }
    return sum;
}

void testSummationAccuracy() {
    printf("=== Summation Accuracy Test ===\n");
    
    float *h_data = (float*)malloc(N * sizeof(float));
    
    // Create challenging data: alternating large and small values
    for (int i = 0; i < N; i++) {
        if (i % 2 == 0) {
            h_data[i] = 1e6f;
        } else {
            h_data[i] = 1e-6f;
        }
    }
    
    // CPU reference
    double refSum = cpuReferenceSum(h_data, N);
    printf("Reference sum (double): %.10e\n", refSum);
    
    // Device allocation
    float *d_data, *d_result;
    cudaMalloc((void**)&d_data, N * sizeof(float));
    cudaMalloc((void**)&d_result, sizeof(float));
    
    cudaMemcpy(d_data, h_data, N * sizeof(float), cudaMemcpyHostToDevice);
    
    int blockSize = 256;
    int gridSize = 64;
    
    // Naive sum
    cudaMemset(d_result, 0, sizeof(float));
    naiveSum<<<gridSize, blockSize, blockSize * sizeof(float)>>>(
        d_data, d_result, N);
    cudaDeviceSynchronize();
    
    float naiveResult;
    cudaMemcpy(&naiveResult, d_result, sizeof(float), cudaMemcpyDeviceToHost);
    printf("Naive GPU sum:      %.10e\n", naiveResult);
    printf("Naive error:        %.10e\n", fabs(naiveResult - refSum));
    
    // Kahan sum
    cudaMemset(d_result, 0, sizeof(float));
    kahanSum<<<gridSize, blockSize, blockSize * sizeof(float)>>>(
        d_data, d_result, N);
    cudaDeviceSynchronize();
    
    float kahanResult;
    cudaMemcpy(&kahanResult, d_result, sizeof(float), cudaMemcpyDeviceToHost);
    printf("Kahan GPU sum:      %.10e\n", kahanResult);
    printf("Kahan error:        %.10e\n", fabs(kahanResult - refSum));
    
    printf("Kahan improvement:  %.2fx\n", 
           fabs(naiveResult - refSum) / fabs(kahanResult - refSum));
    
    cudaFree(d_data);
    cudaFree(d_result);
    free(h_data);
}

int main() {
    testSummationAccuracy();
    return 0;
}

4.3 混合精度矩阵乘法

cpp

#include <cuda_runtime.h>
#include <cuda_fp16.h>
#include <stdio.h>

#define TILE_SIZE 16

// FP32 matrix multiplication
__global__ void matMulFP32(const float* A, const float* B, float* C,
                            int M, int N, int K) {
    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    
    if (row < M && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < K; k++) {
            sum += A[row * K + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

// Mixed precision: FP16 computation, FP32 accumulation
__global__ void matMulMixed(const __half* A, const __half* B, float* C,
                             int M, int N, int K) {
    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    
    if (row < M && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < K; k++) {
            // Convert to float, multiply, accumulate in float
            float a = __half2float(A[row * K + k]);
            float b = __half2float(B[k * N + col]);
            sum += a * b;
        }
        C[row * N + col] = sum;
    }
}

// Using half2 for vectorized operations
__global__ void matMulHalf2(const __half2* A, const __half2* B, float* C,
                             int M, int N, int K2) {  // K2 = K/2
    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    
    if (row < M && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < K2; k++) {
            __half2 a = A[row * K2 + k];
            __half2 b = B[k * N + col];
            
            // Multiply two half pairs and accumulate
            __half2 prod = __hmul2(a, b);
            sum += __half2float(prod.x) + __half2float(prod.y);
        }
        C[row * N + col] = sum;
    }
}

void benchmarkMixedPrecision(int M, int N, int K) {
    printf("=== Mixed Precision Matrix Multiplication (%dx%dx%d) ===\n", 
           M, N, K);
    
    size_t sizeA = M * K * sizeof(float);
    size_t sizeB = K * N * sizeof(float);
    size_t sizeC = M * N * sizeof(float);
    
    float *h_A = (float*)malloc(sizeA);
    float *h_B = (float*)malloc(sizeB);
    float *h_C_fp32 = (float*)malloc(sizeC);
    float *h_C_mixed = (float*)malloc(sizeC);
    
    // Initialize
    for (int i = 0; i < M * K; i++) h_A[i] = (float)rand() / RAND_MAX;
    for (int i = 0; i < K * N; i++) h_B[i] = (float)rand() / RAND_MAX;
    
    // FP32 on GPU
    float *d_A_fp32, *d_B_fp32, *d_C_fp32;
    cudaMalloc((void**)&d_A_fp32, sizeA);
    cudaMalloc((void**)&d_B_fp32, sizeB);
    cudaMalloc((void**)&d_C_fp32, sizeC);
    
    cudaMemcpy(d_A_fp32, h_A, sizeA, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B_fp32, h_B, sizeB, cudaMemcpyHostToDevice);
    
    dim3 blockSize(TILE_SIZE, TILE_SIZE);
    dim3 gridSize((N + TILE_SIZE - 1) / TILE_SIZE, 
                  (M + TILE_SIZE - 1) / TILE_SIZE);
    
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // FP32 benchmark
    cudaEventRecord(start);
    for (int i = 0; i < 20; i++) {
        matMulFP32<<<gridSize, blockSize>>>(d_A_fp32, d_B_fp32, d_C_fp32, 
                                             M, N, K);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float msFP32 = 0;
    cudaEventElapsedTime(&msFP32, start, stop);
    
    cudaMemcpy(h_C_fp32, d_C_fp32, sizeC, cudaMemcpyDeviceToHost);
    
    // Mixed precision
    __half *d_A_half, *d_B_half;
    cudaMalloc((void**)&d_A_half, M * K * sizeof(__half));
    cudaMalloc((void**)&d_B_half, K * N * sizeof(__half));
    
    // Convert to half
    __half *h_A_half = (__half*)malloc(M * K * sizeof(__half));
    __half *h_B_half = (__half*)malloc(K * N * sizeof(__half));
    
    for (int i = 0; i < M * K; i++) {
        h_A_half[i] = __float2half(h_A[i]);
    }
    for (int i = 0; i < K * N; i++) {
        h_B_half[i] = __float2half(h_B[i]);
    }
    
    cudaMemcpy(d_A_half, h_A_half, M * K * sizeof(__half), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B_half, h_B_half, K * N * sizeof(__half), cudaMemcpyHostToDevice);
    
    float *d_C_mixed;
    cudaMalloc((void**)&d_C_mixed, sizeC);
    
    // Mixed precision benchmark
    cudaEventRecord(start);
    for (int i = 0; i < 20; i++) {
        matMulMixed<<<gridSize, blockSize>>>(d_A_half, d_B_half, d_C_mixed,
                                              M, N, K);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float msMixed = 0;
    cudaEventElapsedTime(&msMixed, start, stop);
    
    cudaMemcpy(h_C_mixed, d_C_mixed, sizeC, cudaMemcpyDeviceToHost);
    
    printf("FP32 time:    %.3f ms\n", msFP32 / 20);
    printf("Mixed time:   %.3f ms\n", msMixed / 20);
    
    // Check accuracy
    float maxError = 0.0f;
    for (int i = 0; i < M * N; i++) {
        float error = fabs(h_C_mixed[i] - h_C_fp32[i]);
        if (error > maxError) maxError = error;
    }
    printf("Max error:    %.2e\n", maxError);
    
    // Cleanup
    cudaFree(d_A_fp32);
    cudaFree(d_B_fp32);
    cudaFree(d_C_fp32);
    cudaFree(d_A_half);
    cudaFree(d_B_half);
    cudaFree(d_C_mixed);
    free(h_A);
    free(h_B);
    free(h_C_fp32);
    free(h_C_mixed);
    free(h_A_half);
    free(h_B_half);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}

int main() {
    benchmarkMixedPrecision(512, 512, 512);
    return 0;
}

五、深入思考

5.1 历史背景

IEEE 754标准于1985年首次发布，统一了之前各厂商不一致的浮点表示。在此之前，不同计算机可能使用不同的浮点格式，使得数值软件的可移植性极差。IEEE 754的成功在于它在硬件实现的复杂性和数值准确性之间取得了平衡。

GPU早期的浮点支持非常有限，甚至不支持IEEE 754合规的单精度运算。NVIDIA Tesla架构（2006）首次实现了完整的IEEE 754单精度支持。Fermi架构（2010）添加了双精度支持，但性能受限。Volta（2017）引入Tensor Core，开创了混合精度计算的新时代。

近年来，新的浮点格式不断涌现：BF16（Brain Float 16）由Google提出，保留了FP32的指数范围但减少了尾数位；FP8正在成为深度学习推理的新标准；TF32是NVIDIA的专有格式，旨在简化混合精度编程。

5.2 前沿方向

Posit算术：Posit是一种替代浮点表示，由John Gustafson提出。Posit使用可变长度的指数字段，可以在相同的位数下提供更高的动态范围和精度。Posit在机器学习应用中显示出潜力，但硬件支持仍在早期阶段。

随机舍入：传统的确定性舍入可能引入系统性偏差。随机舍入以概率方式舍入，统计上无偏。这对于长期运行的模拟和深度学习训练有益。

任意精度算术：对于需要极高精度的应用（如密码学、高精度物理模拟），GPU上的任意精度算术库正在发展中。这些库使用多个浮点数表示一个高精度数。

量子数值计算：量子计算机的数值表示与经典计算机完全不同。量子比特可以处于叠加态，数值计算范式需要重新思考。

5.3 延伸阅读

IEEE 754-2019 Standard for Floating-Point Arithmetic

David Goldberg, "What Every Computer Scientist Should Know About Floating-Point Arithmetic", 1991

NVIDIA CUDA C Programming Guide - Floating-Point Numbers

"Accuracy and Stability of Numerical Algorithms", Nicholas J. Higham

六、实验笔记

实验1：Kahan求和精度提升量化

目的：测量Kahan求和相比朴素求和的精度提升。

步骤：

生成具有挑战性的数据集（交替大小数、大范围值）

分别用朴素求和和Kahan求和计算总和

与双精度参考值比较，计算相对误差

改变数据分布，观察误差变化

预期：对于均匀分布的数据，两者差异不大；对于大范围数据，Kahan求和可将误差降低1-2个数量级。

实验2：混合精度训练的收敛性测试

目的：验证混合精度是否影响模型收敛。

步骤：

实现一个简单的神经网络（如MNIST分类器）

分别用纯FP32和混合精度训练

比较训练曲线和最终准确率

预期：混合精度训练应达到与FP32相近的准确率，但训练时间缩短。可能需要调整学习率和loss scaling。

七、本章小结

本章系统介绍了GPU编程中的数值考虑，从IEEE 754浮点数标准的基础知识出发，深入分析了舍入误差的来源、传播和影响。我们理解了为什么浮点数不能精确表示许多十进制小数，以及这种表示误差如何在运算中累积和放大。

数值稳定性是高质量数值软件的核心要求。通过方差计算的例子，我们看到了看似等价的算法在数值行为上可能有天壤之别。Welford算法和Kahan求和等技术通过巧妙的数学变换，显著减少了舍入误差的影响。

在GPU上，不同精度的浮点运算有着巨大的性能差异。单精度通常是默认选择，但在科学计算中可能需要双精度。半精度和混合精度计算在深度学习中日益普及，可以在保持模型精度的同时大幅提升吞吐量。

CUDA提供了丰富的工具来处理数值问题：快速数学函数用于性能敏感的场景，特殊值检测函数用于鲁棒性，half类型用于混合精度。程序员需要根据应用的精度要求和性能目标，做出明智的选择。

通过三个代码示例，我们实践了浮点数精度演示、Kahan求和对比和混合精度矩阵乘法，将理论知识转化为实际的编程技能。

八、习题

浮点表示：将十进制数0.25、0.1、-3.5转换为IEEE 754单精度二进制表示。说明转换过程。

机器Epsilon：编写一个CUDA kernel计算机器epsilon。提示：从1.0开始，不断除以2，直到1.0 + x == 1.0。

抵消误差：给定a = 1.00001和b = 1.00000，计算a - b的相对误差。如果使用双精度计算，相对误差会如何变化？

稳定性分析：以下计算方差的公式哪个更稳定？为什么？

- Var(X) = E[X²] - (E[X])²

- Var(X) = E[(X - μ)²]

Kahan求和实现：修改第4章的归约kernel，使每个线程使用Kahan求和。测量精度提升。

混合精度转换：编写一个函数，将FP32数组转换为FP16，再转回FP32，计算最大绝对误差。误差主要来自什么？

快速函数精度：比较sinf(x)和__sinf(x)在x = 0.1, 1.0, 10.0处的结果差异。最大相对误差是多少？

非规格化数性能：编写一个测试，比较正常数和非常小的非规格化数的运算速度。性能差异有多大？

条件数分析：矩阵的条件数定义为||A||·||A⁻¹||。高条件数矩阵的病态程度如何影响数值计算的稳定性？举例说明。

综合设计：设计一个数值稳定的并行标准差计算kernel。说明如何处理均值计算、平方差累加，以及如何减少舍入误差。

第6章 数值考虑 - 浮点数的陷阱 (Numerical Considerations - Pitfalls of Floating-Point Arithmetic)