Annotation of OpenXM_contrib/gmp/mpq/set_d.c, Revision 1.1.1.1
1.1 maekawa 1: /* mpq_set_d(mpq_t q, double d) -- Set q to d without rounding.
2:
3: Copyright (C) 2000 Free Software Foundation, Inc.
4:
5: This file is part of the GNU MP Library.
6:
7: The GNU MP Library is free software; you can redistribute it and/or modify
8: it under the terms of the GNU Lesser General Public License as published by
9: the Free Software Foundation; either version 2.1 of the License, or (at your
10: option) any later version.
11:
12: The GNU MP Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15: License for more details.
16:
17: You should have received a copy of the GNU Lesser General Public License
18: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include "gmp.h"
23: #include "gmp-impl.h"
24: #include "longlong.h"
25:
26: #if BITS_PER_MP_LIMB != 32 && BITS_PER_MP_LIMB != 64
27: This file does not handle any the used BITS_PER_MP_LIMB.
28: #endif
29:
30: void
31: #if __STDC__
32: mpq_set_d (mpq_ptr dest, double d)
33: #else
34: mpq_set_d (dest, d)
35: mpq_ptr dest;
36: double d;
37: #endif
38: {
39: int negative;
40: mp_exp_t exp;
41: mp_limb_t tp[3];
42: mp_ptr np, dp;
43: mp_size_t nn, dn;
44: int c;
45:
46: negative = d < 0;
47: d = ABS (d);
48:
49: exp = __gmp_extract_double (tp, d);
50:
51: /* There are two main version of the conversion. The `then' arm handles
52: things that have a fractional part, while the `else' part handles
53: only integers. */
54: #if BITS_PER_MP_LIMB == 32
55: if (exp <= 1 || (exp == 2 && tp[0] != 0))
56: #else
57: if (exp <= 1)
58: #endif
59: {
60: if (d == 0.0)
61: {
62: SIZ(&(dest->_mp_num)) = 0;
63: SIZ(&(dest->_mp_den)) = 1;
64: PTR(&(dest->_mp_den))[0] = 1;
65: return;
66: }
67:
68: dn = -exp;
69: if (dest->_mp_num._mp_alloc < 3)
70: _mpz_realloc (&(dest->_mp_num), 3);
71: np = PTR(&(dest->_mp_num));
72: #if BITS_PER_MP_LIMB == 32
73: if ((tp[0] | tp[1]) == 0)
74: np[0] = tp[2], nn = 1;
75: else if (tp[0] == 0)
76: np[1] = tp[2], np[0] = tp[1], nn = 2;
77: else
78: np[2] = tp[2], np[1] = tp[1], np[0] = tp[0], nn = 3;
79: #else
80: if (tp[0] == 0)
81: np[0] = tp[1], nn = 1;
82: else
83: np[1] = tp[1], np[0] = tp[0], nn = 2;
84: #endif
85: dn += nn + 1;
86: if (dest->_mp_den._mp_alloc < dn)
87: _mpz_realloc (&(dest->_mp_den), dn);
88: dp = PTR(&(dest->_mp_den));
89: MPN_ZERO (dp, dn - 1);
90: dp[dn - 1] = 1;
91: count_trailing_zeros (c, np[0] | dp[0]);
92: if (c != 0)
93: {
94: mpn_rshift (np, np, nn, c);
95: nn -= np[nn - 1] == 0;
96: mpn_rshift (dp, dp, dn, c);
97: dn -= dp[dn - 1] == 0;
98: }
99: SIZ(&(dest->_mp_den)) = dn;
100: SIZ(&(dest->_mp_num)) = negative ? -nn : nn;
101: }
102: else
103: {
104: nn = exp;
105: if (dest->_mp_num._mp_alloc < nn)
106: _mpz_realloc (&(dest->_mp_num), nn);
107: np = PTR(&(dest->_mp_num));
108: #if BITS_PER_MP_LIMB == 32
109: switch (nn)
110: {
111: default:
112: MPN_ZERO (np, nn - 3);
113: np += nn - 3;
114: /* fall through */
115: case 3:
116: np[2] = tp[2], np[1] = tp[1], np[0] = tp[0];
117: break;
118: case 2:
119: np[1] = tp[2], np[0] = tp[1];
120: break;
121: }
122: #else
123: switch (nn)
124: {
125: default:
126: MPN_ZERO (np, nn - 2);
127: np += nn - 2;
128: /* fall through */
129: case 2:
130: np[1] = tp[1], np[0] = tp[0];
131: break;
132: }
133: #endif
134: dp = PTR(&(dest->_mp_den));
135: dp[0] = 1;
136: SIZ(&(dest->_mp_den)) = 1;
137: SIZ(&(dest->_mp_num)) = negative ? -nn : nn;
138: }
139: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>