1 /* mpf_sub -- Subtract two floats.
2 
3 Copyright 1993-1996, 1999-2002, 2004, 2005, 2011, 2014 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.h"
32 #include "gmp-impl.h"
33 
34 void
mpf_sub(mpf_ptr r,mpf_srcptr u,mpf_srcptr v)35 mpf_sub (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
36 {
37   mp_srcptr up, vp;
38   mp_ptr rp, tp;
39   mp_size_t usize, vsize, rsize;
40   mp_size_t prec;
41   mp_exp_t exp;
42   mp_size_t ediff;
43   int negate;
44   TMP_DECL;
45 
46   usize = SIZ (u);
47   vsize = SIZ (v);
48 
49   /* Handle special cases that don't work in generic code below.  */
50   if (usize == 0)
51     {
52       mpf_neg (r, v);
53       return;
54     }
55   if (vsize == 0)
56     {
57       if (r != u)
58         mpf_set (r, u);
59       return;
60     }
61 
62   /* If signs of U and V are different, perform addition.  */
63   if ((usize ^ vsize) < 0)
64     {
65       __mpf_struct v_negated;
66       v_negated._mp_size = -vsize;
67       v_negated._mp_exp = EXP (v);
68       v_negated._mp_d = PTR (v);
69       mpf_add (r, u, &v_negated);
70       return;
71     }
72 
73   TMP_MARK;
74 
75   /* Signs are now known to be the same.  */
76   negate = usize < 0;
77 
78   /* Make U be the operand with the largest exponent.  */
79   if (EXP (u) < EXP (v))
80     {
81       mpf_srcptr t;
82       t = u; u = v; v = t;
83       negate ^= 1;
84       usize = SIZ (u);
85       vsize = SIZ (v);
86     }
87 
88   usize = ABS (usize);
89   vsize = ABS (vsize);
90   up = PTR (u);
91   vp = PTR (v);
92   rp = PTR (r);
93   prec = PREC (r) + 1;
94   exp = EXP (u);
95   ediff = exp - EXP (v);
96 
97   /* If ediff is 0 or 1, we might have a situation where the operands are
98      extremely close.  We need to scan the operands from the most significant
99      end ignore the initial parts that are equal.  */
100   if (ediff <= 1)
101     {
102       if (ediff == 0)
103 	{
104 	  /* Skip leading limbs in U and V that are equal.  */
105 	      /* This loop normally exits immediately.  Optimize for that.  */
106 	      while (up[usize - 1] == vp[vsize - 1])
107 		{
108 		  usize--;
109 		  vsize--;
110 		  exp--;
111 
112 		  if (usize == 0)
113 		    {
114                       /* u cancels high limbs of v, result is rest of v */
115 		      negate ^= 1;
116                     cancellation:
117                       /* strip high zeros before truncating to prec */
118                       while (vsize != 0 && vp[vsize - 1] == 0)
119                         {
120                           vsize--;
121                           exp--;
122                         }
123 		      if (vsize > prec)
124 			{
125 			  vp += vsize - prec;
126 			  vsize = prec;
127 			}
128                       MPN_COPY_INCR (rp, vp, vsize);
129                       rsize = vsize;
130                       goto done;
131 		    }
132 		  if (vsize == 0)
133 		    {
134                       vp = up;
135                       vsize = usize;
136                       goto cancellation;
137 		    }
138 		}
139 
140 	  if (up[usize - 1] < vp[vsize - 1])
141 	    {
142 	      /* For simplicity, swap U and V.  Note that since the loop above
143 		 wouldn't have exited unless up[usize - 1] and vp[vsize - 1]
144 		 were non-equal, this if-statement catches all cases where U
145 		 is smaller than V.  */
146 	      MPN_SRCPTR_SWAP (up,usize, vp,vsize);
147 	      negate ^= 1;
148 	      /* negating ediff not necessary since it is 0.  */
149 	    }
150 
151 	  /* Check for
152 	     x+1 00000000 ...
153 	      x  ffffffff ... */
154 	  if (up[usize - 1] != vp[vsize - 1] + 1)
155 	    goto general_case;
156 	  usize--;
157 	  vsize--;
158 	  exp--;
159 	}
160       else /* ediff == 1 */
161 	{
162 	  /* Check for
163 	     1 00000000 ...
164 	     0 ffffffff ... */
165 
166 	  if (up[usize - 1] != 1 || vp[vsize - 1] != GMP_NUMB_MAX
167 	      || (usize >= 2 && up[usize - 2] != 0))
168 	    goto general_case;
169 
170 	  usize--;
171 	  exp--;
172 	}
173 
174       /* Skip sequences of 00000000/ffffffff */
175       while (vsize != 0 && usize != 0 && up[usize - 1] == 0
176 	     && vp[vsize - 1] == GMP_NUMB_MAX)
177 	{
178 	  usize--;
179 	  vsize--;
180 	  exp--;
181 	}
182 
183       if (usize == 0)
184 	{
185 	  while (vsize != 0 && vp[vsize - 1] == GMP_NUMB_MAX)
186 	    {
187 	      vsize--;
188 	      exp--;
189 	    }
190 	}
191       else if (usize > prec - 1)
192 	{
193 	  up += usize - (prec - 1);
194 	  usize = prec - 1;
195 	}
196       if (vsize > prec - 1)
197 	{
198 	  vp += vsize - (prec - 1);
199 	  vsize = prec - 1;
200 	}
201 
202       tp = TMP_ALLOC_LIMBS (prec);
203       {
204 	mp_limb_t cy_limb;
205 	if (vsize == 0)
206 	  {
207 	    MPN_COPY (tp, up, usize);
208 	    tp[usize] = 1;
209 	    rsize = usize + 1;
210 	    exp++;
211 	    goto normalized;
212 	  }
213 	if (usize == 0)
214 	  {
215 	    cy_limb = mpn_neg (tp, vp, vsize);
216 	    rsize = vsize;
217 	  }
218 	else if (usize >= vsize)
219 	  {
220 	    /* uuuu     */
221 	    /* vv       */
222 	    mp_size_t size;
223 	    size = usize - vsize;
224 	    MPN_COPY (tp, up, size);
225 	    cy_limb = mpn_sub_n (tp + size, up + size, vp, vsize);
226 	    rsize = usize;
227 	  }
228 	else /* (usize < vsize) */
229 	  {
230 	    /* uuuu     */
231 	    /* vvvvvvv  */
232 	    mp_size_t size;
233 	    size = vsize - usize;
234 	    cy_limb = mpn_neg (tp, vp, size);
235 	    cy_limb = mpn_sub_nc (tp + size, up, vp + size, usize, cy_limb);
236 	    rsize = vsize;
237 	  }
238 	if (cy_limb == 0)
239 	  {
240 	    tp[rsize] = 1;
241 	    rsize++;
242 	    exp++;
243 	    goto normalized;
244 	  }
245 	goto normalize;
246       }
247     }
248 
249 general_case:
250   /* If U extends beyond PREC, ignore the part that does.  */
251   if (usize > prec)
252     {
253       up += usize - prec;
254       usize = prec;
255     }
256 
257   /* If V extends beyond PREC, ignore the part that does.
258      Note that this may make vsize negative.  */
259   if (vsize + ediff > prec)
260     {
261       vp += vsize + ediff - prec;
262       vsize = prec - ediff;
263     }
264 
265   if (ediff >= prec)
266     {
267       /* V completely cancelled.  */
268       if (rp != up)
269 	MPN_COPY (rp, up, usize);
270       rsize = usize;
271     }
272   else
273     {
274       /* Allocate temp space for the result.  Allocate
275 	 just vsize + ediff later???  */
276       tp = TMP_ALLOC_LIMBS (prec);
277 
278       /* Locate the least significant non-zero limb in (the needed
279 	 parts of) U and V, to simplify the code below.  */
280       for (;;)
281 	{
282 	  if (vsize == 0)
283 	    {
284 	      MPN_COPY (rp, up, usize);
285 	      rsize = usize;
286 	      goto done;
287 	    }
288 	  if (vp[0] != 0)
289 	    break;
290 	  vp++, vsize--;
291 	}
292       for (;;)
293 	{
294 	  if (usize == 0)
295 	    {
296 	      MPN_COPY (rp, vp, vsize);
297 	      rsize = vsize;
298 	      negate ^= 1;
299 	      goto done;
300 	    }
301 	  if (up[0] != 0)
302 	    break;
303 	  up++, usize--;
304 	}
305 
306       /* uuuu     |  uuuu     |  uuuu     |  uuuu     |  uuuu    */
307       /* vvvvvvv  |  vv       |    vvvvv  |    v      |       vv */
308 
309       if (usize > ediff)
310 	{
311 	  /* U and V partially overlaps.  */
312 	  if (ediff == 0)
313 	    {
314 	      /* Have to compare the leading limbs of u and v
315 		 to determine whether to compute u - v or v - u.  */
316 	      if (usize >= vsize)
317 		{
318 		  /* uuuu     */
319 		  /* vv       */
320 		  mp_size_t size;
321 		  size = usize - vsize;
322 		  MPN_COPY (tp, up, size);
323 		  mpn_sub_n (tp + size, up + size, vp, vsize);
324 		  rsize = usize;
325 		}
326 	      else /* (usize < vsize) */
327 		{
328 		  /* uuuu     */
329 		  /* vvvvvvv  */
330 		  mp_size_t size;
331 		  size = vsize - usize;
332 		  ASSERT_CARRY (mpn_neg (tp, vp, size));
333 		  mpn_sub_nc (tp + size, up, vp + size, usize, CNST_LIMB (1));
334 		  rsize = vsize;
335 		}
336 	    }
337 	  else
338 	    {
339 	      if (vsize + ediff <= usize)
340 		{
341 		  /* uuuu     */
342 		  /*   v      */
343 		  mp_size_t size;
344 		  size = usize - ediff - vsize;
345 		  MPN_COPY (tp, up, size);
346 		  mpn_sub (tp + size, up + size, usize - size, vp, vsize);
347 		  rsize = usize;
348 		}
349 	      else
350 		{
351 		  /* uuuu     */
352 		  /*   vvvvv  */
353 		  mp_size_t size;
354 		  rsize = vsize + ediff;
355 		  size = rsize - usize;
356 		  ASSERT_CARRY (mpn_neg (tp, vp, size));
357 		  mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
358 		  /* Should we use sub_nc then sub_1? */
359 		  MPN_DECR_U (tp + size, usize, CNST_LIMB (1));
360 		}
361 	    }
362 	}
363       else
364 	{
365 	  /* uuuu     */
366 	  /*      vv  */
367 	  mp_size_t size, i;
368 	  size = vsize + ediff - usize;
369 	  ASSERT_CARRY (mpn_neg (tp, vp, vsize));
370 	  for (i = vsize; i < size; i++)
371 	    tp[i] = GMP_NUMB_MAX;
372 	  mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1);
373 	  rsize = size + usize;
374 	}
375 
376     normalize:
377       /* Full normalize.  Optimize later.  */
378       while (rsize != 0 && tp[rsize - 1] == 0)
379 	{
380 	  rsize--;
381 	  exp--;
382 	}
383     normalized:
384       MPN_COPY (rp, tp, rsize);
385     }
386 
387  done:
388   TMP_FREE;
389   if (rsize == 0) {
390     SIZ (r) = 0;
391     EXP (r) = 0;
392   } else {
393     SIZ (r) = negate ? -rsize : rsize;
394     EXP (r) = exp;
395   }
396 }
397