add mp_get_double, mp_set_double

[minad]

basic tests in demo.c

[minad]

@sjaeckel as discussed in #120

[minad]

this resolves #32 and #3

[czurnieden]

Are you sure that mp_set _double is endianess independent?

I think that NaN and +/-infinity as an input to mp_set _double should result in an error instead of MP_OKAY, especially with the fact in mind that Libtommath has no "error-number" like NaN.

The union trick is AFAIK a problem in C++, but I am not sure.

The double in question might not be an IEEE-754 binary64 and sizeof(double) might not even be 8.
(check with #ifdef __STDC_IEC_559__ for IEEE-754 compliance)

isinf[inite]() and isnan() are in newer standards (check with #if __STDC_VERSION__ >= 199901L)

A negative de-biased exponent means the double is smaller than unity already (the mantissa m is in the range 1 <= m < 2 assuming IEEE-754 compliance). It is an extra check but on the other side saves an unnecessary call to mp_div_2d.
You could also check the biased exponent if it is zero (+-0 or denormal number).
If the biased exponent is 2^w-1 it stands for NaN and +/- infinity. (Tests like e.g.: d = 1.0/0.0 can be problematic and a test like d != d can be optimized away, there are still a lot of "stupid" compilers out there.)

I think the error checking (overflow only in this case, that is +/- infinity) should be done inside mp_get_double and it should return that error (MP_RANGE) instead of the double (int mp_get_double(const mp_int *a, double *d)). I think a simple test if the highest bit in the mp_int is at >=1024 will do.

There is also this "hint": #3 (comment)
So it might be a good idea to implement some methods to include the floating point code into the library, for example with some macro checkpoints (#ifdef LTM_USE_FLOATS ...) and an entry in the makefiles to define that macro at compile time (methods vary between the different compilers).

Yes, floating points are a pain in the behind, I know.

[minad]

• I agree with returning an error for NaN and infinity.
• The code depends on endianness and on IEEE754 specifics. Supporting big endian would be ok, but I have no interest in supporting other floating point formats.
• The union trick is at least standarized in C11 afaik. We could also use memcpy instead.
• Yes, I am aware of isfinite. I did not use it since I don't know tommaths policy on additional includes.
• Using an additional floating point macro might make sense, but tommath already has this very fine grained macro system - so I am not sure if it is needed.

[minad]

@czurnieden Ah, know I understand your point about "the hint", you mean defining an additional function class like LTM_ALL_WITHOUT_FLOAT and LTM_REALLY_ALL or something like that. Is there a decision on how this should be done?

[minad]

Now that I think about it - maybe tommath should stay floating point free for embedded purposes and because of those subtleties concerning their representation? It might make sense to split this PR?

[sjaeckel]

It might make sense to split this PR?

please yes! I'd say 3 PR's

1. mp_get_double() (should be more or less a noop review-wise)
2. mp_set_double() (could stay this PR as there are already comments regarding this function)
3. two's complement functions

Already one comment regarding 3. how about changing their naming to mp_XX_tc? @hikari-no-yume you asked for this kind of function first IIRC, any input?

[minad]

@sjaeckel get_double/set_double is still in this pr due to the test.

[sjaeckel]

1. not sure though if what @czurnieden provided in Feature request: conversion from mp_int to float or double precision #3 (comment) ff. would be a tiny bit more portable?
2. Also I'm pretty surprised by the difference in LOC of your mp_get_double() vs. the one of @czurnieden
3. Regarding the macro I'd propose LTM_WITHOUT_FLOAT as IMO it should be an opt-out

[sjaeckel]

these changes should go in the other PR

[sjaeckel]

how about

diff --git a/demo/demo.c b/demo/demo.c
index e2c8aed..728d132 100644
--- a/demo/demo.c
+++ b/demo/demo.c
@@ -264,6 +264,17 @@ int main(void)
    // test mp_get_double/mp_set_double
-   printf("\n\nTesting: mp_get_double");
+   q = (unsigned long long)rand();
+   r = q;
+   printf("\n\nTesting: mp_get_double\n");
    for (i = 0; i < 1000; ++i) {
-      int tmp = rand();
-      double dbl = (double)tmp * rand() + 1;
+      double dbl = (double)q + (rand() * rand()) + 1;
+      q = (unsigned long long)dbl;
+      if ((i%72) == 0) {
+         putchar('\r');
+      }
+      if (r < q) {
+         putchar('-');
+      } else {
+         putchar('+');
+      }
+      r = q;
       mp_set_double(&a, dbl);

so dbl grows with each iteration and we'll have a wider range covered!?

[minad]

Yes, this is only very basic here. Fractions should also be tested.

[sjaeckel]

i --> 0 ❤️

🤣

[minad]

:D

[sjaeckel]

can you please create a separate commit for changes in the code and updates of the makefiles, callgraph & tommath_class.h ? that makes cherry-picking a lot easier

[czurnieden]

My old code from #3 is still a bit lacking. Slightly more elaborated versions:

mp_get_double.c

#include <tommath_private.h>
#ifdef BN_MP_GET_DOUBLE_C

#ifndef LTM_WITHOUT_FLOAT

#ifndef __STDC_IEC_559__
#error "mp_get_double needs IEEE-754/IEC-60559 compliant floating point arithmetic."
#endif

/* To avoid inclusion of float.h */
#ifndef DBL_MANT_DIG
#define DBL_MANT_DIG 53
#endif

static double st_pow(double d, int e)
{
   double t;

   if (e == 0) {
      return d;
   }
   t = 1.0;
   while (e > 1) {
      if (e % 2 == 0) {
         d *= d;
         e /= 2;
      } else {
         t *= d;
         d *= d;
         e = (e - 1)/2;
      }
   }
   return d * t;
}

/* Convert to double, assumes IEEE-754 conforming double. From code by
   gerdr (Gerhard R.) https://github.com/gerdr */
int mp_get_double(const mp_int *a, double *d)
{
   int digits_needed, i, limit;
   double multiplier;

   if (sizeof(double) != 8) {
      return MP_VAL;
   }

   if (mp_iszero(a)) {
      *d = 0.0;
      return MP_OKAY;
   }

   if (mp_count_bits(a) > 2047) {
      /* or return +/- infinity? */
      return MP_VAL;
   }

   /* digits_needed * DIGIT_BIT >= 64 bit */
   digits_needed = ((DBL_MANT_DIG + DIGIT_BIT) / DIGIT_BIT) + 1;

   multiplier = (double)(MP_MASK + 1);

   *d = 0.0;

   /* Could be assumed, couldn't it? */
   mp_clamp(a);

   i = a->used;

   limit = (i <= digits_needed) ? 0 : i - digits_needed;

   while (i-- > limit) {
      *d += a->dp[i];
      *d *= multiplier;
   }

   if (a->sign == MP_NEG) {
      *d *= -1.0;
   }

   *d *= st_pow(2.0, i * DIGIT_BIT);

   return MP_OKAY;
}

#endif
#endif

And mp_set_double.c

#include <tommath_private.h>
#ifdef BN_MP_SET_DOUBLE_C

#ifndef LTM_WITHOUT_FLOAT

#ifndef __STDC_IEC_559__
#error "mp_set_double needs IEEE-754/IEC-60559 compliant floating point arithmetic."
#endif

typedef union {
   struct {
      uint32_t a; /* low in little endian */
      uint32_t b; /* high in little endian */
   } w;
   double u;
} dblintern;

static int is_little_endian()
{
   int i = 1;
   char *c = (char *)&i;

   if (c[0] == 1) {
      return 1;
   } else {
      return 0;
   }
}
/* reduced frexp (neither checks nor balances!) */
static double mp_frexp(double d, int *e)
{
   int little_endian;
   uint32_t high;

   dblintern v;

   little_endian = is_little_endian();

   if (little_endian) {
      v.u = d;
      high = v.w.b;
   } else {
      v.u = d;
      high = v.w.a;
   }
   /* get exponent */
   /* shift right 20 bits */
   *e = high / 1048576;
   /* extract exponent from sign+exponent */
   *e %= 2048;
   /*  subtract bias */
   *e -= 1022;
   /* get fraction */
   high = high % 1048576;
   /* set exponent to bias */
   high += 1071644672;
   /* reassemble number */
   if (little_endian) {
      v.w.b = high;
   } else {
      v.w.a = high;
   }
   return v.u;
}

/* integer part of a double */
int mp_set_double(const mp_int *c, double d)
{
   int expnt, res, sign;
   double frac;

   if (sizeof(double) != 8) {
      return MP_VAL;
   }

   /* is the input usable? */
   frac = mp_frexp(d, &expnt);
   /* (s,q)NaN */
   if ((expnt + 1022) == 2047 && frac != 0) {
      return MP_VAL;
   }
   /* +/- infinity */
   if ((expnt + 1022) == 2047 && frac == 0) {
      return MP_VAL;
   }

   sign = (d < 0) ? MP_NEG : MP_ZPOS;
   mp_zero(c);
   d = (d < 0) ? d * -1.0 : d;

   if (d < 1.0 && d > 0.0) {
      return MP_OKAY;
   }

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

   /*
      This can be done directly if double is an IEEE-745 binary64
      but needs a branch for endiannes in that case (just copy the
      one from mp_frexp()).
   */
   for(;;){
      frac *= 2.0;
      if (frac >= 1.0) {
         if ((res = mp_add_d(c, 1, c)) != MP_OKAY) {
            return res;
         }
         frac -= 1.0;
      }
   }
   if (expnt > 0) {
      if ((res = mp_mul_2d(c, expnt, c)) != MP_OKAY) {
         return res;
      }
   }

   c->sign = sign;
   return MP_OKAY;
}

#endif
#endif

It is quite a lot but I avoided the inclusion of math.h and float.h.
Obligatory caveat: not tested yet and uses the union-trick that seems to be problematic in C++.

[minad]

@czurnieden What do you suggest? Instead of using isfinite in my code I can also check the nan/infinity bitpattern as in your code. I could also (with a little effort) ensure that my set_double function works on big endian. However I don't have access to any big endian machine and I don't want to fiddle with emulators just to test this one thing. In another project I tested my function quick-check style and found no differences to other big integer implementations. However I don't have a conclusive test suite available which is appropriate to be included here.

[czurnieden]

@minad
I tried to minimize the use of external functions, especially from the math library, that's all. isfinite had been put into the C-standard since 9899:1999, the macro to test for it is #if __STDC_VERSION__ >=199901L (the use of normal operations instead of bit-juggling was not intentional, it was originally meant to test the optimizations of a NIH-C-compiler and was the first implementation of frexp I found on my HD to C&P into that code. Sorry for that. )

Any time you do bit-juggling you get an endiannes problem. You can either do the usual ignore&document or make a branch (small and big is almost always sufficient). Some of the modern processors in wide use can switch at runtime, so a compile-time test is insufficient (see e.g.: https://community.arm.com/processors/f/discussions/3313/how-and-where-do-i-change-endianness-in-cortex-a9 for a short example or https://stackoverflow.com/questions/6212951/endianness-of-android-ndk for a slightly deeper discussion).

To avoid any bit-juggling and get rid of the problem with endianness , you would need the frexp from the C-library (also there since C-99; C-89 works, too. Only for double but that's what is needed here) and use a loop to read the bits of the fraction instead (I just saw that I deleted the shift after the addition which was in the original code. Why, oh, why did I change it? *grrr*).

You can do the frexp without bit-juggling but it is much slower, even with the optimization.

double frexp(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
  };
  /* can be done with the other method once you have the exponent
  if (isinf(x)) {
    *eptr = 0;
    return x;
  }
  if (isnan(x)) {
    *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;
}

That way you don't have any problems with endiannes but it is a lot of code, assumes IEEE-754 compatibility, and is slow, too.

For the test: you need tests for the extremes and some random values. The extremes here are NaN, +/- infinity, 0, denormal numbers (one example is enough), DBL_MAX, -DBL_MAX, two numbers, one positive and one negative that are one bit smaller than 1 and one bit bigger than -1 respectively (e.g.: with double smallerthanone = 1.0 - DBL_EPSILON/FLT_RADIX), and a number that is different depending on (wrong/correct) endiannes. That together with a handful of random values will do as a test for the implementation of mp_set_double.
You don't need all of those for a test for mp_get_double, of course.

(You should, in theory, also check that FLT_RADIX == 2. The value of FLT_EVAL_METHOD is not needed here but if you do a test you need to check if that value is anything other than -1, 0, 1, 2. (both macros in float.h))

Did I already mentioned that working with floating points is a real pain in the behind? ;-)

[minad]

@czurnieden Thank you. Yes, floating point is subtle. However my functions are tested for the special cases. I just didn't bother to support non-le and non-ieee754 architectures. Simply because it doesn't make sense for my project, which uses libtommath as an alternative to gmp. If the use case comes up to support ppc, sparc or endian switching arm, I can adapt. Those arms are usually running in le mode (at least debian arm assumes the le mode afaik).

The question is what makes sense for tommath. I think supporting fixed be and le, and ieee754 should be enough. It was a historical mistake that C left too many things unspecified, like CHAR_BIT==7?! So i would focus on the main use cases. Bit fiddling is ok if IEEE754 is assumed. FLT_RADIX!=2 is the same nonsense as CHAR_BIT. Are there any such architectures around? I would not support weird endianness switching at runtime. Avoiding libc is preferable.

But then given all those complications and since there is no other floating point code in tommath, I would just keep it outside of tommath and continue to use my own code in my project. Maybe the functions could also be added as unsupported extension to tommath for people who need it.

@sjaeckel What do you think?

[czurnieden]

@minad

I think supporting fixed be and le, and ieee754 should be enough.

OK for me.

It was a historical mistake that C left too many things unspecified, like CHAR_BIT==7?

There were a lot of architectures at the time of its birth, some of them really weird. You won't find CHAR_BIT==7 in modern machines but 16 and 32 are used quite regularly in DSPs.

FLT_RADIX!=2 is the same nonsense as CHAR_BIT

Yes, it is not really necessary but it is one extra independent point to check for IEEE-754 compliance and it costs nothing to check for. (The other value for FLT_RADIX would be ten but that gets a different type in IEEE-754, so FLT_RADIX is always two in IEEE-754 compliant systems)

Oh, I nearly forgot: one of the systems with FLT_RADIX != 2 is the IBM 370 (still supported by Redhat, if I'm not mistaken), of course! See a short excerpt of its float.h here: https://gcc.gnu.org/ml/gcc/2002-09/msg00038.html

I would not support weird endianness switching at runtime.

Its not weird, just uncomfortable ;-)
But serious: its completely ok to ignore it but it needs to be documented in that case!

Avoiding libc is preferable.

The Libc is already used in Libtommath, but we should avoid Libmath if possible, because it is an extra library in some systems. and it saves program-memory if not loaded.

But then given all those complications and since there is no other floating point code in tommath,

I have FFT-multiplication in my fork (in my master branch) so I would support a branch of libtommath that uses floating point types. Based on develop maybe?

Maybe the functions could also be added as unsupported extension to tommath for people who need it.

No, nothing unsupported, please. It shouldn't be a lot of work to support those two functions, just keep the branch up to date and adapt to API-changes but those changes happen once in a blue moon. The last one was the addition of the const qualifier throughout Libtommath where applicable.

I also started to fill the blanks in libtomfloat some time ago, so some of the functions in libtomfloat that manipulate mp_int can go to that branch of libtommath.

[czurnieden]

@minad
I took the liberty to try my teeth at mp_set_double and put the result in the branch https://github.com/czurnieden/libtommath/tree/withfloats

• A handful of values get tested at demo/demo.c lines 266ff - I might have forgot some.
• Should be endianness agnostic (big and little endian only) at runtime and does not rely on any "third-party" functions.
• Needs an IEEE-754 compliant double in binary64 format as input because it does some bit-juggling.
• Code is extensively commented, but not a single line in the actual documentation has been added.

Please check if it is to your taste.

[minad]

@czurnieden Sure! Great if you take a stab at it. Some comments:

• Endianness won't switch at runtime or the compiler will optimize your check away in the end I think. I would prefer a macro check, eliminating the helper function is_little_endian. However an endian-at-runtime-switching-aware compiler might not optimize away the check ;) Just to make sure, it would be great if you check the generated assembly.
• Your mp_get_double function seems overly complicated to my version. Why? I don't see the need for st_pow.
• I would like to have so random test additional to the edge cases, in particular checking round tripping between get/set. Maybe another few cases are interesting, testing integer numbers larger than 2^64, testing integers larger than the limit where doubles cannot represent integers precisely anymore (2^53?)

[czurnieden]

@minad

I would prefer a macro check, eliminating the helper function is_little_endian.

Makes the code itself simpler but moves the problem to compile time, so it needs a note in the documentation and a line or two in the makefiles or tommath.h which should not be a large problem.
That's why I put that branch in my own fork for us to be able to play with it.

Your mp_get_double function seems overly complicated to my version. Why?

I just took the code gerdr (Gerhard R.) https://github.com/gerdr offered last time.
You can simplify it, of course, but only slightly because some of teh methods used reduce the rounding errors.

I don't see the need for st_pow

To avoid the use of a function from the mathlib.
I changed it to direct manipulation of the exponent which eliminates any rounding errors in the old function. It could be put inline now to avoid the need of that extra function call.

I would like to have so random test additional to the edge cases,

Feel free to add some ;-)

testing integers larger than the limit where doubles cannot represent integers precisely anymore

I added one of it already: 123e123. That number can be represented exactly with a bigint but not with a binary64 double. And it is also a case where a round-trip would only work in one direction double -> mp_int -> double (which should always work) but not mp_int -> double -> mp_int.

Tests for mp_get_double will be done later this day (it's two o'clock in the morning here). Yes, that means that I don't even know yet if my version of mp_get_double works at all.

[sjaeckel]

@minad can you please rebase on top of current develop and force-push?

[sjaeckel]

Regarding the last comments/questions:

• Endianness ...

ltm was endianness unaware until now (it simply didn't need to know about it), but has run-time support for both, little&big, in mp_import() resp. mp_export()
IMO it should stay this way, as there exist still architectures that allow switching the endianness AFAIK

To avoid the use of a function from the mathlib.

I'm fully pro staying independent of libm, especially since ltm is a lot in use on embedded devices where there's normally no need for libm.
On the other side one can argue that if you're dealing with doubles you're normally already linking against libm or not!? ... I'd make it easy to turn-off those functions with one switch for now, something like LTM_NO_DOUBLE

I'd happily accept a second PR with the implementation of @czurnieden added with a compile-time switch like LTM_DOUBLE_WITHOUT_LIBM ... or something like that ;)

Tests

yep, they need some more coverage of edge-/corner-cases, I think then we could be nearly ready to merge after the build turned green :)

[sjaeckel]

IMO it should stay this way, as there exist still architectures that allow switching the endianness AFAIK

TBH after re-reading my comment I'm not sure if this would require a run-time switch ...

[minad]

ltm was endianness unaware until now (it simply didn't need to know about it), but has run-time support for both, little&big, in mp_import() resp. mp_export()
IMO it should stay this way, as there exist still architectures that allow switching the endianness AFAIK

Agree. But switching at runtime should not be supported. It is sufficient if the code compiles for both little and big.

On the other side one can argue that if you're dealing with doubles you're normally already linking against libm or not!? ... I'd make it easy to turn-off those functions with one switch for now, something like LTM_NO_DOUBLE

libm is not needed for mp_set_double since no double operations are performed. libm might be needed for mp_get_double if the platform does not support double multiplication etc natively.

yep, they need some more coverage of edge-/corner-cases, I think then we could be nearly ready to merge after the build turned green :)

Agree.

[sjaeckel]

TBH after re-reading my comment I'm not sure if this would require a run-time switch ...

+

Agree. But switching at runtime should not be supported. It is sufficient if the code compiles for both little and big.

= 👍

[minad]

Hmm, concerning endianness - I am not sure where endianness plays a role when I extract the double bits, since uint64_t and double should have the same endianness. Right now we test for little endian and throw a compile error on other platforms. Could someone help who has access to a big endian system (sparc...)? @czurnieden?

[sjaeckel]

Could someone help who has access to a big endian system

this reminds me of https://gist.github.com/sjaeckel/94120aeb2dacfa693b1dd360de451b89 and that I wanted to turn this into a real blogpost & build-script ...or something like that... ;)

[sjaeckel]

that's why I proposed to introduce LTM_NO_DOUBLE

alternatively you can also handle the API definition in tommath.h the same way :)

[minad]

Sure I could use an emulator but this will cost time which I could spend in better ways... ;)

[minad]

that's why I proposed to introduce LTM_NO_DOUBLE
alternatively you can also handle the API definition in tommath.h the same way :)

This will only make the header less clean. Right now if only the prototype exists, linking will fail. And compilation of the double function will print a warning.

[minad]

@karel-m Which floating format do those machines use? We rely on IEEE754, maybe they are just missing the predefined macros?

@sjaeckel Is it ready to be merged, what do you think?

[sjaeckel]

Neither AIX/Power (IBM xlc compiler) nor HP_UX/PA-RISC (HP cc compiler) have __STDC_IEC_559__ but both boxes are pretty old.

@karel-m so none of them defines __STDC_IEC_559__ but only the HP-UX box gives the following??

bn_mp_set_double.c:
cpp: "bn_mp_set_double.c", line 48: warning 2029: "mp_set_double implementation is only available on plat

@minad as soon as we've figured this out it should be ready to merge

[sjaeckel]

@minad I think the error checks in https://github.com/czurnieden/libtommath/blob/withfloats/bn_mp_get_double.c#L110 are a good idea to incorporate now!

[minad]

@sjaeckel I guess on AIX karel did not get the message since he had the ifdef removed there completely. What kind of error checks do you mean? sizeof(double) != 8 should be a static_assert, but I don't think it is necessary. The check for overflow should really set the double to +-inf instead of returning MP_VAL. My current implementation returns inf in that case as desired. There is no reason for an additional separate check to handle the corner case.

[sjaeckel]

OK, for me it's ready then, @karel-m and/or @czurnieden do you have any objections?

[karel-m]

Which floating format do those machines use?

I do not know.

so none of them defines STDC_IEC_559 but only the HP-UX box gives the following?

Both seem to be missing the define and both show the warning.

IBM xlc compiler doesn't like (except C++ comments) also the following:

   if (mp_set_double(&a, +1.0/0.0) != MP_VAL) {
   ..
   if (mp_set_double(&a, -1.0/0.0) != MP_VAL) {
   ..
   if (mp_set_double(&a, +0.0/0.0) != MP_VAL) {
   ..
   if (mp_set_double(&a, -0.0/0.0) != MP_VAL) {

The warnings:

"demo/demo.c", line 419.31: 1506-232 (I) Divisor for modulus or division operator cannot be zero.
"demo/demo.c", line 423.31: 1506-232 (I) Divisor for modulus or division operator cannot be zero.
"demo/demo.c", line 427.31: 1506-232 (I) Divisor for modulus or division operator cannot be zero.
"demo/demo.c", line 431.31: 1506-232 (I) Divisor for modulus or division operator cannot be zero.

[karel-m]

OK, for me it's ready then, @karel-m and/or @czurnieden do you have any objections?

No objections from me

[sjaeckel]

IBM xlc compiler doesn't like (except C++ comments) also the following:

shouldn't we then fix that before merging?

[sjaeckel]

except C++ comments

I'll fix those C++ comments after this PR is merged as they're allover demo.c ...

[sjaeckel]

"demo/demo.c", line 419.31: 1506-232 (I) Divisor for modulus or division operator cannot be zero. "demo/demo.c", line 423.31: 1506-232 (I) Divisor for modulus or division operator cannot be zero. "demo/demo.c", line 427.31: 1506-232 (I) Divisor for modulus or division operator cannot be zero. "demo/demo.c", line 431.31: 1506-232 (I) Divisor for modulus or division operator cannot be zero.

but those lines are within an #ifdef __STDC_IEC_559__ so it seems that the compiler defines it or did you remove the ifdef?!

[czurnieden]

The divisions by zero in demo.c can be skipped if some…uhm…oversensitive compilers try to be overly helpful but the respective macros like INFINITY and NAN are defined in math.h and these (standard-conforming) workarounds are necessary if we do not want to include math.h.

OK, for me it's ready then, @karel-m and/or @czurnieden do you have any objections?

As far as I was able to follow the discussion: no objections.

[karel-m]

or did you remove the ifdef?

yes, after removing the ifdef

[karel-m]

BTW this is Windows 10/x64 + mingw-w64 gcc-7.1.0

bn_mp_set_double.c:48:4: warning: #warning "mp_set_double implementation is only available on platforms with IEEE754 floating point format" [-Wcpp]
 #  warning "mp_set_double implementation is only available on platforms with IEEE754 floating point format"
    ^~~~~~~

[karel-m]

And the same warning with cygwin/64bit + gcc-7.3 (Target: x86_64-pc-cygwin)

[czurnieden]

@karel-m would you mind to check the value of the macro __GCC_IEC_559 ?

[minad]

@karel-m @czurnieden It seems STDC_IEC_559 is not defined on many platforms despite using IEEE754 format. The reason might be that the platforms are not fully compliant. I guess I warning has to suffice then.

[karel-m]

@czurnieden changing the ifdef to #if defined(__STDC_IEC_559__) || defined(__GCC_IEC_559) fixes mingw-w64 ang cygwin64 builds

[minad]

@karel-m Thank you. I added #if defined(__STDC_IEC_559__) || defined(__GCC_IEC_559)

[sjaeckel]

This is ready now?

[minad]

Yes!