Return to div_ui.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

Annotation of OpenXM_contrib/gmp/mpfr/div_ui.c, Revision 1.1.1.1

1.1       maekawa     1: /* mpfr_div_ui -- divide a floating-point number by a machine integer
                      2:
                      3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
                      4:
                      5: This file is part of the MPFR Library.
                      6:
                      7: The MPFR 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 MPFR 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 MPFR 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 <stdio.h>
                     23: #include "gmp.h"
                     24: #include "gmp-impl.h"
                     25: #include "longlong.h"
                     26: #include "mpfr.h"
                     27:
                     28: /* #define DEBUG */
                     29:
                     30: /* returns 0 if result exact, non-zero otherwise */
                     31: int
                     32: #ifdef __STDC__
                     33: mpfr_div_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long u, unsigned char rnd_mode)
                     34: #else
                     35: mpfr_div_ui(y, x, u, rnd_mode)
                     36:      mpfr_ptr y;
                     37:      mpfr_srcptr x;
                     38:      unsigned long u;
                     39:      unsigned char rnd_mode;
                     40: #endif
                     41: {
                     42:   int xn, yn, dif, sh, i; mp_limb_t *xp, *yp, *tmp, c, d;
                     43:   TMP_DECL(marker);
                     44:
                     45:   if (FLAG_NAN(x)) { SET_NAN(y); return 1; }
                     46:   if (u==0) { printf("infinity\n"); return 1; }
                     47:
                     48:   TMP_MARK(marker);
                     49:   xn = (PREC(x)-1)/BITS_PER_MP_LIMB + 1;
                     50:   yn = (PREC(y)-1)/BITS_PER_MP_LIMB + 1;
                     51:
                     52:   xp = MANT(x);
                     53:   yp = MANT(y);
                     54:   EXP(y) = EXP(x);
                     55:   if (SIGN(x)!=SIGN(y)) CHANGE_SIGN(y);
                     56:
                     57:   dif = yn+1-xn;
                     58: #ifdef DEBUG
                     59:   printf("dif=%d u=%lu xn=%d\n",dif,u,xn);
                     60:   printf("x="); mpfr_print_raw(x); putchar('\n');
                     61: #endif
                     62:
                     63:   /* we need to store yn+1 = xn + dif limbs of the quotient */
                     64:   if (ABSSIZE(y)>=yn+1) tmp=yp;
                     65:   else tmp=TMP_ALLOC((yn+1)*BYTES_PER_MP_LIMB);
                     66:
                     67:   c = (mp_limb_t) u;
                     68:   if (dif>=0) {
                     69:     /* patch for bug in mpn_divrem_1 */
                     70: #if (UDIV_NEEDS_NORMALIZATION==1)
                     71:     count_leading_zeros(sh, c);
                     72:     c <<= sh;
                     73:     EXP(y) += sh;
                     74: #endif
                     75:     c = mpn_divrem_1(tmp, dif, xp, xn, c);
                     76:   }
                     77:   else /* dif < 0 i.e. xn > yn */
                     78:     c = mpn_divrem_1(tmp, 0, xp-dif, yn, c);
                     79:
                     80:   if (tmp[yn]==0) { tmp--; sh=0; EXP(y) -= mp_bits_per_limb; }
                     81:   /* shift left to normalize */
                     82:   count_leading_zeros(sh, tmp[yn]);
                     83:   if (sh) {
                     84:     mpn_lshift(yp, tmp+1, yn, sh);
                     85:     yp[0] += tmp[0] >> (BITS_PER_MP_LIMB-sh);
                     86:     EXP(y) -= sh;
                     87:   }
                     88:   else MPN_COPY(yp, tmp+1, yn);
                     89: #ifdef DEBUG
                     90: printf("y="); mpfr_print_raw(y); putchar('\n');
                     91: #endif
                     92:
                     93:   sh = yn*BITS_PER_MP_LIMB - PREC(y);
                     94:   /* it remains sh bits in less significant limb of y */
                     95:
                     96:   d = *yp & (((mp_limb_t)1 << sh) - 1);
                     97:   *yp ^= d; /* set to zero lowest sh bits */
                     98:
                     99:   TMP_FREE(marker);
                    100:   if ((c | d)==0) {
                    101:     for (i=0; i<-dif && xp[i]==0; i++);
                    102:     if (i>=-dif) return 0; /* result is exact */
                    103:   }
                    104:
                    105:   switch (rnd_mode) {
                    106:   case GMP_RNDZ:
                    107:     return 1; /* result is inexact */
                    108:   case GMP_RNDU:
                    109:     if (SIGN(y)>0) mpfr_add_one_ulp(y);
                    110:     return 1; /* result is inexact */
                    111:   case GMP_RNDD:
                    112:     if (SIGN(y)<0) mpfr_add_one_ulp(y);
                    113:     return 1; /* result is inexact */
                    114:   case GMP_RNDN:
                    115:     if (d < ((mp_limb_t)1 << (sh-1))) return 1;
                    116:     else if (d > ((mp_limb_t)1 << (sh-1))) {
                    117:       mpfr_add_one_ulp(y);
                    118:     }
                    119:     else { /* d = (mp_limb_t)1 << (sh-1) */
                    120:       if (c) mpfr_add_one_ulp(y);
                    121:       else {
                    122:        for (i=0; i<-dif && xp[i]==0; i++);
                    123:        if (i<-dif) mpfr_add_one_ulp(y);
                    124:        else { /* exactly in the middle */
                    125:          if (*yp & ((mp_limb_t)1 << sh)) mpfr_add_one_ulp(y);
                    126:        }
                    127:       }
                    128:     }
                    129:     return 1;
                    130:   }
                    131:   return 0; /* to prevent warning from gcc */
                    132: }
                    133:
                    134:
                    135: