Austin Schuh | dace2a6 | 2020-08-18 10:56:48 -0700 | [diff] [blame] | 1 | /* mpq_set_d(mpq_t q, double d) -- Set q to d without rounding. |
| 2 | |
| 3 | Copyright 2000, 2002, 2003, 2012, 2014, 2018 Free Software Foundation, |
| 4 | Inc. |
| 5 | |
| 6 | This file is part of the GNU MP Library. |
| 7 | |
| 8 | The GNU MP Library is free software; you can redistribute it and/or modify |
| 9 | it under the terms of either: |
| 10 | |
| 11 | * the GNU Lesser General Public License as published by the Free |
| 12 | Software Foundation; either version 3 of the License, or (at your |
| 13 | option) any later version. |
| 14 | |
| 15 | or |
| 16 | |
| 17 | * the GNU General Public License as published by the Free Software |
| 18 | Foundation; either version 2 of the License, or (at your option) any |
| 19 | later version. |
| 20 | |
| 21 | or both in parallel, as here. |
| 22 | |
| 23 | The GNU MP Library is distributed in the hope that it will be useful, but |
| 24 | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
| 25 | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
| 26 | for more details. |
| 27 | |
| 28 | You should have received copies of the GNU General Public License and the |
| 29 | GNU Lesser General Public License along with the GNU MP Library. If not, |
| 30 | see https://www.gnu.org/licenses/. */ |
| 31 | |
| 32 | #include "config.h" |
| 33 | |
| 34 | #if HAVE_FLOAT_H |
| 35 | #include <float.h> /* for DBL_MAX */ |
| 36 | #endif |
| 37 | |
| 38 | #include "gmp-impl.h" |
| 39 | #include "longlong.h" |
| 40 | |
| 41 | #if LIMBS_PER_DOUBLE > 4 |
| 42 | choke me |
| 43 | #endif |
| 44 | |
| 45 | void |
| 46 | mpq_set_d (mpq_ptr dest, double d) |
| 47 | { |
| 48 | int negative; |
| 49 | mp_exp_t exp; |
| 50 | mp_limb_t tp[LIMBS_PER_DOUBLE]; |
| 51 | mp_ptr np, dp; |
| 52 | mp_size_t nn, dn; |
| 53 | int c; |
| 54 | |
| 55 | DOUBLE_NAN_INF_ACTION (d, |
| 56 | __gmp_invalid_operation (), |
| 57 | __gmp_invalid_operation ()); |
| 58 | |
| 59 | negative = d < 0; |
| 60 | d = ABS (d); |
| 61 | |
| 62 | exp = __gmp_extract_double (tp, d); |
| 63 | |
| 64 | /* There are two main version of the conversion. The `then' arm handles |
| 65 | numbers with a fractional part, while the `else' arm handles integers. */ |
| 66 | #if LIMBS_PER_DOUBLE == 4 |
| 67 | if (exp <= 1 || (exp == 2 && (tp[0] | tp[1]) != 0)) |
| 68 | #endif |
| 69 | #if LIMBS_PER_DOUBLE == 3 |
| 70 | if (exp <= 1 || (exp == 2 && tp[0] != 0)) |
| 71 | #endif |
| 72 | #if LIMBS_PER_DOUBLE == 2 |
| 73 | if (exp <= 1) |
| 74 | #endif |
| 75 | { |
| 76 | if (d == 0.0) |
| 77 | { |
| 78 | SIZ(NUM(dest)) = 0; |
| 79 | SIZ(DEN(dest)) = 1; |
| 80 | MPZ_NEWALLOC (DEN(dest), 1)[0] = 1; |
| 81 | return; |
| 82 | } |
| 83 | |
| 84 | #if LIMBS_PER_DOUBLE == 4 |
| 85 | np = MPZ_NEWALLOC (NUM(dest), 4); |
| 86 | if ((tp[0] | tp[1] | tp[2]) == 0) |
| 87 | np[0] = tp[3], nn = 1; |
| 88 | else if ((tp[0] | tp[1]) == 0) |
| 89 | np[1] = tp[3], np[0] = tp[2], nn = 2; |
| 90 | else if (tp[0] == 0) |
| 91 | np[2] = tp[3], np[1] = tp[2], np[0] = tp[1], nn = 3; |
| 92 | else |
| 93 | np[3] = tp[3], np[2] = tp[2], np[1] = tp[1], np[0] = tp[0], nn = 4; |
| 94 | #endif |
| 95 | #if LIMBS_PER_DOUBLE == 3 |
| 96 | np = MPZ_NEWALLOC (NUM(dest), 3); |
| 97 | if ((tp[0] | tp[1]) == 0) |
| 98 | np[0] = tp[2], nn = 1; |
| 99 | else if (tp[0] == 0) |
| 100 | np[1] = tp[2], np[0] = tp[1], nn = 2; |
| 101 | else |
| 102 | np[2] = tp[2], np[1] = tp[1], np[0] = tp[0], nn = 3; |
| 103 | #endif |
| 104 | #if LIMBS_PER_DOUBLE == 2 |
| 105 | np = MPZ_NEWALLOC (NUM(dest), 2); |
| 106 | if (tp[0] == 0) |
| 107 | np[0] = tp[1], nn = 1; |
| 108 | else |
| 109 | np[1] = tp[1], np[0] = tp[0], nn = 2; |
| 110 | #endif |
| 111 | dn = nn + 1 - exp; |
| 112 | ASSERT (dn > 0); /* -exp >= -1; nn >= 1*/ |
| 113 | dp = MPZ_NEWALLOC (DEN(dest), dn); |
| 114 | MPN_ZERO (dp, dn - 1); |
| 115 | dp[dn - 1] = 1; |
| 116 | count_trailing_zeros (c, np[0] | dp[0]); |
| 117 | if (c != 0) |
| 118 | { |
| 119 | mpn_rshift (np, np, nn, c); |
| 120 | nn -= np[nn - 1] == 0; |
| 121 | --dn; |
| 122 | dp[dn - 1] = CNST_LIMB(1) << (GMP_LIMB_BITS - c); |
| 123 | } |
| 124 | SIZ(DEN(dest)) = dn; |
| 125 | } |
| 126 | else |
| 127 | { |
| 128 | nn = exp; |
| 129 | np = MPZ_NEWALLOC (NUM(dest), nn); |
| 130 | switch (nn) |
| 131 | { |
| 132 | default: |
| 133 | MPN_ZERO (np, nn - LIMBS_PER_DOUBLE); |
| 134 | np += nn - LIMBS_PER_DOUBLE; |
| 135 | /* fall through */ |
| 136 | #if LIMBS_PER_DOUBLE == 2 |
| 137 | case 2: |
| 138 | np[1] = tp[1], np[0] = tp[0]; |
| 139 | break; |
| 140 | #endif |
| 141 | #if LIMBS_PER_DOUBLE == 3 |
| 142 | case 3: |
| 143 | np[2] = tp[2], np[1] = tp[1], np[0] = tp[0]; |
| 144 | break; |
| 145 | case 2: |
| 146 | np[1] = tp[2], np[0] = tp[1]; |
| 147 | break; |
| 148 | #endif |
| 149 | #if LIMBS_PER_DOUBLE == 4 |
| 150 | case 4: |
| 151 | np[3] = tp[3], np[2] = tp[2], np[1] = tp[1], np[0] = tp[0]; |
| 152 | break; |
| 153 | case 3: |
| 154 | np[2] = tp[3], np[1] = tp[2], np[0] = tp[1]; |
| 155 | break; |
| 156 | case 2: |
| 157 | np[1] = tp[3], np[0] = tp[2]; |
| 158 | break; |
| 159 | #endif |
| 160 | } |
| 161 | MPZ_NEWALLOC (DEN(dest), 1)[0] = 1; |
| 162 | SIZ(DEN(dest)) = 1; |
| 163 | } |
| 164 | SIZ(NUM(dest)) = negative ? -nn : nn; |
| 165 | } |