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

帖子 #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

帖子 #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

帖子 #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

sqrt_vector_2019-02-03.cpp

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

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

减法运算效率测试（分为 vec_minus 和 s_minus 函数）

Code: [全选] [展开/折叠] [Download] (vec_minus.cpp)

// 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 位数的情况

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

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

BasecaseMultiply.png

2019-02-01 初步版本

Code: [全选] [展开/折叠] [Download] (BasecaseMultiply.cpp)

// 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）

优化版，“静态”容器

Code: [全选] [展开/折叠] [Download] (BasecaseMultiply_Opt.cpp)

// 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 · 帖子 #6 由 **rubyish** » 2019年02月03日 00:24

這個難度太高

综合型编程论坛

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

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

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

性能分析

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

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

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

在线用户