xf.li | bdd93d5 | 2023-05-12 07:10:14 -0700 | [diff] [blame^] | 1 | /* Copyright (C) 1995-2016 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, see |
| 16 | <http://www.gnu.org/licenses/>. */ |
| 17 | |
| 18 | #include "gmp.h" |
| 19 | #include "gmp-impl.h" |
| 20 | #include "longlong.h" |
| 21 | #include <ieee754.h> |
| 22 | #include <float.h> |
| 23 | #include <stdlib.h> |
| 24 | |
| 25 | /* Convert a `long double' in IEEE854 standard double-precision format to a |
| 26 | multi-precision integer representing the significand scaled up by its |
| 27 | number of bits (64 for long double) and an integral power of two |
| 28 | (MPN frexpl). */ |
| 29 | |
| 30 | mp_size_t |
| 31 | __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, |
| 32 | int *expt, int *is_neg, |
| 33 | long double value) |
| 34 | { |
| 35 | union ieee854_long_double u; |
| 36 | u.d = value; |
| 37 | |
| 38 | *is_neg = u.ieee.negative; |
| 39 | *expt = (int) u.ieee.exponent - IEEE854_LONG_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 32 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] = ((mp_limb_t) 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 | |
| 54 | if (u.ieee.exponent == 0) |
| 55 | { |
| 56 | /* A biased exponent of zero is a special case. |
| 57 | Either it is a zero or it is a denormal number. */ |
| 58 | if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2. */ |
| 59 | /* It's zero. */ |
| 60 | *expt = 0; |
| 61 | else |
| 62 | { |
| 63 | /* It is a denormal number, meaning it has no implicit leading |
| 64 | one bit, and its exponent is in fact the format minimum. */ |
| 65 | int cnt; |
| 66 | |
| 67 | /* One problem with Intel's 80-bit format is that the explicit |
| 68 | leading one in the normalized representation has to be zero |
| 69 | for denormalized number. If it is one, the number is according |
| 70 | to Intel's specification an invalid number. We make the |
| 71 | representation unique by explicitly clearing this bit. */ |
| 72 | res_ptr[N - 1] &= ~((mp_limb_t) 1 << ((LDBL_MANT_DIG - 1) % BITS_PER_MP_LIMB)); |
| 73 | |
| 74 | if (res_ptr[N - 1] != 0) |
| 75 | { |
| 76 | count_leading_zeros (cnt, res_ptr[N - 1]); |
| 77 | if (cnt != 0) |
| 78 | { |
| 79 | #if N == 2 |
| 80 | res_ptr[N - 1] = res_ptr[N - 1] << cnt |
| 81 | | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); |
| 82 | res_ptr[0] <<= cnt; |
| 83 | #else |
| 84 | res_ptr[N - 1] <<= cnt; |
| 85 | #endif |
| 86 | } |
| 87 | *expt = LDBL_MIN_EXP - 1 - cnt; |
| 88 | } |
| 89 | else if (res_ptr[0] != 0) |
| 90 | { |
| 91 | count_leading_zeros (cnt, res_ptr[0]); |
| 92 | res_ptr[N - 1] = res_ptr[0] << cnt; |
| 93 | res_ptr[0] = 0; |
| 94 | *expt = LDBL_MIN_EXP - 1 - BITS_PER_MP_LIMB - cnt; |
| 95 | } |
| 96 | else |
| 97 | { |
| 98 | /* This is the special case of the pseudo denormal number |
| 99 | with only the implicit leading bit set. The value is |
| 100 | in fact a normal number and so we have to treat this |
| 101 | case differently. */ |
| 102 | #if N == 2 |
| 103 | res_ptr[N - 1] = 0x80000000ul; |
| 104 | #else |
| 105 | res_ptr[0] = 0x8000000000000000ul; |
| 106 | #endif |
| 107 | *expt = LDBL_MIN_EXP - 1; |
| 108 | } |
| 109 | } |
| 110 | } |
| 111 | else if (u.ieee.exponent < 0x7fff |
| 112 | #if N == 2 |
| 113 | && res_ptr[0] == 0 |
| 114 | #endif |
| 115 | && res_ptr[N - 1] == 0) |
| 116 | /* Pseudo zero. */ |
| 117 | *expt = 0; |
| 118 | |
| 119 | return N; |
| 120 | } |