github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpq/get_d.c (about)

     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  }