The biggest native integer libtommath allowed to set directly seems to be[1] an `unsigned long`

in the function `mp_set_int`

. The biggest native integer used, on the other side, is hidden behind `mp_word`

which is the type able to hold twice the size of `mp_digit`

and can be larger than an `unsigned long`

.

I need for my calculator some ways to work with native numbers without much ado where *ado* means a lot of conditionals, preprocessor directives, complicated data structures and all that mess. One of the ways to avoid it is to use the digits of the large integers directly if the large integer has only one. An example? Ok, an example.

Take the partial harmonic series, for example

If you calculate it with the help of the binary-splitting algorithm, especially than, a lot of the numbers involved are in the range of native integers and hold only on digit of the big numbers. The initialization of the big numbers in libtommath set them to 8 digits at least (the responsible variable is `MP_PREC`

in `tommath.h`

) and consumes costly heap memory to do so.

Fredrik Johansson proposed in a blogpost to postpone the reducing of the fraction to the very end. It is not much but it is something so let’s follow his advice and do so using my rational library (as mostly: without any error checking for less code-cluttering).

static mp_rat h_temp; mp_rat *_harmonics(unsigned long a, unsigned long b) { unsigned long m; mp_rat ta, tb; mp_int p, q, r, s; mp_word ps, qr; int e; if (b - a == 1) { mpq_set_int(&h_temp, (long) 1, (long) a); return &h_temp; } m = (a + b) >> 1; mpq_init_multi(&ta, &tb, NULL); // This is not necessarily necessary mp_init_multi(&p, &q, &r, &s, NULL); mpq_exch(_harmonics(a, m), &ta); mpq_exch(_harmonics(m, b), &tb); mp_exch(&ta.numerator, &p); mp_exch(&ta.denominator, &q); mp_exch(&tb.numerator, &r); mp_exch(&tb.denominator, &s); if ((&p)->used == 1 && (&s)->used == 1) { ps = (&p)->dp[0] * (mp_word) (&s)->dp[0]; mp_set_word(&ta.numerator, ps); } else { mp_mul(&p, &s, &ta.numerator); } if ((&q)->used == 1 && (&r)->used == 1) { qr = (&q)->dp[0] * (mp_word) (&r)->dp[0]; mp_set_word(&tb.numerator, qr); } else { mp_mul(&q, &r, &tb.numerator); } mp_add(&ta.numerator, &tb.numerator, &h_temp.numerator); mp_mul(&q, &s, &h_temp.denominator); mp_clear_multi(&p, &q, &r, &s, NULL); mpq_clear_multi(&ta, &tb, NULL); return &h_temp; } int harmonics(unsigned long n, mp_rat * c) { mpq_init(&h_temp); mpq_exch(_harmonics(1, n + 1), c); mpq_reduce(c); mpq_clear(&h_temp); return 1; }

The library `libtommath`

is not very friendly if used in such a way but I cannot find out a better way to implement the binary splitting algorithm. This implementation of the partial harmonic series for example, is much slower than Fredriks implementation with `gmpy`

(but faster than the native `Python`

one, at least 😉 ). It takes about 0.67 seconds for but already 193.61 seconds—yes, over 3 minutes!— for . That is definitely too much.

Funny thing: the normal algorithm is much faster, just 40 seconds for but also slower for smaller values, like about 1.09 seconds for whith the cut-off point at about on my machine. And it is *asymptotically* slower, I measured some points in between to find out. It is really time to implement fast multiplication in full in `libtommath`

.

Some of the problems with the first algorithm might have their reason in the nonexistent balancing of the multiplicands. There is a balancing function in the pull-queue but it seems to be a port from Ruby which makes it impossible to accept because of Ruby’s license (given that the submitter is not the original programmer of the Ruby code which I haven’t checked.)

static mp_rat h_temp; mp_rat * _harmonics(unsigned long a,unsigned long b){ unsigned long m; mp_rat ta,tb; if(b-a==1){ mpq_set_int(&h_temp,(long)1,(long)a); return &h_temp; } m = (a+b)>>1; mpq_init_multi(&ta,&tb,NULL); mpq_exch(_harmonics(a,m),&ta); mpq_exch(_harmonics(m,b),&tb); mpq_add(&ta,&tb,&h_temp); mpq_clear_multi(&ta,&tb,NULL); return &h_temp; } // use the same caller as above

However, it was just meant to be used as an example for `mp_set_word`

which I’ll present here:

#if (MP_PREC > 1) int mp_set_word(mp_int *c,mp_word w){ mp_zero(c); if(w == 0){ return MP_OKAY; } do{ c->dp[c->used++] = (mp_digit)w&MP_MASK; }while( (w >>= DIGIT_BIT) > 0 && c->used < MP_PREC); if( w != 0 ){ return MP_VAL; } return MP_OKAY; } #else #warning variable "MP_PREC" should be at least 2 #endif

The preprocessor mess is necessary even if the constant `MP_PREC`

should be at least 8 (eight) but, as I can tell you from some very bad experience, one never knows.

[1] The reason for the stylistically awkward subjunctive is: I mean it. I really could have overseen an already implemented function doing exactly what I wanted in the first place and hence is not a case of NIH. This time 😉