[Gmp-commit] /var/hg/gmp: New function mpn_div_qr_2.
mercurial at gmplib.org
mercurial at gmplib.org
Sun Mar 20 18:51:58 CET 2011
details: /var/hg/gmp/rev/ca5189aa27ce
changeset: 14072:ca5189aa27ce
user: Niels M?ller <nisse at lysator.liu.se>
date: Sun Mar 20 18:51:36 2011 +0100
description:
New function mpn_div_qr_2.
diffstat:
ChangeLog | 8 ++
configure.in | 1 +
gmp-h.in | 3 +
mpn/generic/div_qr_2.c | 153 +++++++++++++++++++++++++++++++++++++++++++++++++
tests/mpn/t-div.c | 46 ++++++++++++++
5 files changed, 211 insertions(+), 0 deletions(-)
diffs (252 lines):
diff -r 79f0ae1c3b9b -r ca5189aa27ce ChangeLog
--- a/ChangeLog Sun Mar 20 11:57:33 2011 +0100
+++ b/ChangeLog Sun Mar 20 18:51:36 2011 +0100
@@ -1,3 +1,11 @@
+2011-03-20 Niels Möller <nisse at lysator.liu.se>
+
+ * tests/mpn/t-div.c (main): Added tests for mpn_divrem_2 and
+ mpn_div_qr_2.
+
+ * mpn/generic/div_qr_2.c (mpn_div_qr_2): New file and function.
+ Intended to eventually replace divrem_2.
+
2011-03-16 Marc Glisse <marc.glisse at inria.fr>
* gmpxx.h (__gmp_set_expr): Remove broken declarations.
diff -r 79f0ae1c3b9b -r ca5189aa27ce configure.in
--- a/configure.in Sun Mar 20 11:57:33 2011 +0100
+++ b/configure.in Sun Mar 20 18:51:36 2011 +0100
@@ -2530,6 +2530,7 @@
toom_interpolate_5pts toom_interpolate_6pts toom_interpolate_7pts \
toom_interpolate_8pts toom_interpolate_12pts toom_interpolate_16pts \
invertappr invert binvert mulmod_bnm1 sqrmod_bnm1 \
+ div_qr_2 \
sbpi1_div_q sbpi1_div_qr sbpi1_divappr_q \
dcpi1_div_q dcpi1_div_qr dcpi1_divappr_q \
mu_div_qr mu_divappr_q mu_div_q \
diff -r 79f0ae1c3b9b -r ca5189aa27ce gmp-h.in
--- a/gmp-h.in Sun Mar 20 11:57:33 2011 +0100
+++ b/gmp-h.in Sun Mar 20 18:51:36 2011 +0100
@@ -1535,6 +1535,9 @@
#define mpn_divrem_2 __MPN(divrem_2)
__GMP_DECLSPEC mp_limb_t mpn_divrem_2 __GMP_PROTO ((mp_ptr, mp_size_t, mp_ptr, mp_size_t, mp_srcptr));
+#define mpn_div_qr_2 __MPN(div_qr_2)
+__GMP_DECLSPEC mp_limb_t mpn_div_qr_2 __GMP_PROTO ((mp_ptr, mp_ptr, mp_size_t, mp_srcptr));
+
#define mpn_gcd __MPN(gcd)
__GMP_DECLSPEC mp_size_t mpn_gcd __GMP_PROTO ((mp_ptr, mp_ptr, mp_size_t, mp_ptr, mp_size_t));
diff -r 79f0ae1c3b9b -r ca5189aa27ce mpn/generic/div_qr_2.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpn/generic/div_qr_2.c Sun Mar 20 18:51:36 2011 +0100
@@ -0,0 +1,153 @@
+/* mpn_div_qr_2 -- Divide natural numbers, producing both remainder and
+ quotient. The divisor is two limbs.
+
+ Contributed to the GNU project by Torbjorn Granlund and Niels Möller
+
+ THIS FILE CONTAINS INTERNAL FUNCTIONS WITH MUTABLE INTERFACES. IT IS
+ ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS
+ ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP
+ RELEASE.
+
+
+Copyright 1993, 1994, 1995, 1996, 1999, 2000, 2001, 2002, 2011 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 Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, 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 Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+
+static mp_limb_t
+mpn_div_qr_2_norm (mp_ptr qp, mp_ptr np, mp_size_t nn,
+ mp_limb_t d1, mp_limb_t d0, mp_limb_t di)
+{
+ mp_limb_t qh;
+ mp_size_t i;
+ mp_limb_t r1, r0;
+
+ ASSERT (nn >= 2);
+ ASSERT (d1 & GMP_NUMB_HIGHBIT);
+
+ np += nn - 2;
+ r1 = np[1];
+ r0 = np[0];
+
+ qh = 0;
+ if (r1 >= d1 && (r1 > d1 || r0 >= d0))
+ {
+#if GMP_NAIL_BITS == 0
+ sub_ddmmss (r1, r0, r1, r0, d1, d0);
+#else
+ r0 = r0 - d0;
+ r1 = r1 - d1 - (r0 >> GMP_LIMB_BITS - 1);
+ r0 &= GMP_NUMB_MASK;
+#endif
+ qh = 1;
+ }
+
+ for (i = nn - 2 - 1; i >= 0; i--)
+ {
+ mp_limb_t n0, q;
+ n0 = np[-1];
+ udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di);
+ np--;
+ qp[i] = q;
+ }
+
+ np[1] = r1;
+ np[0] = r0;
+
+ return qh;
+}
+
+static mp_limb_t
+mpn_div_qr_2_unnorm (mp_ptr qp, mp_ptr np, mp_size_t nn,
+ mp_limb_t d1, mp_limb_t d0, int shift, mp_limb_t di)
+{
+ mp_limb_t qh;
+ mp_limb_t r2, r1, r0;
+ mp_size_t i;
+
+ ASSERT (nn >= 2);
+ ASSERT (d1 & GMP_NUMB_HIGHBIT);
+ ASSERT (shift > 0);
+
+ r2 = np[nn-1] >> (GMP_LIMB_BITS - shift);
+ r1 = (np[nn-1] << shift) | (np[nn-2] >> (GMP_LIMB_BITS - shift));
+ r0 = np[nn-2] << shift;
+
+ udiv_qr_3by2 (qh, r2, r1, r2, r1, r0, d1, d0, di);
+
+ for (i = nn - 2 - 1; i >= 0; i--)
+ {
+ mp_limb_t q;
+ r0 = np[i];
+ r1 |= r0 >> (GMP_LIMB_BITS - shift);
+ r0 <<= shift;
+ udiv_qr_3by2 (q, r2, r1, r2, r1, r0, d1, d0, di);
+ qp[i] = q;
+ }
+
+ np[0] = (r1 >> shift) | (r2 << (GMP_LIMB_BITS - shift));
+ np[1] = r2 >> shift;
+ return qh;
+}
+
+/* Divide num {np,nn} by den {dp,2} and write the nn-2 least
+ significant quotient limbs at qp and the 2 long remainder at np.
+ Return the most significant limb of the quotient.
+
+ Preconditions:
+ 1. qp must either not overlap with the input operands at all, or
+ qp >= np + 2 must hold true. (This means that it's possible to put
+ the quotient in the high part of {np,nn}, right above the remainder.
+ 2. nn >= 2. */
+
+mp_limb_t
+mpn_div_qr_2 (mp_ptr qp, mp_ptr np, mp_size_t nn,
+ mp_srcptr dp)
+{
+ mp_limb_t d1;
+ mp_limb_t d0;
+ gmp_pi1_t dinv;
+
+ ASSERT (nn >= 2);
+ ASSERT (! MPN_OVERLAP_P (qp, nn-2, np, nn) || qp >= np + 2);
+ ASSERT_MPN (np, nn);
+ ASSERT_MPN (dp, 2);
+
+ d1 = dp[1]; d0 = dp[0];
+
+ ASSERT (d1 > 0);
+
+ if (UNLIKELY (d1 & GMP_NUMB_HIGHBIT))
+ {
+ gmp_pi1_t dinv;
+ invert_pi1 (dinv, d1, d0);
+ return mpn_div_qr_2_norm (qp, np, nn, d1, d0, dinv.inv32);
+ }
+ else
+ {
+ int shift;
+ count_leading_zeros (shift, d1);
+ d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
+ d0 <<= shift;
+ invert_pi1 (dinv, d1, d0);
+ return mpn_div_qr_2_unnorm (qp, np, nn, d1, d0, shift, dinv.inv32);
+ }
+}
diff -r 79f0ae1c3b9b -r ca5189aa27ce tests/mpn/t-div.c
--- a/tests/mpn/t-div.c Sun Mar 20 11:57:33 2011 +0100
+++ b/tests/mpn/t-div.c Sun Mar 20 18:51:36 2011 +0100
@@ -410,6 +410,52 @@
check_one (qp, NULL, np, nn, dp, dn, "mpn_div_q", 0);
}
+ if (dn >= 2 && nn >= 2)
+ {
+ mp_limb_t qh;
+
+ /* mpn_divrem_2 */
+ MPN_COPY (rp, np, nn);
+ qp[nn - 2] = qp[nn-1] = qran1;
+
+ qh = mpn_divrem_2 (qp, 0, rp, nn, dp + dn - 2);
+ ASSERT_ALWAYS (qp[nn - 2] == qran1);
+ ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1);
+ qp[nn - 2] = qh;
+
+ check_one (qp, rp, np, nn, dp + dn - 2, 2, "mpn_divrem_2", 0);
+
+ /* Missing: divrem_2 with fraction limbs. */
+
+ /* mpn_div_qr_2 (normalized) */
+ MPN_COPY (rp, np, nn);
+ qp[nn - 2] = qran1;
+
+ qh = mpn_div_qr_2 (qp, rp, nn, dp + dn - 2);
+ ASSERT_ALWAYS (qp[nn - 2] == qran1);
+ ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1);
+ qp[nn - 2] = qh;
+
+ check_one (qp, rp, np, nn, dp + dn - 2, 2, "mpn_div_qr_2 (normalized)", 0);
+
+ /* mpn_div_qr_2 (unnormalized) */
+ dp[dn - 1] &= ~GMP_NUMB_HIGHBIT;
+ if (dp[dn - 1] == 0)
+ continue;
+
+ MPN_COPY (rp, np, nn);
+ qp[nn - 2] = qran1;
+
+ qh = mpn_div_qr_2 (qp, rp, nn, dp + dn - 2);
+ ASSERT_ALWAYS (qp[nn - 2] == qran1);
+ ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1);
+ qp[nn - 2] = qh;
+
+ check_one (qp, rp, np, nn, dp + dn - 2, 2, "mpn_div_qr_2 (unnormalized)", 0);
+
+ qp[nn - dn + 1] = qran1;
+ }
+
/* Finally, test mpn_div_q without msb set. */
dp[dn - 1] &= ~GMP_NUMB_HIGHBIT;
if (dp[dn - 1] == 0)
More information about the gmp-commit
mailing list