Feature request: conversion from mp_int to float or double precision

[moritz]

It would be nice if I could convert a mp_int into a float or double, preserving the magnitude and some of the first digits, though of course the general case will be lossy. If the mp_int is too large, it could return Inf or -Inf.

We use mp_int for storing big integers in a Perl 6 compiler, and it we'll need to offer such functionality to our uses eventually.

[moritz]

Here is a first shot at such a function. I'm neither a good C hacker nor very well versed with libtommath, so beware of bugs -- the few cases I tested seemed to work though.

#include "tommath.h"
#include <math.h>

double mp_get_double(mp_int *a) {
    double d    = 0.0;
    double sign = SIGN(a) == MP_NEG ? -1.0 : 1.0;
    if (USED(a) == 0)
        return d;
    if (USED(a) == 1)
        return sign * (double) mp_get_int(a);

    int i;
    for (i = USED(a) - 1; DIGIT(a, i) == 0 && i > 0; i--) {
        /* do nothing */
    }
    d = (double) DIGIT(a, i);
    i--;
    if (i == -1) {
        return sign * d;
    }
    d *= pow(2.0, DIGIT_BIT);
    d += (double) DIGIT(a, i);

    d *= pow(2.0, DIGIT_BIT * i);
    return sign * d;
}

[gerdr]

Depending on the value of DIGIT_BIT, moritz' code loses precision. Barring off-by-one errors, the following code should work in general:

#include <tommath.h>
#include <math.h>

#define PRECISION 53

double mp_get_double(mp_int *a)
{
    static const int NEED_DIGITS = (PRECISION + 2 * DIGIT_BIT - 2) / DIGIT_BIT;
    static const double DIGIT_MULTI = (mp_digit)1 << DIGIT_BIT;

    int i, limit;
    double d = 0.0;

    mp_clamp(a);
    i = USED(a);
    limit = i <= NEED_DIGITS ? 0 : i - NEED_DIGITS;

    while (i-- > limit) {
        d += DIGIT(a, i);
        d *= DIGIT_MULTI;
    }

    if(SIGN(a) == MP_NEG)
        d *= -1.0;

    d *= pow(2.0, i * DIGIT_BIT);
    return d;
}

[hikari-no-yume]

I'd find this useful for my patch to add native bigints to PHP's Zend Engine.

[hikari-no-yume]

The reverse, exporting an mp_int to a double would also be handy, to save me from implementing it myself. Especially since I don't want a dependency on LibTomMath's internals... abstraction is my friend.

[czurnieden]

The reverse [...] would also be handy,

That would be mp_set_double?

#include <math.h>
#include <tommath.h>
int mp_set_double(mp_int * c, double d){
  int exp, res, sign;
  double frac;

  sign = (d < 0)?MP_NEG:MP_ZPOS;
  frac = frexp(abs(d), &exp);

  if(exp <= 0){
    return MP_OKAY;
  }
  if(exp == 1 && frac < .5){
    return MP_OKAY;
  }

  mp_zero(c);

  if(frac == 0){
    c->sign = sign;
    return MP_OKAY;
  }

  while(exp-- >= 0){
    frac *= 2.0;
    if(frac >= 1.0){
      if ((res = mp_add_d(c, 1, c)) != MP_OKAY) {
        return res;
      }
      frac -= 1.0;
    }
    if(exp > 0){
      if ((res = mp_mul_2d(c, 1, c)) != MP_OKAY) { 
        return res;
      }
    }
  }
  c->sign = sign;                
  return MP_OKAY;
}

Try it out with e.g.:

#include <stdio.h>
#include <stdlib.h>
int main(int argc, char **argv){
  char *eptr, retstr;
  double v;
  mp_int r;

  v = strtod(argv[1],&eptr);
  mp_init(&r);
  mp_set_double(&r,v);
  mp_toradix(&r, &retstr, 10);

  printf("%f = %s\n",v,&retstr);

  exit(EXIT_SUCCESS);
}

The problem might be the rounding mode — it doesn't respect it (doesn't even read it ;-) ) it always rounds to zero (truncates). That's the reason the code is listed here instead of being committed regularly: any code handling floating point values must respect rounding modes.

PS: did the FFT work?

[czurnieden]

It seems as if nearbyint() would do here, so

No, it doesn't. Silly me! Sorry.

[czurnieden]

Nothing in TV, so here wethey are ;-)

https://github.com/czurnieden/libtommath/blob/master/bn_mp_get_double.c
https://github.com/czurnieden/libtommath/blob/master/bn_mp_set_double.c

Please check.

[hikari-no-yume]

@czurnieden You do a == HUGE_VAL check. On systems where it's available (C99), perhaps better to do == INFINITY?

[czurnieden]

You do a == HUGE_VAL check. On systems where it's available (C99), perhaps better to do == INFINITY?

I just bracketed all c99 stuff out before I came back here ;-)
But one more or less ...

BTW: does the FFT work? I cannot repair something when I do not know that it is broken (a variation of "works for me") save how it is broken.

[sjaeckel]

I just had an idea how we could include that functionality in ltm without adding a dependency to libm ...

We provide a macro with the implementation, so someone who wants to use it puts that macro in his code and then he can use the implementation...

I only proposed this so someone else can propose something better ;-)
I revived this because I think that it's useful and we should just do it somehow (but without adding a dependency to libm)

[czurnieden]

mp_get_double (convert mp_int to double) does not need the libm but mp_set_double (double to mp_int) does use frexp (isnan, isinf, and nearbyint if available). It is not that complicated to write a frexp if you know the endinaness in advance, otherwise it is quite tedious (vid. e.g.: the implementation in the coreutils in coreutils-8.24/lib/frexp.c ), so it's a rewrite without frexp et al.?

[czurnieden]

Done (for big and little endian only), to be found under the same adress. Pleeeeease check before use.

[hikari-no-yume]

What's wrong with depending on libm? math.h is part of the C standard library.

[czurnieden]

@TazeTSchnitzel LibTomMath gets used quite widely, especially on embedded stuff. Those PICs do not necessarily have a FPU--some don't even support floats, they have to do it all with integers. So it might be already slow to use floating point math. On top of that: they don't have a lot of memory available, loading a library for one or two of its functions is a large overhead that might not be acceptable or it is even impossible to load in the first place because of lack of memory.

And while I'm writing that I get increasingly uncomfortable with the bit-pushing in my implementation of frexp(3). I think I will do it the long way instead. But not today ;-)

[hikari-no-yume]

I can understand the embedded use-case. Hmm.

Perhaps we could have an option to not build the floating-point parts of the library. Or maybe we could use libm if available... oh I don't know.