contrib/gmp/mpf/sub.c

   1 /* mpf_sub -- Subtract two floats.
   2
   3 Copyright 1993, 1994, 1995, 1996, 1999, 2000, 2001, 2002, 2004, 2005 Free
   4 Software Foundation, Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   8 The GNU MP Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MP Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  20
  21 #include "gmp.h"
  22 #include "gmp-impl.h"
  23
  24 void
  25 mpf_sub (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
  26 {
  27   mp_srcptr up, vp;
  28   mp_ptr rp, tp;
  29   mp_size_t usize, vsize, rsize;
  30   mp_size_t prec;
  31   mp_exp_t exp;
  32   mp_size_t ediff;
  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       mpf_neg (r, v);
  43       return;
  44     }
  45   if (vsize == 0)
  46     {
  47       if (r != u)
  48         mpf_set (r, u);
  49       return;
  50     }
  51
  52   /* If signs of U and V are different, perform addition.  */
  53   if ((usize ^ vsize) < 0)
  54     {
  55       __mpf_struct v_negated;
  56       v_negated._mp_size = -vsize;
  57       v_negated._mp_exp = v->_mp_exp;
  58       v_negated._mp_d = v->_mp_d;
  59       mpf_add (r, u, &v_negated);
  60       return;
  61     }
  62
  63   TMP_MARK;
  64
  65   /* Signs are now known to be the same.  */
  66   negate = usize < 0;
  67
  68   /* Make U be the operand with the largest exponent.  */
  69   if (u->_mp_exp < v->_mp_exp)
  70     {
  71       mpf_srcptr t;
  72       t = u; u = v; v = t;
  73       negate ^= 1;
  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 + 1;
  84   exp = u->_mp_exp;
  85   ediff = u->_mp_exp - v->_mp_exp;
  86
  87   /* If ediff is 0 or 1, we might have a situation where the operands are
  88      extremely close.  We need to scan the operands from the most significant
  89      end ignore the initial parts that are equal.  */
  90   if (ediff <= 1)
  91     {
  92       if (ediff == 0)
  93         {
  94           /* Skip leading limbs in U and V that are equal.  */
  95           if (up[usize - 1] == vp[vsize - 1])
  96             {
  97               /* This loop normally exits immediately.  Optimize for that.  */
  98               do
  99                 {
 100                   usize--;
 101                   vsize--;
 102                   exp--;
 103
 104                   if (usize == 0)
 105                     {
 106                       /* u cancels high limbs of v, result is rest of v */
 107                       negate ^= 1;
 108                     cancellation:
 109                       /* strip high zeros before truncating to prec */
 110                       while (vsize != 0 && vp[vsize - 1] == 0)
 111                         {
 112                           vsize--;
 113                           exp--;
 114                         }
 115                       if (vsize > prec)
 116                         {
 117                           vp += vsize - prec;
 118                           vsize = prec;
 119                         }
 120                       MPN_COPY_INCR (rp, vp, vsize);
 121                       rsize = vsize;
 122                       goto done;
 123                     }
 124                   if (vsize == 0)
 125                     {
 126                       vp = up;
 127                       vsize = usize;
 128                       goto cancellation;
 129                     }
 130                 }
 131               while (up[usize - 1] == vp[vsize - 1]);
 132             }
 133
 134           if (up[usize - 1] < vp[vsize - 1])
 135             {
 136               /* For simplicity, swap U and V.  Note that since the loop above
 137                  wouldn't have exited unless up[usize - 1] and vp[vsize - 1]
 138                  were non-equal, this if-statement catches all cases where U
 139                  is smaller than V.  */
 140               MPN_SRCPTR_SWAP (up,usize, vp,vsize);
 141               negate ^= 1;
 142               /* negating ediff not necessary since it is 0.  */
 143             }
 144
 145           /* Check for
 146              x+1 00000000 ...
 147               x  ffffffff ... */
 148           if (up[usize - 1] != vp[vsize - 1] + 1)
 149             goto general_case;
 150           usize--;
 151           vsize--;
 152           exp--;
 153         }
 154       else /* ediff == 1 */
 155         {
 156           /* Check for
 157              1 00000000 ...
 158              0 ffffffff ... */
 159
 160           if (up[usize - 1] != 1 || vp[vsize - 1] != GMP_NUMB_MAX
 161               || (usize >= 2 && up[usize - 2] != 0))
 162             goto general_case;
 163
 164           usize--;
 165           exp--;
 166         }
 167
 168       /* Skip sequences of 00000000/ffffffff */
 169       while (vsize != 0 && usize != 0 && up[usize - 1] == 0
 170              && vp[vsize - 1] == GMP_NUMB_MAX)
 171         {
 172           usize--;
 173           vsize--;
 174           exp--;
 175         }
 176
 177       if (usize == 0)
 178         {
 179           while (vsize != 0 && vp[vsize - 1] == GMP_NUMB_MAX)
 180             {
 181               vsize--;
 182               exp--;
 183             }
 184         }
 185
 186       if (usize > prec - 1)
 187         {
 188           up += usize - (prec - 1);
 189           usize = prec - 1;
 190         }
 191       if (vsize > prec - 1)
 192         {
 193           vp += vsize - (prec - 1);
 194           vsize = prec - 1;
 195         }
 196
 197       tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
 198       {
 199         mp_limb_t cy_limb;
 200         if (vsize == 0)
 201           {
 202             mp_size_t size, i;
 203             size = usize;
 204             for (i = 0; i < size; i++)
 205               tp[i] = up[i];
 206             tp[size] = 1;
 207             rsize = size + 1;
 208             exp++;
 209             goto normalize;
 210           }
 211         if (usize == 0)
 212           {
 213             mp_size_t size, i;
 214             size = vsize;
 215             for (i = 0; i < size; i++)
 216               tp[i] = ~vp[i] & GMP_NUMB_MASK;
 217             cy_limb = 1 - mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1);
 218             rsize = vsize;
 219             if (cy_limb == 0)
 220               {
 221                 tp[rsize] = 1;
 222                 rsize++;
 223                 exp++;
 224               }
 225             goto normalize;
 226           }
 227         if (usize >= vsize)
 228           {
 229             /* uuuu     */
 230             /* vv       */
 231             mp_size_t size;
 232             size = usize - vsize;
 233             MPN_COPY (tp, up, size);
 234             cy_limb = mpn_sub_n (tp + size, up + size, vp, vsize);
 235             rsize = usize;
 236           }
 237         else /* (usize < vsize) */
 238           {
 239             /* uuuu     */
 240             /* vvvvvvv  */
 241             mp_size_t size, i;
 242             size = vsize - usize;
 243             for (i = 0; i < size; i++)
 244               tp[i] = ~vp[i] & GMP_NUMB_MASK;
 245             cy_limb = mpn_sub_n (tp + size, up, vp + size, usize);
 246             cy_limb+= mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
 247             cy_limb-= mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1);
 248             rsize = vsize;
 249           }
 250         if (cy_limb == 0)
 251           {
 252             tp[rsize] = 1;
 253             rsize++;
 254             exp++;
 255           }
 256         goto normalize;
 257       }
 258     }
 259
 260 general_case:
 261   /* If U extends beyond PREC, ignore the part that does.  */
 262   if (usize > prec)
 263     {
 264       up += usize - prec;
 265       usize = prec;
 266     }
 267
 268   /* If V extends beyond PREC, ignore the part that does.
 269      Note that this may make vsize negative.  */
 270   if (vsize + ediff > prec)
 271     {
 272       vp += vsize + ediff - prec;
 273       vsize = prec - ediff;
 274     }
 275
 276   /* Allocate temp space for the result.  Allocate
 277      just vsize + ediff later???  */
 278   tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
 279
 280   if (ediff >= prec)
 281     {
 282       /* V completely cancelled.  */
 283       if (tp != up)
 284         MPN_COPY (rp, up, usize);
 285       rsize = usize;
 286     }
 287   else
 288     {
 289       /* Locate the least significant non-zero limb in (the needed
 290          parts of) U and V, to simplify the code below.  */
 291       for (;;)
 292         {
 293           if (vsize == 0)
 294             {
 295               MPN_COPY (rp, up, usize);
 296               rsize = usize;
 297               goto done;
 298             }
 299           if (vp[0] != 0)
 300             break;
 301           vp++, vsize--;
 302         }
 303       for (;;)
 304         {
 305           if (usize == 0)
 306             {
 307               MPN_COPY (rp, vp, vsize);
 308               rsize = vsize;
 309               negate ^= 1;
 310               goto done;
 311             }
 312           if (up[0] != 0)
 313             break;
 314           up++, usize--;
 315         }
 316
 317       /* uuuu     |  uuuu     |  uuuu     |  uuuu     |  uuuu    */
 318       /* vvvvvvv  |  vv       |    vvvvv  |    v      |       vv */
 319
 320       if (usize > ediff)
 321         {
 322           /* U and V partially overlaps.  */
 323           if (ediff == 0)
 324             {
 325               /* Have to compare the leading limbs of u and v
 326                  to determine whether to compute u - v or v - u.  */
 327               if (usize >= vsize)
 328                 {
 329                   /* uuuu     */
 330                   /* vv       */
 331                   mp_size_t size;
 332                   size = usize - vsize;
 333                   MPN_COPY (tp, up, size);
 334                   mpn_sub_n (tp + size, up + size, vp, vsize);
 335                   rsize = usize;
 336                 }
 337               else /* (usize < vsize) */
 338                 {
 339                   /* uuuu     */
 340                   /* vvvvvvv  */
 341                   mp_size_t size, i;
 342                   size = vsize - usize;
 343                   tp[0] = -vp[0] & GMP_NUMB_MASK;
 344                   for (i = 1; i < size; i++)
 345                     tp[i] = ~vp[i] & GMP_NUMB_MASK;
 346                   mpn_sub_n (tp + size, up, vp + size, usize);
 347                   mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
 348                   rsize = vsize;
 349                 }
 350             }
 351           else
 352             {
 353               if (vsize + ediff <= usize)
 354                 {
 355                   /* uuuu     */
 356                   /*   v      */
 357                   mp_size_t size;
 358                   size = usize - ediff - vsize;
 359                   MPN_COPY (tp, up, size);
 360                   mpn_sub (tp + size, up + size, usize - size, vp, vsize);
 361                   rsize = usize;
 362                 }
 363               else
 364                 {
 365                   /* uuuu     */
 366                   /*   vvvvv  */
 367                   mp_size_t size, i;
 368                   size = vsize + ediff - usize;
 369                   tp[0] = -vp[0] & GMP_NUMB_MASK;
 370                   for (i = 1; i < size; i++)
 371                     tp[i] = ~vp[i] & GMP_NUMB_MASK;
 372                   mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
 373                   mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
 374                   rsize = vsize + ediff;
 375                 }
 376             }
 377         }
 378       else
 379         {
 380           /* uuuu     */
 381           /*      vv  */
 382           mp_size_t size, i;
 383           size = vsize + ediff - usize;
 384           tp[0] = -vp[0] & GMP_NUMB_MASK;
 385           for (i = 1; i < vsize; i++)
 386             tp[i] = ~vp[i] & GMP_NUMB_MASK;
 387           for (i = vsize; i < size; i++)
 388             tp[i] = GMP_NUMB_MAX;
 389           mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1);
 390           rsize = size + usize;
 391         }
 392
 393     normalize:
 394       /* Full normalize.  Optimize later.  */
 395       while (rsize != 0 && tp[rsize - 1] == 0)
 396         {
 397           rsize--;
 398           exp--;
 399         }
 400       MPN_COPY (rp, tp, rsize);
 401     }
 402
 403  done:
 404   r->_mp_size = negate ? -rsize : rsize;
 405   if (rsize == 0)
 406     exp = 0;
 407   r->_mp_exp = exp;
 408   TMP_FREE;
 409 }