contrib/libgmp/mpf/sub.c

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