source/libs/gmp/gmp-src/mpf/sub.c

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