contrib/gmp/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 the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19
  20 #include "gmp.h"
  21 #include "gmp-impl.h"
  22
  23 void
  24 mpf_add (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
  25 {
  26   mp_srcptr up, vp;
  27   mp_ptr rp, tp;
  28   mp_size_t usize, vsize, rsize;
  29   mp_size_t prec;
  30   mp_exp_t uexp;
  31   mp_size_t ediff;
  32   mp_limb_t cy;
  33   int negate;
  34   TMP_DECL;
  35
  36   usize = u->_mp_size;
  37   vsize = v->_mp_size;
  38
  39   /* Handle special cases that don't work in generic code below.  */
  40   if (usize == 0)
  41     {
  42     set_r_v_maybe:
  43       if (r != v)
  44         mpf_set (r, v);
  45       return;
  46     }
  47   if (vsize == 0)
  48     {
  49       v = u;
  50       goto set_r_v_maybe;
  51     }
  52
  53   /* If signs of U and V are different, perform subtraction.  */
  54   if ((usize ^ vsize) < 0)
  55     {
  56       __mpf_struct v_negated;
  57       v_negated._mp_size = -vsize;
  58       v_negated._mp_exp = v->_mp_exp;
  59       v_negated._mp_d = v->_mp_d;
  60       mpf_sub (r, u, &v_negated);
  61       return;
  62     }
  63
  64   TMP_MARK;
  65
  66   /* Signs are now known to be the same.  */
  67   negate = usize < 0;
  68
  69   /* Make U be the operand with the largest exponent.  */
  70   if (u->_mp_exp < v->_mp_exp)
  71     {
  72       mpf_srcptr t;
  73       t = u; u = v; v = t;
  74       usize = u->_mp_size;
  75       vsize = v->_mp_size;
  76     }
  77
  78   usize = ABS (usize);
  79   vsize = ABS (vsize);
  80   up = u->_mp_d;
  81   vp = v->_mp_d;
  82   rp = r->_mp_d;
  83   prec = r->_mp_prec;
  84   uexp = u->_mp_exp;
  85   ediff = u->_mp_exp - v->_mp_exp;
  86
  87   /* If U extends beyond PREC, ignore the part that does.  */
  88   if (usize > prec)
  89     {
  90       up += usize - prec;
  91       usize = prec;
  92     }
  93
  94   /* If V extends beyond PREC, ignore the part that does.
  95      Note that this may make vsize negative.  */
  96   if (vsize + ediff > prec)
  97     {
  98       vp += vsize + ediff - prec;
  99       vsize = prec - ediff;
 100     }
 101
 102 #if 0
 103   /* Locate the least significant non-zero limb in (the needed parts
 104      of) U and V, to simplify the code below.  */
 105   while (up[0] == 0)
 106     up++, usize--;
 107   while (vp[0] == 0)
 108     vp++, vsize--;
 109 #endif
 110
 111   /* Allocate temp space for the result.  Allocate
 112      just vsize + ediff later???  */
 113   tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
 114
 115   if (ediff >= prec)
 116     {
 117       /* V completely cancelled.  */
 118       if (rp != up)
 119         MPN_COPY_INCR (rp, up, usize);
 120       rsize = usize;
 121     }
 122   else
 123     {
 124       /* uuuu     |  uuuu     |  uuuu     |  uuuu     |  uuuu    */
 125       /* vvvvvvv  |  vv       |    vvvvv  |    v      |       vv */
 126
 127       if (usize > ediff)
 128         {
 129           /* U and V partially overlaps.  */
 130           if (vsize + ediff <= usize)
 131             {
 132               /* uuuu     */
 133               /*   v      */
 134               mp_size_t size;
 135               size = usize - ediff - vsize;
 136               MPN_COPY (tp, up, size);
 137               cy = mpn_add (tp + size, up + size, usize - size, vp, vsize);
 138               rsize = usize;
 139             }
 140           else
 141             {
 142               /* uuuu     */
 143               /*   vvvvv  */
 144               mp_size_t size;
 145               size = vsize + ediff - usize;
 146               MPN_COPY (tp, vp, size);
 147               cy = mpn_add (tp + size, up, usize, vp + size, usize - ediff);
 148               rsize = vsize + ediff;
 149             }
 150         }
 151       else
 152         {
 153           /* uuuu     */
 154           /*      vv  */
 155           mp_size_t size;
 156           size = vsize + ediff - usize;
 157           MPN_COPY (tp, vp, vsize);
 158           MPN_ZERO (tp + vsize, ediff - usize);
 159           MPN_COPY (tp + size, up, usize);
 160           cy = 0;
 161           rsize = size + usize;
 162         }
 163
 164       MPN_COPY (rp, tp, rsize);
 165       rp[rsize] = cy;
 166       rsize += cy;
 167       uexp += cy;
 168     }
 169
 170   r->_mp_size = negate ? -rsize : rsize;
 171   r->_mp_exp = uexp;
 172   TMP_FREE;
 173 }