1 /* Copyright (C) 1995,1996,1997,1998,1999,2002,2003 2 Free Software Foundation, Inc. 3 This file is part of the GNU C Library. 4 5 The GNU C Library is free software; you can redistribute it and/or 6 modify it under the terms of the GNU Lesser General Public 7 License as published by the Free Software Foundation; either 8 version 2.1 of the License, or (at your option) any later version. 9 10 The GNU C Library is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 13 Lesser General Public License for more details. 14 15 You should have received a copy of the GNU Lesser General Public 16 License along with the GNU C Library; if not, write to the Free 17 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 18 02111-1307 USA. */ 19 20 #include <float.h> 21 #include <math.h> 22 #include <stdlib.h> 23 #include "gmp-impl.h" 24 25 /* Convert a `__float128' in IEEE854 quad-precision format to a 26 multi-precision integer representing the significand scaled up by its 27 number of bits (113 for long double) and an integral power of two 28 (MPN frexpl). */ 29 30 mp_size_t 31 mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size, 32 int *expt, int *is_neg, 33 __float128 value) 34 { 35 ieee854_float128 u; 36 u.value = value; 37 38 *is_neg = u.ieee.negative; 39 *expt = (int) u.ieee.exponent - IEEE854_FLOAT128_BIAS; 40 41 #if BITS_PER_MP_LIMB == 32 42 res_ptr[0] = u.ieee.mantissa3; /* Low-order 32 bits of fraction. */ 43 res_ptr[1] = u.ieee.mantissa2; 44 res_ptr[2] = u.ieee.mantissa1; 45 res_ptr[3] = u.ieee.mantissa0; /* High-order 32 bits. */ 46 #define N 4 47 #elif BITS_PER_MP_LIMB == 64 48 res_ptr[0] = ((mp_limb_t) u.ieee.mantissa2 << 32) | u.ieee.mantissa3; 49 res_ptr[1] = ((mp_limb_t) u.ieee.mantissa0 << 32) | u.ieee.mantissa1; 50 #define N 2 51 #else 52 #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" 53 #endif 54 /* The format does not fill the last limb. There are some zeros. */ 55 #define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \ 56 - (FLT128_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB))) 57 58 if (u.ieee.exponent == 0) 59 { 60 /* A biased exponent of zero is a special case. 61 Either it is a zero or it is a denormal number. */ 62 if (res_ptr[0] == 0 && res_ptr[1] == 0 63 && res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4. */ 64 /* It's zero. */ 65 *expt = 0; 66 else 67 { 68 /* It is a denormal number, meaning it has no implicit leading 69 one bit, and its exponent is in fact the format minimum. */ 70 int cnt; 71 72 #if N == 2 73 if (res_ptr[N - 1] != 0) 74 { 75 count_leading_zeros (cnt, res_ptr[N - 1]); 76 cnt -= NUM_LEADING_ZEROS; 77 res_ptr[N - 1] = res_ptr[N - 1] << cnt 78 | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); 79 res_ptr[0] <<= cnt; 80 *expt = FLT128_MIN_EXP - 1 - cnt; 81 } 82 else 83 { 84 count_leading_zeros (cnt, res_ptr[0]); 85 if (cnt >= NUM_LEADING_ZEROS) 86 { 87 res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS); 88 res_ptr[0] = 0; 89 } 90 else 91 { 92 res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt); 93 res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt); 94 } 95 *expt = FLT128_MIN_EXP - 1 96 - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt; 97 } 98 #else 99 int j, k, l; 100 101 for (j = N - 1; j > 0; j--) 102 if (res_ptr[j] != 0) 103 break; 104 105 count_leading_zeros (cnt, res_ptr[j]); 106 cnt -= NUM_LEADING_ZEROS; 107 l = N - 1 - j; 108 if (cnt < 0) 109 { 110 cnt += BITS_PER_MP_LIMB; 111 l--; 112 } 113 if (!cnt) 114 for (k = N - 1; k >= l; k--) 115 res_ptr[k] = res_ptr[k-l]; 116 else 117 { 118 for (k = N - 1; k > l; k--) 119 res_ptr[k] = res_ptr[k-l] << cnt 120 | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt); 121 res_ptr[k--] = res_ptr[0] << cnt; 122 } 123 124 for (; k >= 0; k--) 125 res_ptr[k] = 0; 126 *expt = FLT128_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt; 127 #endif 128 } 129 } 130 else 131 /* Add the implicit leading one bit for a normalized number. */ 132 res_ptr[N - 1] |= (mp_limb_t) 1 << (FLT128_MANT_DIG - 1 133 - ((N - 1) * BITS_PER_MP_LIMB)); 134 135 return N; 136 } 137