[日常刷题] C语言：在 平方和 算例中用浮点型还是整形？

原文（自己的Draft）标题：一个double引发的long （漫长）思考

• 前言：先说说刷题。我的好朋友去了某知名高薪大企业，跟我闲聊第一句话就是“你日常刷题么？”他说，这个大企，除了工作任务，还要日常刷题。我说：“我要去这个大企，是不是现在就得刷起来？”，他说，去了再刷也不迟。我内心：虽然从不刷题，但我日常随便写写代码也算是刷题吧，练的是…手速，也是程序员引以为豪常拿出来攀比的一个技能呢，哈哈。希望读者也能日常动手。

问题

这个问题来自于C语言圣经《C Primer Plus》(第六版) 例8-7：

“使用输入检测为一个进行算术运算的函数提供整数，该函数计算特定范围内所有整数的平方和。程序限制了范围的上限是10000000，下限是-10000000”

我的实现

我的源代码可以参考github。
说明：SOS是平方和Sum Of Square的缩写。

答案

其中有一段：

static long get_sum_square(long begin, long end)
{
	long i = 0;
	long sum_square = 0;
	for (i=begin ; i<=end ; i++) {
		sum_square += i*i;
	}

	return sum_square;
}

和书上的例子不一样：

double sum_squares(long a, long b)
{
	double total = 0;
	long i;

	for (i=a ; i<=b ; i++)
		total += (double)i * (double)i;

	return total;
}

探究

看上面两段代码，其本质区别是返回的值的类型，我用long形，答案用double型。那么二者返回的结果会一样么？都是对的么？
测试代码：将两个函数都放到代码中，并同时输出进行比较：

int main(void)
{
	long begin = 0, end = 0;
	do {
		printf("Please input the range :\n");
		//Get the range [begin,end]
		printf("low limit:");
		begin = get_long();
		printf("high limit:");
		end = get_long();
		//Judge the range available.
	} while(is_limit_bad(begin,end));

	//Calculate the sum and output
	long sum_square = get_sum_square(begin, end);
	double sum_square_double = sum_squares(begin, end);
	printf("\r\nlong sos = %ld\r\ndouble sos = %f\r\n", sum_square, sum_square_double);
	printf("Same? %d\r\n",(sum_square == sum_square_double));

	getchar();
	getchar();

	return 0;
}

测试用例：

low	high	long sos	double sos	Same?	Target
1	2	5	5.000000	1	通过一个可以口算的测试用例来确认当前写的代码是可行且基本正确的
-1	-2			RangeError	通过一个不明显的错误说明Error handling起效的
0	10000000	-762584128	333, 333, 383, 333, 688, 000, 000.00	1	直接上半个范围测

根据结果，很明显long sos是错的，但double sos不知道是不是正确的。这是一个瓶颈性的问题，算大数的时候，由于口算无法算出来，没有办法判断正确性。于是寻求网络，还真被我找到了SOS的计算器

测试得真知：

low	high	long sos	double sos	Same?	Sos Calculator	Target	Result
1	100	338350	338350	1	338350	测试正常计算结果，以自然数1开始	结果正确，该计算器可用
20	25	3055	3055.000000	1	3055	测试SOS calculator是否和我的程序相匹配，下限为任意值，而不一定是自然数1	恰好满足需求
1	10000000	-762584128	333, 333, 383, 333, 688, 000, 000.00	0	333, 333, 383, 333, 335, 000, 000	大数stress 测试	这下情况不妙，无论是long的结果，还是double的结果，似乎都不正确

这时，我甚至怀疑SOS Calculator算出的结果。好在这个网页给了下面的公式：

哦呵呵，胜造七级浮屠啊，公式就是不论你数多大，都得按我这个道理。所以，公式算出来的就是正确答案啊，于是将程序中的main函数再进行修改：

long sum_square = get_sum_square(begin, end);
double sum_square_double = sum_squares(begin, end);
printf("\r\nlong sos = %ld\r\ndouble sos = %f\r\n", sum_square, sum_square_double);
printf("Same? %d\r\n",(sum_square == sum_square_double));
printf("Right sos = %ld", right_sum_square(begin, end));

其中 right_sum_square ()的实现如下：

static long right_sum_square(long begin, long end)
{
	return end * (end+1) * (2*end+1) / 6 - begin * (begin-1) * (2*begin-1) / 6;
}

但悲剧了，公式给我的结果竟是这：

丧心病狂的我决定带进去手算一遍，得到的结果是：33333338333333500000，和网上计算器的结果一样，看来还是网上计算器靠谱。
long 型出现了负数，八成是因为越界了。这样看，这可能是两个问题混在一起了，一个是越界，另一个是long的结果和double结果的不同。下面分别进行讨论。

long 越界

为了便于探究和区分，将函数改名get_sum_square 修改为 get_sum_square_long，将函数名sum_squares修改为get_sum_square_double，这样，如果认为long类型会越界，咱就搞个long long类型，放到main()函数中：

printf("sum_square_double = %f\r\n", sum_square_double);
printf("sum_square_long(%d) = %ld\r\n", sizeof(long), sum_square_long);
printf("sum_square_long_long(%d) = %I64d\r\n", sizeof(long long int), get_sum_square_long_long(begin, end));

其中：

static long long get_sum_square_long_long(long begin, long end)
{
	long long i=0;
	long long sum_square = 0;
	for (i=begin ; i<=end ; i++) {
		sum_square += i*i;
	}

	return sum_square;
}

我用Eclipse没能使用%lld输出64位的整型数，%I64d参考了 printf和scanf处理long long int型数据。

悲剧，看起来还是越界了，用电脑计算器直接将结果打进去，发现打64位后，就没有办法再往里输入了，少了个0，怎么也输入不进去：

看来是“真的”越界了，这也有拓展到一个信息，在64位的计算机里的Calculator限制了输入内容要小于64位最大可以表示的数（9223372036854775807），共19位十进制数字:

这样看来确实是double的结果相对更准确呢。
简单分析一下，整形的存储和浮点型不一样。整形值越大，这个数的二进制位数就越多，相当于往左移，当左移到第64位的时候，移不动了，剩下的就舍弃了，相当于100和1000的差别。而浮点型，是存储了一个小数与2的指数，这个小数是对原数不断除2 的结果，后面放不下的舍弃了，会把这后面舍弃的位数加到前面2的阶乘里，舍弃的也只是小数后面影响最小的一些数。参考：关于C++ double浮点数精度丢失的分析。

long 和 double 结果不同的探究

根据上面的测试，可以猜测，[low,high]范围越大就越有可能出现结果不同。那么我就找到第一个“出轨”的位置，以方便探究，在main中加入下面代码段：

    double sum_square_double = get_sum_square_double(begin, end);
	long long sum_square_long_long = get_sum_square_long_long(begin, end);
	long i = 0;
	for (i=begin ; i<=end ; i++) {
		sum_square_double = get_sum_square_double(begin, i);
		sum_square_long_long = get_sum_square_long_long(begin, i);
		if (sum_square_double != sum_square_long_long) {
			printf("i = %ld\r\n", i);
			printf("sum_square_double = %f\r\n", sum_square_double);
			printf("sum_square_long_long = %I64d\r\n", sum_square_long_long);
			break;
		}
	}

Come on，我就是个疯子，因为计算机在轮询计算没到满足条件的时候没有打印，让我没有办法确定是不是进入了死循环，所以我就把它们全打印出来了，也就是没有了if (sum_square_double != sum_square_long_long)这句判断，得到如下结果（另外，这里建议将等号对齐，我没有对齐，找起来有点慢）：

可以确定，问题出在第300081项上。根据公式有：

SOS(300081)
= 300081*(300081+1)(3000812+1) /6
= 9,007,336,992,830,441‬。

所以是浮点型的结果错了。
接下来，涉及到的是浮点型的表示和浮点数的加法，参考这篇文章：
浮点数的加减法运算。

看过文章后，动手！将两个加数都写成在计算机中的浮点数表示的形式：

300081*300081
= 90,048,606,561（Dec）
= 14 F750 B161‬（Hex）
=+ (2^36) * (1.4F750B161)

SOS(300080)
= 9007246944223880（Dec）
= 20 000B 1A83 C288‬（Hex）
=+(2^53) * (1.000058D41E1440)

按照上面给的参考内容有关浮点数的加法，第一步是对阶，之后相加，即：

2^53 * (1.000058D41E1440) + 2^53 * (0.0000A7BA858B01)
= 2^53 * (1.000058D41E1440 + 0.0000A7BA858B01)
= 2^53 * 1.0001008EA39F41
= 9,007,336,992,830,440‬ （计算出来只是为了说明咱们笔算的是没错的）

上面这个结果2^53 * 1.0001008EA39F41，在小数点后有56 bits，double最多能放52 bits，所以最后4位会被舍弃，也就是 2^53 * 1.0001008EA39F4 = 9,007,336,992,830,440‬。诺，这就是问题所在了。

结论

• 在浮点数加法中，如果有一个数的有效位数为52bit，再加一个新的数就很容易越界，而产生舍入误差。
• 在整形型和浮点型型所占用位数相同的时候，求平方和，或者大数加法时，推荐使用整形。
• 但是在结果的范围有可能造成整形越界的时候，就要使用浮点型，虽然结果的精度不够，但是不会产生巨大误差。
• 这有点类似于“粗调”和“细调”两个旋钮。

反思

孟子曰：“尽信书不如无书”。 书上还真没写错，但是探究一下发现学问蛮深的。所以孔子又说：“学而不思则罔”，只是把书上代码抄一遍，就真的只是练了练打字速度了（首尾呼应，有没有！）。