Kahan

bxz

2023-04-23

Floating-point Calculation

Kahan前复习一些浮点数的运算

A=Ma*2^Ea

B=Mb*2^Eb

假设(Ea>=Eb)

A±B=(Ma±Mb*2^Eb-Ea)*2^Ea

A*B=(Ma*Mb)*2^Ea+Eb

A/B=(Ma/Mb)*2^Ea-Eb

加减法先求阶差，小的阶码往大的转换，所以小的阶码的有效数需要右移两阶码差的绝对值，如果右移的过程中出现移出的情况，会把移出的位保存到临时的附加位上，有效数相加，结果规格化

运算过程中可能出现以下几种情况

除以0
阶码上溢，单精度> 127 ±inf
阶码下溢出，单精度< -126 ± 0
有效数溢出例如11.01 需要右归，需要把可能移出的位保存到临时附加位，这里最多右移一位，因为两有效数之和不可能达到4，所以有效数整数部分最多2位(11.)，保留一个1(1.)，最多右移一位
非规格化有效数，如0.11 需要左归
有效数全零，需要把阶码也置全零

因为有了附加位精度会得到提升，但是也需要依靠附加位做舍入

IEEE 754给出了四种rounding方式

就近舍入(默认)
- 舍入到最接近，在一样接近的情况下偶数优先(二进制里以0结尾的)
向+inf舍入
向-inf舍入
向0舍入

Kahan

先看段代码

#define TIMES 40000000
float sum = 0;
float test;
test= 1234567890;
printf("%f\n",test);
for (size_t i = 0; i < TIMES; i++){
    sum += 0.1;
}
grxer@Ubuntu22-wsl ~/s/test> gcc -O2 -o test test.c
grxer@Ubuntu22-wsl ~/s/test> ./test
1234567936.000000
2097152.000000

结果并没有我们想的那样，误差非常大

造成这种现象的原因主要有

0.1在内的好多数在我们浮点数里的无法精确表示 0 01111011 10011001100110011001101 有效数1001无效循环，类似与十进制的三分之一，相加的时候由于需要小的阶码往大的转换，小的有效数需要右移(我们的附加位有限)，两数相差小时，精度损失，当两数相差过大时，甚至存在大数吃小数的情况(右移成0.0的情况)，同时由于我们单精度浮点数有效数是23位，外加一个隐藏的1 可以表示十进制有效位数log(2²⁴)约等于7位

误差补偿算法Kahan

for (size_t i = 0; i < TIMES; i++){
    y = 0.1 - c;
    t = sum2 + y;
    c = (t - sum2) - y;
    sum2 = t;
}
sum2-=c；

主要思想就是计算出每次累计的舍入误差，添加到下一次计算上

c是对丢失的低位进行运算补偿的变量

当我们的sum2和y相差过大时，t会得到一个已经产生误差的值，我们用这个值去减原来的sum2(他们的两数相差不是很大，结果比较准确)然后减去这次的加数(同样两数相差不大，结果准确)就可以计算出这次计算的误差值，在下一次加法时去用这个误差修正结果

#include <stdio.h>
#define TIMES 400000000
int main() {
    float sum = 0;
    float sum2 = 0, y = 0, t = 0, c = 0;
    for (size_t i = 0; i < TIMES; i++) {
        sum += 0.1;
    }
    for (size_t i = 0; i < TIMES; i++) {
        y = 0.1 - c;
        t = sum2 + y;
        c = (t - sum2) - y;
        sum2 = t;
    }
    sum2 -= c;
    printf("%f\n", sum);
    printf("%f", sum2);
}

grxer@Ubuntu22-wsl ~/s/test> gcc -O2 -o test test.c
grxer@Ubuntu22-wsl ~/s/test> ./test
2097152.000000
40000000.000000⏎

同时也可以看到这次我在TIMES比第一次后面多加了一个0，第一种累计法结果还是2097152.000000,也印证了相差过大时大数字会把小数字完全吃掉，本例中0.1阶码01111011=123-127=-4，2097152阶码10010100=148-127=21

0.1有效数需要右移25位，但是他的有效数算上隐藏1最多24位，推测一下应该有一位的附加位作保护位

9.21日补档GRS

IEE754标准保证浮点数舍入误差在0.5upl内

ulp(units in the last place):最后一位上的单位值或称最小精度单位

为了达到0.5upl误差，在upl后多加了保护位G（guard bit),舍入位R（round bit）和黏着位S（sticky bit）来提升精度

舍入位右侧有非零是粘合位才置1

https://zhuanlan.zhihu.com/p/356960443 ：舍入与异常