1 /* Copyright (C) 1993, 1994, 1995, 1996, 1997 Free Software Foundation, Inc. 2 This file is part of the GNU C Library. 3 4 The GNU C Library is free software; you can redistribute it and/or 5 modify it under the terms of the GNU Lesser General Public 6 License as published by the Free Software Foundation; either 7 version 2.1 of the License, or (at your option) any later version. 8 9 The GNU C Library is distributed in the hope that it will be useful, 10 but WITHOUT ANY WARRANTY; without even the implied warranty of 11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12 Lesser General Public License for more details. 13 14 You should have received a copy of the GNU Lesser General Public 15 License along with the GNU C Library; if not, write to the Free 16 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 17 02111-1307 USA. */ 18 19 #include "gmp.h" 20 #include "gmp-impl.h" 21 #include "longlong.h" 22 #include <ieee754.h> 23 #include <float.h> 24 #include <stdlib.h> 25 26 /* Convert a `double' in IEEE754 standard double-precision format to a 27 multi-precision integer representing the significand scaled up by its 28 number of bits (52 for double) and an integral power of two (MPN frexp). */ 29 30 mp_size_t 31 __mpn_extract_double (mp_ptr res_ptr, mp_size_t size, 32 int *expt, int *is_neg, 33 double value) 34 { 35 union ieee754_double u; 36 u.d = value; 37 38 *is_neg = u.ieee.negative; 39 *expt = (int) u.ieee.exponent - IEEE754_DOUBLE_BIAS; 40 41 #if BITS_PER_MP_LIMB == 32 42 res_ptr[0] = u.ieee.mantissa1; /* Low-order 32 bits of fraction. */ 43 res_ptr[1] = u.ieee.mantissa0; /* High-order 20 bits. */ 44 #define N 2 45 #elif BITS_PER_MP_LIMB == 64 46 /* Hopefully the compiler will combine the two bitfield extracts 47 and this composition into just the original quadword extract. */ 48 res_ptr[0] = ((unsigned long int) u.ieee.mantissa0 << 32) | u.ieee.mantissa1; 49 #define N 1 50 #else 51 #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" 52 #endif 53 /* The format does not fill the last limb. There are some zeros. */ 54 #define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \ 55 - (DBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB))) 56 57 if (u.ieee.exponent == 0) 58 { 59 /* A biased exponent of zero is a special case. 60 Either it is a zero or it is a denormal number. */ 61 if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2. */ 62 /* It's zero. */ 63 *expt = 0; 64 else 65 { 66 /* It is a denormal number, meaning it has no implicit leading 67 one bit, and its exponent is in fact the format minimum. */ 68 int cnt; 69 70 if (res_ptr[N - 1] != 0) 71 { 72 count_leading_zeros (cnt, res_ptr[N - 1]); 73 cnt -= NUM_LEADING_ZEROS; 74 #if N == 2 75 res_ptr[N - 1] = res_ptr[1] << cnt 76 | (N - 1) 77 * (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); 78 res_ptr[0] <<= cnt; 79 #else 80 res_ptr[N - 1] <<= cnt; 81 #endif 82 *expt = DBL_MIN_EXP - 1 - cnt; 83 } 84 else 85 { 86 count_leading_zeros (cnt, res_ptr[0]); 87 if (cnt >= NUM_LEADING_ZEROS) 88 { 89 res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS); 90 res_ptr[0] = 0; 91 } 92 else 93 { 94 res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt); 95 res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt); 96 } 97 *expt = DBL_MIN_EXP - 1 98 - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt; 99 } 100 } 101 } 102 else 103 /* Add the implicit leading one bit for a normalized number. */ 104 res_ptr[N - 1] |= 1L << (DBL_MANT_DIG - 1 - ((N - 1) * BITS_PER_MP_LIMB)); 105 106 return N; 107 } 108