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-impl.h"
32
33 void
34 mpf_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 }