--- /dev/null
+/* double mpq_get_d (mpq_t src) -- mpq to double, rounding towards zero.
+
+Copyright 1995, 1996, 2001, 2002, 2003, 2004, 2005 Free Software Foundation,
+Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include <stdio.h> /* for NULL */
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+
+/* All that's needed is to get the high 53 bits of the quotient num/den,
+ rounded towards zero. More than 53 bits is fine, any excess is ignored
+ by mpn_get_d.
+
+ N_QLIMBS is how many quotient limbs we need to satisfy the mantissa of a
+ double, assuming the highest of those limbs is non-zero. The target
+ qsize for mpn_tdiv_qr is then 1 more than this, since that function may
+ give a zero in the high limb (and non-zero in the second highest).
+
+ The use of 8*sizeof(double) in N_QLIMBS is an overestimate of the
+ mantissa bits, but it gets the same result as the true value (53 or 48 or
+ whatever) when rounded up to a multiple of GMP_NUMB_BITS, for non-nails.
+
+ Enhancements:
+
+ Use the true mantissa size in the N_QLIMBS formala, to save a divide step
+ in nails.
+
+ Examine the high limbs of num and den to see if the highest 1 bit of the
+ quotient will fall high enough that just N_QLIMBS-1 limbs is enough to
+ get the necessary bits, thereby saving a division step.
+
+ Bit shift either num or den to arrange for the above condition on the
+ high 1 bit of the quotient, to save a division step always. A shift to
+ save a division step is definitely worthwhile with mpn_tdiv_qr, though we
+ may want to reassess this on big num/den when a quotient-only division
+ exists.
+
+ Maybe we could estimate the final exponent using nsize-dsize (and
+ possibly the high limbs of num and den), so as to detect overflow and
+ return infinity or zero quickly. Overflow is never very helpful to an
+ application, and can therefore probably be regarded as abnormal, but we
+ may still like to optimize it if the conditions are easy. (This would
+ only be for float formats we know, unknown formats are not important and
+ can be left to mpn_get_d.)
+
+ Future:
+
+ If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of
+ padding n with zeros in temporary space.
+
+ If/when a quotient-only division exists it can be used here immediately.
+ remp is only to satisfy mpn_tdiv_qr, the remainder is not used.
+
+ Alternatives:
+
+ An alternative algorithm, that may be faster:
+ 0. Let n be somewhat larger than the number of significant bits in a double.
+ 1. Extract the most significant n bits of the denominator, and an equal
+ number of bits from the numerator.
+ 2. Interpret the extracted numbers as integers, call them a and b
+ respectively, and develop n bits of the fractions ((a + 1) / b) and
+ (a / (b + 1)) using mpn_divrem.
+ 3. If the computed values are identical UP TO THE POSITION WE CARE ABOUT,
+ we are done. If they are different, repeat the algorithm from step 1,
+ but first let n = n * 2.
+ 4. If we end up using all bits from the numerator and denominator, fall
+ back to a plain division.
+ 5. Just to make life harder, The computation of a + 1 and b + 1 above
+ might give carry-out... Needs special handling. It might work to
+ subtract 1 in both cases instead.
+
+ Not certain if this approach would be faster than a quotient-only
+ division. Presumably such optimizations are the sort of thing we would
+ like to have helping everywhere that uses a quotient-only division. */
+
+double
+mpq_get_d (const MP_RAT *src)
+{
+ double res;
+ mp_srcptr np, dp;
+ mp_ptr remp, tp;
+ mp_size_t nsize = src->_mp_num._mp_size;
+ mp_size_t dsize = src->_mp_den._mp_size;
+ mp_size_t qsize, prospective_qsize, zeros, chop, tsize;
+ mp_size_t sign_quotient = nsize;
+ long exp;
+#define N_QLIMBS (1 + (sizeof (double) + BYTES_PER_MP_LIMB-1) / BYTES_PER_MP_LIMB)
+ mp_limb_t qarr[N_QLIMBS + 1];
+ mp_ptr qp = qarr;
+ TMP_DECL;
+
+ ASSERT (dsize > 0); /* canonical src */
+
+ /* mpn_get_d below requires a non-zero operand */
+ if (UNLIKELY (nsize == 0))
+ return 0.0;
+
+ TMP_MARK;
+ nsize = ABS (nsize);
+ dsize = ABS (dsize);
+ np = src->_mp_num._mp_d;
+ dp = src->_mp_den._mp_d;
+
+ prospective_qsize = nsize - dsize + 1; /* from using given n,d */
+ qsize = N_QLIMBS + 1; /* desired qsize */
+
+ zeros = qsize - prospective_qsize; /* padding n to get qsize */
+ exp = (long) -zeros * GMP_NUMB_BITS; /* relative to low of qp */
+
+ chop = MAX (-zeros, 0); /* negative zeros means shorten n */
+ np += chop;
+ nsize -= chop;
+ zeros += chop; /* now zeros >= 0 */
+
+ tsize = nsize + zeros; /* size for possible copy of n */
+
+ if (WANT_TMP_DEBUG)
+ {
+ /* separate blocks, for malloc debugging */
+ remp = TMP_ALLOC_LIMBS (dsize);
+ tp = (zeros > 0 ? TMP_ALLOC_LIMBS (tsize) : NULL);
+ }
+ else
+ {
+ /* one block with conditionalized size, for efficiency */
+ remp = TMP_ALLOC_LIMBS (dsize + (zeros > 0 ? tsize : 0));
+ tp = remp + dsize;
+ }
+
+ /* zero extend n into temporary space, if necessary */
+ if (zeros > 0)
+ {
+ MPN_ZERO (tp, zeros);
+ MPN_COPY (tp+zeros, np, nsize);
+ np = tp;
+ nsize = tsize;
+ }
+
+ ASSERT (qsize == nsize - dsize + 1);
+ mpn_tdiv_qr (qp, remp, (mp_size_t) 0, np, nsize, dp, dsize);
+
+ /* strip possible zero high limb */
+ qsize -= (qp[qsize-1] == 0);
+
+ res = mpn_get_d (qp, qsize, sign_quotient, exp);
+ TMP_FREE;
+ return res;
+}
projects / msp430-gcc.git / blobdiff