Big Integers Study

Intro

Prolog still stands apart in languages landscape, and SWI-Prolog has grown big, features rich, since Jan is boosting it in practical SW engineering field.

Since I’d like to use `library(clpfd)` for some simple tasks, and it depends on unlimited integers arithmetic, I need efficient unlimited integer arithmetic.

After the announce that the licence transition, from GPL to BSD, would require abandoning the awesome GMP, I started to evaluate if alternatives in my reach were available.

Big integers

A small benchmark, of little practical interest, but requiring unlimited precision integer arithmetic, is the evaluation of `choose(N, K)`, best known as binomial coefficient, with large N.

In case of `choose(50000, 50)`, the result is this big number `284958500315333290867708487072990268397101930544468658476216100935982755508148971449700622210078705183923286686402942943816349032142836981589618876226813174803825580124000`.
I will show it abbreviated in output of the following tests, like `28495...24000`.

The benchmark originated from Jason, as a study about Python’ native foreign interfaces.

```def factorial(n, stop=0):
o = 1
while n > stop:
o *= n
n -= 1
return o

def choose(n, k):
return factorial(n, stop=k) / factorial(n - k)

if __name__ == '__main__':
print choose(50000, 50)```

To run with Python 3.4, a small change is required – just added parenthesis to print

```def factorial(n, stop=0):
o = 1
while n > stop:
o *= n
n -= 1
return o

def choose(n, k):
return factorial(n, stop=k) // factorial(n - k)

if __name__ == '__main__':
print(choose(50000, 50))```

Building C++

I do prefer to use QtCreator when developing C++, so there are .pro files you can use to see dependencies and build details.

Building can be simple as opening the .pro file, configuring the project for release, and build.

After done, right click on the .pro file in Projects, and select `'Open Terminal Here'`.

Then paste the bash command line – for instance `\$ time ../build-choose-gmp-Desktop_Qt_5_7_0_GCC_64bit-Release/choose-gmp`

Timings

Here is a table summarizing my tests:

source interface performance runtime
SWI-prolog GMP 6.2 sec
```                    \$ time swipl -O -g 'choose,halt' choose.pl
28495...24000
real        0m6.272s
user        0m6.189s
sys 0m0.064s
```
python 2.7 native 6 sec
```                    \$ time python choose.py
28495...24000
real        0m6.010s
user        0m5.886s
sys 0m0.112s
```
python 3.4 native 6.6 sec
```                        \$ time python3 choose-3.py
28495...24000
real    0m6.691s
user    0m6.625s
sys     0m0.056s
```
node (js V8) OpenSSL BN 15.5 sec
```                        \$ npm install bignum
\$ time node choose.js
<BigNum 28495...24000>
real    0m15.501s
user    0m12.608s
sys     0m2.333s
```
C++ GMP 2.2 sec
```                    \$ time ../build-choose-gmp-Desktop_Qt_5_7_0_GCC_64bit-Release/choose-gmp
28495...24000
real        0m2.231s
user        0m2.095s
sys 0m0.003s
```
C++ OpenSSL BN 2.5 sec
```                    \$ time ../build-choose-openssl-Desktop_Qt_5_7_0_GCC_64bit-Release/choose-openssl
28495...24000
real        0m2.534s
user        0m2.504s
sys 0m0.028s
```
C++ KNST – native 26.9 sec
```                    \$ time ../build-choose-knst-Desktop_Qt_5_7_0_GCC_64bit-Release/choose-knst
28495...24000
real        0m26.933s
user        0m26.711s
sys 0m0.200s
```

SWI-Prolog

```/** <module> choose
*
*  benchmark BigNum GMP implementations
*  from http://jasonstitt.com/c-extension-n-choose-k
*/

:- module(choose, [choose/0]).

choose :-
choose(50000, 50, C),
writeln(C).

choose(N, K, C) :-
factorial(N, K, A),
M is N - K,
factorial(M, B),
C is A / B.

factorial(N, F) :-
factorial(N, 1, F).

factorial(N, Stop, F) :-
factorial(N, Stop, N, F).

factorial(N, Stop, A, F) :-
N > Stop + 1,
M is N - 1,
A1 is A * M,
!,  % without the cut, I get a stack overflow with -O flag
factorial(M, Stop, A1, F).
factorial(N, Stop, A, A) :-
N =< Stop + 1.

/* this is about 30% slower ~ 8.5 sec
factorial(N, Stop, A, F) :-
(   N > Stop + 1
->  M is N - 1,
A1 is A * M,
factorial(M, Stop, A1, F)
;   A = F
).
*/```

SWI-Prolog code requires absolutely a tail recursive loop, otherwise the performance is very slow.

I started with a naive translation, and the execution time was about 95 sec – due to inability to run GC and a lot of big numbers allocation, I guess.

A first attempt to optimize (i.e. `\$ time swipl -O -g 'choose,halt' choose.pl`) caused a stack overflow.

Javascript

```var bignum = require('bignum');

function factorial(n, stop) {
var o = bignum(1)
while (n > stop) {
o = o.mul(n)
n -= 1
}
return o
}
//console.log(factorial(6, 0))

function choose(n, k) {
return factorial(n, k).div(factorial(n - k, 0))
}
console.log(choose(50000, 50))```

So far, I’ve been unable to install GMP support in Node. `npm install bigint` fails compiling the code,
I think it refers to this repo, indeed rather old.

Also bignum installation is weird, it just works for a session.
When I restart the system, I must reinstall. Anyway, it’s used in `choose.js`

C++

Knst refers to this github project.

```#include <iostream>
#include <BN.h>

using namespace std;

BN factorial(BN n, BN stop = BN(0)) {
BN one(1), o = one;
while (n > stop) {
o = o * n;
n = n - one;
}
return o;
}

BN choose(BN n, BN k) {
return factorial(n, k) / factorial(n - k);
}

int main() {
cout << to_string(choose(BN(50000), BN(50))) << endl;
}```

OpenSSL BN is practically on par with GMP, but the interface is notably lower level, and requires optimized calls.

```#include <iostream>
#include <stdexcept>

// https://www.openssl.org/docs/man1.0.1/crypto/bn.html
#include <openssl/bn.h>

using namespace std;
inline void err(string msg) { throw runtime_error(msg); }

BIGNUM* factorial(BN_CTX* ctx, int n, int stop) {

BIGNUM *o = 0;
if (!BN_dec2bn(&o, "1"))
err("o: cannot BN_dec2bn");

BIGNUM *N = 0;
if (!BN_dec2bn(&N, to_string(n).c_str()))
err("N: cannot BN_dec2bn");

BIGNUM *Stop = 0;
if (!BN_dec2bn(&Stop, to_string(stop).c_str()))
err("Stop: cannot BN_dec2bn");

while (BN_cmp(N, Stop) > 0) {
BN_mul(o, o, N, ctx);
BN_sub_word(N, 1);
}
BN_free(N);
BN_free(Stop);

return o;
}

BIGNUM* choose(BN_CTX *ctx, int n, int k) {
BIGNUM* a = factorial(ctx, n, k), *b = factorial(ctx, n - k, 0), *c = BN_new();
if (!BN_div(c, 0, a, b, ctx))
err("cannot BN_div");
BN_free(a);
BN_free(b);
return c;
}

int main_old() {
try {
BN_CTX *ctx = BN_CTX_new();
if (!ctx)
err("cannot BN_CTX_new");

BIGNUM *n_k = choose(ctx, 50000, 50);
char *s = BN_bn2dec(n_k);
BN_free(n_k);

cout << s << endl;
OPENSSL_free(s);

BN_CTX_free(ctx);
return 0;
}
catch(exception &e) {
cout << "exception " << e.what() << endl;
return 1;
}
}

BIGNUM* factorial(int n, int stop) {

BIGNUM *o = 0;
if (!BN_dec2bn(&o, "1"))
err("o: cannot BN_dec2bn");

BIGNUM *N = 0;
if (!BN_dec2bn(&N, to_string(n).c_str()))
err("N: cannot BN_dec2bn");

BIGNUM *Stop = 0;
if (!BN_dec2bn(&Stop, to_string(stop).c_str()))
err("Stop: cannot BN_dec2bn");

while (n > stop) {
BN_mul_word(o, n);
--n;
}

return o;
}
BIGNUM* choose(int n, int k) {
BN_CTX *ctx = BN_CTX_new();
if (!ctx)
err("cannot BN_CTX_new");

BIGNUM* a = factorial(n, k), *b = factorial(n - k, 0), *c = BN_new();
if (!BN_div(c, nullptr, a, b, ctx))
err("cannot BN_div");
BN_free(a);
BN_free(b);
BN_CTX_free(ctx);
return c;
}

int main() {
try {
BIGNUM *n_k = choose(50000, 50);
char *s = BN_bn2dec(n_k);
BN_free(n_k);

cout << s << endl;
OPENSSL_free(s);
return 0;
}
catch(exception &e) {
cout << "exception " << e.what() << endl;
return 1;
}
}```

You can see the not optimized calls following `main_old`, the performance was worse – about 4.5 secs.

In GMP, the code is really simple, thanks to a good C++ interface, but I cannot say if it could be further optimized, maybe applying the same trick used in OpenSSL BN.

```#include <iostream>
#include <gmpxx.h>

using namespace std;

mpz_class factorial(int n, int stop = 0) {
mpz_class o = 1;
while (n > stop) {
o *= n;
n -= 1;
}
return o;
}

mpz_class choose(int n, int k) {
return factorial(n, k) / factorial(n - k);
}

int main() {
cout << choose(50000, 50);
}```