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-impl.h"
      33  
      34  
      35  /* Not done:
      36  
      37     No attempt is made to identify an overlap u==v.  The result will be
      38     correct (1.0), but a full actual division is done whereas of course
      39     x/x==1 needs no work.  Such a call is not a sensible thing to make, and
      40     it's left to an application to notice and optimize if it might arise
      41     somehow through pointer aliasing or whatever.
      42  
      43     Enhancements:
      44  
      45     The high quotient limb is non-zero when high{up,vsize} >= {vp,vsize}.  We
      46     could make that comparison and use qsize==prec instead of qsize==prec+1,
      47     to save one limb in the division.
      48  
      49     If r==u but the size is enough bigger than prec that there won't be an
      50     overlap between quotient and dividend in mpn_div_q, then we can avoid
      51     copying up,usize.  This would only arise from a prec reduced with
      52     mpf_set_prec_raw and will be pretty unusual, but might be worthwhile if
      53     it could be worked into the copy_u decision cleanly.  */
      54  
      55  void
      56  mpf_div (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
      57  {
      58    mp_srcptr up, vp;
      59    mp_ptr rp, tp, new_vp;
      60    mp_size_t usize, vsize, rsize, prospective_rsize, tsize, zeros;
      61    mp_size_t sign_quotient, prec, high_zero, chop;
      62    mp_exp_t rexp;
      63    int copy_u;
      64    TMP_DECL;
      65  
      66    usize = SIZ(u);
      67    vsize = SIZ(v);
      68  
      69    if (UNLIKELY (vsize == 0))
      70      DIVIDE_BY_ZERO;
      71  
      72    if (usize == 0)
      73      {
      74        SIZ(r) = 0;
      75        EXP(r) = 0;
      76        return;
      77      }
      78  
      79    sign_quotient = usize ^ vsize;
      80    usize = ABS (usize);
      81    vsize = ABS (vsize);
      82    prec = PREC(r);
      83  
      84    TMP_MARK;
      85    rexp = EXP(u) - EXP(v) + 1;
      86  
      87    rp = PTR(r);
      88    up = PTR(u);
      89    vp = PTR(v);
      90  
      91    prospective_rsize = usize - vsize + 1; /* quot from using given u,v sizes */
      92    rsize = prec + 1;			 /* desired quot */
      93  
      94    zeros = rsize - prospective_rsize;	 /* padding u to give rsize */
      95    copy_u = (zeros > 0 || rp == up);	 /* copy u if overlap or padding */
      96  
      97    chop = MAX (-zeros, 0);		 /* negative zeros means shorten u */
      98    up += chop;
      99    usize -= chop;
     100    zeros += chop;			 /* now zeros >= 0 */
     101  
     102    tsize = usize + zeros;		 /* size for possible copy of u */
     103  
     104    /* copy and possibly extend u if necessary */
     105    if (copy_u)
     106      {
     107        tp = TMP_ALLOC_LIMBS (tsize + 1);	/* +1 for mpn_div_q's scratch needs */
     108        MPN_ZERO (tp, zeros);
     109        MPN_COPY (tp+zeros, up, usize);
     110        up = tp;
     111        usize = tsize;
     112      }
     113    else
     114      {
     115        tp = TMP_ALLOC_LIMBS (usize + 1);
     116      }
     117  
     118    /* ensure divisor doesn't overlap quotient */
     119    if (rp == vp)
     120      {
     121        new_vp = TMP_ALLOC_LIMBS (vsize);
     122        MPN_COPY (new_vp, vp, vsize);
     123        vp = new_vp;
     124      }
     125  
     126    ASSERT (usize-vsize+1 == rsize);
     127    mpn_div_q (rp, up, usize, vp, vsize, tp);
     128  
     129    /* strip possible zero high limb */
     130    high_zero = (rp[rsize-1] == 0);
     131    rsize -= high_zero;
     132    rexp -= high_zero;
     133  
     134    SIZ(r) = sign_quotient >= 0 ? rsize : -rsize;
     135    EXP(r) = rexp;
     136    TMP_FREE;
     137  }