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