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