X-Git-Url: https://oss.titaniummirror.com/gitweb?a=blobdiff_plain;f=gmp%2Fmpq%2Fget_d.c;fp=gmp%2Fmpq%2Fget_d.c;h=c8f91ca4d8b97cd38d8f8e7489c672083f5d90f6;hb=6fed43773c9b0ce596dca5686f37ac3fc0fa11c0;hp=0000000000000000000000000000000000000000;hpb=27b11d56b743098deb193d510b337ba22dc52e5c;p=msp430-gcc.git

diff --git a/gmp/mpq/get_d.c b/gmp/mpq/get_d.c
new file mode 100644
index 00000000..c8f91ca4
--- /dev/null
+++ b/gmp/mpq/get_d.c
@@ -0,0 +1,165 @@
+/* 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;
+}