source/libs/gmp/gmp-src/mpq/get_d.c

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