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 } 174