Click here to Skip to main content
13,624,914 members
Click here to Skip to main content
Add your own
alternative version

Stats

6.9K views
104 downloads
7 bookmarked
Posted 10 Aug 2016
Licenced CPOL

N: An experimental multithreaded bignum library for Windows

, 10 Aug 2016
Rate this:
Please Sign up or sign in to vote.
A quick implementation of my ideas

Introduction

Yes this is not optimal, windows only, not fully complete. Why bother?

Dunno, just because both you and me are interested in low level stuff. Keep reading!

Background

Numbers are nice, but they can only be as large as 2^64 in a form of unsigned long long. Sometimes you need math between much larger numbers (for example, in cryptography). This class allows math between any size big number, as long as your available memory can handle it.

Multicore Library

For multicore operations, the library uses my own tpoollib, which uses the new Vista+ API for threadpools.

The N class internals

N stores all its digits in a vector<D> (D being a typedef for int at the moment), as a decimal array, a base value (default 10) to know what base system the digits are, and a bool flag to keep the sign value. Construction is done via either an integer, or a string:

N n(5LL);
N n(-12438LL);
N n("384747274747282882818373"); // If a non digit is encountered, parsing stops
N n("-274727372737");
N n("111",2); // n = 7

N has also a member function ParseBase:

N n;
n.ParseBase("FFFF",16);  // n = 65535

N representation and access

In the debug version of the library, after each operation that changes N a special function is called to aid you in debugging.

N::operator[](size_t) gives the specific digit from the right, for example if n = 12345 then n[0] is 5. If the specified digit exceeds the number, the function does not crash, but returns 0 instead, allowing itself to be used when adding or multiplying.

Function s(unsigned long long base = b); returns a string representation of N, in the specified base system.

Function tobase(unsigned long long nb) returns a new N which stores the same value in a different base system.

Function upperpart() and lowerpart() return a subpart of the number:

N a("12345");
auto a2 = a.upperpart(2); // a2 = 12
auto a3 = a.lowerpart(3); // a3 = 345
a.Set(16LL);
auto str = a.s();      // str = "16"
str = a.tobase(2).s(); // str = "10000"
str = a.tobase(16).s(); // str = "10"
str = a.tobase(10).s(); / str = "16"

Single Add

Adding two same-signed Ns (via operators or function calls or subtracting different-sign Ns) eventually calls
static N w_add(const N& n1, const N& n2);

This function loops through all digits of the numbers and adds them with carry:

static N w_add(const N& n1, const N& n2)
    {
    if (n1.b != n2.b)
        return w_add(n1, n2.tobase(n1.b));
    if (n1.m != n2.m)
        {
        if (n1.m)
            return w_subx(n2, n1.negative());
        return w_subx(n1, n2.negative());
        }
    if (n1.n.empty()) return n2;
    if (n2.n.empty()) return n1;

    N n;
    n.ChangeInternalBase(n1.b);
    n.n.reserve(max(n1.NumDigits(), n2.NumDigits()));
    D carry = 0;

    if (n1.m && n2.m)
        n.m = true;

    size_t j = 0;
    for (size_t i = 0; i < n1.NumDigits() || i < n2.NumDigits() ; i++)
        {
        j = i;
        D sum = n1[i] + n2[i] + carry;
        carry = 0;
        if (sum >= n1.b)
            {
            carry = 1;
            sum -= n1.b;
            }
        n.n.push_back(sum);
        }
    n.n.push_back(carry);
    n.RemoveZeroes();
    }

Single Sub

Subtraction eventually calls function w_subx():

static N w_subx(const N& n1, const N& n2)
    {
    if (n1.m != n2.m)
        return w_add(n1, n2.negative());

    if (n1.absolute() < n2.absolute())
        return w_subx(n2, n1).negative();

    if (n2.IsZero())
        return n1;
    if (n1.IsZero())
        return n2.negative();

    N n;
    n.ChangeInternalBase(n1.b);
    n.m = n1.m;
    int carry = 0;

    for (size_t i = 0; i < n1.NumDigits() || i < n2.NumDigits(); i++)
        {
        signed long long sum = n1[i] - n2[i] + carry;
        carry = 0;
        if (sum < 0)
            {
            sum = n1.b + sum;
            carry = -1;
            }
        n.n.push_back(sum);
        }
    n.n.push_back(carry);
    n.RemoveZeroes();
    return n;

    }

Logical operations

w_logical() is called for operations like AND, OR and XOR. It always converts the numbers to binary. Note that operator ^ by default executes w_pow(), unless both operands are in binary format.

    static N w_logical(const N& n1, const N& n2,int x)
        {
        if (n1.b != 2)
            return w_logical(n1.tobase(2),n2,x);
        if (n2.b != 2)
            return w_logical(n1, n2.tobase(2),x);

        N n;
        n.ChangeInternalBase(2);
        n.n.reserve(max(n1.NumDigits(), n2.NumDigits()));

        for (size_t i = 0; i < n1.NumDigits() || i < n2.NumDigits(); i++)
            {
            D sum = 0;
            if (x == 0) sum = n1[i] & n2[i];
            if (x == 1) sum = n1[i] | n2[i];
            if (x == 2) sum = n1[i] ^ n2[i];
            n.n.push_back(sum);
            }
        n.RemoveZeroes();
        return n;
        }


N a("10",2);
N b("11",2);
auto r = a | b; // r = 5

In addition, cshl(), cshr(), shl2(), shl2(), negative() and absolute() perform the relevant logical operations. shl() and shr() perform self-shl/shr, where the other functions generate a new N.

Multiplication

The multiplication function creates a vector<N> with all multiplications, then adds them together:

static N w_mul(const N& n1, const N& n2)
    {
    if (n1.b != n2.b)
        return w_mul(n1, n2.tobase(n1.b));
    N n;
    n.ChangeInternalBase(n1.b);
    vector<N> addiz;
    for (size_t i = 0; i < n1.n.size(); i++)
        {
        D d1 = n1.n[i];
        N addi;
        addi.n.reserve(i + n2.n.size());
        for (size_t j = 0; j < i; j++)
            addi.n.push_back(0);
        D carry = 0;
        for (size_t y = 0; y < n2.n.size(); y++)
            {
            D d2 = n2.n[y];
            D dm = (d1*d2) + carry;
            carry = 0;
            carry = dm / n1.b;
            dm %= n1.b;
            addi.n.push_back(dm);
            }
        addi.n.push_back(carry);
        addi.RemoveZeroes();
        addiz.push_back(addi);
        }
    for (auto& a : addiz)
        n += a;
    if (n1.m != n2.m)
        n.m = true;
    return n;
    }

Division

The division function (used by division stuff and operators like /,%,/=,%= etc) returns a tuple of the result and the remainder:

static tuple<N, N> w_div(const N& n1, const N& n2,bool NoChangeBase = false)
    {
    if (n1.b != n2.b && NoChangeBase == false)
        return w_div(n1.b, n2.tobase(n1.b));
    if (n2 > n1)
        {
        N res = n1;
        return std::make_tuple<N, N>(0LL, std::forward<N>(res));
        }
    if (n2 == n1)
        return std::make_tuple<N, N>(1LL, 0LL);

    N rem = n1;
    N res;
    res.ChangeInternalBase(n1.b);

    for (;;)
        {
        auto nd2 = n2.NumDigits();
        auto upper = rem.upperpart(nd2);
        if (upper < n2)
            {
            nd2++;
            upper = rem.upperpart(nd2);
            if (upper < n2)
                {
                // End...
                return std::make_tuple<N, N>(forward<N>(res), forward<N>(rem));
                }
            }

        // This needs optimization
        unsigned long long js = 9;
        N m1;
        for (; js >= 1; js--)
            {
            m1 = w_mul(n2, js);
            if (m1 < upper)
                break;
            }

        res.n.insert(res.n.begin(),js);
        upper -= m1;
        upper.shl(rem.NumDigits() - nd2);
        upper += rem.lowerpart(rem.NumDigits() - nd2);
        rem = upper;
        }
    }

Power

Nothing much complicated.

    static N w_pow(const N& n1, const N& n2)
        {
        N z = n1;
        if (n2 == 0ll)
            return 1ll;
        if (n2 == 1ll)
            return n1;
        if (n1 == 1ll)
            return 1ll;

        for (N j = 1ll; j < n2; ++j)
            z *= n1;
        return z;
        }

N a("-5");
auto a2 = N::w_pow(a, 3LL); // a2 = -125

Multicore Operations

A multiple core adding accepts a vector of Ns, and adds them in pairs, then calls itself with the result until there are only 2 numbers to be added:

static N w_add2(vector<N>& n, tpoollib::tpool<>& pool)
    {
    if (n.size() == 0)
        return 0LL;
    if (n.size() == 1)
        return n[0];
    if (n.size() == 2)
        {
        N res = n[0];
        res += n[1];
        return res;
        }
    struct Z
        {
        N* n1;
        N* n2;
        N* res;
        };
    vector<Z> z(n.size());
    vector<N> res(n.size() / 2);

    for (size_t i = 0; i < n.size(); i += 2)
        {
        if (i == (n.size() - 1))
            break; // odd number of additions

        auto a = [](PTP_CALLBACK_INSTANCE, PVOID j, PTP_WORK)
            {
            Z* z = (Z*)j;
            *z->res = w_add(*z->n1,*z->n2);
            };
        Z& zz = z[i];
        zz.n1 = &n[i];
        zz.n2 = &n[i + 1];
        zz.res = &res[i / 2];

        auto wo = pool.CreateItem<PTP_WORK, PTP_WORK_CALLBACK>(a, (PVOID)&zz);
        pool.RunItem(wo);
        }
    pool.Join();
    if (n.size() % 2)
        res.push_back(n[n.size() - 1]);
    return w_add2(res,pool);
    }

It uses a temporary struct Z to pass it to the threadpool function.

A multiple core multiplication creates a vector of multiplications, then adds them with w_add2():

static N w_mul2(const N& n1, const N& n2, tpoollib::tpool<>& pool)
    {
    size_t muls = n1.NumDigits() * n2.NumDigits();
    vector<N> a;
    a.reserve(muls);
    for (size_t i = 0; i < n1.NumDigits(); i++)
        {
        for (size_t ii = 0; ii < n2.NumDigits(); ii++)
            {
            N rr;
            D d1 = n1[i];
            D d2 = n2[ii];
            unsigned long long r = d1 * d2;
            rr = r;
            rr.shl(ii + i);
            a.push_back(rr);
            }
        }
    return w_add2(a,pool);
    }

This can take optimizations, for example to find the sums within a multicore for loop itself.

A multiple core power just uses w_mul2():

static N w_pow2(const N& n1, const N& n2,tpoollib::tpool<>& pool)
    {
    N z = n1;
    if (n2 == 0ll)
        return 1ll;
    if (n2 == 1ll)
        return n1;
    if (n1 == 1ll)
        return 1ll;

    for (N j = 1ll; j < n2; ++j)
        z = w_mul2(z, n1,pool);
    return z;
    }

This can take optimizations as well.

The code

All included in one header file, feel free to enhance!

History

11 - 8 - 2016 : First release

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Michael Chourdakis
Engineer
Greece Greece
I'm working in C++, PHP , Java, Windows, iOS and Android.

I 've a PhD in Digital Signal Processing and I specialize in Pro Audio applications.

My home page: http://www.michaelchourdakis.com

You may also be interested in...

Pro

Comments and Discussions

 
QuestionFormat Pin
Nelek10-Aug-16 23:54
protectorNelek10-Aug-16 23:54 
AnswerRe: Format Pin
Michael Chourdakis11-Aug-16 0:41
mvpMichael Chourdakis11-Aug-16 0:41 
GeneralRe: Format Pin
Nelek11-Aug-16 0:51
protectorNelek11-Aug-16 0:51 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

Permalink | Advertise | Privacy | Cookies | Terms of Use | Mobile
Web04-2016 | 2.8.180712.1 | Last Updated 11 Aug 2016
Article Copyright 2016 by Michael Chourdakis
Everything else Copyright © CodeProject, 1999-2018
Layout: fixed | fluid