1/* double mpq_get_d (mpq_t src) -- mpq to double, rounding towards zero.
2
3Copyright 1995, 1996, 2001, 2002, 2003, 2004, 2005 Free Software Foundation,
4Inc.
5
6This file is part of the GNU MP Library.
7
8The GNU MP Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MP Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
20
21#include <stdio.h>  /* for NULL */
22#include "gmp.h"
23#include "gmp-impl.h"
24#include "longlong.h"
25
26
27/* All that's needed is to get the high 53 bits of the quotient num/den,
28   rounded towards zero.  More than 53 bits is fine, any excess is ignored
29   by mpn_get_d.
30
31   N_QLIMBS is how many quotient limbs we need to satisfy the mantissa of a
32   double, assuming the highest of those limbs is non-zero.  The target
33   qsize for mpn_tdiv_qr is then 1 more than this, since that function may
34   give a zero in the high limb (and non-zero in the second highest).
35
36   The use of 8*sizeof(double) in N_QLIMBS is an overestimate of the
37   mantissa bits, but it gets the same result as the true value (53 or 48 or
38   whatever) when rounded up to a multiple of GMP_NUMB_BITS, for non-nails.
39
40   Enhancements:
41
42   Use the true mantissa size in the N_QLIMBS formula, to save a divide step
43   in nails.
44
45   Examine the high limbs of num and den to see if the highest 1 bit of the
46   quotient will fall high enough that just N_QLIMBS-1 limbs is enough to
47   get the necessary bits, thereby saving a division step.
48
49   Bit shift either num or den to arrange for the above condition on the
50   high 1 bit of the quotient, to save a division step always.  A shift to
51   save a division step is definitely worthwhile with mpn_tdiv_qr, though we
52   may want to reassess this on big num/den when a quotient-only division
53   exists.
54
55   Maybe we could estimate the final exponent using nsize-dsize (and
56   possibly the high limbs of num and den), so as to detect overflow and
57   return infinity or zero quickly.  Overflow is never very helpful to an
58   application, and can therefore probably be regarded as abnormal, but we
59   may still like to optimize it if the conditions are easy.  (This would
60   only be for float formats we know, unknown formats are not important and
61   can be left to mpn_get_d.)
62
63   Future:
64
65   If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of
66   padding n with zeros in temporary space.
67
68   If/when a quotient-only division exists it can be used here immediately.
69   remp is only to satisfy mpn_tdiv_qr, the remainder is not used.
70
71   Alternatives:
72
73   An alternative algorithm, that may be faster:
74   0. Let n be somewhat larger than the number of significant bits in a double.
75   1. Extract the most significant n bits of the denominator, and an equal
76      number of bits from the numerator.
77   2. Interpret the extracted numbers as integers, call them a and b
78      respectively, and develop n bits of the fractions ((a + 1) / b) and
79      (a / (b + 1)) using mpn_divrem.
80   3. If the computed values are identical UP TO THE POSITION WE CARE ABOUT,
81      we are done.  If they are different, repeat the algorithm from step 1,
82      but first let n = n * 2.
83   4. If we end up using all bits from the numerator and denominator, fall
84      back to a plain division.
85   5. Just to make life harder, The computation of a + 1 and b + 1 above
86      might give carry-out...  Needs special handling.  It might work to
87      subtract 1 in both cases instead.
88
89   Not certain if this approach would be faster than a quotient-only
90   division.  Presumably such optimizations are the sort of thing we would
91   like to have helping everywhere that uses a quotient-only division. */
92
93double
94mpq_get_d (const MP_RAT *src)
95{
96  double res;
97  mp_srcptr np, dp;
98  mp_ptr remp, tp;
99  mp_size_t nsize = src->_mp_num._mp_size;
100  mp_size_t dsize = src->_mp_den._mp_size;
101  mp_size_t qsize, prospective_qsize, zeros, chop, tsize;
102  mp_size_t sign_quotient = nsize;
103  long exp;
104#define N_QLIMBS (1 + (sizeof (double) + BYTES_PER_MP_LIMB-1) / BYTES_PER_MP_LIMB)
105  mp_limb_t qarr[N_QLIMBS + 1];
106  mp_ptr qp = qarr;
107  TMP_DECL;
108
109  ASSERT (dsize > 0);    /* canonical src */
110
111  /* mpn_get_d below requires a non-zero operand */
112  if (UNLIKELY (nsize == 0))
113    return 0.0;
114
115  TMP_MARK;
116  nsize = ABS (nsize);
117  dsize = ABS (dsize);
118  np = src->_mp_num._mp_d;
119  dp = src->_mp_den._mp_d;
120
121  prospective_qsize = nsize - dsize + 1;   /* from using given n,d */
122  qsize = N_QLIMBS + 1;                    /* desired qsize */
123
124  zeros = qsize - prospective_qsize;       /* padding n to get qsize */
125  exp = (long) -zeros * GMP_NUMB_BITS;     /* relative to low of qp */
126
127  chop = MAX (-zeros, 0);                  /* negative zeros means shorten n */
128  np += chop;
129  nsize -= chop;
130  zeros += chop;                           /* now zeros >= 0 */
131
132  tsize = nsize + zeros;                   /* size for possible copy of n */
133
134  if (WANT_TMP_DEBUG)
135    {
136      /* separate blocks, for malloc debugging */
137      remp = TMP_ALLOC_LIMBS (dsize);
138      tp = (zeros > 0 ? TMP_ALLOC_LIMBS (tsize) : NULL);
139    }
140  else
141    {
142      /* one block with conditionalized size, for efficiency */
143      remp = TMP_ALLOC_LIMBS (dsize + (zeros > 0 ? tsize : 0));
144      tp = remp + dsize;
145    }
146
147  /* zero extend n into temporary space, if necessary */
148  if (zeros > 0)
149    {
150      MPN_ZERO (tp, zeros);
151      MPN_COPY (tp+zeros, np, nsize);
152      np = tp;
153      nsize = tsize;
154    }
155
156  ASSERT (qsize == nsize - dsize + 1);
157  mpn_tdiv_qr (qp, remp, (mp_size_t) 0, np, nsize, dp, dsize);
158
159  /* strip possible zero high limb */
160  qsize -= (qp[qsize-1] == 0);
161
162  res = mpn_get_d (qp, qsize, sign_quotient, exp);
163  TMP_FREE;
164  return res;
165}
166