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

Diff for /OpenXM_contrib/gmp/mpfr/Attic/div_ui.c between version 1.1.1.1 and 1.1.1.2

version 1.1.1.1, 2000/09/09 14:12:19 version 1.1.1.2, 2003/08/25 16:06:06
Line 1 
Line 1 
 /* mpfr_div_ui -- divide a floating-point number by a machine integer  /* mpfr_div_ui -- divide a floating-point number by a machine integer
   
 Copyright (C) 1999 PolKA project, Inria Lorraine and Loria  Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
   
 This file is part of the MPFR Library.  This file is part of the MPFR Library.
   
 The MPFR Library is free software; you can redistribute it and/or modify  The MPFR Library is free software; you can redistribute it and/or modify
 it under the terms of the GNU Library General Public License as published by  it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2 of the License, or (at your  the Free Software Foundation; either version 2.1 of the License, or (at your
 option) any later version.  option) any later version.
   
 The MPFR Library is distributed in the hope that it will be useful, but  The MPFR Library is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 License for more details.  License for more details.
   
 You should have received a copy of the GNU Library General Public License  You should have received a copy of the GNU Lesser General Public License
 along with the MPFR Library; see the file COPYING.LIB.  If not, write to  along with the MPFR Library; see the file COPYING.LIB.  If not, write to
 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,  the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */  MA 02111-1307, USA. */
Line 22  MA 02111-1307, USA. */
Line 22  MA 02111-1307, USA. */
 #include <stdio.h>  #include <stdio.h>
 #include "gmp.h"  #include "gmp.h"
 #include "gmp-impl.h"  #include "gmp-impl.h"
 #include "longlong.h"  #include "longlong.h"
 #include "mpfr.h"  #include "mpfr.h"
   #include "mpfr-impl.h"
   
 /* #define DEBUG */  
   
 /* returns 0 if result exact, non-zero otherwise */  /* returns 0 if result exact, non-zero otherwise */
 int  int
 #ifdef __STDC__  mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode)
 mpfr_div_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long u, unsigned char rnd_mode)  
 #else  
 mpfr_div_ui(y, x, u, rnd_mode)  
      mpfr_ptr y;  
      mpfr_srcptr x;  
      unsigned long u;  
      unsigned char rnd_mode;  
 #endif  
 {  {
   int xn, yn, dif, sh, i; mp_limb_t *xp, *yp, *tmp, c, d;    long int xn, yn, dif, sh, i;
     mp_limb_t *xp, *yp, *tmp, c, d;
     int inexact, middle = 1;
   TMP_DECL(marker);    TMP_DECL(marker);
   
   if (FLAG_NAN(x)) { SET_NAN(y); return 1; }    if (MPFR_IS_NAN(x))
   if (u==0) { printf("infinity\n"); return 1; }      {
         MPFR_SET_NAN(y);
         MPFR_RET_NAN;
       }
   
     MPFR_CLEAR_NAN(y); /* clear NaN flag */
   
     if (MPFR_IS_INF(x))
       {
         MPFR_SET_INF(y);
         MPFR_SET_SAME_SIGN(y, x);
         MPFR_RET(0);
       }
   
     if (u == 0)
       {
         if (MPFR_IS_ZERO(x)) /* 0/0 is NaN */
           {
             MPFR_SET_NAN(y);
             MPFR_RET_NAN;
           }
         else /* x/0 is Inf */
           {
             MPFR_SET_INF(y);
             MPFR_SET_SAME_SIGN(y, x);
             MPFR_RET(0);
           }
       }
   
     MPFR_CLEAR_INF(y);
     MPFR_SET_SAME_SIGN(y, x);
   
     if (MPFR_IS_ZERO(x))
       {
         MPFR_SET_ZERO(y);
         MPFR_RET(0);
       }
   
   TMP_MARK(marker);    TMP_MARK(marker);
   xn = (PREC(x)-1)/BITS_PER_MP_LIMB + 1;    xn = (MPFR_PREC(x) - 1)/BITS_PER_MP_LIMB + 1;
   yn = (PREC(y)-1)/BITS_PER_MP_LIMB + 1;    yn = (MPFR_PREC(y) - 1)/BITS_PER_MP_LIMB + 1;
   
   xp = MANT(x);    xp = MPFR_MANT(x);
   yp = MANT(y);    yp = MPFR_MANT(y);
   EXP(y) = EXP(x);    MPFR_EXP(y) = MPFR_EXP(x);
   if (SIGN(x)!=SIGN(y)) CHANGE_SIGN(y);  
   
   dif = yn+1-xn;    dif = yn + 1 - xn;
 #ifdef DEBUG  
   printf("dif=%d u=%lu xn=%d\n",dif,u,xn);  
   printf("x="); mpfr_print_raw(x); putchar('\n');  
 #endif  
   
   /* we need to store yn+1 = xn + dif limbs of the quotient */    /* we need to store yn+1 = xn + dif limbs of the quotient */
   if (ABSSIZE(y)>=yn+1) tmp=yp;    /* don't use tmp=yp since the mpn_lshift call below requires yp >= tmp+1 */
   else tmp=TMP_ALLOC((yn+1)*BYTES_PER_MP_LIMB);    tmp = TMP_ALLOC((yn + 1) * BYTES_PER_MP_LIMB);
   
   c = (mp_limb_t) u;    c = (mp_limb_t) u;
   if (dif>=0) {    if (dif >= 0)
     /* patch for bug in mpn_divrem_1 */      c = mpn_divrem_1 (tmp, dif, xp, xn, c); /* used all the dividend */
 #if (UDIV_NEEDS_NORMALIZATION==1)    else /* dif < 0 i.e. xn > yn, don't use the (-dif) low limbs from x */
     count_leading_zeros(sh, c);      c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, c);
     c <<= sh;  
     EXP(y) += sh;  
 #endif  
     c = mpn_divrem_1(tmp, dif, xp, xn, c);  
   }  
   else /* dif < 0 i.e. xn > yn */  
     c = mpn_divrem_1(tmp, 0, xp-dif, yn, c);  
   
   if (tmp[yn]==0) { tmp--; sh=0; EXP(y) -= mp_bits_per_limb; }    inexact = (c != 0);
     if (rnd_mode == GMP_RNDN)
       {
         if (2 * c < (mp_limb_t) u)
           middle = -1;
         else if (2 * c > (mp_limb_t) u)
           middle = 1;
         else
           middle = 0; /* exactly in the middle */
       }
     for (i=0; ((inexact == 0) || (middle == 0)) && (i < -dif); i++)
       if (xp[i])
         inexact = middle = 1; /* larger than middle */
   
     if (tmp[yn] == 0) /* high limb is zero */
       {
         tmp--;
         sh = 0;
         MPFR_EXP(y) -= BITS_PER_MP_LIMB;
       }
   
     /* now we have yn limbs starting from tmp[1], with tmp[yn]<>0 */
   
   /* shift left to normalize */    /* shift left to normalize */
   count_leading_zeros(sh, tmp[yn]);    count_leading_zeros(sh, tmp[yn]);
   if (sh) {    if (sh)
     mpn_lshift(yp, tmp+1, yn, sh);      {
     yp[0] += tmp[0] >> (BITS_PER_MP_LIMB-sh);        mpn_lshift (yp, tmp + 1, yn, sh);
     EXP(y) -= sh;        yp[0] += tmp[0] >> (BITS_PER_MP_LIMB - sh);
   }        middle = middle || ((tmp[0] << sh) != 0);
   else MPN_COPY(yp, tmp+1, yn);        inexact = inexact || ((tmp[0] << sh) != 0);
 #ifdef DEBUG        MPFR_EXP(y) -= sh;
 printf("y="); mpfr_print_raw(y); putchar('\n');      }
 #endif    else
       MPN_COPY(yp, tmp + 1, yn);
   
   sh = yn*BITS_PER_MP_LIMB - PREC(y);    sh = yn * BITS_PER_MP_LIMB - MPFR_PREC(y);
   /* it remains sh bits in less significant limb of y */    /* it remains sh bits in less significant limb of y */
   
   d = *yp & (((mp_limb_t)1 << sh) - 1);    d = *yp & ((MP_LIMB_T_ONE << sh) - MP_LIMB_T_ONE);
   *yp ^= d; /* set to zero lowest sh bits */    *yp ^= d; /* set to zero lowest sh bits */
   
   TMP_FREE(marker);    TMP_FREE(marker);
   if ((c | d)==0) {    if ((d == 0) && (inexact == 0))
     for (i=0; i<-dif && xp[i]==0; i++);      return 0; /* result is exact */
     if (i>=-dif) return 0; /* result is exact */  
   }  
   
   switch (rnd_mode) {    switch (rnd_mode)
   case GMP_RNDZ:      {
     return 1; /* result is inexact */      case GMP_RNDZ:
   case GMP_RNDU:        MPFR_RET(-MPFR_SIGN(x)); /* result is inexact */
     if (SIGN(y)>0) mpfr_add_one_ulp(y);  
     return 1; /* result is inexact */      case GMP_RNDU:
   case GMP_RNDD:        if (MPFR_SIGN(y) > 0)
     if (SIGN(y)<0) mpfr_add_one_ulp(y);          mpfr_add_one_ulp (y, rnd_mode);
     return 1; /* result is inexact */        MPFR_RET(1); /* result is inexact */
   case GMP_RNDN:  
     if (d < ((mp_limb_t)1 << (sh-1))) return 1;      case GMP_RNDD:
     else if (d > ((mp_limb_t)1 << (sh-1))) {        if (MPFR_SIGN(y) < 0)
       mpfr_add_one_ulp(y);          mpfr_add_one_ulp (y, rnd_mode);
     }        MPFR_RET(-1); /* result is inexact */
     else { /* d = (mp_limb_t)1 << (sh-1) */  
       if (c) mpfr_add_one_ulp(y);      case GMP_RNDN:
       else {        if (sh && d < (MP_LIMB_T_ONE << (sh - 1)))
         for (i=0; i<-dif && xp[i]==0; i++);          MPFR_RET(-MPFR_SIGN(x));
         if (i<-dif) mpfr_add_one_ulp(y);        else if (sh && d > (MP_LIMB_T_ONE << (sh - 1)))
         else { /* exactly in the middle */          {
           if (*yp & ((mp_limb_t)1 << sh)) mpfr_add_one_ulp(y);            mpfr_add_one_ulp (y, rnd_mode);
             MPFR_RET(MPFR_SIGN(x));
         }          }
       else /* sh = 0 or d = 1 << (sh-1) */
         {
           /* we are in the middle if:
              (a) sh > 0 and inexact == 0
              (b) sh=0 and middle=1
            */
           if ((sh && inexact) || (!sh && (middle > 0)) || (*yp & (MP_LIMB_T_ONE << sh)))
             {
               mpfr_add_one_ulp (y, rnd_mode);
               MPFR_RET(MPFR_SIGN(x));
             }
               else
                 MPFR_RET(-MPFR_SIGN(x));
       }        }
     }      }
     return 1;    MPFR_RET(inexact); /* should never go here */
   }  
   return 0; /* to prevent warning from gcc */  
 }  }
   
   
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.1.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>