Return to dc_divrem_n.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpn / generic

Annotation of OpenXM_contrib/gmp/mpn/generic/dc_divrem_n.c, Revision 1.1

1.1     ! ohara       1: /* mpn_dc_divrem_n and auxilliary routines.
        !             2:
        !             3:    THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE
        !             4:    INTERFACES.  IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.
        !             5:    IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A
        !             6:    FUTURE GNU MP RELEASE.
        !             7:
        !             8:
        !             9: Copyright 2000, 2001, 2002 Free Software Foundation, Inc.
        !            10: Contributed by Paul Zimmermann.
        !            11:
        !            12: This file is part of the GNU MP Library.
        !            13:
        !            14: The GNU MP Library is free software; you can redistribute it and/or modify
        !            15: it under the terms of the GNU Lesser General Public License as published by
        !            16: the Free Software Foundation; either version 2.1 of the License, or (at your
        !            17: option) any later version.
        !            18:
        !            19: The GNU MP Library is distributed in the hope that it will be useful, but
        !            20: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
        !            21: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
        !            22: License for more details.
        !            23:
        !            24: You should have received a copy of the GNU Lesser General Public License
        !            25: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
        !            26: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
        !            27: MA 02111-1307, USA. */
        !            28:
        !            29: #include "gmp.h"
        !            30: #include "gmp-impl.h"
        !            31:
        !            32: /*
        !            33: [1] Fast Recursive Division, by Christoph Burnikel and Joachim Ziegler,
        !            34:     Technical report MPI-I-98-1-022, october 1998.
        !            35:     http://www.mpi-sb.mpg.de/~ziegler/TechRep.ps.gz
        !            36: */
        !            37:
        !            38: static mp_limb_t mpn_dc_div_3_halves_by_2
        !            39:   _PROTO ((mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n));
        !            40:
        !            41:
        !            42: #if 0
        !            43: static
        !            44: unused_mpn_divrem (qp, qxn, np, nn, dp, dn)
        !            45:      mp_ptr qp;
        !            46:      mp_size_t qxn;
        !            47:      mp_ptr np;
        !            48:      mp_size_t nn;
        !            49:      mp_srcptr dp;
        !            50:      mp_size_t dn;
        !            51: {
        !            52:   /* This might be useful: */
        !            53:   if (qxn != 0)
        !            54:     {
        !            55:       mp_limb_t c;
        !            56:       mp_ptr tp = alloca ((nn + qxn) * BYTES_PER_MP_LIMB);
        !            57:       MPN_COPY (tp + qxn - nn, np, nn);
        !            58:       MPN_ZERO (tp, qxn);
        !            59:       c = mpn_divrem (qp, 0L, tp, nn + qxn, dp, dn);
        !            60:       /* Maybe copy proper part of tp to np?  Documentation is unclear about
        !            61:         the returned np value when qxn != 0 */
        !            62:       return c;
        !            63:     }
        !            64: }
        !            65: #endif
        !            66:
        !            67:
        !            68: /* mpn_dc_divrem_n - Implements algorithm of page 8 in [1]: divides (np,2n)
        !            69:    by (dp,n) and puts the quotient in (qp,n), the remainder in (np,n).
        !            70:    Returns most significant limb of the quotient, which is 0 or 1.
        !            71:    Requires that the most significant bit of the divisor is set.  */
        !            72:
        !            73: mp_limb_t
        !            74: mpn_dc_divrem_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n)
        !            75: {
        !            76:   mp_limb_t qhl, cc;
        !            77:
        !            78:   if (n % 2 != 0)
        !            79:     {
        !            80:       qhl = mpn_dc_divrem_n (qp + 1, np + 2, dp + 1, n - 1);
        !            81:       cc = mpn_submul_1 (np + 1, qp + 1, n - 1, dp[0]);
        !            82:       cc = mpn_sub_1 (np + n, np + n, 1, cc);
        !            83:       if (qhl) cc += mpn_sub_1 (np + n, np + n, 1, dp[0]);
        !            84:       while (cc)
        !            85:        {
        !            86:          qhl -= mpn_sub_1 (qp + 1, qp + 1, n - 1, (mp_limb_t) 1);
        !            87:          cc -= mpn_add_n (np + 1, np + 1, dp, n);
        !            88:        }
        !            89:       qhl += mpn_add_1 (qp + 1, qp + 1, n - 1,
        !            90:                        mpn_sb_divrem_mn (qp, np, n + 1, dp, n));
        !            91:     }
        !            92:   else
        !            93:     {
        !            94:       mp_size_t n2 = n/2;
        !            95:       qhl = mpn_dc_div_3_halves_by_2 (qp + n2, np + n2, dp, n2);
        !            96:       qhl += mpn_add_1 (qp + n2, qp + n2, n2,
        !            97:                        mpn_dc_div_3_halves_by_2 (qp, np, dp, n2));
        !            98:     }
        !            99:   return qhl;
        !           100: }
        !           101:
        !           102:
        !           103: /* divides (np, 3n) by (dp, 2n) and puts the quotient in (qp, n),
        !           104:    the remainder in (np, 2n) */
        !           105:
        !           106: static mp_limb_t
        !           107: mpn_dc_div_3_halves_by_2 (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n)
        !           108: {
        !           109:   mp_size_t twon = n + n;
        !           110:   mp_limb_t qhl, cc;
        !           111:   mp_ptr tmp;
        !           112:   TMP_DECL (marker);
        !           113:
        !           114:   TMP_MARK (marker);
        !           115:   if (n < DIV_DC_THRESHOLD)
        !           116:     qhl = mpn_sb_divrem_mn (qp, np + n, twon, dp + n, n);
        !           117:   else
        !           118:     qhl = mpn_dc_divrem_n (qp, np + n, dp + n, n);
        !           119:   tmp = (mp_ptr) TMP_ALLOC (twon * BYTES_PER_MP_LIMB);
        !           120:   mpn_mul_n (tmp, qp, dp, n);
        !           121:   cc = mpn_sub_n (np, np, tmp, twon);
        !           122:   TMP_FREE (marker);
        !           123:   if (qhl) cc += mpn_sub_n (np + n, np + n, dp, n);
        !           124:   while (cc)
        !           125:     {
        !           126:       qhl -= mpn_sub_1 (qp, qp, n, (mp_limb_t) 1);
        !           127:       cc -= mpn_add_n (np, np, dp, twon);
        !           128:     }
        !           129:   return qhl;
        !           130: }