contrib/libgmp/mpf/ui_sub.c

   1 /* mpf_ui_sub -- Subtract a float from an unsigned long int.
   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_ui_sub (mpf_ptr r, unsigned long int u, mpf_srcptr v)
  28 #else
  29 mpf_ui_sub (r, u, v)
  30      mpf_ptr r;
  31      unsigned long int 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 uexp;
  40   mp_size_t ediff;
  41   int negate;
  42   mp_limb_t ulimb;
  43   TMP_DECL (marker);
  44
  45   vsize = v->_mp_size;
  46
  47   /* Handle special cases that don't work in generic code below.  */
  48   if (u == 0)
  49     {
  50       mpf_neg (r, v);
  51       return;
  52     }
  53   if (vsize == 0)
  54     {
  55       mpf_set_ui (r, u);
  56       return;
  57     }
  58
  59   /* If signs of U and V are different, perform addition.  */
  60   if (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_ui (r, &v_negated, u);
  67       return;
  68     }
  69
  70   TMP_MARK (marker);
  71
  72   /* Signs are now known to be the same.  */
  73
  74   ulimb = u;
  75   /* Make U be the operand with the largest exponent.  */
  76   if (1 < v->_mp_exp)
  77     {
  78       negate = 1;
  79       usize = ABS (vsize);
  80       vsize = 1;
  81       up = v->_mp_d;
  82       vp = &ulimb;
  83       rp = r->_mp_d;
  84       prec = r->_mp_prec + 1;
  85       uexp = v->_mp_exp;
  86       ediff = uexp - 1;
  87     }
  88   else
  89     {
  90       negate = 0;
  91       usize = 1;
  92       vsize = ABS (vsize);
  93       up = &ulimb;
  94       vp = v->_mp_d;
  95       rp = r->_mp_d;
  96       prec = r->_mp_prec;
  97       uexp = 1;
  98       ediff = 1 - v->_mp_exp;
  99     }
 100
 101   /* Ignore leading limbs in U and V that are equal.  Doing
 102      this helps increase the precision of the result.  */
 103   if (ediff == 0)
 104     {
 105       /* This loop normally exits immediately.  Optimize for that.  */
 106       for (;;)
 107         {
 108           usize--;
 109           vsize--;
 110           if (up[usize] != vp[vsize])
 111             break;
 112           uexp--;
 113           if (usize == 0)
 114             goto Lu0;
 115           if (vsize == 0)
 116             goto Lv0;
 117         }
 118       usize++;
 119       vsize++;
 120       /* Note that either operand (but not both operands) might now have
 121          leading zero limbs.  It matters only that U is unnormalized if
 122          vsize is now zero, and vice versa.  And it is only in that case
 123          that we have to adjust uexp.  */
 124       if (vsize == 0)
 125       Lv0:
 126         while (usize != 0 && up[usize - 1] == 0)
 127           usize--, uexp--;
 128       if (usize == 0)
 129       Lu0:
 130         while (vsize != 0 && vp[vsize - 1] == 0)
 131           vsize--, uexp--;
 132     }
 133
 134   /* If U extends beyond PREC, ignore the part that does.  */
 135   if (usize > prec)
 136     {
 137       up += usize - prec;
 138       usize = prec;
 139     }
 140
 141   /* If V extends beyond PREC, ignore the part that does.
 142      Note that this may make vsize negative.  */
 143   if (vsize + ediff > prec)
 144     {
 145       vp += vsize + ediff - prec;
 146       vsize = prec - ediff;
 147     }
 148
 149   /* Allocate temp space for the result.  Allocate
 150      just vsize + ediff later???  */
 151   tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
 152
 153   if (ediff >= prec)
 154     {
 155       /* V completely cancelled.  */
 156       if (tp != up)
 157         MPN_COPY (rp, up, usize);
 158       rsize = usize;
 159     }
 160   else
 161     {
 162       /* Locate the least significant non-zero limb in (the needed
 163          parts of) U and V, to simplify the code below.  */
 164       for (;;)
 165         {
 166           if (vsize == 0)
 167             {
 168               MPN_COPY (rp, up, usize);
 169               rsize = usize;
 170               goto done;
 171             }
 172           if (vp[0] != 0)
 173             break;
 174           vp++, vsize--;
 175         }
 176       for (;;)
 177         {
 178           if (usize == 0)
 179             {
 180               MPN_COPY (rp, vp, vsize);
 181               rsize = vsize;
 182               negate ^= 1;
 183               goto done;
 184             }
 185           if (up[0] != 0)
 186             break;
 187           up++, usize--;
 188         }
 189
 190       /* uuuu     |  uuuu     |  uuuu     |  uuuu     |  uuuu    */
 191       /* vvvvvvv  |  vv       |    vvvvv  |    v      |       vv */
 192
 193       if (usize > ediff)
 194         {
 195           /* U and V partially overlaps.  */
 196           if (ediff == 0)
 197             {
 198               /* Have to compare the leading limbs of u and v
 199                  to determine whether to compute u - v or v - u.  */
 200               if (usize > vsize)
 201                 {
 202                   /* uuuu     */
 203                   /* vv       */
 204                   int cmp;
 205                   cmp = mpn_cmp (up + usize - vsize, vp, vsize);
 206                   if (cmp >= 0)
 207                     {
 208                       mp_size_t size;
 209                       size = usize - vsize;
 210                       MPN_COPY (tp, up, size);
 211                       mpn_sub_n (tp + size, up + size, vp, vsize);
 212                       rsize = usize;
 213                     }
 214                   else
 215                     {
 216                       /* vv       */  /* Swap U and V. */
 217                       /* uuuu     */
 218                       mp_size_t size, i;
 219                       size = usize - vsize;
 220                       tp[0] = -up[0];
 221                       for (i = 1; i < size; i++)
 222                         tp[i] = ~up[i];
 223                       mpn_sub_n (tp + size, vp, up + size, vsize);
 224                       mpn_sub_1 (tp + size, tp + size, vsize, (mp_limb_t) 1);
 225                       negate ^= 1;
 226                       rsize = usize;
 227                     }
 228                 }
 229               else if (usize < vsize)
 230                 {
 231                   /* uuuu     */
 232                   /* vvvvvvv  */
 233                   int cmp;
 234                   cmp = mpn_cmp (up, vp + vsize - usize, usize);
 235                   if (cmp > 0)
 236                     {
 237                       mp_size_t size, i;
 238                       size = vsize - usize;
 239                       tp[0] = -vp[0];
 240                       for (i = 1; i < size; i++)
 241                         tp[i] = ~vp[i];
 242                       mpn_sub_n (tp + size, up, vp + size, usize);
 243                       mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
 244                       rsize = vsize;
 245                     }
 246                   else
 247                     {
 248                       /* vvvvvvv  */  /* Swap U and V. */
 249                       /* uuuu     */
 250                       /* This is the only place we can get 0.0.  */
 251                       mp_size_t size;
 252                       size = vsize - usize;
 253                       MPN_COPY (tp, vp, size);
 254                       mpn_sub_n (tp + size, vp + size, up, usize);
 255                       negate ^= 1;
 256                       rsize = vsize;
 257                     }
 258                 }
 259               else
 260                 {
 261                   /* uuuu     */
 262                   /* vvvv     */
 263                   int cmp;
 264                   cmp = mpn_cmp (up, vp + vsize - usize, usize);
 265                   if (cmp > 0)
 266                     {
 267                       mpn_sub_n (tp, up, vp, usize);
 268                       rsize = usize;
 269                     }
 270                   else
 271                     {
 272                       mpn_sub_n (tp, vp, up, usize);
 273                       negate ^= 1;
 274                       rsize = usize;
 275                       /* can give zero */
 276                     }
 277                 }
 278             }
 279           else
 280             {
 281               if (vsize + ediff <= usize)
 282                 {
 283                   /* uuuu     */
 284                   /*   v      */
 285                   mp_size_t size;
 286                   size = usize - ediff - vsize;
 287                   MPN_COPY (tp, up, size);
 288                   mpn_sub (tp + size, up + size, usize - size, vp, vsize);
 289                   rsize = usize;
 290                 }
 291               else
 292                 {
 293                   /* uuuu     */
 294                   /*   vvvvv  */
 295                   mp_size_t size, i;
 296                   size = vsize + ediff - usize;
 297                   tp[0] = -vp[0];
 298                   for (i = 1; i < size; i++)
 299                     tp[i] = ~vp[i];
 300                   mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
 301                   mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
 302                   rsize = vsize + ediff;
 303                 }
 304             }
 305         }
 306       else
 307         {
 308           /* uuuu     */
 309           /*      vv  */
 310           mp_size_t size, i;
 311           size = vsize + ediff - usize;
 312           tp[0] = -vp[0];
 313           for (i = 1; i < vsize; i++)
 314             tp[i] = ~vp[i];
 315           for (i = vsize; i < size; i++)
 316             tp[i] = ~(mp_limb_t) 0;
 317           mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1);
 318           rsize = size + usize;
 319         }
 320
 321       /* Full normalize.  Optimize later.  */
 322       while (rsize != 0 && tp[rsize - 1] == 0)
 323         {
 324           rsize--;
 325           uexp--;
 326         }
 327       MPN_COPY (rp, tp, rsize);
 328     }
 329
 330  done:
 331   r->_mp_size = negate ? -rsize : rsize;
 332   r->_mp_exp = uexp;
 333   TMP_FREE (marker);
 334 }