快速求(sqrt(2)-1)的幂展开形式

gxqcn

我用类似快速乘幂的算法，消除了递归，效率大幅提升：

1. r = 1000 耗时：0.000062 s 位数：765
2. r = 10000 耗时：0.001045 s 位数：7655
3. r = 100000 耗时：0.028692 s 位数：76555
4. r = 1000000 耗时：0.362074 s 位数：765551
5. r = 10000000 耗时：3.179194 s 位数：7655514
6. 测试环境：Windows XP SP2，AMD Athlon 64 Processor 3200+，1GB DDR - 200MHz

复制代码

代码如下：

1. #include < iostream.h >
3. #include < HugeCalc.h > // 公共接口
4. #include < HugeInt.h > // 10进制系统
5. #include < HugeIntX.h > // 16进制系统
7. #pragma message( "automatic link to .../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
8. #pragma comment( lib, "HugeCalc.lib" )
10. int main(int argc, char* argv[])
11. {
12. cout << "Call " << HugeCalc::GetVer() << endl;
14. if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
15. {
16. cout << "\r\n警告：您未通过 HugeCalc.dll 的许可认证！\r\n" \
17. << "\r\n解决方案可选下列方案之一：" \
18. << "\r\n 一、请将本程序移动到“/CopyrightByGuoXianqiang/[../]”目录下运行；" \
19. << "\r\n 二、或请在 HugeCalc.chm 中进行注册（一劳永逸）。" \
20. << endl;
21. }
22. else
23. {
24. CHugeInt a, b, k;
25. UINT32 i, r;
27. cout << "*** [sqrt(2)-1]^r = sqrt(k) - sqrt(k-1) ***" << endl;
28. cout << "input r, output k ...(if r==0, then exit)" << endl;
30. for ( ; ; )
31. {
32. cout << endl << "r = ";
33. cin >> r;
35. if ( 0 == r ) break;
37. HugeCalc::EnableTimer();
38. HugeCalc::ResetTimer();
40. // a + b*sqrt(2)
41. a = 1;
42. b = -1;
44. i = 1;
45. while ( i <= r )
46. {
47. i <<= 1;
48. }
49. i >>= 1;
51. while ( 1 != i )
52. {
53. k.Mul( a, b ) <<= 1;
54. a.Pow(2) += ( b.Pow(2) <<= 1 );
55. b.Swap(k);
57. if ( 0 != ( r & ( i >>= 1 )))
58. {
59. ( k = b ).Negate();
60. b -= a;
61. a.Swap( k ) -= b;
62. }
63. }
65. if ( 0 == ( 1 & r ))
66. {
67. k.Swap(a).Pow(2);
68. }
69. else
70. {
71. k.Swap(b).Pow(2) <<= 1;
72. }
74. HugeCalc::EnableTimer( FALSE );
76. #if 1 // 设“0”可关闭输出具体数值
77. cout << "k = " << k.GetStr(FS_NORMAL) << endl;
78. #else
79. cout << "k have " << k.GetDigits() << " digits!" << endl;
80. #endif
81. cout << "computation took " << HugeCalc::GetTimerStr() << endl;
82. }
83. }
85. cout << endl;
86. system( "pause" );
88. return 0;
89. }

复制代码

编译好的程序压缩包：

(sqrt(2)-1)^r_fast.zip (306.51 KB, 下载次数: 13)

2008-1-24 15:22 上传

点击文件名下载附件
快速计算(sqrt(2)-1)^r

（注：非 HugeCalc 注册用户也可正常运行）

mathe

呵呵，程序运行结果输出太花时间了， 我在我的笔记本上运行 1,000,000计算时间是1.322009 s，不过输出好像花了更长的时间。 10,000,000的情况竟然也能够算下去，计算时间是16.754020s.

gxqcn

原帖由 mathe 于 2008-1-24 15:49 发表 呵呵，程序运行结果输出太花时间了， 我在我的笔记本上运行 1,000,000计算时间是1.322009 s，不过输出好像花了更长的时间。 10,000,000的情况竟然也能够算下去，计算时间是16.754020s.

这不是我的程序吧？（在我印象中，你的笔记本速度比我的台式机速度还要好） 我的输出也会很快的。

mathe

是你的程序。 如果不重定向，输出应该没有办法提高速度吧，比如r=10,000,000 显示的计算时间只有1.3秒左右，可是整个输出过程，我就看见窗口中白花花一片数字在向上闪去，闪了1分多钟， 你是不是运行时将输出重定向了？

gxqcn

对照 14# 与 12# 的描述，感觉在12# 中的“r”值均少打了个“0”， 所以我才会有所疑惑。 我没有重定向输出，但由于不存在复杂的进制转换，所以 HugeCalc 输出会很高效， 如果不想输出具体数值，可把我在 11# 中源代码的“#if 1”改成“#if 0”即可。

mathe

不是进制转换问题，只是纯粹屏幕输入输出问题。12#的确少了个0 通常我都在Linux上面写代码的，所以没有去编译你的代码，而是直接在笔记本上运行你编译后的程序。

gxqcn

屏幕输出确实比较耗时。 我重定向输出到文件测试，所耗时间也只是秒级。

gxqcn

下面对该问题进行理论推导并尝试推广。 为叙述方便，记 sqrt(2)=j，

数 (a + b*j) 中 a 可称为“有理部”，b 可称为“无理部”。 则 r=0 时，(j-1)^r = (j+1)^(-r) = 1 ± 0*j 当 r>0 时，(j-1)^r = (j+1)^(-r) = a_r - b_r*j = sqrt( k_r ) - sqrt( k_r - 1 ) 当 r<0 时，(j-1)^r = (j+1)^(-r) = a_r + b_r*j = sqrt( k_r ) + sqrt( k_r - 1 ) 楼主所要求的 k = max{ a_r^2, 2*b_r^2 }，更一般地： 当 r 为偶数时，k_r = a_r^2；当 r 为奇数时 k_r = 2*b_r^2 从以上等式可看出，(j-1)^r、(j+1)^r、(j-1)^(-r)、(j+1)^(-r) 的 k_r 值彼此相等。 （因为它们的“有理部”、“无理部”绝对值彼此相等） 有了这个结论就好办了，甚至可为下一步程序优化作理论指导， 因为大数同号加法更简单，所以我们尽可能采用“有理部”与“无理部”同号之情形。 (1+j)*(a + b*j) = (a+2*b) + (a+b)*j => b_r+1 = a_r + b_r, a_r+1 = b_r+1 + b_r (a + b*j)^2 = (a^2 + 2*b^2) + 2*a*b*j => a_2r = a_r^2 + 2*b_r^2, b_2r = 2*a_r*b_r 注意到 a_r^2 - 2*b_r^2 = ±1，所以 a_2r 的计算可以少算一次平方。 更好的结论则为 a_2r = 2*k_r - 1 于是乎，就有了下面的优化算法：

1. //////////////////////////////////////////////////////////////////////////
2. //
3. // 目的：快速计算 (sqrt(2) - 1)^r
4. // 设计：郭先强 ( gxqcn@163.com; HugeCalc@Gmail.com )
5. // 日期：2008-01-24
6. // 备注：运行编译好的程序，后面若带参数则输出到文件 Output.txt 中
7. //
8. // Web: http://www.emath.ac.cn/
9. // BBS: http://bbs.emath.ac.cn/
10. //
11. //////////////////////////////////////////////////////////////////////////
13. // Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
14. // Win32 Debug: Debug Multithreaded DLL
15. // Win32 Release: Multithreaded DLL
17. #include < iostream.h >
18. #include < fstream.h >
20. #include < HugeCalc.h > // 公共接口
21. #include < HugeInt.h > // 10进制系统
22. #include < HugeIntX.h > // 16进制系统
24. #pragma message( "automatic link to .../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
25. #pragma comment( lib, "HugeCalc.lib" )
27. int main(int argc, char* argv[])
28. {
29. cout << "Call " << HugeCalc::GetVer() << endl << endl;
31. if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
32. {
33. cout << "\r\n警告：您未通过 HugeCalc.dll 的许可认证！\r\n" \
34. << "\r\n解决方案可选下列方案之一：" \
35. << "\r\n 一、请将本程序移动到“/CopyrightByGuoXianqiang/[../]”目录下运行；" \
36. << "\r\n 二、或请在 HugeCalc.chm 中进行注册（一劳永逸）。" \
37. << endl;
38. }
39. else
40. {
41. #if 1 // 用10进制内核
42. CHugeInt a, b, k;
43. #else // 用16进制内核
44. CHugeIntX a, b, k;
45. #endif
46. SINT32 r; // r 可正可负
48. BOOL bSave2File = ( argc > 1 ); // 后面若有参数则存盘
50. cout << "*** [sqrt(2)-1]^r = sqrt(k) - sqrt(k-1) ***" << endl;
51. cout << "input r, output k ...(if r==0, then exit)" << endl;
53. for ( ; ; )
54. {
55. cout << endl << "r = ";
56. cin >> r;
58. if ( 0 == r )
59. {
60. cout << "k = 1" << endl;
61. break;
62. }
64. if ( r < 0 ) r = -r;
66. HugeCalc::EnableTimer();
67. HugeCalc::ResetTimer();
69. // a + b*sqrt(2)
70. a = 1;
71. b = 1;
72. k = 2;
74. UINT32 i = 1UL << 31;
75. while ( i > (UINT32)r )
76. {
77. i >>= 1;
78. }
80. while ( 0 != ( r & ( i - 1 )))
81. {
82. ( b *= a ) <<= 1;
83. --( a.Swap( k ) <<= 1 );
85. if ( 0 != ( r & ( i >>= 1 )))
86. {
87. k = b;
88. b += a;
89. a.Swap( k ) += b;
91. ( k = b ).Pow( 2 ) <<= 1;
92. }
93. else
94. {
95. ( k = a ).Pow( 2 );
96. }
97. }
99. // 变量 a,b 不再使用，清空其所占资源
100. a = b = 0;
102. while ( 1 != i )
103. {
104. ( --( k <<= 1 )).Pow( 2 );
105. i >>= 1;
106. }
108. HugeCalc::EnableTimer( FALSE );
110. if ( bSave2File )
111. {
112. ofstream outputFile( "Output.txt", ios::ate );
113. if ( outputFile )
114. {
115. outputFile << endl;
116. outputFile << "r = " << r << endl;
117. outputFile << "k = " << k.GetStr( FS_NORMAL ) << endl;
118. outputFile << "*** k have " << k.GetDigits() << " digits! ***" << endl;
119. outputFile << "computation took " << HugeCalc::GetTimerStr() << endl;
120. outputFile.close();
121. }
122. else
123. {
124. cerr << "File could not be opened" << endl;
125. bSave2File = FALSE;
126. }
127. }
129. if ( !bSave2File )
130. {
131. cout << "k = " << k.GetStr( FS_NORMAL ) << endl;
132. }
134. cout << "*** k have " << k.GetDigits() << " digits! ***" << endl;
135. cout << "computation took " << HugeCalc::GetTimerStr() << endl;
136. }
137. }
139. cout << endl;
140. system( "pause" );
142. return 0;
143. }

复制代码

编译好的程序压缩包：

r_k.zip (340.27 KB, 下载次数: 8)

2008-1-25 06:25 上传

点击文件名下载附件
快速计算 (sqrt(2)-1)^r 终极版

（无论注册 HugeCalc 与否均可正常运行；可存盘输出；内含源代码）

gxqcn

1. 调用 HugeCalc 的 10 进制内核版本：
2. r = 1000 耗时：0.000048 s 位数：765
3. r = 10000 耗时：0.000632 s 位数：7655
4. r = 100000 耗时：0.012068 s 位数：76555
5. r = 1000000 耗时：0.236285 s 位数：765551
6. r = 10000000 耗时：1.586116 s 位数：7655514
7. r = 100000000 耗时：19.850545 s 位数：76555137
9. 调用 HugeCalc 的 16 进制内核版本：
10. r = 1000 耗时：0.000067 s 位数：636
11. r = 10000 耗时：0.000379 s 位数：6358
12. r = 100000 耗时：0.021730 s 位数：63578
13. r = 1000000 耗时：0.227757 s 位数：635777
14. r = 10000000 耗时：1.378083 s 位数：6357767
15. r = 100000000 耗时：19.163445 s 位数：63577665
17. 测试环境：Windows XP SP2，AMD Athlon 64 Processor 3200+，1GB DDR - 200MHz

复制代码

mathe

对于特殊情况应该还有优化机会，不过机会不是很大了 现在对于(a+bj)做一次平方，需要对系数做一次平方，一次乘法，两次移位。（一次简单减法可以忽略） 而(a+bj)乘上1+j或(1-j)需要两次加减法运算。 不知道现在通常情况现在乘法需要时间是加法的多少倍（当然同数据大小相关），如果我们对于这个关系有一个比较的表达式，那么应该还有少数优化机会。 比如如果计算(j-1)^1023，有可能写成(j-1)^1024*(j+1)还更加快一些。