Adventures of a Programmer: Parser Writing Peril XIII

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

{\displaystyle {  \mathrm{H}_n = \sum_{k=1}^n \frac{1}{k}  }}

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_exch(_harmonics(1, n + 1), c);
  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 \mathop{\mathrm{H}}_{10\,000} but already 193.61 seconds—yes, over 3 minutes!— for \mathop{\mathrm{H}}_{100\,000} . That is definitely too much.
Funny thing: the normal algorithm is much faster, just 40 seconds for \mathop{\mathrm{H}}_{100\,000} but also slower for smaller values, like about 1.09 seconds for \mathop{\mathrm{H}}_{10\,000} whith the cut-off point at about \mathop{\mathrm{H}}_{21\,000} 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;
   return &h_temp;
  m = (a+b)>>1;


  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){
  if(w == 0){
    return MP_OKAY;
    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;
#warning variable "MP_PREC" should be at least 2

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 😉


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s