数学研发网设为首页收藏本站

数学研发论坛

 找回密码
 欢迎注册
楼主: northwolves

[擂台] 快速求(sqrt(2)-1)的幂展开形式

[复制链接]
发表于 2008-1-24 15:22:16 | 显示全部楼层

我的可以计算更大的 :)

我用类似快速乘幂的算法,消除了递归,效率大幅提升:
  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. #include < HugeCalc.h > // 公共接口
  3. #include < HugeInt.h >  // 10进制系统
  4. #include < HugeIntX.h > // 16进制系统

  5. #pragma message( "automatic link to .../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
  6. #pragma comment( lib, "HugeCalc.lib" )

  7. int main(int argc, char* argv[])
  8. {
  9.     cout << "Call " << HugeCalc::GetVer() << endl;

  10.     if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
  11.     {
  12.         cout << "\r\n警告:您未通过 HugeCalc.dll 的许可认证!\r\n" \
  13.             << "\r\n解决方案可选下列方案之一:" \
  14.             << "\r\n    一、请将本程序移动到“/CopyrightByGuoXianqiang/[../]”目录下运行;" \
  15.             << "\r\n    二、或请在 HugeCalc.chm 中进行注册(一劳永逸)。" \
  16.             << endl;
  17.     }
  18.     else
  19.     {
  20.         CHugeInt a, b, k;
  21.         UINT32 i, r;

  22.         cout << "*** [sqrt(2)-1]^r = sqrt(k) - sqrt(k-1) ***"  << endl;
  23.         cout << "input r, output k ...(if r==0, then exit)" << endl;
  24.         
  25.         for ( ; ; )
  26.         {
  27.             cout << endl << "r = ";
  28.             cin >> r;

  29.             if ( 0 == r ) break;

  30.             HugeCalc::EnableTimer();
  31.             HugeCalc::ResetTimer();

  32.             // a + b*sqrt(2)
  33.             a = 1;
  34.             b = -1;

  35.             i = 1;
  36.             while ( i <= r )
  37.             {
  38.                 i <<= 1;
  39.             }
  40.             i >>= 1;

  41.             while ( 1 != i )
  42.             {
  43.                 k.Mul( a, b ) <<= 1;
  44.                 a.Pow(2) += ( b.Pow(2) <<= 1 );
  45.                 b.Swap(k);

  46.                 if ( 0 != ( r & ( i >>= 1 )))
  47.                 {
  48.                     ( k = b ).Negate();
  49.                     b -= a;
  50.                     a.Swap( k ) -= b;
  51.                 }               
  52.             }

  53.             if ( 0 == ( 1 & r ))
  54.             {
  55.                 k.Swap(a).Pow(2);
  56.             }
  57.             else
  58.             {
  59.                 k.Swap(b).Pow(2) <<= 1;
  60.             }

  61.             HugeCalc::EnableTimer( FALSE );

  62. #if 1   // 设“0”可关闭输出具体数值
  63.             cout << "k = " << k.GetStr(FS_NORMAL) << endl;
  64. #else
  65.             cout << "k have " << k.GetDigits() << " digits!" << endl;
  66. #endif
  67.             cout << "computation took " << HugeCalc::GetTimerStr() << endl;
  68.         }
  69.     }

  70.     cout << endl;
  71.     system( "pause" );

  72.     return 0;
  73. }
复制代码
编译好的程序压缩包:
推荐
(sqrt(2)-1)^r_fast.zip (306.51 KB, 下载次数: 13)

评分

参与人数 1威望 +1 贡献 +1 鲜花 +1 收起 理由
mathe + 1 + 1 + 1 呵呵,程序很快,也给郭大管理员献献花

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-1-24 15:49:11 | 显示全部楼层
呵呵,程序运行结果输出太花时间了,
我在我的笔记本上运行
1,000,000计算时间是1.322009 s,不过输出好像花了更长的时间。
10,000,000的情况竟然也能够算下去,计算时间是16.754020s.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-1-24 15:57:21 | 显示全部楼层
原帖由 mathe 于 2008-1-24 15:49 发表
呵呵,程序运行结果输出太花时间了,
我在我的笔记本上运行
1,000,000计算时间是1.322009 s,不过输出好像花了更长的时间。
10,000,000的情况竟然也能够算下去,计算时间是16.754020s.


这不是我的程序吧?(在我印象中,你的笔记本速度比我的台式机速度还要好)
我的输出也会很快的。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-1-24 16:13:04 | 显示全部楼层
是你的程序。
如果不重定向,输出应该没有办法提高速度吧,比如r=10,000,000
显示的计算时间只有1.3秒左右,可是整个输出过程,我就看见窗口中白花花一片数字在向上闪去,闪了1分多钟,
你是不是运行时将输出重定向了?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-1-24 16:24:38 | 显示全部楼层
对照 14# 与 12# 的描述,感觉在12# 中的“r”值均少打了个“0”,
所以我才会有所疑惑。

我没有重定向输出,但由于不存在复杂的进制转换,所以 HugeCalc 输出会很高效,
如果不想输出具体数值,可把我在 11# 中源代码的“#if 1”改成“#if 0”即可。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-1-24 16:31:13 | 显示全部楼层
不是进制转换问题,只是纯粹屏幕输入输出问题。12#的确少了个0
通常我都在Linux上面写代码的,所以没有去编译你的代码,而是直接在笔记本上运行你编译后的程序。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-1-24 16:36:39 | 显示全部楼层
屏幕输出确实比较耗时。

我重定向输出到文件测试,所耗时间也只是秒级。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-1-24 20:37:24 | 显示全部楼层

“终极”版本 :)

下面对该问题进行理论推导并尝试推广。

为叙述方便,记 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) = ar - br*j  = sqrt( kr ) - sqrt( kr - 1 )
当 r<0 时,(j-1)^r = (j+1)^(-r) = ar + br*j = sqrt( kr ) + sqrt( kr - 1 )

楼主所要求的 k = max{ ar^2, 2*br^2 },更一般地:
当 r 为偶数时,kr = ar^2;当 r 为奇数时 kr = 2*br^2

从以上等式可看出,(j-1)^r、(j+1)^r、(j-1)^(-r)、(j+1)^(-r) 的 kr 值彼此相等
(因为它们的“有理部”、“无理部”绝对值彼此相等)

有了这个结论就好办了,甚至可为下一步程序优化作理论指导,
因为大数同号加法更简单,所以我们尽可能采用“有理部”与“无理部”同号之情形。

(1+j)*(a + b*j) = (a+2*b) + (a+b)*j    => br+1 = ar + br, ar+1 = br+1 + br
(a + b*j)^2 = (a^2 + 2*b^2) + 2*a*b*j => a2r = ar^2 + 2*br^2, b2r = 2*ar*br
注意到 ar^2 - 2*br^2 = ±1,所以 a2r 的计算可以少算一次平方。
更好的结论则为 a2r = 2*kr - 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. // Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
  13. //      Win32 Debug:    Debug Multithreaded DLL
  14. //      Win32 Release:  Multithreaded DLL

  15. #include < iostream.h >
  16. #include < fstream.h >

  17. #include < HugeCalc.h > // 公共接口
  18. #include < HugeInt.h >  // 10进制系统
  19. #include < HugeIntX.h > // 16进制系统

  20. #pragma message( "automatic link to .../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
  21. #pragma comment( lib, "HugeCalc.lib" )

  22. int main(int argc, char* argv[])
  23. {
  24.     cout << "Call " << HugeCalc::GetVer() << endl << endl;

  25.     if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
  26.     {
  27.         cout << "\r\n警告:您未通过 HugeCalc.dll 的许可认证!\r\n" \
  28.             << "\r\n解决方案可选下列方案之一:" \
  29.             << "\r\n    一、请将本程序移动到“/CopyrightByGuoXianqiang/[../]”目录下运行;" \
  30.             << "\r\n    二、或请在 HugeCalc.chm 中进行注册(一劳永逸)。" \
  31.             << endl;
  32.     }
  33.     else
  34.     {
  35. #if 1    // 用10进制内核
  36.         CHugeInt    a, b, k;
  37. #else    // 用16进制内核
  38.         CHugeIntX   a, b, k;
  39. #endif
  40.         SINT32 r;   // r 可正可负

  41.         BOOL bSave2File = ( argc > 1 );     // 后面若有参数则存盘

  42.         cout << "*** [sqrt(2)-1]^r = sqrt(k) - sqrt(k-1) ***" << endl;
  43.         cout << "input r, output k ...(if r==0, then exit)" << endl;

  44.         for ( ; ; )
  45.         {
  46.             cout << endl << "r = ";
  47.             cin >> r;

  48.             if ( 0 == r )
  49.             {
  50.                 cout << "k = 1" << endl;
  51.                 break;
  52.             }

  53.             if ( r < 0 ) r = -r;

  54.             HugeCalc::EnableTimer();
  55.             HugeCalc::ResetTimer();

  56.             // a + b*sqrt(2)
  57.             a = 1;
  58.             b = 1;
  59.             k = 2;

  60.             UINT32 i = 1UL << 31;
  61.             while ( i > (UINT32)r )
  62.             {
  63.                 i >>= 1;
  64.             }

  65.             while ( 0 != ( r & ( i - 1 )))
  66.             {
  67.                 ( b *= a ) <<= 1;
  68.                 --( a.Swap( k ) <<= 1 );

  69.                 if ( 0 != ( r & ( i >>= 1 )))
  70.                 {
  71.                     k = b;
  72.                     b += a;
  73.                     a.Swap( k ) += b;

  74.                     ( k = b ).Pow( 2 ) <<= 1;
  75.                 }
  76.                 else
  77.                 {
  78.                     ( k = a ).Pow( 2 );
  79.                 }
  80.             }

  81.             // 变量 a,b 不再使用,清空其所占资源
  82.             a = b = 0;

  83.             while ( 1 != i )
  84.             {
  85.                 ( --( k <<= 1 )).Pow( 2 );
  86.                 i >>= 1;
  87.             }

  88.             HugeCalc::EnableTimer( FALSE );

  89.             if ( bSave2File )
  90.             {
  91.                 ofstream outputFile( "Output.txt", ios::ate );
  92.                 if ( outputFile )
  93.                 {
  94.                     outputFile << endl;
  95.                     outputFile << "r = " << r << endl;
  96.                     outputFile << "k = " << k.GetStr( FS_NORMAL ) << endl;
  97.                     outputFile << "*** k have " << k.GetDigits() << " digits! ***" << endl;
  98.                     outputFile << "computation took " << HugeCalc::GetTimerStr() << endl;
  99.                     outputFile.close();
  100.                 }
  101.                 else
  102.                 {
  103.                     cerr << "File could not be opened" << endl;
  104.                     bSave2File = FALSE;
  105.                 }
  106.             }

  107.             if ( !bSave2File )
  108.             {
  109.                 cout << "k = " << k.GetStr( FS_NORMAL ) << endl;
  110.             }

  111.             cout << "*** k have " << k.GetDigits() << " digits! ***" << endl;
  112.             cout << "computation took " << HugeCalc::GetTimerStr() << endl;
  113.         }
  114.     }

  115.     cout << endl;
  116.     system( "pause" );

  117.     return 0;
  118. }
复制代码
编译好的程序压缩包:
r_k.zip (340.27 KB, 下载次数: 8)
无论注册 HugeCalc 与否均可正常运行;可存盘输出;内含源代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-1-25 07:58:42 | 显示全部楼层

“终极版”的测试报告

  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. 调用 HugeCalc 的 16 进制内核版本:
  9. r = 1000        耗时:0.000067 s    位数:636
  10. r = 10000       耗时:0.000379 s    位数:6358
  11. r = 100000      耗时:0.021730 s    位数:63578
  12. r = 1000000     耗时:0.227757 s    位数:635777
  13. r = 10000000    耗时:1.378083 s    位数:6357767
  14. r = 100000000   耗时:19.163445 s   位数:63577665

  15. 测试环境:Windows XP SP2,AMD Athlon 64 Processor 3200+,1GB DDR - 200MHz
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-1-25 08:18:38 | 显示全部楼层
对于特殊情况应该还有优化机会,不过机会不是很大了
现在对于(a+bj)做一次平方,需要对系数做一次平方,一次乘法,两次移位。(一次简单减法可以忽略)
而(a+bj)乘上1+j或(1-j)需要两次加减法运算。
不知道现在通常情况现在乘法需要时间是加法的多少倍(当然同数据大小相关),如果我们对于这个关系有一个比较的表达式,那么应该还有少数优化机会。
比如如果计算(j-1)^1023,有可能写成(j-1)^1024*(j+1)还更加快一些。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|Archiver|数学研发网 ( 苏ICP备07505100号  

GMT+8, 2017-5-1 10:38 , Processed in 0.210810 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表