mpz_limbs interface
Niels Möller
nisse at lysator.liu.se
Thu Feb 6 08:56:59 UTC 2014
nisse at lysator.liu.se (Niels Möller) writes:
> For mpn_set_d, I think it would make some sense to have it return a
> base-2 exponent, and write the mantissa to a few limbs. Number of limbs
> would be a constant, part of the ABI, similar to LIMBS_PER_DOUBLE but
> renamed for external use.
>
> mp_bitcnt_t
> mpn_set_d (mp_limb_t rp[LIMBS_PER_BOUBLE], d);
Below is a patch to do this (and return value is long, not mp_bitcnt_t,
since it needs to be signed).
What do you think?
Regards,
/Niels
diff -Nprc gmp.2462ec1bc65c/gmp-h.in gmp/gmp-h.in
*** gmp.2462ec1bc65c/gmp-h.in 2014-02-06 09:29:33.665957246 +0100
--- gmp/gmp-h.in 2014-02-06 09:29:33.669957220 +0100
*************** __GMP_DECLSPEC mp_bitcnt_t mpn_scan1 (mp
*** 1584,1589 ****
--- 1584,1593 ----
#define mpn_set_str __MPN(set_str)
__GMP_DECLSPEC mp_size_t mpn_set_str (mp_ptr, const unsigned char *, size_t, int);
+ #define MPN_SET_D_SIZE ((53 + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS)
+ #define mpn_set_d __MPN(set_d)
+ __GMP_DECLSPEC long mpn_set_d (mp_ptr, double);
+
#define mpn_sizeinbase __MPN(sizeinbase)
__GMP_DECLSPEC size_t mpn_sizeinbase (mp_srcptr, mp_size_t, int);
diff -Nprc gmp.2462ec1bc65c/mpn/generic/set_d.c gmp/mpn/generic/set_d.c
*** gmp.2462ec1bc65c/mpn/generic/set_d.c 1970-01-01 01:00:00.000000000 +0100
--- gmp/mpn/generic/set_d.c 2014-02-06 09:29:33.669957220 +0100
***************
*** 0 ****
--- 1,160 ----
+ /* mpn_set_d convert from double to array of mp_limb_t.
+
+ Copyright 1996, 1999-2002, 2006, 2012, 2014 Free Software Foundation, Inc.
+
+ This file is part of the GNU MP Library.
+
+ The GNU MP Library is free software; you can redistribute it and/or modify
+ it under the terms of either:
+
+ * the GNU Lesser General Public License as published by the Free
+ Software Foundation; either version 3 of the License, or (at your
+ option) any later version.
+
+ or
+
+ * the GNU General Public License as published by the Free Software
+ Foundation; either version 2 of the License, or (at your option) any
+ later version.
+
+ or both in parallel, as here.
+
+ The GNU MP Library is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received copies of the GNU General Public License and the
+ GNU Lesser General Public License along with the GNU MP Library. If not,
+ see https://www.gnu.org/licenses/. */
+
+ #include "gmp.h"
+ #include "gmp-impl.h"
+
+ #ifdef XDEBUG
+ #undef _GMP_IEEE_FLOATS
+ #endif
+
+ #ifndef _GMP_IEEE_FLOATS
+ #define _GMP_IEEE_FLOATS 0
+ #endif
+
+ /* Stores the mantissa as MPN_SET_D_SIZE limbs at rp, and returns the
+ exponent. Defined so that after e = mpn_set_d (rp, double d),
+
+ d = {rp, MPN_SET_D_SIZE} 2^e
+
+ Always returns with rp[MPN_SET_D_SIZE - 1] > 0, but besides this
+ requirement, the precise normalization is unspecified.
+ */
+ long
+ mpn_set_d (mp_ptr rp, double d)
+ {
+ long exp;
+ ASSERT (d > 0); /* Exclude negatives */
+ ASSERT (d == d); /* Exclude NaN */
+ ASSERT (d != 0.5*d); /* Exclude infinities */
+
+ MPN_ZERO (rp, MPN_SET_D_SIZE);
+
+ #if _GMP_IEEE_FLOATS
+ {
+ #if defined (__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 8
+ /* Work around alpha-specific bug in GCC 2.8.x. */
+ volatile
+ #endif
+ union ieee_double_extract x;
+ mp_size_t i;
+ uint_least32_t manl, manh;
+ mp_bitcnt_t shift;
+ x.d = d;
+ exp = x.s.exp;
+ manl = x.s.manl;
+ manh = x.s.manh;
+
+ if (exp > 0)
+ /* Set implicit one bit */
+ manh |= 1UL << 20;
+ else
+ /* Denormalized case, not optimized */
+ while ( (manh & (1UL << 22)) == 0)
+ {
+ manh = (manh << 1) | (manl >> 31);
+ manl = (manl << 1) & 0xffffffff;
+ exp--;
+ }
+
+ /* Assign mantissa */
+ #if GMP_NUMB_BITS >= 32
+ rp[0] = manl;
+ #else
+ for (i = 0; manl > 0; manl >>= GMP_NUMB_BITS, i++)
+ rp[i] = manl & GMP_NUMB_MASK;
+ #endif
+
+ shift = 32 % GMP_NUMB_BITS;
+ i = 32 / GMP_NUMB_BITS;
+
+ rp[i++] |= ( (mp_limb_t) manh << shift) & GMP_NUMB_MASK;
+ if (GMP_NUMB_BITS - shift < 20)
+ {
+ manh >>= (GMP_NUMB_BITS - shift);
+ #if GMP_NUMB_BITS >= 32
+ rp[i] = manh;
+ #else
+ for (; manh > 0; manh >>= GMP_NUMB_BITS, i++)
+ rp[i] = manh & GMP_NUMB_MASK;
+ #endif
+ }
+
+ exp -= (1023 + 52);
+ }
+ #else /* Non-ieee-floats */
+ {
+ mp_size_t i;
+ mp_limb_t x;
+
+ exp = 0;
+ /* Normalize to range 0.5 <= d < 1.0. Arbitrarily, start with 32-bit steps. */
+ if (d < 1.0)
+ do
+ {
+ d *= 4294967296.0; /* 2^32 */
+ exp -= 32;
+ }
+ while (d < 1.0);
+ else if (d >= 4294967296.0)
+ do
+ {
+ d *= 1.0 / 4294967296.0;
+ exp += 32;
+ }
+ while (d >= 4294967296.0);
+ while (d >= 1.0)
+ {
+ d *= 0.5;
+ exp++;
+ }
+ /* Bits for the high limb. */
+ i = MPN_SET_D_SIZE - 1;
+
+ ASSERT (GMP_NUMB_BITS * i < 53);
+ x = (mp_limb_t) (2.0 * (double) (CNST_LIMB(1) << (53 - GMP_NUMB_BITS * i - 1)) * d);
+ rp[i] = x;
+ d -= x;
+ ASSERT (d < 1.0);
+
+ while (i > 0)
+ {
+ x = (mp_limb_t) (2.0 * (double) (CNST_LIMB(1) << (GMP_NUMB_BITS - 1)) * d);
+ rp[--i] = x;
+ d -= x;
+ ASSERT (d < 1.0);
+ }
+ exp -= 53;
+ }
+ #endif
+
+ ASSERT (rp[MPN_SET_D_SIZE - 1] > 0);
+ return exp;
+ }
diff -Nprc gmp.2462ec1bc65c/mpz/set_d.c gmp/mpz/set_d.c
*** gmp.2462ec1bc65c/mpz/set_d.c 2014-02-06 09:29:33.657957297 +0100
--- gmp/mpz/set_d.c 2014-02-06 09:29:33.665957246 +0100
*************** see https://www.gnu.org/licenses/. */
*** 38,117 ****
#include "gmp-impl.h"
- /* We used to have a special case for d < MP_BASE_AS_DOUBLE, just casting
- double -> limb. Unfortunately gcc 3.3 on powerpc970-apple-darwin6.8.5
- got this wrong. (It assumed __fixunsdfdi returned its result in a single
- 64-bit register, where instead that function followed the calling
- conventions and gave the result in two parts r3 and r4.) Hence the use
- of __gmp_extract_double in all cases. */
-
void
mpz_set_d (mpz_ptr r, double d)
{
int negative;
! mp_limb_t tp[LIMBS_PER_DOUBLE];
mp_ptr rp;
mp_size_t rn;
DOUBLE_NAN_INF_ACTION (d,
__gmp_invalid_operation (),
__gmp_invalid_operation ());
!
negative = d < 0;
d = ABS (d);
! rn = __gmp_extract_double (tp, d);
!
! if (ALLOC(r) < rn)
! _mpz_realloc (r, rn);
!
! if (rn <= 0)
! rn = 0;
!
! rp = PTR (r);
!
! switch (rn)
{
! default:
! MPN_ZERO (rp, rn - LIMBS_PER_DOUBLE);
! rp += rn - LIMBS_PER_DOUBLE;
! /* fall through */
! #if LIMBS_PER_DOUBLE == 2
! case 2:
! rp[1] = tp[1], rp[0] = tp[0];
! break;
! case 1:
! rp[0] = tp[1];
! break;
! #endif
! #if LIMBS_PER_DOUBLE == 3
! case 3:
! rp[2] = tp[2], rp[1] = tp[1], rp[0] = tp[0];
! break;
! case 2:
! rp[1] = tp[2], rp[0] = tp[1];
! break;
! case 1:
! rp[0] = tp[2];
! break;
! #endif
! #if LIMBS_PER_DOUBLE == 4
! case 4:
! rp[3] = tp[3], rp[2] = tp[2], rp[1] = tp[1], rp[0] = tp[0];
! break;
! case 3:
! rp[2] = tp[3], rp[1] = tp[2], rp[0] = tp[1];
! break;
! case 2:
! rp[1] = tp[3], rp[0] = tp[2];
! break;
! case 1:
! rp[0] = tp[3];
! break;
! #endif
! case 0:
! break;
}
! SIZ(r) = negative ? -rn : rn;
}
--- 38,73 ----
#include "gmp-impl.h"
void
mpz_set_d (mpz_ptr r, double d)
{
int negative;
! mp_limb_t tp[MPN_SET_D_SIZE];
! mpz_t t = MPZ_ROINIT_N (tp, MPN_SET_D_SIZE);
mp_ptr rp;
mp_size_t rn;
+ long exp;
DOUBLE_NAN_INF_ACTION (d,
__gmp_invalid_operation (),
__gmp_invalid_operation ());
!
negative = d < 0;
d = ABS (d);
! if (d < 1.0)
{
! SIZ(r) = 0;
! return;
}
! exp = mpn_set_d (tp, d);
!
! if (exp >= 0)
! mpz_mul_2exp (r, t, exp);
! else
! mpz_fdiv_q_2exp (r, t, -exp);
!
! if (negative)
! SIZ (r) = -SIZ (r);
}
--
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.
More information about the gmp-devel
mailing list