ni...@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. _______________________________________________ gmp-devel mailing list gmp-devel@gmplib.org https://gmplib.org/mailman/listinfo/gmp-devel
