N个数平方和的平方根,公式如下:
我们都知道,计算机浮点数的表示,越靠近0点,所表示的浮点数精度越高,如下图。
那么,如果按照上图的公式直接计算,所得的结果精度会丢失得非常严重。
算法设计
为了避免丢失精度,我们就有了设计算法的核心思想:在计算平方和时,将每个加数,都转换为0到1之间的小数。
在计算出所有元素的平方和之后,再乘以放大倍数,这样就可以减少计算过程中的精度损失。用下图所示的公式表示。
图中A即为放大倍数。那么我们如何获取放大倍数A呢?
方法一:
在开始计算之前,先遍历所有元素,取最大绝对值,作为放大倍数。
但这样会多一次循环,参与计算的元素越多,这一次循环的代价越大。
方法二:
我们从已经参与过计算的元素中,取最大绝对值,作为放大倍数。
其中i为下标,初始值 i = 0。放大倍数A的初始值取0。
那么按照我们的核心思想开始设计,就有如下的设计过程:
我们先假设几个参与计算的元素:x0到x5,各个元素之间的大小关系如下图:
为了保证每个加数都是0到1之间的小数,MAX作为分子还是分母,要根据MAX和xi之间的大小来决定:
当MAX小于xi时,作为分子;否则作为分母。这样我们就有了如下的计算过程,下图中的MAX值用粗体表示。
于是我们得到了下图所示的计算结果:
显然,如果我们能在计算过程中,将 i = 1, 2时的分母,用新的MAX替换掉,并补充项,就可以得到我们想要的结果。
那么如何替换呢?下边的算法会告诉我们答案。
高精度算法的逻辑及公式如下图:
用表示所有已处理操作数中绝对值最大的数。
初始值 i = 0,MAX = 0,sum = 0。
其中,if语句下的--乘法和加1--操作,起到了替换分母的作用和补充项的作用。
验证代码
验证如下,编译指令 gcc squareSum.c -lm。
// squareSum.c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define N 100
#define CYCLE 10000000
int main(){
double res;
int i_sum = 0;
double d_sum = 0.0;
int max = 0;
double tmp;
time_t start, end;
int x[N];
int i = 0;
srand((unsigned)time(NULL));
for (; i < N; i++){
x[i] = rand() % 1000;
printf("x[%d]: %d\n", i+1, x[i]);
}
start = time(NULL);
for (int j = 0; j < CYCLE; j++){
for (i = 0; i < N; i++){
i_sum += x[i] * x[i];
}
res = sqrt(i_sum);
i_sum = 0;
}
end = time(NULL);
printf("squart root of i_sum: %f\n", res);
printf("difftime: %f\n", difftime(end, start));
start = time(NULL);
for (int j = 0; j < CYCLE; j++){
for (i = 0; i < N; i++){
if (max < x[i]){
tmp = (double)max / x[i];
d_sum = tmp * tmp * d_sum + 1;
max = x[i];
}
else{
tmp = (double)x[i] / max;
d_sum += tmp * tmp;
}
}
res = sqrt(d_sum) * max;
d_sum = 0;
}
end = time(NULL);
printf("squart root of d_sum: %f\n", res);
printf("difftime: %f\n", difftime(end, start));
return 0;
}
执行结果如下图。
当然,这个程序计算执行时间的方式有问题,没有考虑整型运算和浮点运算的区别,但基本也能看出,高精度算法在时间上没有任何优势。感兴趣的朋友可以全部替换成浮点数验证一下。
至于算法的严密证明,还请大佬们指点。