快速求(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 >
   
2. 
   
3. #include < HugeCalc.h > // 公共接口
   
4. #include < HugeInt.h >  // 10进制系统
   
5. #include < HugeIntX.h > // 16进制系统
   
6. 
   
7. #pragma message( "automatic link to .../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
   
8. #pragma comment( lib, "HugeCalc.lib" )
   
9. 
   
10. int main(int argc, char* argv[])
    
11. {
    
12.     cout << "Call " << HugeCalc::GetVer() << endl;
    
13. 
    
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;
    
26. 
    
27.         cout << "*** [sqrt(2)-1]^r = sqrt(k) - sqrt(k-1) ***"  << endl;
    
28.         cout << "input r, output k ...(if r==0, then exit)" << endl;
    
29.         
    
30.         for ( ; ; )
    
31.         {
    
32.             cout << endl << "r = ";
    
33.             cin >> r;
    
34. 
    
35.             if ( 0 == r ) break;
    
36. 
    
37.             HugeCalc::EnableTimer();
    
38.             HugeCalc::ResetTimer();
    
39. 
    
40.             // a + b*sqrt(2)
    
41.             a = 1;
    
42.             b = -1;
    
43. 
    
44.             i = 1;
    
45.             while ( i <= r )
    
46.             {
    
47.                 i <<= 1;
    
48.             }
    
49.             i >>= 1;
    
50. 
    
51.             while ( 1 != i )
    
52.             {
    
53.                 k.Mul( a, b ) <<= 1;
    
54.                 a.Pow(2) += ( b.Pow(2) <<= 1 );
    
55.                 b.Swap(k);
    
56. 
    
57.                 if ( 0 != ( r & ( i >>= 1 )))
    
58.                 {
    
59.                     ( k = b ).Negate();
    
60.                     b -= a;
    
61.                     a.Swap( k ) -= b;
    
62.                 }               
    
63.             }
    
64. 
    
65.             if ( 0 == ( 1 & r ))
    
66.             {
    
67.                 k.Swap(a).Pow(2);
    
68.             }
    
69.             else
    
70.             {
    
71.                 k.Swap(b).Pow(2) <<= 1;
    
72.             }
    
73. 
    
74.             HugeCalc::EnableTimer( FALSE );
    
75. 
    
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.     }
    
84. 
    
85.     cout << endl;
    
86.     system( "pause" );
    
87. 
    
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. //////////////////////////////////////////////////////////////////////////
    
12. 
    
13. // Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
    
14. //      Win32 Debug:    Debug Multithreaded DLL
    
15. //      Win32 Release:  Multithreaded DLL
    
16. 
    
17. #include < iostream.h >
    
18. #include < fstream.h >
    
19. 
    
20. #include < HugeCalc.h > // 公共接口
    
21. #include < HugeInt.h >  // 10进制系统
    
22. #include < HugeIntX.h > // 16进制系统
    
23. 
    
24. #pragma message( "automatic link to .../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
    
25. #pragma comment( lib, "HugeCalc.lib" )
    
26. 
    
27. int main(int argc, char* argv[])
    
28. {
    
29.     cout << "Call " << HugeCalc::GetVer() << endl << endl;
    
30. 
    
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 可正可负
    
47. 
    
48.         BOOL bSave2File = ( argc > 1 );     // 后面若有参数则存盘
    
49. 
    
50.         cout << "*** [sqrt(2)-1]^r = sqrt(k) - sqrt(k-1) ***" << endl;
    
51.         cout << "input r, output k ...(if r==0, then exit)" << endl;
    
52. 
    
53.         for ( ; ; )
    
54.         {
    
55.             cout << endl << "r = ";
    
56.             cin >> r;
    
57. 
    
58.             if ( 0 == r )
    
59.             {
    
60.                 cout << "k = 1" << endl;
    
61.                 break;
    
62.             }
    
63. 
    
64.             if ( r < 0 ) r = -r;
    
65. 
    
66.             HugeCalc::EnableTimer();
    
67.             HugeCalc::ResetTimer();
    
68. 
    
69.             // a + b*sqrt(2)
    
70.             a = 1;
    
71.             b = 1;
    
72.             k = 2;
    
73. 
    
74.             UINT32 i = 1UL << 31;
    
75.             while ( i > (UINT32)r )
    
76.             {
    
77.                 i >>= 1;
    
78.             }
    
79. 
    
80.             while ( 0 != ( r & ( i - 1 )))
    
81.             {
    
82.                 ( b *= a ) <<= 1;
    
83.                 --( a.Swap( k ) <<= 1 );
    
84. 
    
85.                 if ( 0 != ( r & ( i >>= 1 )))
    
86.                 {
    
87.                     k = b;
    
88.                     b += a;
    
89.                     a.Swap( k ) += b;
    
90. 
    
91.                     ( k = b ).Pow( 2 ) <<= 1;
    
92.                 }
    
93.                 else
    
94.                 {
    
95.                     ( k = a ).Pow( 2 );
    
96.                 }
    
97.             }
    
98. 
    
99.             // 变量 a,b 不再使用，清空其所占资源
    
100.             a = b = 0;
     
101. 
     
102.             while ( 1 != i )
     
103.             {
     
104.                 ( --( k <<= 1 )).Pow( 2 );
     
105.                 i >>= 1;
     
106.             }
     
107. 
     
108.             HugeCalc::EnableTimer( FALSE );
     
109. 
     
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.             }
     
128. 
     
129.             if ( !bSave2File )
     
130.             {
     
131.                 cout << "k = " << k.GetStr( FS_NORMAL ) << endl;
     
132.             }
     
133. 
     
134.             cout << "*** k have " << k.GetDigits() << " digits! ***" << endl;
     
135.             cout << "computation took " << HugeCalc::GetTimerStr() << endl;
     
136.         }
     
137.     }
     
138. 
     
139.     cout << endl;
     
140.     system( "pause" );
     
141. 
     
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
   
8. 
   
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
    
16. 
    
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)还更加快一些。