version 1.1.1.1, 2000/01/10 15:35:23 |
version 1.1.1.3, 2003/08/25 16:06:20 |
|
|
/* mpn/bdivmod.c: mpn_bdivmod for computing U/V mod 2^d. |
/* mpn/bdivmod.c: mpn_bdivmod for computing U/V mod 2^d. |
|
|
Copyright (C) 1991, 1993, 1994, 1995, 1996 Free Software Foundation, Inc. |
Copyright 1991, 1993, 1994, 1995, 1996, 1999, 2000, 2001, 2002 Free Software |
|
Foundation, Inc. |
|
|
This file is part of the GNU MP Library. |
This file is part of the GNU MP Library. |
|
|
The GNU MP Library is free software; you can redistribute it and/or modify |
The GNU MP Library is free software; you can redistribute it and/or modify |
it under the terms of the GNU Library General Public License as published by |
it under the terms of the GNU Lesser General Public License as published by |
the Free Software Foundation; either version 2 of the License, or (at your |
the Free Software Foundation; either version 2.1 of the License, or (at your |
option) any later version. |
option) any later version. |
|
|
The GNU MP Library is distributed in the hope that it will be useful, but |
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 |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |
License for more details. |
License for more details. |
|
|
You should have received a copy of the GNU Library General Public License |
You should have received a copy of the GNU Lesser General Public License |
along with the GNU MP Library; see the file COPYING.LIB. If not, write to |
along with the GNU MP Library; see the file COPYING.LIB. If not, write to |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
MA 02111-1307, USA. */ |
MA 02111-1307, USA. */ |
Line 51 MA 02111-1307, USA. */ |
|
Line 52 MA 02111-1307, USA. */ |
|
#include "gmp-impl.h" |
#include "gmp-impl.h" |
#include "longlong.h" |
#include "longlong.h" |
|
|
|
|
mp_limb_t |
mp_limb_t |
#if __STDC__ |
|
mpn_bdivmod (mp_ptr qp, mp_ptr up, mp_size_t usize, |
mpn_bdivmod (mp_ptr qp, mp_ptr up, mp_size_t usize, |
mp_srcptr vp, mp_size_t vsize, unsigned long int d) |
mp_srcptr vp, mp_size_t vsize, unsigned long int d) |
#else |
|
mpn_bdivmod (qp, up, usize, vp, vsize, d) |
|
mp_ptr qp; |
|
mp_ptr up; |
|
mp_size_t usize; |
|
mp_srcptr vp; |
|
mp_size_t vsize; |
|
unsigned long int d; |
|
#endif |
|
{ |
{ |
/* Cache for v_inv is used to make mpn_accelgcd faster. */ |
mp_limb_t v_inv; |
static mp_limb_t previous_low_vlimb = 0; |
|
static mp_limb_t v_inv; /* 1/V mod 2^BITS_PER_MP_LIMB. */ |
|
|
|
if (vp[0] != previous_low_vlimb) /* Cache miss; compute v_inv. */ |
ASSERT (usize >= 1); |
{ |
ASSERT (vsize >= 1); |
mp_limb_t v = previous_low_vlimb = vp[0]; |
ASSERT (usize * GMP_NUMB_BITS >= d); |
mp_limb_t make_zero = 1; |
ASSERT (! MPN_OVERLAP_P (up, usize, vp, vsize)); |
mp_limb_t two_i = 1; |
ASSERT (! MPN_OVERLAP_P (qp, d/GMP_NUMB_BITS, vp, vsize)); |
v_inv = 0; |
ASSERT (MPN_SAME_OR_INCR2_P (qp, d/GMP_NUMB_BITS, up, usize)); |
do |
ASSERT_MPN (up, usize); |
{ |
ASSERT_MPN (vp, vsize); |
while ((two_i & make_zero) == 0) |
|
two_i <<= 1, v <<= 1; |
|
v_inv += two_i; |
|
make_zero -= v; |
|
} |
|
while (make_zero); |
|
} |
|
|
|
/* Need faster computation for some common cases in mpn_accelgcd. */ |
/* 1/V mod 2^GMP_NUMB_BITS. */ |
|
modlimb_invert (v_inv, vp[0]); |
|
|
|
/* Fast code for two cases previously used by the accel part of mpn_gcd. |
|
(Could probably remove this now it's inlined there.) */ |
if (usize == 2 && vsize == 2 && |
if (usize == 2 && vsize == 2 && |
(d == BITS_PER_MP_LIMB || d == 2*BITS_PER_MP_LIMB)) |
(d == GMP_NUMB_BITS || d == 2*GMP_NUMB_BITS)) |
{ |
{ |
mp_limb_t hi, lo; |
mp_limb_t hi, lo; |
mp_limb_t q = up[0] * v_inv; |
mp_limb_t q = (up[0] * v_inv) & GMP_NUMB_MASK; |
umul_ppmm (hi, lo, q, vp[0]); |
umul_ppmm (hi, lo, q, vp[0] << GMP_NAIL_BITS); |
up[0] = 0, up[1] -= hi + q*vp[1], qp[0] = q; |
up[0] = 0; |
if (d == 2*BITS_PER_MP_LIMB) |
up[1] -= hi + q*vp[1]; |
q = up[1] * v_inv, up[1] = 0, qp[1] = q; |
qp[0] = q; |
|
if (d == 2*GMP_NUMB_BITS) |
|
{ |
|
q = (up[1] * v_inv) & GMP_NUMB_MASK; |
|
up[1] = 0; |
|
qp[1] = q; |
|
} |
return 0; |
return 0; |
} |
} |
|
|
/* Main loop. */ |
/* Main loop. */ |
while (d >= BITS_PER_MP_LIMB) |
while (d >= GMP_NUMB_BITS) |
{ |
{ |
mp_limb_t q = up[0] * v_inv; |
mp_limb_t q = (up[0] * v_inv) & GMP_NUMB_MASK; |
mp_limb_t b = mpn_submul_1 (up, vp, MIN (usize, vsize), q); |
mp_limb_t b = mpn_submul_1 (up, vp, MIN (usize, vsize), q); |
if (usize > vsize) |
if (usize > vsize) |
mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b); |
mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b); |
d -= BITS_PER_MP_LIMB; |
d -= GMP_NUMB_BITS; |
up += 1, usize -= 1; |
up += 1, usize -= 1; |
*qp++ = q; |
*qp++ = q; |
} |
} |
Line 114 mpn_bdivmod (qp, up, usize, vp, vsize, d) |
|
Line 107 mpn_bdivmod (qp, up, usize, vp, vsize, d) |
|
{ |
{ |
mp_limb_t b; |
mp_limb_t b; |
mp_limb_t q = (up[0] * v_inv) & (((mp_limb_t)1<<d) - 1); |
mp_limb_t q = (up[0] * v_inv) & (((mp_limb_t)1<<d) - 1); |
switch (q) |
if (q <= 1) |
{ |
{ |
case 0: return 0; |
if (q == 0) |
case 1: b = mpn_sub_n (up, up, vp, MIN (usize, vsize)); break; |
return 0; |
default: b = mpn_submul_1 (up, vp, MIN (usize, vsize), q); break; |
else |
|
b = mpn_sub_n (up, up, vp, MIN (usize, vsize)); |
} |
} |
|
else |
|
b = mpn_submul_1 (up, vp, MIN (usize, vsize), q); |
|
|
if (usize > vsize) |
if (usize > vsize) |
mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b); |
mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b); |
return q; |
return q; |