N个数平方和的平方根--高精度算法

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;
}

执行结果如下图。

当然，这个程序计算执行时间的方式有问题，没有考虑整型运算和浮点运算的区别，但基本也能看出，高精度算法在时间上没有任何优势。感兴趣的朋友可以全部替换成浮点数验证一下。

至于算法的严密证明，还请大佬们指点。