xref: /netbsd-src/external/lgpl3/gmp/dist/mpq/set_d.c (revision 72c7faa4dbb41dbb0238d6b4a109da0d4b236dd4)
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
mpq_set_d(mpq_ptr dest,double d)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 }
166