| 1 | /* Split a double into fraction and mantissa. |
| 2 | Copyright (C) 2007-2019 Free Software Foundation, Inc. |
| 3 | |
| 4 | This program is free software: you can redistribute it and/or modify |
| 5 | it under the terms of the GNU General Public License as published by |
| 6 | the Free Software Foundation; either version 3 of the License, or |
| 7 | (at your option) any later version. |
| 8 | |
| 9 | This program 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 |
| 12 | GNU General Public License for more details. |
| 13 | |
| 14 | You should have received a copy of the GNU General Public License |
| 15 | along with this program. If not, see <https://www.gnu.org/licenses/>. */ |
| 16 | |
| 17 | /* Written by Paolo Bonzini <bonzini@gnu.org>, 2003, and |
| 18 | Bruno Haible <bruno@clisp.org>, 2007. */ |
| 19 | |
| 20 | #if ! defined USE_LONG_DOUBLE |
| 21 | # include <config.h> |
| 22 | #endif |
| 23 | |
| 24 | /* Specification. */ |
| 25 | #include <math.h> |
| 26 | |
| 27 | #include <float.h> |
| 28 | #ifdef USE_LONG_DOUBLE |
| 29 | # include "isnanl-nolibm.h" |
| 30 | # include "fpucw.h" |
| 31 | #else |
| 32 | # include "isnand-nolibm.h" |
| 33 | #endif |
| 34 | |
| 35 | /* This file assumes FLT_RADIX = 2. If FLT_RADIX is a power of 2 greater |
| 36 | than 2, or not even a power of 2, some rounding errors can occur, so that |
| 37 | then the returned mantissa is only guaranteed to be <= 1.0, not < 1.0. */ |
| 38 | |
| 39 | #ifdef USE_LONG_DOUBLE |
| 40 | # define FUNC frexpl |
| 41 | # define DOUBLE long double |
| 42 | # define ISNAN isnanl |
| 43 | # define DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING |
| 44 | # define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING () |
| 45 | # define END_ROUNDING() END_LONG_DOUBLE_ROUNDING () |
| 46 | # define L_(literal) literal##L |
| 47 | #else |
| 48 | # define FUNC frexp |
| 49 | # define DOUBLE double |
| 50 | # define ISNAN isnand |
| 51 | # define DECL_ROUNDING |
| 52 | # define BEGIN_ROUNDING() |
| 53 | # define END_ROUNDING() |
| 54 | # define L_(literal) literal |
| 55 | #endif |
| 56 | |
| 57 | DOUBLE |
| 58 | FUNC (DOUBLE x, int *expptr) |
| 59 | { |
| 60 | int sign; |
| 61 | int exponent; |
| 62 | DECL_ROUNDING |
| 63 | |
| 64 | /* Test for NaN, infinity, and zero. */ |
| 65 | if (ISNAN (x) || x + x == x) |
| 66 | { |
| 67 | *expptr = 0; |
| 68 | return x; |
| 69 | } |
| 70 | |
| 71 | sign = 0; |
| 72 | if (x < 0) |
| 73 | { |
| 74 | x = - x; |
| 75 | sign = -1; |
| 76 | } |
| 77 | |
| 78 | BEGIN_ROUNDING (); |
| 79 | |
| 80 | { |
| 81 | /* Since the exponent is an 'int', it fits in 64 bits. Therefore the |
| 82 | loops are executed no more than 64 times. */ |
| 83 | DOUBLE pow2[64]; /* pow2[i] = 2^2^i */ |
| 84 | DOUBLE powh[64]; /* powh[i] = 2^-2^i */ |
| 85 | int i; |
| 86 | |
| 87 | exponent = 0; |
| 88 | if (x >= L_(1.0)) |
| 89 | { |
| 90 | /* A positive exponent. */ |
| 91 | DOUBLE pow2_i; /* = pow2[i] */ |
| 92 | DOUBLE powh_i; /* = powh[i] */ |
| 93 | |
| 94 | /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, |
| 95 | x * 2^exponent = argument, x >= 1.0. */ |
| 96 | for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); |
| 97 | ; |
| 98 | i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) |
| 99 | { |
| 100 | if (x >= pow2_i) |
| 101 | { |
| 102 | exponent += (1 << i); |
| 103 | x *= powh_i; |
| 104 | } |
| 105 | else |
| 106 | break; |
| 107 | |
| 108 | pow2[i] = pow2_i; |
| 109 | powh[i] = powh_i; |
| 110 | } |
| 111 | /* Avoid making x too small, as it could become a denormalized |
| 112 | number and thus lose precision. */ |
| 113 | while (i > 0 && x < pow2[i - 1]) |
| 114 | { |
| 115 | i--; |
| 116 | powh_i = powh[i]; |
| 117 | } |
| 118 | exponent += (1 << i); |
| 119 | x *= powh_i; |
| 120 | /* Here 2^-2^i <= x < 1.0. */ |
| 121 | } |
| 122 | else |
| 123 | { |
| 124 | /* A negative or zero exponent. */ |
| 125 | DOUBLE pow2_i; /* = pow2[i] */ |
| 126 | DOUBLE powh_i; /* = powh[i] */ |
| 127 | |
| 128 | /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, |
| 129 | x * 2^exponent = argument, x < 1.0. */ |
| 130 | for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); |
| 131 | ; |
| 132 | i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) |
| 133 | { |
| 134 | if (x < powh_i) |
| 135 | { |
| 136 | exponent -= (1 << i); |
| 137 | x *= pow2_i; |
| 138 | } |
| 139 | else |
| 140 | break; |
| 141 | |
| 142 | pow2[i] = pow2_i; |
| 143 | powh[i] = powh_i; |
| 144 | } |
| 145 | /* Here 2^-2^i <= x < 1.0. */ |
| 146 | } |
| 147 | |
| 148 | /* Invariants: x * 2^exponent = argument, and 2^-2^i <= x < 1.0. */ |
| 149 | while (i > 0) |
| 150 | { |
| 151 | i--; |
| 152 | if (x < powh[i]) |
| 153 | { |
| 154 | exponent -= (1 << i); |
| 155 | x *= pow2[i]; |
| 156 | } |
| 157 | } |
| 158 | /* Here 0.5 <= x < 1.0. */ |
| 159 | } |
| 160 | |
| 161 | if (sign < 0) |
| 162 | x = - x; |
| 163 | |
| 164 | END_ROUNDING (); |
| 165 | |
| 166 | *expptr = exponent; |
| 167 | return x; |
| 168 | } |