Return to divexact.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/divexact.c, Revision 1.1

1.1     ! maekawa     1: /* mpz_divexact -- finds quotient when known that quot * den == num && den != 0.
        !             2:
        !             3: Copyright (C) 1991, 1993, 1994, 1995 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: /*  Ken Weber (kweber@mat.ufrgs.br, kweber@mcs.kent.edu)
        !            23:
        !            24:     Funding for this work has been partially provided by Conselho Nacional
        !            25:     de Desenvolvimento Cienti'fico e Tecnolo'gico (CNPq) do Brazil, Grant
        !            26:     301314194-2, and was done while I was a visiting reseacher in the Instituto
        !            27:     de Matema'tica at Universidade Federal do Rio Grande do Sul (UFRGS).
        !            28:
        !            29:     References:
        !            30:         T. Jebelean, An algorithm for exact division, Journal of Symbolic
        !            31:         Computation, v. 15, 1993, pp. 169-180.  */
        !            32:
        !            33: #include "gmp.h"
        !            34: #include "gmp-impl.h"
        !            35: #include "longlong.h"
        !            36:
        !            37: void
        !            38: #if __STDC__
        !            39: mpz_divexact (mpz_ptr quot, mpz_srcptr num, mpz_srcptr den)
        !            40: #else
        !            41: mpz_divexact (quot, num, den)
        !            42:      mpz_ptr quot;
        !            43:      mpz_srcptr num;
        !            44:      mpz_srcptr den;
        !            45: #endif
        !            46: {
        !            47:   mp_ptr qp, tp;
        !            48:   mp_size_t qsize, tsize;
        !            49:
        !            50:   mp_srcptr np = num->_mp_d;
        !            51:   mp_srcptr dp = den->_mp_d;
        !            52:   mp_size_t nsize = ABS (num->_mp_size);
        !            53:   mp_size_t dsize = ABS (den->_mp_size);
        !            54:   TMP_DECL (marker);
        !            55:
        !            56:   /*  Generate divide-by-zero error if dsize == 0.  */
        !            57:   if (dsize == 0)
        !            58:     {
        !            59:       quot->_mp_size = 1 / dsize;
        !            60:       return;
        !            61:     }
        !            62:
        !            63:   if (nsize == 0)
        !            64:     {
        !            65:       quot->_mp_size = 0;
        !            66:       return;
        !            67:     }
        !            68:
        !            69:   qsize = nsize - dsize + 1;
        !            70:   if (quot->_mp_alloc < qsize)
        !            71:     _mpz_realloc (quot, qsize);
        !            72:   qp = quot->_mp_d;
        !            73:
        !            74:   TMP_MARK (marker);
        !            75:
        !            76:   /*  QUOT <-- NUM/2^r, T <-- DEN/2^r where = r number of twos in DEN.  */
        !            77:   while (dp[0] == 0)
        !            78:     np += 1, nsize -= 1, dp += 1, dsize -= 1;
        !            79:   tsize = MIN (qsize, dsize);
        !            80:   if (dp[0] & 1)
        !            81:     {
        !            82:       if (qp != dp)
        !            83:        MPN_COPY (qp, np, qsize);
        !            84:       if (qp == dp)            /*  QUOT and DEN overlap.  */
        !            85:        {
        !            86:          tp = (mp_ptr) TMP_ALLOC (sizeof (mp_limb_t) * tsize);
        !            87:          MPN_COPY (tp, dp, tsize);
        !            88:        }
        !            89:       else
        !            90:        tp = (mp_ptr) dp;
        !            91:     }
        !            92:   else
        !            93:     {
        !            94:       unsigned long int r;
        !            95:       tp = (mp_ptr) TMP_ALLOC (sizeof (mp_limb_t) * tsize);
        !            96:       count_trailing_zeros (r, dp[0]);
        !            97:       mpn_rshift (tp, dp, tsize, r);
        !            98:       if (dsize > tsize)
        !            99:        tp[tsize-1] |= dp[tsize] << (BITS_PER_MP_LIMB - r);
        !           100:       mpn_rshift (qp, np, qsize, r);
        !           101:       if (nsize > qsize)
        !           102:        qp[qsize-1] |= np[qsize] << (BITS_PER_MP_LIMB - r);
        !           103:     }
        !           104:
        !           105:   /*  Now QUOT <-- QUOT/T.  */
        !           106:   mpn_bdivmod (qp, qp, qsize, tp, tsize, qsize * BITS_PER_MP_LIMB);
        !           107:   MPN_NORMALIZE (qp, qsize);
        !           108:
        !           109:   quot->_mp_size = (num->_mp_size < 0) == (den->_mp_size < 0) ? qsize : -qsize;
        !           110:
        !           111:   TMP_FREE (marker);
        !           112: }