语言不限。可以使用容器库等辅助操作,可以内联汇编,但是不允许直接使用大数相关的库。
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
- 帖子: 492
- 注册时间: 2016年07月19日 12:14
- 拥有现金: 锁定
- 储蓄: 锁定
- Has thanked: 58 times
- Been thanked: 96 times
- 联系:
Re: 不依赖外部库实现大数加、减、乘、除、开根
先说开根吧,Wikipedia 上有很齐全的方案收集
手算开根法
Methods of computing square roots
https://en.wikipedia.org/wiki/Methods_o ... uare_roots
手算开根法
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
- 帖子: 492
- 注册时间: 2016年07月19日 12:14
- 拥有现金: 锁定
- 储蓄: 锁定
- Has thanked: 58 times
- Been thanked: 96 times
- 联系:
性能分析
照着手算的方案,撸了一段C++版的,用VS2015 做了性能分析,好像开销主要在于字符串数组的[]取址操作,
待优化
2019-02-03
手算法 C++ vector 容器版,BASE=10^8,i7 CPU 4GHz, 1W位0.11秒,可能某些数字会触发问题吧?不管了过年了
待优化
2019-02-03
手算法 C++ vector 容器版,BASE=10^8,i7 CPU 4GHz, 1W位0.11秒,可能某些数字会触发问题吧?不管了过年了

// SQRT - Paper and Pencil method, C++ implementation (Only for Positive Integer Number)
// 523066680/vicyang
// 2019-02
// vector solution, BASE = 10^8, not finished yet
您没有权限查看这个主题的附件。
- 523066680
- Administrator
- 帖子: 492
- 注册时间: 2016年07月19日 12:14
- 拥有现金: 锁定
- 储蓄: 锁定
- Has thanked: 58 times
- Been thanked: 96 times
- 联系:
大数减法、加法 C++实现 Re: 不依赖外部库实现大数加、减、乘、除、开根
使用向量容器存储 32位int值,每一段int存储8位数字( 按10^8 为进制处理)
减法运算效率测试(分为 vec_minus 和 s_minus 函数)
1W位数字减法,2W次循环,向量(10^8 进制) 和 string 逐位操作,效率对比
《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 位数的情况
减法运算效率测试(分为 vec_minus 和 s_minus 函数)
- // bignum minus
- // 523066680/vicyang
- // 2019-01
- #include <iostream>
- #include <chrono>
- #include <vector>
- #include <cstdio>
- using namespace std;
- vector<int> vec_minus(const vector<int> &a, const vector<int> &b);
- string s_minus(const string &a, const string &b);
- int s_cmp( const string &a, const string &b );
- string vec2str( const vector<int> &vec );
- void check(const vector<int> &va, const vector<int> &vb)
- {
- string a = vec2str(va);
- string b = vec2str(vb);
- string cmd = "perl -Mbignum -e \"print " + a + "-" + b + "\"";
- system( cmd.c_str() );
- cout << " <- check " << endl;
- }
- int main(int argc, char *argv[] )
- {
- auto stage0 = chrono::system_clock::now();
- chrono::duration<double> diff;
- string a(10000, '8');
- string b(10000, '9');
- string c;
- // s_minus 耗时测试 10000位数,2W次
- for (int i = 0; i < 20000; i++) s_minus(a, b);
- auto stage1 = chrono::system_clock::now();
- diff = stage1-stage0;
- cout << " s_minus, Time used: " << diff.count() << endl;
- // 小范围的"向量"减法测试
- vector<int> va{552, 443, 12345678, 87654321};
- vector<int> vb{ 93456, 92432100};
- //check(va, vb);
- vector<int> vc = vec_minus(va,vb);
- cout << vec2str(vc) << endl;
- // 结果应为1的测试
- va = {1, 0, 0};
- vb = {99999999,99999999};
- //check(va, vb);
- vc = vec_minus(va, vb);
- cout << vec2str(vc) << endl;
- // vec_minus 耗时测试 10000位数,2W次
- stage1 = chrono::system_clock::now();
- va.assign( 1250, 99999999 );
- vb.assign( 1250, 12345678 );
- for (int i = 0; i < 20000; i++) vec_minus(va, vb);
- auto stage2 = chrono::system_clock::now();
- diff = stage2-stage1;
- cout << "vec_minus, Time used: " << diff.count() << endl;
- return 0;
- }
- // 测试vector方案
- vector<int> vec_minus(const vector<int> &a, const vector<int> &b)
- {
- static int ia; // iter
- const int base = 100000000;
- register const int *pa = a.data();
- register const int *pb = b.data();
- vector<int> c( a.size() );
- register int *pc = c.data();
- int t, cut=0, ib=b.size()-1, zero=0;
- for (ia = a.size()-1; ia >= 0; ia-- )
- {
- t = ib >= 0 ? (pa[ia]) - (pb[ib--]) + cut
- : (pa[ia]) + cut;
- t < 0 ? t += base, cut = -1 : cut = 0;
- zero = t == 0 ? zero+1 : 0; // 此判断须独立,t有可能+base后才为0
- pc[ia] = t;
- }
- c.erase(c.begin(), c.begin()+zero);
- return c;
- }
- // 大数减法 字符串操作, 暂时假设 a >= b
- string s_minus(const string &a, const string &b)
- {
- static int ia;
- // 获取指针以绕过[]重载,更快地取址
- register const char* pa=a.data();
- register const char* pb=b.data();
- int cmp = s_cmp(a, b);
- if (cmp == 0) return "0";
- string s( a.length(), '0');
- register const char* ps=s.data();
- int t, cut=0, ib=b.length()-1, zero=0;
- for (ia = a.length()-1; ia >= 0; ia-- )
- {
- t = ib >= 0 ? (pa[ia]-'0') - (pb[ib--]-'0') - cut
- : (pa[ia]-'0') - cut;
- cut = t < 0 ? 1 : 0;
- s[ia] = t + '0' + cut*10;
- ps[ia] == '0' ? zero++ : zero=0;
- }
- if (zero > 0) s.erase(0, zero);
- return s;
- }
- string vec2str( const vector<int> &vec )
- {
- const int base = 100000000;
- string s("");
- s += to_string( vec[0] );
- for ( int it = 1; it < vec.size(); it++ )
- s += to_string(vec[it]+base).substr(1,8);
- return s;
- }
- int s_cmp(const string &a, const string &b )
- {
- if ( a.length() > b.length() ) return 1;
- else if ( a.length() < b.length() ) return -1;
- else return a.compare(b);
- }
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
- 帖子: 492
- 注册时间: 2016年07月19日 12:14
- 拥有现金: 锁定
- 储蓄: 锁定
- Has thanked: 58 times
- Been thanked: 96 times
- 联系:
大数乘法 - BasecaseMultiply C++ 实现 Re: 不依赖外部库实现大数加、减、乘、除、开根
《MCA》中的BasecaseMultiply 函数实现,采用10^8为进制。
2019-02-01 初步版本
批量测试,1W位*1W位,100次,本机耗时 4.2秒,仍有很大优化空间
(按字符串处理的版本,1W位乘法10次已经耗时13秒。)
其中 容器相关的 []操作符重载、push_back、insert操作占较大比例。
使用"指针"代替容器[]操作符,使用 emplace_back 代替 push_back 可以获得轻微的提高。
考虑提前预判新容器的长度并使用静态大小的容器,又或者换成C数组操作(不要 insert 和 push_back,提供offset和length)
优化版,“静态”容器
对比初版,耗时为1.66s,时间缩短一半以上。
下一次尝试 KaratsubaMultiply (此算法仅针对两个乘数长度相同的情况),要等过年后了,放松一下。
2019-02-01 初步版本
- // BasecaseMultiply - C++ implementation
- // 523066680/vicyang
- // 2019-02
- #include <iostream>
- #include <string>
- #include <vector>
- #include <chrono>
- using namespace std;
- typedef unsigned long long ULL;
- string vec2str( const vector<ULL> &vec );
- vector<ULL> vec_plus(const vector<ULL> &a, const vector<ULL> &b);
- vector<ULL> BasecaseMultiply( const vector<ULL>& a, const vector<ULL>& b);
- vector<ULL> vec_mp_single( const vector<ULL>& a, int b);
- void shift( vector<ULL>& vec, int n );
- const unsigned long long BASE = 100000000;
- const int MAXLEN = 8;
- int main(int argc, char *argv[] )
- {
- vector<ULL> a{777, 0, 12345678};
- vector<ULL> b{6, 99999999, 99999999};
- vector<ULL> c;
- c = BasecaseMultiply( a, b );
- cout << vec2str(c);
- return 0;
- }
- vector<ULL> vec_mp_single( const vector<ULL>& a, int b)
- {
- vector<ULL> c( a.size() );
- if ( b == 0 ) { return vector<ULL>{0}; }
- ULL pool = 0, v;
- for ( int i = a.size()-1; i >=0 ; i-- ) {
- v = a[i] * b + pool;
- c[i] = v % BASE, pool = v / BASE;
- }
- if (pool > 0) c.insert( c.begin(), pool );
- return c;
- }
- // 参数:向量引用、 前移的段数(也许需要考虑vec为0的情况)
- void shift( vector<ULL>& vec, int n )
- {
- for (int i = 1; i <= n; i++) vec.push_back(0);
- }
- // 假设 b.size() >= a.size()
- vector<ULL> BasecaseMultiply( const vector<ULL>& a, const vector<ULL>& b)
- {
- vector<ULL> c;
- vector<ULL> t;
- int bi = b.size() - 1, indent = 1;
- c = vec_mp_single( a, b[bi--] );
- while ( bi >= 0 )
- {
- if ( b[bi] > 0 )
- {
- t = vec_mp_single(a, b[bi]);
- shift(t, indent);
- c = vec_plus(t, c);
- }
- bi--, indent++;
- }
- return c;
- }
- // 测试vector作为参数
- vector<ULL> vec_plus(const vector<ULL> &a, const vector<ULL> &b)
- {
- static int ia; // iter
- vector<ULL> c( a.size() );
- int t, pool=0, ib=b.size()-1;
- int v, r;
- for (ia = a.size()-1; ia >= 0; ia-- )
- {
- t = ib >= 0 ? (a[ia]) + (b[ib--]) + pool
- : (a[ia]) + pool;
- c[ia] = t % BASE, pool = t/BASE;
- }
- if ( pool > 0 ) c.insert(c.begin(), pool);
- return c;
- }
- string vec2str( const vector<ULL> &vec )
- {
- string s("");
- s += to_string( vec[0] );
- for ( int it = 1; it < vec.size(); it++ )
- s += to_string(vec[it]+BASE).substr(1, MAXLEN);
- return s;
- }
批量测试,1W位*1W位,100次,本机耗时 4.2秒,仍有很大优化空间
(按字符串处理的版本,1W位乘法10次已经耗时13秒。)
其中 容器相关的 []操作符重载、push_back、insert操作占较大比例。
使用"指针"代替容器[]操作符,使用 emplace_back 代替 push_back 可以获得轻微的提高。
考虑提前预判新容器的长度并使用静态大小的容器,又或者换成C数组操作(不要 insert 和 push_back,提供offset和length)
优化版,“静态”容器
- // BasecaseMultiply - C++ Static Container
- // 523066680/vicyang
- // 2019-01
- #include <iostream>
- #include <string>
- #include <vector>
- #include <chrono>
- using namespace std;
- using namespace std::chrono;
- typedef unsigned long long ULL;
- string vec2str( const vector<ULL> &vec );
- vector<ULL> BasecaseMultiply( const vector<ULL>& a, const vector<ULL>& b);
- int vec_accum(vector<ULL> &a, vector<ULL> &b, int bbegin);
- int vec_mp_single( const vector<ULL>& a, ULL b, vector<ULL>& c, ULL indent);
- void time_used(system_clock::time_point& time_a);
- const ULL BASE = 100000000;
- const int MAXLEN = 8;
- int main(int argc, char *argv[] )
- {
- system_clock::time_point start;
- vector<ULL> a{777, 0, 12345678};
- vector<ULL> b{6, 0, 1, 99999999, 99999999};
- vector<ULL> c;
- c = BasecaseMultiply( a, b );
- if (c[0] == 0ULL ) c.erase( c.begin() );
- cout << vec2str(c) << endl;
- start = system_clock::now();
- a.assign( 1250, 99999999 );
- b.assign( 1250, 11111111 );
- for (int i=0; i<100; i++) c = BasecaseMultiply( a, b );
- //cout << vec2str(c) << endl;
- time_used( start );
- return 0;
- }
- // 假设 b.size() >= a.size()
- vector<ULL> BasecaseMultiply( const vector<ULL>& a, const vector<ULL>& b)
- {
- //使用固定尺寸的容器
- vector<ULL> c(a.size() + b.size(), 0);
- vector<ULL> t(a.size() + b.size(), 0);
- register const ULL* pa = a.data();
- register const ULL* pb = b.data();
- int tpos = 1, tbegin;
- int bi = b.size() - 1, indent = 1;
- //初始化C值,末尾offset为0
- vec_mp_single( a, pb[bi--], c, 0 );
- while ( bi >= 0 )
- {
- //如果对应的b==0则不计入
- if ( pb[bi] > 0 )
- { //参数:vector_a, int_b, vector_t, t的末尾offset;返回t的实际起始下标
- tbegin = vec_mp_single(a, pb[bi], t, tpos);
- vec_accum(c, t, tbegin);
- }
- bi--, indent++, tpos++;
- }
- return c;
- }
- int vec_mp_single( const vector<ULL>& a, ULL b, vector<ULL>& c, ULL indent)
- {
- register const ULL *pa = a.data();
- register ULL *pc = c.data();
- if ( b == 0 ) { c[0] = 0; return 1; }
- //实际末位置
- int ci = c.size() - 1 - indent;
- ULL pool = 0, v;
- for ( int i = a.size()-1; i >=0 ; i-- ) {
- v = pa[i] * b + pool;
- pc[ci--] = v % BASE, pool = v/BASE;
- }
- if (pool > 0) { pc[ci--] = pool; }
- //返回实际起始位置
- return ci+1;
- }
- // 在原来的基础上叠加数值
- int vec_accum(vector<ULL> &a, vector<ULL> &b, int bbegin)
- {
- static int ia; // iter
- register ULL* pa = a.data();
- register ULL* pb = b.data();
- ULL t, pool=0;
- register int ib = b.size()-1;
- for (ia = a.size()-1; ia >= 0; ia-- )
- {
- t = ib >= 0 ? (pa[ia]) + (pb[ib--]) + pool
- : (pa[ia]) + pool;
- pa[ia] = t % BASE, pool = t/BASE;
- if ( pool == 0 and ib < bbegin ) break;
- }
- if ( pool > 0ULL ) pa[ia] = pool;
- return ia;
- }
- string vec2str( const vector<ULL> &vec )
- {
- string s("");
- s += to_string( vec[0] );
- for (unsigned int it = 1; it < vec.size(); it++ )
- s += to_string(vec[it]+BASE).substr(1, MAXLEN);
- return s;
- }
- void time_used(system_clock::time_point& time_a) {
- duration<double> diff;
- diff = chrono::system_clock::now() - time_a;
- cout << "Time used: " << diff.count() << endl;
- time_a = chrono::system_clock::now();
- }
对比初版,耗时为1.66s,时间缩短一半以上。
下一次尝试 KaratsubaMultiply (此算法仅针对两个乘数长度相同的情况),要等过年后了,放松一下。
您没有权限查看这个主题的附件。
在线用户
用户浏览此论坛: 没有注册用户 和 0 访客