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     ! 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: