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

回复
头像
523066680
Administrator
Administrator
帖子: 573
注册时间: 2016年07月19日 12:14
联系:

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

帖子 523066680 »

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

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
帖子: 573
注册时间: 2016年07月19日 12:14
联系:

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

帖子 523066680 »

先说开根吧,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
帖子: 573
注册时间: 2016年07月19日 12:14
联系:

性能分析

帖子 523066680 »

照着手算的方案,撸了一段C++版的,用VS2015 做了性能分析,好像开销主要在于字符串数组的[]取址操作,
待优化
sqrt_decimal.cpp
(5.4 KiB) 已下载 91 次
prof.png
(28.38 KiB) 已下载 239 次
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
(5.52 KiB) 已下载 75 次
头像
523066680
Administrator
Administrator
帖子: 573
注册时间: 2016年07月19日 12:14
联系:

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

帖子 523066680 »

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

减法运算效率测试(分为 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
Administrator
帖子: 573
注册时间: 2016年07月19日 12:14
联系:

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

帖子 523066680 »

《MCA》中的BasecaseMultiply 函数实现,采用10^8为进制。
BasecaseMultiply.png
(22.16 KiB) 已下载 36 次
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 (此算法仅针对两个乘数长度相同的情况),要等过年后了,放松一下。
头像
rubyish
渐入佳境
渐入佳境
帖子: 52
注册时间: 2018年04月23日 09:58
联系:

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

帖子 rubyish »

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

在线用户

正浏览此版面之用户: 没有注册用户 和 1 访客