[BACK]Return to ui_sub.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpf

Annotation of OpenXM_contrib/gmp/mpf/ui_sub.c, Revision 1.1.1.2

1.1       maekawa     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
1.1.1.2 ! maekawa     8: it under the terms of the GNU Lesser General Public License as published by
        !             9: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1       maekawa    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
1.1.1.2 ! maekawa    14: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
1.1       maekawa    15: License for more details.
                     16:
1.1.1.2 ! maekawa    17: You should have received a copy of the GNU Lesser General Public License
1.1       maekawa    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: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>