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

[分享] 重大消息-一种新的任意精度对数算法研制成功

[复制链接]
发表于 2018-1-5 11:07:27 | 显示全部楼层
一发不可收拾,我用下面的代码能够说明时间统计的稳定性以及误差的传递,大家可以看看
  1. #include <iostream>
  2. #include <chrono>

  3. void hello()
  4. {
  5.     for( int i = 0,j=0; i < 1e8; ++i ) j++;
  6. }

  7. int main()
  8. {
  9.     for(int k=0;k<10;++k){
  10.         auto t1 = std::chrono::high_resolution_clock::now();
  11.         hello();
  12.         auto t2 = std::chrono::high_resolution_clock::now();

  13.         auto duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();

  14.         std::cout << duration/1000. <<"ms"<< std::endl;

  15.     }
  16.     return 0;
  17. }

复制代码

运行结果是
  1. 222.795ms
  2. 223.371ms
  3. 218.557ms
  4. 221.014ms
  5. 217.021ms
  6. 221.605ms
  7. 221.48ms
  8. 219.964ms
  9. 221.761ms
  10. 220.376ms
复制代码

可见,波动范围是2ms, 统计时间的最后一个不可靠度的位置是1ms, 加大被测试的函数里的循环次数1e8变成10e9, 误差发生了传递,或者说被放大了,最后一个不可靠的位置是10ms,所以,真正严谨的benchmarking还得统计一下耗时的稳定性,即你的统计数据的有效位是哪一位.

  1. 2236.23ms
  2. 2242.03ms
  3. 2209.1ms
  4. 2220.01ms
  5. 2210.53ms
  6. 2212.43ms
  7. 2205.2ms
  8. 2214.93ms
  9. 2221.3ms
  10. 2239.92ms
复制代码

点评

这很重要的  发表于 2018-1-5 22:07
统计出 benchmarking耗时的稳定性,又能如何?  发表于 2018-1-5 18:16
这个已经非常稳定了。如果你经常测试一些运算,包括GMP的性能,你会发现,性能的波动比这大的多。另外,如果你测一测maple的的log函数的性能,在计算精度达10万位时,性能的波动不是千分之几,而是百分之千(10倍)  发表于 2018-1-5 17:22
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2018-1-5 18:11:17 | 显示全部楼层
对于CPU,虽然我们可以从文档中得到每条指令的时钟周期(如http://www.agner.org/optimize/instruction_tables.pdf),但是,你如果做性能测试时会发现,这些数据最多只有一点儿参考价值而已,总的运行时间和指令的时钟周期的累加值差别很大。你永远无法计算出实际运算时间。甚至每次的运行时间都不一样。这不仅和线程切换,中断等外部因素有关,还包括大量的未知因素。过去,我的性能测试总是测试多次(如5次或者7次),然后取最快的那个结果(纵然这样,结果依然是不稳定的,波动范围在10%以内,已经是非常好的了)。但是,为了公平起见(如Maple多次测试的时间差别极大),我改为多次测试的均值。

你如果运行一下GMP源码build出来的那个speed 程序,你会发现,每次得到的结果都不一样。在通常情况下,最大值,最小值在平均值+-%5上下波动。如“speed -C -E -s 1000  mpn_mul_basecase”
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2018-1-5 18:22:10 | 显示全部楼层
mathe 发表于 2018-1-5 09:26
能够直接在C里面调用mathematica的内部接口,或者将你自己的算法编写成能够让mathematica调用的接口,然后 ...

不错的想法,需要研究下mathematica的外部接口怎么做。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-1-5 22:12:16 | 显示全部楼层
liangbch 发表于 2018-1-5 18:11
对于CPU,虽然我们可以从文档中得到每条指令的时钟周期(如http://www.agner.org/optimize/instruction_tab ...


nope, 个人认为这是很要命的。如果不事先声明数据的有效程度[方差,稳定性],我完全可以很安全的说 这些数据是没有什么实际的参考意义的。虽然从感情上,我相信你的人品,相信你有实力在某些算法细节上胜过Mathematica.  

确实,网上很多各种benchmarking根本就没给出对应的测试数据的统计方差,我平时也不太当回事。不知怎么的,我竟然在这里较真起来了。莫非只有我一个人是这么对待数据的真实性吗?
=======================================================================

数据的方差大也不要紧,但就是不能没有。 方差大了,咱们还可以另外解释,我们可以说平均性能怎么样,最好性能怎么样,最坏性能怎么样。 是吧
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-1-6 00:26:51 | 显示全部楼层
仿照例子,安装目录下【/usr/local/Wolfram/Mathematica/11.2/SystemFiles/Links/WSTP/DeveloperKit/Linux-x86-64】,写了一个测试C/C++测试程序。基本跑通需求。
存在的问题就是到底怎么设置时间的check point才算合理,不是特别的清楚。我先把代码贴出来,有兴趣的自己可以试试。暂时不打算继续较真下去了,到此为止。(反正Mathematica是输了还是赢了,跟我没任何关系, ^_^)

编译的命令如下:
  1. g++ main.cxx -I. -L. -lWSTP64i4 -lm -lpthread -lrt -ldl -luuid -std=c++11
  2. LD_LIBRARY_PATH=. ./a.out -linkname "math -wstp"
复制代码


  1. #include <iostream>
  2. #include <chrono>
  3. #include <cmath>
  4. #include "wstp.h"

  5. static WSENV ep = (WSENV)0;
  6. static WSLINK lp = (WSLINK)0;

  7. static void error( WSLINK lp)
  8. {
  9.     if( WSError( lp)){
  10.         fprintf( stderr, "Error detected by WSTP: %s.\n",
  11.                  WSErrorMessage(lp));
  12.     }else{
  13.         fprintf( stderr, "Error detected by this program.\n");
  14.     }
  15.     exit(3);
  16. }


  17. static void deinit( void)
  18. {
  19.     if( ep) WSDeinitialize( ep);
  20. }


  21. static void closelink( void)
  22. {
  23.     if( lp) WSClose( lp);
  24. }

  25. static void init_and_openlink( int argc, char* argv[])
  26. {
  27.     int err;

  28.     ep =  WSInitialize( (WSParametersPointer)0);
  29.     if( ep == (WSENV)0) exit(1);
  30.     atexit( deinit);

  31.     lp = WSOpenArgv( ep, argv, argv + argc, &err);
  32.     if(lp == (WSLINK)0) exit(2);
  33.     atexit( closelink);
  34. }


  35. static void read_and_print_expression( WSLINK lp)
  36. {
  37.     const char *s;
  38.     int n;
  39.     int i, len;
  40.     double r;
  41.     static int indent;

  42.     switch( WSGetNext( lp)) {
  43.     case WSTKSYM:
  44.         std::cout<< __LINE__ <<std::endl;
  45.         WSGetSymbol( lp, &s);
  46.         printf( "%s ", s);
  47.         WSReleaseSymbol( lp, s);
  48.         break;
  49.     case WSTKSTR:
  50.         std::cout<< __LINE__ <<std::endl;
  51.         WSGetString( lp, &s);
  52.         printf( ""%s" ", s);
  53.         WSReleaseString( lp, s);
  54.         break;
  55.     case WSTKINT:
  56.         std::cout<< __LINE__ <<std::endl;
  57.         WSGetInteger( lp, &n);
  58.         printf( "%d ", n);
  59.         break;
  60.     case WSTKREAL:
  61. //        std::cout<< __LINE__ <<std::endl;
  62.         WSGetReal( lp, &r);
  63.         printf( "%g ", r);
  64.         break;
  65.     case WSTKFUNC:
  66.         std::cout<< __LINE__ <<std::endl;
  67.         indent += 3;
  68.         printf( "\n %*.*s", indent, indent, "");
  69.         if( WSGetArgCount( lp, &len) == 0){
  70.             error( lp);
  71.         }else{
  72.             read_and_print_expression( lp);
  73.             printf( "[");
  74.             for( i = 1; i <= len; ++i){
  75.                 read_and_print_expression( lp);
  76.                 if( i != len) printf( ", ");
  77.             }
  78.             printf( "]");
  79.         }
  80.         indent -= 3;
  81.         break;
  82.     case WSTKERROR:
  83.         std::cout<< __LINE__ <<std::endl;
  84.     default:
  85.         error( lp);
  86.     }
  87. }


  88. int main(int argc, char* argv[])
  89. {
  90.     init_and_openlink( argc, argv);

  91.     int pkt,x,n;
  92.     std::cout<< "put x: "<<std::endl;
  93.     std::cin>>x;
  94.     std::cout<< "put n: "<<std::endl;
  95.     std::cin>>n;
  96.     std::string cmd = "N[Log["+std::to_string(x)+"],"+std::to_string(n)+"]";
  97.     std::cout<< "Expression: "<<cmd<<std::endl;

  98. //    WSPutFunction( lp, "N", 2);
  99. //    WSPutFunction( lp, "Log", 1);
  100. //    WSPutInteger( lp, x);
  101. //    WSPutInteger( lp, n);
  102. //    WSEndPacket( lp);

  103.     WSPutFunction(lp, "ToExpression", 1);
  104.     WSPutString(lp, cmd.c_str());
  105.     WSEndPacket(lp);

  106.     while( (pkt = WSNextPacket(lp), pkt) && pkt != RETURNPKT) {
  107.         WSNewPacket( lp);
  108.         if (WSError( lp)) error( lp);
  109.     }

  110.     auto t1 = std::chrono::high_resolution_clock::now();
  111.     read_and_print_expression( lp);
  112.     auto t2 = std::chrono::high_resolution_clock::now();
  113.     std::cout<< "cost:"<<std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()/1000.<<"ms"<<std::endl;

  114.     WSPutFunction( lp, "Exit", 0L);

  115.     return 0;

  116. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-3-17 23:47:15 | 显示全部楼层
不知论文进度怎样了
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2018-3-20 09:55:53 | 显示全部楼层
这个事情暂时放下,也没有重启的时间表
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-4-8 17:37:44 | 显示全部楼层
我提供个测试例子吧,统一使用一个很大的常数,比如10000内的质数乘到一起,结果减1,再平移小数点变成整数部分小于10的小数,预先存储成文件,然后读进程序来运算,然后再测试求时间,测试时候从读完这个小数后计时,这样就既能求的比最小测量误差大很多的时间了,又避免了可能的预先储存结果的问题
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2022-11-15 14:17:40 | 显示全部楼层
我就好奇,干啥工作需要把对数算上万位?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2022-12-13 11:14:56 | 显示全部楼层
nyy 发表于 2022-11-15 14:17
我就好奇,干啥工作需要把对数算上万位?

对技术的追求是无止尽的. 圆周率已经计算到几万亿位了,但就应用而言,几十位圆周率也足够了,不是吗?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-3-28 21:50 , Processed in 0.054743 second(s), 16 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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