source/libs/gmp/gmp-src/mpf/add.c

   1 /* mpf_add -- Add two floats.
   2
   3 Copyright 1993, 1994, 1996, 2000, 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 "gmp.h"
  32 #include "gmp-impl.h"
  33
  34 void
  35 mpf_add (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
  36 {
  37   mp_srcptr up, vp;
  38   mp_ptr rp, tp;
  39   mp_size_t usize, vsize, rsize;
  40   mp_size_t prec;
  41   mp_exp_t uexp;
  42   mp_size_t ediff;
  43   mp_limb_t cy;
  44   int negate;
  45   TMP_DECL;
  46
  47   usize = u->_mp_size;
  48   vsize = v->_mp_size;
  49
  50   /* Handle special cases that don't work in generic code below.  */
  51   if (usize == 0)
  52     {
  53     set_r_v_maybe:
  54       if (r != v)
  55         mpf_set (r, v);
  56       return;
  57     }
  58   if (vsize == 0)
  59     {
  60       v = u;
  61       goto set_r_v_maybe;
  62     }
  63
  64   /* If signs of U and V are different, perform subtraction.  */
  65   if ((usize ^ vsize) < 0)
  66     {
  67       __mpf_struct v_negated;
  68       v_negated._mp_size = -vsize;
  69       v_negated._mp_exp = v->_mp_exp;
  70       v_negated._mp_d = v->_mp_d;
  71       mpf_sub (r, u, &v_negated);
  72       return;
  73     }
  74
  75   TMP_MARK;
  76
  77   /* Signs are now known to be the same.  */
  78   negate = usize < 0;
  79
  80   /* Make U be the operand with the largest exponent.  */
  81   if (u->_mp_exp < v->_mp_exp)
  82     {
  83       mpf_srcptr t;
  84       t = u; u = v; v = t;
  85       usize = u->_mp_size;
  86       vsize = v->_mp_size;
  87     }
  88
  89   usize = ABS (usize);
  90   vsize = ABS (vsize);
  91   up = u->_mp_d;
  92   vp = v->_mp_d;
  93   rp = r->_mp_d;
  94   prec = r->_mp_prec;
  95   uexp = u->_mp_exp;
  96   ediff = u->_mp_exp - v->_mp_exp;
  97
  98   /* If U extends beyond PREC, ignore the part that does.  */
  99   if (usize > prec)
 100     {
 101       up += usize - prec;
 102       usize = prec;
 103     }
 104
 105   /* If V extends beyond PREC, ignore the part that does.
 106      Note that this may make vsize negative.  */
 107   if (vsize + ediff > prec)
 108     {
 109       vp += vsize + ediff - prec;
 110       vsize = prec - ediff;
 111     }
 112
 113 #if 0
 114   /* Locate the least significant non-zero limb in (the needed parts
 115      of) U and V, to simplify the code below.  */
 116   while (up[0] == 0)
 117     up++, usize--;
 118   while (vp[0] == 0)
 119     vp++, vsize--;
 120 #endif
 121
 122   /* Allocate temp space for the result.  Allocate
 123      just vsize + ediff later???  */
 124   tp = TMP_ALLOC_LIMBS (prec);
 125
 126   if (ediff >= prec)
 127     {
 128       /* V completely cancelled.  */
 129       if (rp != up)
 130         MPN_COPY_INCR (rp, up, usize);
 131       rsize = usize;
 132     }
 133   else
 134     {
 135       /* uuuu     |  uuuu     |  uuuu     |  uuuu     |  uuuu    */
 136       /* vvvvvvv  |  vv       |    vvvvv  |    v      |       vv */
 137
 138       if (usize > ediff)
 139         {
 140           /* U and V partially overlaps.  */
 141           if (vsize + ediff <= usize)
 142             {
 143               /* uuuu     */
 144               /*   v      */
 145               mp_size_t size;
 146               size = usize - ediff - vsize;
 147               MPN_COPY (tp, up, size);
 148               cy = mpn_add (tp + size, up + size, usize - size, vp, vsize);
 149               rsize = usize;
 150             }
 151           else
 152             {
 153               /* uuuu     */
 154               /*   vvvvv  */
 155               mp_size_t size;
 156               size = vsize + ediff - usize;
 157               MPN_COPY (tp, vp, size);
 158               cy = mpn_add (tp + size, up, usize, vp + size, usize - ediff);
 159               rsize = vsize + ediff;
 160             }
 161         }
 162       else
 163         {
 164           /* uuuu     */
 165           /*      vv  */
 166           mp_size_t size;
 167           size = vsize + ediff - usize;
 168           MPN_COPY (tp, vp, vsize);
 169           MPN_ZERO (tp + vsize, ediff - usize);
 170           MPN_COPY (tp + size, up, usize);
 171           cy = 0;
 172           rsize = size + usize;
 173         }
 174
 175       MPN_COPY (rp, tp, rsize);
 176       rp[rsize] = cy;
 177       rsize += cy;
 178       uexp += cy;
 179     }
 180
 181   r->_mp_size = negate ? -rsize : rsize;
 182   r->_mp_exp = uexp;
 183   TMP_FREE;
 184 }