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