add d_mul_2exp and replace slow ldexp calls

Author: fredrik-johansson

doc/source/d_vec.rst

@@ -65,7 +65,7 @@ Comparison
     are within ``eps`` of each other, otherwise returns `0`.
 
 
-Addition and subtraction
+Arithmetic
 --------------------------------------------------------------------------------
 
 
@@ -79,6 +79,11 @@ Addition and subtraction
     Sets ``(res, len2)`` to ``(vec1, len2)`` minus ``(vec2, len2)``.
 
 
+.. function:: void _d_vec_mul_2exp(double * res, const double * vec, slong len, int e)
+
+    Sets ``(res, len)`` to ``(vec, len)`` multiplied by `2^e`.
+
+
 Dot product and norm
 --------------------------------------------------------------------------------
 

doc/source/double_extras.rst

@@ -33,6 +33,17 @@ Arithmetic
     ``len`` coefficients. Requires that ``len`` is nonzero.
 
 
+.. function:: double d_mul_2exp_inrange(double x, int i)
+              double d_mul_2exp_inrange2(double x, int i)
+              double d_mul_2exp(double x, int i)
+
+    Returns `x \cdot 2^i`.
+
+    The *inrange* version requires that `2^i` is in the normal exponent
+    range. The *inrange2* version additionally requires that both
+    `x` and `x \cdot 2^i` are in the normal exponent range,
+    and in particular also assumes that `x \ne 0`.
+
 
 Special functions
 --------------------------------------------------------------------------------

src/arb_mat/addmul_rad_mag_fast.c

@@ -18,6 +18,8 @@
 # include <math.h>
 #endif
 
+#include "double_extras.h"
+
 /* Block size for better cache locality. */
 #define BLOCK_SIZE 32
 
@@ -101,7 +103,7 @@ static inline slong _mag_get_exp(const mag_t x)
 static double
 mag_get_d_fixed_si(const mag_t x, slong e)
 {
-    return ldexp(MAG_MAN(x), MAG_EXP(x) - e - MAG_BITS);
+    return d_mul_2exp(MAG_MAN(x), MAG_EXP(x) - e - MAG_BITS);
 }
 
 void

src/arf/get.c

@@ -97,7 +97,7 @@ arf_get_d(const arf_t x, arf_rnd_t rnd)
         else
             v = (double)(tp[1]) + (double)(tp[0]) * ldexp(1,-32);
 
-        v = ldexp(v, ARF_EXP(t) - FLINT_BITS);
+        v = d_mul_2exp(v, ARF_EXP(t) - FLINT_BITS);
 
         if (ARF_SGNBIT(t))
             v = -v;

src/d_vec.h

@@ -12,7 +12,14 @@
 #ifndef D_VEC_H
 #define D_VEC_H
 
+#ifdef D_VEC_INLINES_C
+#define D_VEC_INLINE
+#else
+#define D_VEC_INLINE static inline
+#endif
+
 #include "flint.h"
+#include "double_extras.h"
 
 #ifdef __cplusplus
 extern "C" {
@@ -45,12 +52,32 @@ int _d_vec_is_zero(const double * vec, slong len);
 
 int _d_vec_is_approx_zero(const double * vec, slong len, double eps);
 
-/*  Addition  ****************************************************************/
+/*  Arithmetic  ****************************************************************/
 
 void _d_vec_add(double * res, const double * vec1, const double * vec2, slong len2);
 
 void _d_vec_sub(double * res, const double * vec1, const double * vec2, slong len2);
 
+D_VEC_INLINE
+void _d_vec_mul_2exp(double * res, const double * x, slong len, int e)
+{
+    slong i;
+
+    if (e >= D_MIN_NORMAL_EXPONENT && e <= D_MAX_NORMAL_EXPONENT)
+    {
+        double_uint64_u u;
+        u.i = ((int64_t) (e + D_EXPONENT_BIAS)) << D_EXPONENT_SHIFT;
+
+        for (i = 0; i < len; i++)
+            res[i] = x[i] * u.f;
+    }
+    else
+    {
+        for (i = 0; i < len; i++)
+            res[i] = ldexp(x[i], e);
+    }
+}
+
 /*  Dot product and norm  **************************************/
 
 double _d_vec_dot(const double * vec1, const double * vec2, slong len2);

src/d_vec/dot_heuristic.c

@@ -36,8 +36,8 @@ _d_vec_dot_heuristic(const double *vec1, const double *vec2, slong len2,
     {
         p = frexp(psum, &pexp);
         n = frexp(nsum, &nexp);
-        p = ldexp(1.0, pexp - D_BITS);
-        n = ldexp(1.0, nexp - D_BITS);
+        p = d_mul_2exp(1.0, pexp - D_BITS);
+        n = d_mul_2exp(1.0, nexp - D_BITS);
         d = FLINT_MAX(p, n);
         *err = 2 * len2 * d;
     }

src/double_extras/is_nan.c renamed to src/d_vec/inlines.c

@@ -1,6 +1,5 @@
 /*
-    Copyright (C) 2012 Fredrik Johansson
-    Copyright (C) 2014 Abhinav Baid
+    Copyright (C) 2015 William Hart
 
     This file is part of FLINT.
 
@@ -10,13 +9,7 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
-#include "double_extras.h"
+#define D_VEC_INLINES_C
 
-int
-d_is_nan(double x)
-{
-    if (x != x)
-        return 1;
-    else
-        return 0;
-}
+#include "flint.h"
+#include "d_vec.h"

src/double_extras.h

@@ -18,6 +18,7 @@
 #define DOUBLE_EXTRAS_INLINE static inline
 #endif
 
+#include <stdint.h>
 #include <math.h>
 #include "flint.h"
 
@@ -50,10 +51,56 @@ double d_polyval(const double * poly, int len, double x)
 
 double d_lambertw(double x);
 
-int d_is_nan(double x);
+DOUBLE_EXTRAS_INLINE
+int d_is_nan(double x)
+{
+    return x != x;
+}
 
 double d_log2(double x);
 
+typedef union
+{
+    double f;
+    uint64_t i;
+}
+double_uint64_u;
+
+#define D_MIN_NORMAL_EXPONENT -1022
+#define D_MAX_NORMAL_EXPONENT 1023
+#define D_EXPONENT_BIAS 1023
+#define D_EXPONENT_SHIFT 52
+
+/* Assumes that 2^i is in the normal exponent range. */
+FLINT_FORCE_INLINE double d_mul_2exp_inrange(double x, int i)
+{
+    FLINT_ASSERT(i >= D_MIN_NORMAL_EXPONENT && i <= D_MAX_NORMAL_EXPONENT);
+    double_uint64_u u;
+    u.i = ((int64_t) (i + D_EXPONENT_BIAS)) << D_EXPONENT_SHIFT;
+    return x * u.f;
+}
+
+/* Assumes that 2^i, x and x*2^i are all in the normal exponent range. */
+/* In particular, also assumes x != 0. */
+FLINT_FORCE_INLINE double d_mul_2exp_inrange2(double x, int i)
+{
+    FLINT_ASSERT(i >= D_MIN_NORMAL_EXPONENT && i <= D_MAX_NORMAL_EXPONENT);
+    FLINT_ASSERT(x != 0);
+
+    double_uint64_u u;
+    u.f = x;
+    u.i += ((int64_t) i) << D_EXPONENT_SHIFT;
+    return u.f;
+}
+
+FLINT_FORCE_INLINE double d_mul_2exp(double x, int i)
+{
+    if (i >= D_MIN_NORMAL_EXPONENT && i <= D_MAX_NORMAL_EXPONENT)
+        return d_mul_2exp_inrange(x, i);
+    else
+        return ldexp(x, i);
+}
+
 #ifdef __cplusplus
 }
 #endif

src/double_extras/test/main.c

@@ -17,6 +17,7 @@
 #include "t-is_nan.c"
 #include "t-lambertw.c"
 #include "t-log2.c"
+#include "t-mul_2exp.c"
 #include "t-randtest.c"
 #include "t-randtest_signed.c"
 
@@ -27,6 +28,7 @@ test_struct tests[] =
     TEST_FUNCTION(d_is_nan),
     TEST_FUNCTION(d_lambertw),
     TEST_FUNCTION(d_log2),
+    TEST_FUNCTION(d_mul_2exp),
     TEST_FUNCTION(d_randtest),
     TEST_FUNCTION(d_randtest_signed)
 };

src/double_extras/test/t-mul_2exp.c

@@ -0,0 +1,103 @@
+/*
+    Copyright (C) 2024 Fredrik Johansson
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "test_helpers.h"
+#include <float.h>
+#include "ulong_extras.h"
+#include "double_extras.h"
+#include "d_vec.h"
+
+TEST_FUNCTION_START(d_mul_2exp, state)
+{
+    double x, res1, res2;
+    int e, ok;
+    slong iter;
+
+    for (iter = 0; iter < 10000 * flint_test_multiplier(); iter++)
+    {
+        x = d_randtest_special(state, -1024, 1024);
+        e = n_randint(state, 3000) - 1500;
+
+        res1 = d_mul_2exp(x, e);
+        res2 = ldexp(x, e);
+
+        if (d_is_nan(res2))
+            ok = d_is_nan(res1);
+        else
+            ok = res1 == res2;
+
+        if (!ok)
+            TEST_FUNCTION_FAIL("x = %.20g\n res1 = %.20g\n res2 = %.20g\n", x, res1, res2);
+    }
+
+    for (iter = 0; iter < 10000 * flint_test_multiplier(); iter++)
+    {
+        do {
+            x = d_randtest_signed(state, -256, 256);
+        } while (x == 0.0);
+        e = n_randint(state, 512) - 256;
+
+        res1 = d_mul_2exp_inrange2(x, e);
+        res2 = ldexp(x, e);
+
+        ok = res1 == res2;
+
+        if (!ok)
+            TEST_FUNCTION_FAIL("x = %.20g\n res1 = %.20g\n res2 = %.20g\n", x, res1, res2);
+    }
+
+    for (iter = 0; iter < 10000 * flint_test_multiplier(); iter++)
+    {
+        x = d_randtest_special(state, -1024, 1024);
+        e = n_randint(state, D_MAX_NORMAL_EXPONENT - D_MIN_NORMAL_EXPONENT + 1) + D_MIN_NORMAL_EXPONENT;
+
+        res1 = d_mul_2exp_inrange(x, e);
+        res2 = ldexp(x, e);
+
+        if (d_is_nan(res2))
+            ok = d_is_nan(res1);
+        else
+            ok = res1 == res2;
+
+        if (!ok)
+            TEST_FUNCTION_FAIL("x = %.20g\n res1 = %.20g\n res2 = %.20g\n", x, res1, res2);
+    }
+
+    for (iter = 0; iter < 1000 * flint_test_multiplier(); iter++)
+    {
+        double v[5] = { 0, 0, 0, 0, 0 };
+        double r[10] = { 0, 0, 0, 0, 0 };
+        slong i, n = n_randint(state, 5);
+
+        for (i = 0; i < n; i++)
+            x = d_randtest_special(state, -1024, 1024);
+
+        e = n_randint(state, 3000) - 1500;
+
+        _d_vec_mul_2exp(r, v, n, e);
+
+        for (i = 0; i < n; i++)
+        {
+            res1 = r[i];
+            res2 = ldexp(v[i], e);
+
+            if (d_is_nan(res2))
+                ok = d_is_nan(res1);
+            else
+                ok = res1 == res2;
+
+            if (!ok)
+                TEST_FUNCTION_FAIL("x = %.20g\n res1 = %.20g\n res2 = %.20g\n", x, res1, res2);
+        }
+    }
+
+    TEST_FUNCTION_END(state);
+}