made set_double arch/mem-layout independent

.travis.yml

@@ -2,6 +2,7 @@ language: c
 
 install:
     - sudo apt-get update -qq
+    - sudo apt-get install -qq valgrind
     - sudo apt-get install gcc-multilib
 
 matrix:
@@ -26,6 +27,8 @@ env:
     BUILDOPTIONS="--with-low-mp"
   - |
     BUILDOPTIONS="--with-m64 --with-m32 --with-mx32"
+  - |
+    BUILDOPTIONS="--with-valgrind"
 
 after_failure:
   - cat test_*.log

README.md

@@ -24,6 +24,22 @@ Tests are located in `demo/` and can be built in two flavors.
 * `make test` creates a test binary that is intended to be run against `mtest`. `mtest` can be built with `make mtest` and test execution is done like `./mtest/mtest | ./test`. `mtest` is creating test vectors using an alternative MPI library and `test` is consuming these vectors to verify correct behavior of ltm
 * `make test_standalone` creates a stand-alone test binary that executes several test routines.
 
+[Valgrind](http://valgrind.org/) does not support `long double` but does not reset the limits in `float.h` accordingly.
+Please build `test.c` with the additional compiler flag `-DLTM_MEMCHECK_VALGRIND`, e.g.:
+
+```
+CFLAGS=" -DLTM_MEMCHECK_VALGRIND " make test_standalone
+````
+
+It will skip the tests for `long double` and prints a reminder to compile and run the tests without that flag to assure numerical correctness.
+
+In case your platform doesn't provide Valgrind at all and the auto-detection fails, please build `test.c` with the additional compiler flag `-DLTM_NO_VALGRIND`, e.g.:
+
+```
+CFLAGS=" -DLTM_NO_VALGRIND " make test_standalone
+````
+
+
 ## Building and Installing
 
 Building is straightforward for GNU Linux only, the section "Building LibTomMath" in the documentation in `doc/bn.pdf` has the details.

bn_mp_get_float.c

@@ -0,0 +1,43 @@
+#include "tommath_private.h"
+#ifdef BN_MP_GET_FLOAT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is a library that provides multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library was designed directly after the MPI library by
+ * Michael Fromberger but has been written from scratch with
+ * additional optimizations in place.
+ *
+ * SPDX-License-Identifier: Unlicense
+ */
+#include <float.h>
+#ifdef FLT_MAX
+/* This function is independent of the implementation of the floating point type */
+float mp_get_float(const mp_int *a)
+{
+   int i;
+   float d = 0.0, fac = 1.0;
+   for (i = 0; i < DIGIT_BIT; ++i) {
+      fac *= 2.0;
+   }
+   for (i = a->used; i --> 0;) {
+      d = (d * fac) + (float)DIGIT(a, i);
+   }
+   return (a->sign == MP_NEG) ? -d : d;
+}
+#else
+/* pragma message() not supported by several compilers (in mostly older but still used versions) */
+#  ifdef _MSC_VER
+#    pragma message("The type 'float' does not seem to be supported on your system.")
+#    pragma message("If that is wrong please contact the team at https://github.com/libtom/libtommath")
+#  else
+#    warning "The type 'float' does not seem to be supported on your system."
+#    warning "If that is wrong please contact the team at https://github.com/libtom/libtommath"
+#  endif
+#endif
+#endif
+
+/* ref:         $Format:%D$ */
+/* git commit:  $Format:%H$ */
+/* commit time: $Format:%ai$ */

bn_mp_get_long_double.c

@@ -0,0 +1,43 @@
+#include "tommath_private.h"
+#ifdef BN_MP_GET_LONG_DOUBLE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is a library that provides multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library was designed directly after the MPI library by
+ * Michael Fromberger but has been written from scratch with
+ * additional optimizations in place.
+ *
+ * SPDX-License-Identifier: Unlicense
+ */
+#include <float.h>
+#ifdef LDBL_MAX
+/* This function is independent of the implementation of the floating point type */
+long double mp_get_long_double(const mp_int *a)
+{
+   int i;
+   long double d = 0.0, fac = 1.0;
+   for (i = 0; i < DIGIT_BIT; ++i) {
+      fac *= 2.0;
+   }
+   for (i = a->used; i --> 0;) {
+      d = (d * fac) + (long double)DIGIT(a, i);
+   }
+   return (a->sign == MP_NEG) ? -d : d;
+}
+#else
+/* pragma message() not supported by several compilers (in mostly older but still used versions) */
+#  ifdef _MSC_VER
+#    pragma message("The type 'long double' does not seem to be supported on your system.")
+#    pragma message("If that is wrong please contact the team at https://github.com/libtom/libtommath")
+#  else
+#    warning "The type 'long double' does not seem to be supported on your system."
+#    warning "If that is wrong please contact the team at https://github.com/libtom/libtommath"
+#  endif
+#endif
+#endif
+
+/* ref:         $Format:%D$ */
+/* git commit:  $Format:%H$ */
+/* commit time: $Format:%ai$ */

bn_mp_iseven.c

@@ -9,8 +9,7 @@
  * Michael Fromberger but has been written from scratch with
  * additional optimizations in place.
  *
- * The library is free for all purposes without any express
- * guarantee it works.
+ * SPDX-License-Identifier: Unlicense
  */
 
 int mp_iseven(const mp_int *a)

bn_mp_isodd.c

@@ -9,10 +9,8 @@
  * Michael Fromberger but has been written from scratch with
  * additional optimizations in place.
  *
- * The library is free for all purposes without any express
- * guarantee it works.
+ * SPDX-License-Identifier: Unlicense
  */
-
 int mp_isodd(const mp_int *a)
 {
    return IS_ODD(a) ? MP_YES : MP_NO;

bn_mp_set_double.c

@@ -11,25 +11,29 @@
  *
  * SPDX-License-Identifier: Unlicense
  */
+#include <float.h>
+/* Minimum information needed */
+#if ( (defined DBL_MAX) && (FLT_RADIX == 2) )
 
-#if defined(__STDC_IEC_559__) || defined(__GCC_IEC_559)
+/* We can use a faster method if we have an IEEE compliant machine and a working stdint.h */
+#if ( (defined LTM_NEARLY_IEC_559 ) && (defined UINT64_MAX) )
 int mp_set_double(mp_int *a, double b)
 {
    uint64_t frac;
-   int exp, res;
+   int exp, res = MP_OKAY;
    union {
       double   dbl;
       uint64_t bits;
    } cast;
    cast.dbl = b;
 
-   exp = (int)((unsigned)(cast.bits >> 52) & 0x7FFU);
-   frac = (cast.bits & ((1ULL << 52) - 1ULL)) | (1ULL << 52);
+   exp = (int)((unsigned)(cast.bits >> (DBL_MANT_DIG - 1)) & (unsigned)((2*DBL_MAX_EXP) - 1));
+   frac = (cast.bits & ((1ULL << (DBL_MANT_DIG - 1)) - 1ULL)) | (1ULL << (DBL_MANT_DIG - 1));
 
-   if (exp == 0x7FF) { /* +-inf, NaN */
+   if (exp == ((2*DBL_MAX_EXP) - 1)) { /* +-inf, NaN */
       return MP_VAL;
    }
-   exp -= 1023 + 52;
+   exp -= (DBL_MAX_EXP - 1) + (DBL_MANT_DIG - 1);
 
    res = mp_set_long_long(a, frac);
    if (res != MP_OKAY) {
@@ -40,20 +44,119 @@ int mp_set_double(mp_int *a, double b)
    if (res != MP_OKAY) {
       return res;
    }
-
-   if (((cast.bits >> 63) != 0ULL) && !IS_ZERO(a)) {
+   /* 63 = number of bits without sign-bit in binary64 */
+   if (((cast.bits >> 63) != 0ULL) && (!IS_ZERO(a))) {
       a->sign = MP_NEG;
    }
 
    return MP_OKAY;
 }
+
+#else
+static double s_math_h_less_frexp(double x, int *exp)
+{
+   int exponent = 0, i;
+   double high = 2.0;
+   double low = 0.5;
+
+   if (x >= 1.0) {
+      for (i = 0; x >= high; i++) {
+         exponent += (1 << i);
+         x *= low;
+         low *= low;
+         high *= high;
+      }
+      if (x < 0.5) {
+         while (x < 0.5) {
+            x *= 2.0;
+            exponent--;
+         }
+      } else {
+         while (x >= 1.0) {
+            x /= 2.0;
+            exponent++;
+         }
+      }
+   }
+   *exp = exponent;
+   return x;
+}
+
+int mp_set_double(mp_int *a, double b)
+{
+   int exp = 0, res = MP_OKAY, sign = MP_ZPOS;
+
+   /* Check for NaN */
+   /* TODO: there might be some ompilers who do not return true in case of NaN */
+   if (b != b) {
+      return MP_VAL;
+   }
+
+   if (b < 0.0) {
+      b = b * (-1.0);
+      sign = MP_NEG;
+   }
+
+   /* Check for infinity */
+   if (b > DBL_MAX) {
+      return MP_VAL;
+   }
+
+   mp_zero(a);
+   /* Numbers smaller than 1 truncate to zero */
+   if (b < 1.0) {
+      a->sign = sign;
+      return MP_OKAY;
+   }
+   b = s_math_h_less_frexp(b, &exp);
+
+   while (exp-- >= 0) {
+      b *= 2.0;
+      if (b >= 1.0) {
+         if ((res = mp_add_d(a, 1, a)) != MP_OKAY) {
+            return res;
+         }
+         b -= 1.0;
+      }
+      if (exp >= 0) {
+         if ((res = mp_mul_2d(a, 1, a)) != MP_OKAY) {
+            return res;
+         }
+      }
+      if (b == 0.0) {
+         exp--;
+         break;
+      }
+   }
+   if (res != MP_OKAY) {
+      return res;
+   }
+   res = (exp < 0) ? mp_div_2d(a, -exp, a, NULL) : mp_mul_2d(a, exp, a);
+   if (res != MP_OKAY) {
+      return res;
+   }
+   a->sign = sign;
+   return res;
+}
+#endif
 #else
 /* pragma message() not supported by several compilers (in mostly older but still used versions) */
 #  ifdef _MSC_VER
-#    pragma message("mp_set_double implementation is only available on platforms with IEEE754 floating point format")
+#    pragma message("mp_set_double implementation is only available on platforms with IEEE754 floating point format.")
+#    pragma message("At least DBL_MAX_EXP must be defined and set and, for now, FLT_RADIX must be 2.")
 #  else
 #    warning "mp_set_double implementation is only available on platforms with IEEE754 floating point format"
+#    warning "At least DBL_MAX_EXP must be defined and set and, for now, FLT_RADIX must be 2."
 #  endif
+#if (FLT_RADIX == 16)
+#  ifdef _MSC_VER
+#    pragma message("No radices other than two are supported. IBM's z/OS uses IEEE-754 compliant floats")
+#    pragma message("with the compiler option FLOAT(IEEE)")
+#  else
+#    warning "No radices other than two are supported. IBM's z/OS uses IEEE-754 compliant floats"
+#    warning "with the compiler option FLOAT(IEEE)"
+#  endif
+#endif
 #endif
 #endif