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

     1  /* mpf_div -- Divide two floats.
     2  
     3  Copyright 1993, 1994, 1996, 2000-2002, 2004, 2005, 2010, 2012 Free Software
     4  Foundation, Inc.
     5  
     6  This file is part of the GNU MP Library.
     7  
     8  The GNU MP Library is free software; you can redistribute it and/or modify
     9  it under the terms of either:
    10  
    11    * the GNU Lesser General Public License as published by the Free
    12      Software Foundation; either version 3 of the License, or (at your
    13      option) any later version.
    14  
    15  or
    16  
    17    * the GNU General Public License as published by the Free Software
    18      Foundation; either version 2 of the License, or (at your option) any
    19      later version.
    20  
    21  or both in parallel, as here.
    22  
    23  The GNU MP Library is distributed in the hope that it will be useful, but
    24  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    25  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    26  for more details.
    27  
    28  You should have received copies of the GNU General Public License and the
    29  GNU Lesser General Public License along with the GNU MP Library.  If not,
    30  see https://www.gnu.org/licenses/.  */
    31  
    32  #include "gmp.h"
    33  #include "gmp-impl.h"
    34  
    35  
    36  /* Not done:
    37  
    38     No attempt is made to identify an overlap u==v.  The result will be
    39     correct (1.0), but a full actual division is done whereas of course
    40     x/x==1 needs no work.  Such a call is not a sensible thing to make, and
    41     it's left to an application to notice and optimize if it might arise
    42     somehow through pointer aliasing or whatever.
    43  
    44     Enhancements:
    45  
    46     The high quotient limb is non-zero when high{up,vsize} >= {vp,vsize}.  We
    47     could make that comparison and use qsize==prec instead of qsize==prec+1,
    48     to save one limb in the division.
    49  
    50     If r==u but the size is enough bigger than prec that there won't be an
    51     overlap between quotient and dividend in mpn_div_q, then we can avoid
    52     copying up,usize.  This would only arise from a prec reduced with
    53     mpf_set_prec_raw and will be pretty unusual, but might be worthwhile if
    54     it could be worked into the copy_u decision cleanly.  */
    55  
    56  void
    57  mpf_div (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
    58  {
    59    mp_srcptr up, vp;
    60    mp_ptr rp, tp, new_vp;
    61    mp_size_t usize, vsize, rsize, prospective_rsize, tsize, zeros;
    62    mp_size_t sign_quotient, prec, high_zero, chop;
    63    mp_exp_t rexp;
    64    int copy_u;
    65    TMP_DECL;
    66  
    67    usize = SIZ(u);
    68    vsize = SIZ(v);
    69  
    70    if (UNLIKELY (vsize == 0))
    71      DIVIDE_BY_ZERO;
    72  
    73    if (usize == 0)
    74      {
    75        SIZ(r) = 0;
    76        EXP(r) = 0;
    77        return;
    78      }
    79  
    80    sign_quotient = usize ^ vsize;
    81    usize = ABS (usize);
    82    vsize = ABS (vsize);
    83    prec = PREC(r);
    84  
    85    TMP_MARK;
    86    rexp = EXP(u) - EXP(v) + 1;
    87  
    88    rp = PTR(r);
    89    up = PTR(u);
    90    vp = PTR(v);
    91  
    92    prospective_rsize = usize - vsize + 1; /* quot from using given u,v sizes */
    93    rsize = prec + 1;			 /* desired quot */
    94  
    95    zeros = rsize - prospective_rsize;	 /* padding u to give rsize */
    96    copy_u = (zeros > 0 || rp == up);	 /* copy u if overlap or padding */
    97  
    98    chop = MAX (-zeros, 0);		 /* negative zeros means shorten u */
    99    up += chop;
   100    usize -= chop;
   101    zeros += chop;			 /* now zeros >= 0 */
   102  
   103    tsize = usize + zeros;		 /* size for possible copy of u */
   104  
   105    /* copy and possibly extend u if necessary */
   106    if (copy_u)
   107      {
   108        tp = TMP_ALLOC_LIMBS (tsize + 1);	/* +1 for mpn_div_q's scratch needs */
   109        MPN_ZERO (tp, zeros);
   110        MPN_COPY (tp+zeros, up, usize);
   111        up = tp;
   112        usize = tsize;
   113      }
   114    else
   115      {
   116        tp = TMP_ALLOC_LIMBS (usize + 1);
   117      }
   118  
   119    /* ensure divisor doesn't overlap quotient */
   120    if (rp == vp)
   121      {
   122        new_vp = TMP_ALLOC_LIMBS (vsize);
   123        MPN_COPY (new_vp, vp, vsize);
   124        vp = new_vp;
   125      }
   126  
   127    ASSERT (usize-vsize+1 == rsize);
   128    mpn_div_q (rp, up, usize, vp, vsize, tp);
   129  
   130    /* strip possible zero high limb */
   131    high_zero = (rp[rsize-1] == 0);
   132    rsize -= high_zero;
   133    rexp -= high_zero;
   134  
   135    SIZ(r) = sign_quotient >= 0 ? rsize : -rsize;
   136    EXP(r) = rexp;
   137    TMP_FREE;
   138  }