不依赖外部库实现大数加、减、乘、除、开根

头像
523066680
Administrator
Administrator
帖子: 481
注册时间: 2016年07月19日 12:14
拥有现金: 锁定
储蓄: 锁定
Has thanked: 53 times
Been thanked: 93 times
联系:

不依赖外部库实现大数加、减、乘、除、开根

帖子 #1 523066680 » 2018年12月24日 16:11

语言不限。可以使用容器库等辅助操作,可以内联汇编,但是不允许直接使用大数相关的库。

happy886r 之前已经写过相当完善的计算工具,这里引用一下:
逆波兰计算器revpolish.exe
数学计算工具i

资料整理:
Coding Contest Byte: The Square Root Trick

libgmp 库中的开根算法
"Karatsuba Square Root" algorithm by Paul Zimmermann

GNU MP 库中用到的算法和理论,完整版:
《Modern Computer Arithmetic》
http://maths.anu.edu.au/~brent/pd/mca-cup-0.5.9.pdf

头像
523066680
Administrator
Administrator
帖子: 481
注册时间: 2016年07月19日 12:14
拥有现金: 锁定
储蓄: 锁定
Has thanked: 53 times
Been thanked: 93 times
联系:

Re: 不依赖外部库实现大数加、减、乘、除、开根

帖子 #2 523066680 » 2018年12月25日 14:17

先说开根吧,Wikipedia 上有很齐全的方案收集

手算开根法
Methods of computing square roots
https://en.wikipedia.org/wiki/Methods_o ... uare_roots

1. 4 1 4 2
_______________
\/ 02.00 00 00 00
02 1*1 <= 2 < 2*2 x = 1
01 y = x*x = 1*1 = 1
01 00 24*4 <= 100 < 25*5 x = 4
00 96 y = (20+x)*x = 24*4 = 96
04 00 281*1 <= 400 < 282*2 x = 1
02 81 y = (280+x)*x = 281*1 = 281
01 19 00 2824*4 <= 11900 < 2825*5 x = 4
01 12 96 y = (2820+x)*x = 2824*4 = 11296
06 04 00 28282*2 <= 60400 < 28283*3 x = 2
The desired precision is achieved:
The square root of 2 is about 1.4142

头像
523066680
Administrator
Administrator
帖子: 481
注册时间: 2016年07月19日 12:14
拥有现金: 锁定
储蓄: 锁定
Has thanked: 53 times
Been thanked: 93 times
联系:

性能分析

帖子 #3 523066680 » 2019年01月20日 21:48

照着手算的方案,撸了一段C++版的,用VS2015 做了性能分析,好像开销主要在于字符串数组的[]取址操作,
待优化

sqrt_decimal.cpp


prof.png


2019-02-03
手算法 C++ vector 容器版,BASE=10^8,i7 CPU 4GHz, 1W位0.11秒,可能某些数字会触发问题吧?不管了过年了 :confidence
// SQRT - Paper and Pencil method, C++ implementation (Only for Positive Integer Number)
// 523066680/vicyang
// 2019-02
// vector solution, BASE = 10^8, not finished yet :P

sqrt_vector_2019-02-03.cpp
您没有权限查看这个主题的附件。

头像
523066680
Administrator
Administrator
帖子: 481
注册时间: 2016年07月19日 12:14
拥有现金: 锁定
储蓄: 锁定
Has thanked: 53 times
Been thanked: 93 times
联系:

大数减法、加法 C++实现 Re: 不依赖外部库实现大数加、减、乘、除、开根

帖子 #4 523066680 » 2019年02月01日 15:43

使用向量容器存储 32位int值,每一段int存储8位数字( 按10^8 为进制处理)

减法运算效率测试(分为 vec_minus 和 s_minus 函数)
Code: [全选] [展开/折叠] [Download] (vec_minus.cpp)
  1. // bignum minus
  2. // 523066680/vicyang
  3. // 2019-01
  4.  
  5. #include <iostream>
  6. #include <chrono>
  7. #include <vector>
  8. #include <cstdio>
  9. using namespace std;
  10.  
  11. vector<int> vec_minus(const vector<int> &a, const vector<int> &b);
  12. string s_minus(const string &a, const string &b);
  13. int s_cmp( const string &a, const string &b );
  14. string vec2str( const vector<int> &vec );
  15.  
  16. void check(const vector<int> &va, const vector<int> &vb)
  17. {
  18.     string a = vec2str(va);
  19.     string b = vec2str(vb);
  20.     string cmd = "perl -Mbignum -e \"print " + a + "-" + b + "\"";
  21.     system( cmd.c_str() );
  22.     cout << " <- check " << endl;
  23. }
  24.  
  25. int main(int argc, char *argv[] )
  26. {
  27.     auto stage0 = chrono::system_clock::now();
  28.     chrono::duration<double> diff;
  29.     string a(10000, '8');
  30.     string b(10000, '9');
  31.     string c;
  32.     // s_minus 耗时测试 10000位数,2W次
  33.     for (int i = 0; i < 20000; i++) s_minus(a, b);
  34.     auto stage1 = chrono::system_clock::now();
  35.     diff = stage1-stage0;
  36.     cout << "  s_minus, Time used: " << diff.count() << endl;
  37.  
  38.     // 小范围的"向量"减法测试
  39.     vector<int> va{552, 443, 12345678, 87654321};
  40.     vector<int> vb{ 93456, 92432100};
  41.     //check(va, vb);
  42.     vector<int> vc = vec_minus(va,vb);
  43.     cout << vec2str(vc) << endl;
  44.  
  45.     // 结果应为1的测试
  46.     va = {1, 0, 0};
  47.     vb = {99999999,99999999};
  48.     //check(va, vb);
  49.     vc = vec_minus(va, vb);
  50.     cout << vec2str(vc) << endl;
  51.  
  52.     // vec_minus 耗时测试 10000位数,2W次
  53.     stage1 = chrono::system_clock::now();
  54.     va.assign( 1250, 99999999 );
  55.     vb.assign( 1250, 12345678 );
  56.     for (int i = 0; i < 20000; i++) vec_minus(va, vb);
  57.     auto stage2 = chrono::system_clock::now();
  58.     diff = stage2-stage1;
  59.     cout << "vec_minus, Time used: " << diff.count() << endl;
  60.     return 0;
  61. }
  62.  
  63. // 测试vector方案
  64. vector<int> vec_minus(const vector<int> &a, const vector<int> &b)
  65. {
  66.     static int ia; // iter
  67.     const int base = 100000000;
  68.     register const int *pa = a.data();
  69.     register const int *pb = b.data();
  70.     vector<int> c( a.size() );
  71.     register int *pc = c.data();
  72.     int t, cut=0, ib=b.size()-1, zero=0;
  73.     for (ia = a.size()-1; ia >= 0; ia-- )
  74.     {
  75.         t = ib >= 0 ? (pa[ia]) - (pb[ib--]) + cut
  76.                     : (pa[ia]) + cut;
  77.         t < 0 ? t += base, cut = -1 : cut = 0;
  78.         zero = t == 0 ? zero+1 : 0;  // 此判断须独立,t有可能+base后才为0
  79.         pc[ia] = t;
  80.     }
  81.     c.erase(c.begin(), c.begin()+zero);
  82.     return c;
  83. }
  84.  
  85. // 大数减法 字符串操作, 暂时假设 a >= b
  86. string s_minus(const string &a, const string &b)
  87. {
  88.     static int ia;
  89.     // 获取指针以绕过[]重载,更快地取址
  90.     register const char* pa=a.data();
  91.     register const char* pb=b.data();
  92.     int cmp = s_cmp(a, b);
  93.     if (cmp == 0) return "0";
  94.     string s( a.length(), '0');
  95.     register const char* ps=s.data();
  96.     int t, cut=0, ib=b.length()-1, zero=0;
  97.     for (ia = a.length()-1; ia >= 0; ia-- )
  98.     {
  99.         t = ib >= 0 ? (pa[ia]-'0') - (pb[ib--]-'0') - cut
  100.                     : (pa[ia]-'0') - cut;
  101.         cut = t < 0 ? 1 : 0;
  102.         s[ia] = t + '0' + cut*10;
  103.         ps[ia] == '0' ? zero++ : zero=0;
  104.     }
  105.     if (zero > 0) s.erase(0, zero);
  106.     return s;
  107. }
  108.  
  109. string vec2str( const vector<int> &vec )
  110. {
  111.     const int base = 100000000;
  112.     string s("");
  113.     s += to_string( vec[0] );
  114.     for ( int it = 1; it < vec.size(); it++ )
  115.         s += to_string(vec[it]+base).substr(1,8);
  116.     return s;
  117. }
  118.  
  119. int s_cmp(const string &a, const string &b )
  120. {
  121.     if ( a.length() > b.length() ) return 1;
  122.     else if ( a.length() < b.length() ) return -1;
  123.     else return a.compare(b);
  124. }

1W位数字减法,2W次循环,向量(10^8 进制) 和 string 逐位操作,效率对比

代码: 全选

  s_minus, Time used: 0.860001
vec_minus, Time used: 0.12


《Modern Computer Arithmetic》(后面简称MCA)里面关于采用高进制计算可能产生溢出问题的探讨
对于64位平台 unsigned long long int 支持的最大数字是 2^64-1=18446744073709551615,20位,如果我们充分利用,使用2^64,或者10^19作为进制。
很容易会遇到溢出问题,MCA给出了几种方案:
Let T be the number of different values taken by the data type representing the coefficients ai, bi.
(Clearly, β ≤ T, but equality does not necessarily hold, for example β = 10^9 and T = 2^32.) At step 3, the value of s can be as large as 2β − 1, which is not representable if β = T. Several workarounds are possible:
either use a machine instruction that gives the possible carry of ai + bi,
or use the fact that, if a carry occurs in ai + bi, then the computed sum – if performed modulo T – equals t := ai +bi −T < ai; thus, comparing t and ai will determine if a carry occurred.
A third solution is to keep a bit in reserve, taking β ≤ T/2.

用 T 表示一个存储单元所能表示的数的量,则有 β <= T
(这里 β 表示采用的进制,以及β不一定等于T,例如T=2^32,但采用的进制为10^9)。
考虑加法操作 s=a+b+d (d为1或0,是上一次加法补进的数值),s的最大可能值为 2*β-1,当 β = T 时该公式无法正确计算。
考虑以下方案:
1. 内部编码实现 (水平有限,暂时这么翻译)
2. 通过-T取余数判断,t := a + b - T < a ,对比 t 和 a 断定是否进 1
3. 采用一个保守的进制数 β,令 β ≤ T/2.

偷个懒,暂时采用10^8而不是 10^9, 2^32
如果对相同的进制进行乘法运算,可以使用 unsigned long long int 类型获得最高19位数的运算空间,足以应付 8*2+1 位数的情况

头像
523066680
Administrator
Administrator
帖子: 481
注册时间: 2016年07月19日 12:14
拥有现金: 锁定
储蓄: 锁定
Has thanked: 53 times
Been thanked: 93 times
联系:

大数乘法 - BasecaseMultiply C++ 实现 Re: 不依赖外部库实现大数加、减、乘、除、开根

帖子 #5 523066680 » 2019年02月01日 22:11

《MCA》中的BasecaseMultiply 函数实现,采用10^8为进制。

BasecaseMultiply.png

2019-02-01 初步版本

Code: [全选] [展开/折叠] [Download] (BasecaseMultiply.cpp)
  1. // BasecaseMultiply - C++ implementation
  2. // 523066680/vicyang
  3. // 2019-02
  4.  
  5. #include <iostream>
  6. #include <string>
  7. #include <vector>
  8. #include <chrono>
  9. using namespace std;
  10. typedef unsigned long long ULL;
  11.  
  12. string vec2str( const vector<ULL> &vec );
  13. vector<ULL> vec_plus(const vector<ULL> &a, const vector<ULL> &b);
  14. vector<ULL> BasecaseMultiply( const vector<ULL>& a, const vector<ULL>& b);
  15. vector<ULL> vec_mp_single( const vector<ULL>& a, int b);
  16. void shift( vector<ULL>& vec, int n );
  17.  
  18. const unsigned long long BASE = 100000000;
  19. const int MAXLEN = 8;
  20.  
  21. int main(int argc, char *argv[] )
  22. {
  23.     vector<ULL> a{777, 0, 12345678};
  24.     vector<ULL> b{6, 99999999, 99999999};
  25.     vector<ULL> c;
  26.     c = BasecaseMultiply( a, b );
  27.     cout << vec2str(c);
  28.     return 0;
  29. }
  30.  
  31. vector<ULL> vec_mp_single( const vector<ULL>& a, int b)
  32. {
  33.     vector<ULL> c( a.size() );
  34.     if ( b == 0 ) { return vector<ULL>{0}; }
  35.     ULL pool = 0, v;
  36.     for ( int i = a.size()-1; i >=0 ; i-- ) {
  37.         v = a[i] * b + pool;
  38.         c[i] = v % BASE, pool = v / BASE;
  39.     }
  40.     if (pool > 0) c.insert( c.begin(), pool );
  41.     return c;
  42. }
  43.  
  44. // 参数:向量引用、 前移的段数(也许需要考虑vec为0的情况)
  45. void shift( vector<ULL>& vec, int n )
  46. {
  47.     for (int i = 1; i <= n; i++) vec.push_back(0);
  48. }
  49.  
  50. // 假设 b.size() >= a.size()
  51. vector<ULL> BasecaseMultiply( const vector<ULL>& a, const vector<ULL>& b)
  52. {
  53.     vector<ULL> c;
  54.     vector<ULL> t;
  55.     int bi = b.size() - 1, indent = 1;
  56.     c = vec_mp_single( a, b[bi--] );
  57.     while ( bi >= 0 )
  58.     {
  59.         if ( b[bi] > 0 )
  60.         {
  61.             t = vec_mp_single(a, b[bi]);
  62.             shift(t, indent);
  63.             c = vec_plus(t, c);
  64.         }
  65.         bi--, indent++;
  66.     }
  67.     return c;
  68. }
  69.  
  70. // 测试vector作为参数
  71. vector<ULL> vec_plus(const vector<ULL> &a, const vector<ULL> &b)
  72. {
  73.     static int ia; // iter
  74.     vector<ULL> c( a.size() );
  75.     int t, pool=0, ib=b.size()-1;
  76.     int v, r;
  77.     for (ia = a.size()-1; ia >= 0; ia-- )
  78.     {
  79.         t = ib >= 0 ? (a[ia]) + (b[ib--]) + pool
  80.                     : (a[ia]) + pool;
  81.         c[ia] = t % BASE, pool = t/BASE;
  82.     }
  83.     if ( pool > 0 ) c.insert(c.begin(), pool);
  84.     return c;
  85. }
  86.  
  87. string vec2str( const vector<ULL> &vec )
  88. {
  89.     string s("");
  90.     s += to_string( vec[0] );
  91.     for ( int it = 1; it < vec.size(); it++ )
  92.         s += to_string(vec[it]+BASE).substr(1, MAXLEN);
  93.     return s;
  94. }


批量测试,1W位*1W位,100次,本机耗时 4.2秒,仍有很大优化空间
(按字符串处理的版本,1W位乘法10次已经耗时13秒。)

其中 容器相关的 []操作符重载、push_back、insert操作占较大比例。
使用"指针"代替容器[]操作符,使用 emplace_back 代替 push_back 可以获得轻微的提高。
考虑提前预判新容器的长度并使用静态大小的容器,又或者换成C数组操作(不要 insert 和 push_back,提供offset和length)

优化版,“静态”容器
Code: [全选] [展开/折叠] [Download] (BasecaseMultiply_Opt.cpp)
  1. // BasecaseMultiply - C++ Static Container
  2. // 523066680/vicyang
  3. // 2019-01
  4.  
  5. #include <iostream>
  6. #include <string>
  7. #include <vector>
  8. #include <chrono>
  9. using namespace std;
  10. using namespace std::chrono;
  11. typedef unsigned long long ULL;
  12.  
  13. string vec2str( const vector<ULL> &vec );
  14. vector<ULL> BasecaseMultiply( const vector<ULL>& a, const vector<ULL>& b);
  15. int vec_accum(vector<ULL> &a, vector<ULL> &b, int bbegin);
  16. int vec_mp_single( const vector<ULL>& a, ULL b, vector<ULL>& c, ULL indent);
  17. void time_used(system_clock::time_point& time_a);
  18.  
  19. const ULL BASE = 100000000;
  20. const int MAXLEN = 8;
  21.  
  22. int main(int argc, char *argv[] )
  23. {
  24.     system_clock::time_point start;
  25.     vector<ULL> a{777, 0, 12345678};
  26.     vector<ULL> b{6, 0, 1, 99999999, 99999999};
  27.     vector<ULL> c;
  28.     c = BasecaseMultiply( a, b );
  29.     if (c[0] == 0ULL ) c.erase( c.begin() );
  30.     cout << vec2str(c) << endl;
  31.  
  32.     start = system_clock::now();
  33.     a.assign( 1250, 99999999 );
  34.     b.assign( 1250, 11111111 );
  35.     for (int i=0; i<100; i++) c = BasecaseMultiply( a, b );
  36.     //cout << vec2str(c) << endl;
  37.     time_used( start );
  38.  
  39.     return 0;
  40. }
  41.  
  42. // 假设 b.size() >= a.size()
  43. vector<ULL> BasecaseMultiply( const vector<ULL>& a, const vector<ULL>& b)
  44. {
  45.     //使用固定尺寸的容器
  46.     vector<ULL> c(a.size() + b.size(), 0);
  47.     vector<ULL> t(a.size() + b.size(), 0);
  48.     register const ULL* pa = a.data();
  49.     register const ULL* pb = b.data();
  50.     int tpos = 1, tbegin;
  51.     int bi = b.size() - 1, indent = 1;
  52.     //初始化C值,末尾offset为0
  53.     vec_mp_single( a, pb[bi--], c, 0 );
  54.     while ( bi >= 0 )
  55.     {
  56.         //如果对应的b==0则不计入
  57.         if ( pb[bi] > 0 )
  58.         {   //参数:vector_a, int_b, vector_t, t的末尾offset;返回t的实际起始下标
  59.             tbegin = vec_mp_single(a, pb[bi], t, tpos);
  60.             vec_accum(c, t, tbegin);
  61.         }
  62.         bi--, indent++, tpos++;
  63.     }
  64.     return c;
  65. }
  66.  
  67. int vec_mp_single( const vector<ULL>& a, ULL b, vector<ULL>& c, ULL indent)
  68. {
  69.     register const ULL *pa = a.data();
  70.     register ULL *pc = c.data();
  71.     if ( b == 0 ) { c[0] = 0; return 1; }
  72.     //实际末位置
  73.     int ci = c.size() - 1 - indent;
  74.     ULL pool = 0, v;
  75.     for ( int i = a.size()-1; i >=0 ; i-- ) {
  76.         v = pa[i] * b + pool;
  77.         pc[ci--] = v % BASE, pool = v/BASE;
  78.     }
  79.     if (pool > 0) { pc[ci--] = pool; }
  80.     //返回实际起始位置
  81.     return ci+1;
  82. }
  83.  
  84. // 在原来的基础上叠加数值
  85. int vec_accum(vector<ULL> &a, vector<ULL> &b, int bbegin)
  86. {
  87.     static int ia; // iter
  88.     register ULL* pa = a.data();
  89.     register ULL* pb = b.data();
  90.     ULL t, pool=0;
  91.     register int ib = b.size()-1;
  92.     for (ia = a.size()-1; ia >= 0; ia-- )
  93.     {
  94.         t = ib >= 0 ? (pa[ia]) + (pb[ib--]) + pool
  95.                     : (pa[ia]) + pool;
  96.         pa[ia] = t % BASE, pool = t/BASE;
  97.         if ( pool == 0 and ib < bbegin ) break;
  98.     }
  99.     if ( pool > 0ULL ) pa[ia] = pool;
  100.     return ia;
  101. }
  102.  
  103. string vec2str( const vector<ULL> &vec )
  104. {
  105.     string s("");
  106.     s += to_string( vec[0] );
  107.     for (unsigned int it = 1; it < vec.size(); it++ )
  108.         s += to_string(vec[it]+BASE).substr(1, MAXLEN);
  109.     return s;
  110. }
  111.  
  112. void time_used(system_clock::time_point& time_a) {
  113.     duration<double> diff;
  114.     diff = chrono::system_clock::now() - time_a;
  115.     cout << "Time used: " << diff.count() << endl;
  116.     time_a = chrono::system_clock::now();
  117. }

对比初版,耗时为1.66s,时间缩短一半以上。
下一次尝试 KaratsubaMultiply (此算法仅针对两个乘数长度相同的情况),要等过年后了,放松一下。
您没有权限查看这个主题的附件。

头像
rubyish
崭露头角
崭露头角
帖子: 25
注册时间: 2018年04月23日 09:58
拥有现金: 锁定
Has thanked: 16 times
Been thanked: 13 times
联系:

Re: 不依赖外部库实现大数加、减、乘、除、开根

帖子 #6 rubyish » 2019年02月03日 00:24

這個難度太高 :grimace :grimace
$_


回到 “算法和编码”

在线用户

用户浏览此论坛: 没有注册用户 和 2 访客