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.

[czurnieden]

I did a frexp implementation the long way and without libm. Please run to test. Needs libm but only for the tests because I use libm's frexp to check against. The problem is the check for NaN and Infinity without libm. Some compilers might protest against the way I did it.

Code attached.
Or not?
Well ... sigh ... sorry for that mess.

double frexp_intern(double x, int *eptr)
{
   int sign, exponent;
   int i;

   /*
    * The exponent of an IEEE-754 double (binary64) is an 11-bit large integer
    */
   double ap_2[11] = {
      2.0000000000000000000000000000000000000,
      4.0000000000000000000000000000000000000,
      16.000000000000000000000000000000000000,
      256.00000000000000000000000000000000000,
      65536.000000000000000000000000000000000,
      4294967296.0000000000000000000000000000,
      18446744073709551616.000000000000000000,
      3.4028236692093846346337460743176821146e38,
      1.1579208923731619542357098500868790785e77,
      1.3407807929942597099574024998205846128e154,
      1.7976931348623157e308 // DBL_MAX
   };
   double ap_half[11] = {
      0.50000000000000000000000000000000000000,
      0.25000000000000000000000000000000000000,
      0.062500000000000000000000000000000000000,
      0.0039062500000000000000000000000000000000,
      1.5258789062500000000000000000000000000e-5,
      2.3283064365386962890625000000000000000e-10,
      5.4210108624275221700372640043497085571e-20,
      2.9387358770557187699218413430556141946e-39,
      8.6361685550944446253863518628003995711e-78,
      7.4583407312002067432909653154629338374e-155,
      5.5626846462680034577255817933310101606e-309 // < DBL_MIN
   };


   // TODO: not every compiler might eat this check for Inf and NaN
   // GCC-4.8.4  does
   // TCC 0.9.25 does
   // clang 3.4-1ubuntu3 (based on LLVM 3.4) does
   if (x == 1.0 / 0.0) {
      *eptr = 0;
      return x;
   }
   if (x == 0.0 / 0.0) {
      *eptr = 0;
      return x;
   }

   if (x == 0.0) {
      *eptr = 0;
      return x;
   }

   exponent = 0.0;
   /*
    * Easier to work with positive values
    */
   if (x < 0) {
      x = -x;
      sign = 1;
   } else {
      sign = 0;
   }

   if (x >= 1.0) {
      /*
       * Big steps
       */
      for (i = 0; x >= ap_2[i]; i++) {
         exponent += (1 << i);
         x *= ap_half[i];
      }
      /*
       * Small steps
       */
      if (x < 0.5) {
         while (x < 0.5) {
            x *= 2.0;
            exponent--;
         }
      } else {
         while (x > 1.0) {
            x /= 2.0;
            exponent++;
         }
      }
   } else {
      /*
       * Same as above, but in the opposite direction
       */
      for (i = 0; x < ap_half[i]; i++) {
         exponent -= (1 << i);
         x *= ap_2[i];
      }
      if (x < 0.5) {
         while (x < 0.5) {
            x *= 2.0;
            exponent--;
         }
      } else {
         while (x > 1.0) {
            x /= 2.0;
            exponent++;
         }
      }
   }
   if (sign) {
      x = -x;
   }
   *eptr = exponent;
   return x;
}


#include <stdlib.h>
#include <stdio.h>
#include <time.h>
// for the internal frexp() only
#include <math.h>
int main()
{
   int i, e1, e2;
   double d1, d2, t, f1, f2;

   srand((unsigned int) time(NULL));

   t = 0.0;
   f1 = frexp_intern(t, &e1);
   f2 = frexp(t, &e2);
   if (e1 != e2) {
      printf("0 EXPO FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      exit(EXIT_FAILURE);
   }
   if (f1 != f2) {
      printf("0 MANT NAN:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      exit(EXIT_FAILURE);
   }
   // NaN
   // Will fail MANT-test because NaN != NaN
   t = 0.0 / 0.0;
   f1 = frexp_intern(t, &e1);
   f2 = frexp(t, &e2);
   if (e1 != e2) {
      printf("1 EXPO FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      exit(EXIT_FAILURE);
   }
   if (f1 != f2) {
      printf("1 MANT NAN:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      //exit(EXIT_FAILURE);
   }
   // Infinity
   t = 1e308;
   t *= t;
   f1 = frexp_intern(t, &e1);
   f2 = frexp(t, &e2);
   if (e1 != e2) {
      printf("2 EXPO FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      exit(EXIT_FAILURE);
   }
   if (f1 != f2) {
      printf("2 MANT FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      exit(EXIT_FAILURE);
   }
   // -Infinity
   t = -1e308;
   t *= t;
   f1 = frexp_intern(t, &e1);
   f2 = frexp(t, &e2);
   if (e1 != e2) {
      printf("2 EXPO FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      exit(EXIT_FAILURE);
   }
   if (f1 != f2) {
      printf("2 MANT FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      exit(EXIT_FAILURE);
   }
   // subnormal (denormal) numbers
   t = 4.e-324;
   f1 = frexp_intern(t, &e1);
   f2 = frexp(t, &e2);
   if (e1 != e2) {
      printf("2 EXPO FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      exit(EXIT_FAILURE);
   }
   if (f1 != f2) {
      printf("2 MANT FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
             e2, f1, f2);
      exit(EXIT_FAILURE);
   }

   for (i = 0; i < 1000000; i++) {
      d1 = (double) rand();
      d2 = (double) rand();

      t = d1;
      f1 = frexp_intern(t, &e1);
      f2 = frexp(t, &e2);
      if (e1 != e2) {
         printf("3 EXPO FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
                e2, f1, f2);
         exit(EXIT_FAILURE);
      }
      if (f1 != f2) {
         printf("3 MANT FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
                e2, f1, f2);
         exit(EXIT_FAILURE);
      }

      t = 1 / d1;
      f1 = frexp_intern(t, &e1);
      f2 = frexp(t, &e2);
      if (e1 != e2) {
         printf("4 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
                e2, f1, f2);
         exit(EXIT_FAILURE);
      }
      if (f1 != f2) {
         printf("4 MANT FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
                e2, f1, f2);
         exit(EXIT_FAILURE);
      }

      t = d2 / d1;
      f1 = frexp_intern(t, &e1);
      f2 = frexp(t, &e2);
      if (e1 != e2) {
         printf("5 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
                e2, f1, f2);
         exit(EXIT_FAILURE);
      }
      if (f1 != f2) {
         printf("5 MANT FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
                e2, f1, f2);
         exit(EXIT_FAILURE);
      }

      t = d1 / d2;
      f1 = frexp_intern(t, &e1);
      f2 = frexp(t, &e2);
      if (e1 != e2) {
         printf("6 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
                e2, f1, f2);
         exit(EXIT_FAILURE);
      }
      if (f1 != f2) {
         printf("6 MANT FAIL:t = %g,  e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
                e2, f1, f2);
         exit(EXIT_FAILURE);
      }
   }
   exit(EXIT_SUCCESS);
}

[tomstdenis]

This is not really a 0.43 candidate so let's not merge it into develop at all for now. As for the general comment ... It'd be nice to not have any float code "standard" in LTM. I'd be ok with bundling it with the package but not with making it a core API component.

To be fair though embedded should be using TFM unless they have some specific requirement from LTM.

[sjaeckel]

I also don't mind delaying that one

[minad]

@sjaeckel I guess this can be closed.

[sjaeckel]

closed via #123