/* mpz_div -- divide two integers and produce a quotient. Copyright (C) 1991 Free Software Foundation, Inc. This file is part of the GNU MP Library. The GNU MP Library is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. The GNU MP Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with the GNU MP Library; see the file COPYING. If not, write to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */ #include #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" INTEGER2(primIntegerQuot,mpz_div,mpz_div_need) INTEGER2_C(primIntegerQuotC,mpz_div,mpz_div_need) Int mpz_div_need (MP_INT *num, MP_INT *den) { return 2 + CONINFO_LARGESIZEU(num->sizeTag) + EXTRA - CONINFO_LARGESIZEU(den->sizeTag) /* Result */ + CONINFO_LARGESIZEU(den->sizeTag) + EXTRA /* Normalized denominator */ + CONINFO_LARGESIZEU(num->sizeTag) + EXTRA + 1; /* Normalized remainder */ } #if PROFILE static SInfo divProfInfo = { "Integer","Integer.div","Integer.div"}; #endif NodePtr mpz_div (MP_INT *quot, MP_INT *num, MP_INT *den) { mp_srcptr np, dp; mp_ptr qp, rp; mp_size nsize = CONINFO_LARGESIZES(num->sizeTag); mp_size dsize = CONINFO_LARGESIZES(den->sizeTag); mp_size qsize, rsize; mp_size sign_quotient = CONINFO_LARGEEXTRA(num->sizeTag) ^ CONINFO_LARGEEXTRA(den->sizeTag); unsigned normalization_steps; /* Ensure space is enough for quotient. */ qsize = nsize - dsize + 1; /* qsize cannot be bigger than this. */ if (qsize <= 0) { quot->sizeTag = CONSTRW(0,0); INIT_PROFINFO(quot,&divProfInfo) return 1+EXTRA+(NodePtr)quot; } /* GC CHECK */ qp = quot->d; np = num->d; dp = den->d; #if 0 rp = (mp_ptr) alloca ((size_t)((nsize + 1) * BYTES_PER_MP_LIMB)); #else rp = (mp_ptr)(qsize+1+EXTRA+(NodePtr)quot); /* Remainder after quotient */ #endif count_leading_zeros (normalization_steps, dp[dsize - 1]); /* Normalize the denominator and the numerator. */ if (normalization_steps != 0) { mp_ptr tp; mp_limb ndigit; /* Shift up the denominator setting the most significant bit of the most significant word. Use temporary storage not to clobber the original contents of the denominator. */ #if 0 tp = (mp_ptr) alloca ((size_t)(dsize * BYTES_PER_MP_LIMB)); #else tp = (mp_ptr) (nsize+1+(NodePtr)rp); /* Normalized denominator after remainder */ #endif (void) _mpn_lshift (tp, dp, dsize, normalization_steps); dp = tp; /* Shift up the numerator, possibly introducing a new most significant word. Move the shifted numerator in the remainder meanwhile. */ ndigit = _mpn_lshift (rp, np, nsize, normalization_steps); if (ndigit != 0) { rp[nsize] = ndigit; rsize = nsize + 1; } else rsize = nsize; } else { /* The denominator is already normalized, as required. Copy it to temporary space if it overlaps with the quotient. */ if (dp == qp) { fprintf(stderr,"dp == qp in mpz_div!\n"); exit(-1); } /* Move the numerator to the remainder. */ MPN_COPY (rp, np, nsize); rsize = nsize; } qsize = rsize - dsize + _mpn_div (qp, rp, rsize, dp, dsize); /* Normalize the quotient. We may have at most one leading zero-word, so no loop is needed. */ if (qsize > 0) qsize -= (qp[qsize - 1] == 0); quot->sizeTag = CONSTRW(qsize, sign_quotient?1:0); INIT_PROFINFO(quot,&divProfInfo) return 1+qsize+EXTRA+(NodePtr)quot; }