[BACK]Return to sub.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / gmp-2.0.2 / mpf

Annotation of OpenXM/src/kan96xx/gmp-2.0.2/mpf/sub.c, Revision 1.1.1.1

1.1       maekawa     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: }

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