summaryrefslogtreecommitdiff
path: root/Build/source/libs/gmp/gmp-src/mpn/generic/hgcd2.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/gmp/gmp-src/mpn/generic/hgcd2.c')
-rw-r--r--Build/source/libs/gmp/gmp-src/mpn/generic/hgcd2.c547
1 files changed, 417 insertions, 130 deletions
diff --git a/Build/source/libs/gmp/gmp-src/mpn/generic/hgcd2.c b/Build/source/libs/gmp/gmp-src/mpn/generic/hgcd2.c
index 129637063f9..3fa40122e42 100644
--- a/Build/source/libs/gmp/gmp-src/mpn/generic/hgcd2.c
+++ b/Build/source/libs/gmp/gmp-src/mpn/generic/hgcd2.c
@@ -4,7 +4,8 @@
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 1996, 1998, 2000-2004, 2008, 2012 Free Software Foundation, Inc.
+Copyright 1996, 1998, 2000-2004, 2008, 2012, 2019 Free Software Foundation,
+Inc.
This file is part of the GNU MP Library.
@@ -32,181 +33,466 @@ You should have received copies of the GNU General Public License and the
GNU Lesser General Public License along with the GNU MP Library. If not,
see https://www.gnu.org/licenses/. */
-#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
-#if GMP_NAIL_BITS == 0
+#ifndef HGCD2_DIV1_METHOD
+#define HGCD2_DIV1_METHOD 3
+#endif
+
+#ifndef HGCD2_DIV2_METHOD
+#define HGCD2_DIV2_METHOD 2
+#endif
+
+#if GMP_NAIL_BITS != 0
+#error Nails not implemented
+#endif
+
+#if HAVE_NATIVE_mpn_div_11
+
+#define div1 mpn_div_11
+/* Single-limb division optimized for small quotients.
+ Returned value holds d0 = r, d1 = q. */
+mp_double_limb_t div1 (mp_limb_t, mp_limb_t);
-/* Copied from the old mpn/generic/gcdext.c, and modified slightly to return
- the remainder. */
+#elif HGCD2_DIV1_METHOD == 1
-/* Single-limb division optimized for small quotients. */
-static inline mp_limb_t
-div1 (mp_ptr rp,
- mp_limb_t n0,
- mp_limb_t d0)
+static inline mp_double_limb_t
+div1 (mp_limb_t n0, mp_limb_t d0)
{
- mp_limb_t q = 0;
+ mp_double_limb_t res;
+ res.d1 = n0 / d0;
+ res.d0 = n0 - res.d1 * d0;
+
+ return res;
+}
+
+#elif HGCD2_DIV1_METHOD == 2
+
+static mp_double_limb_t
+div1 (mp_limb_t n0, mp_limb_t d0)
+{
+ mp_double_limb_t res;
+ int ncnt, dcnt, cnt;
+ mp_limb_t q;
+ mp_limb_t mask;
+
+ ASSERT (n0 >= d0);
- if ((mp_limb_signed_t) n0 < 0)
+ count_leading_zeros (ncnt, n0);
+ count_leading_zeros (dcnt, d0);
+ cnt = dcnt - ncnt;
+
+ d0 <<= cnt;
+
+ q = -(mp_limb_t) (n0 >= d0);
+ n0 -= d0 & q;
+ d0 >>= 1;
+ q = -q;
+
+ while (--cnt >= 0)
{
- int cnt;
- for (cnt = 1; (mp_limb_signed_t) d0 >= 0; cnt++)
- {
- d0 = d0 << 1;
- }
+ mask = -(mp_limb_t) (n0 >= d0);
+ n0 -= d0 & mask;
+ d0 >>= 1;
+ q = (q << 1) - mask;
+ }
- q = 0;
- while (cnt)
- {
- q <<= 1;
- if (n0 >= d0)
- {
- n0 = n0 - d0;
- q |= 1;
- }
- d0 = d0 >> 1;
- cnt--;
- }
+ res.d0 = n0;
+ res.d1 = q;
+ return res;
+}
+
+#elif HGCD2_DIV1_METHOD == 3
+
+static inline mp_double_limb_t
+div1 (mp_limb_t n0, mp_limb_t d0)
+{
+ mp_double_limb_t res;
+ if (UNLIKELY ((d0 >> (GMP_LIMB_BITS - 3)) != 0)
+ || UNLIKELY (n0 >= (d0 << 3)))
+ {
+ res.d1 = n0 / d0;
+ res.d0 = n0 - res.d1 * d0;
}
else
{
- int cnt;
- for (cnt = 0; n0 >= d0; cnt++)
- {
- d0 = d0 << 1;
- }
+ mp_limb_t q, mask;
- q = 0;
- while (cnt)
- {
- d0 = d0 >> 1;
- q <<= 1;
- if (n0 >= d0)
- {
- n0 = n0 - d0;
- q |= 1;
- }
- cnt--;
- }
+ d0 <<= 2;
+
+ mask = -(mp_limb_t) (n0 >= d0);
+ n0 -= d0 & mask;
+ q = 4 & mask;
+
+ d0 >>= 1;
+ mask = -(mp_limb_t) (n0 >= d0);
+ n0 -= d0 & mask;
+ q += 2 & mask;
+
+ d0 >>= 1;
+ mask = -(mp_limb_t) (n0 >= d0);
+ n0 -= d0 & mask;
+ q -= mask;
+
+ res.d0 = n0;
+ res.d1 = q;
}
- *rp = n0;
- return q;
+ return res;
+}
+
+#elif HGCD2_DIV1_METHOD == 4
+
+/* Table quotients. We extract the NBITS most significant bits of the
+ numerator limb, and the corresponding bits from the divisor limb, and use
+ these to form an index into the table. This method is probably only useful
+ for short pipelines with slow multiplication.
+
+ Possible improvements:
+
+ * Perhaps extract the highest NBITS of the divisor instead of the same bits
+ as from the numerator. That would require another count_leading_zeros,
+ and a post-multiply shift of the quotient.
+
+ * Compress tables? Their values are tiny, and there are lots of zero
+ entries (which are never used).
+
+ * Round the table entries more cleverly?
+*/
+
+#ifndef NBITS
+#define NBITS 5
+#endif
+
+#if NBITS == 5
+/* This needs full division about 13.2% of the time. */
+static const unsigned char tab[512] = {
+17, 9, 5,4,3,2,2,2,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+18, 9, 6,4,3,2,2,2,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+19,10, 6,4,3,3,2,2,2,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
+20,10, 6,5,3,3,2,2,2,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,
+21,11, 7,5,4,3,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,
+22,11, 7,5,4,3,3,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,
+23,12, 7,5,4,3,3,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,
+24,12, 8,6,4,3,3,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,
+25,13, 8,6,5,4,3,3,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,
+26,13, 8,6,5,4,3,3,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,
+27,14, 9,6,5,4,3,3,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,
+28,14, 9,7,5,4,3,3,3,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,
+29,15,10,7,5,4,4,3,3,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,
+30,15,10,7,6,5,4,3,3,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
+31,16,10,7,6,5,4,3,3,3,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,
+32,16,11,8,6,5,4,3,3,3,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
+};
+#elif NBITS == 6
+/* This needs full division about 9.8% of the time. */
+static const unsigned char tab[2048] = {
+33,17,11, 8, 6, 5,4,4,3,3,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+34,17,11, 8, 6, 5,4,4,3,3,3,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+35,18,12, 9, 7, 5,5,4,3,3,3,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+36,18,12, 9, 7, 6,5,4,3,3,3,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+37,19,13, 9, 7, 6,5,4,4,3,3,3,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+38,19,13, 9, 7, 6,5,4,4,3,3,3,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+39,20,13,10, 7, 6,5,4,4,3,3,3,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+40,20,14,10, 8, 6,5,5,4,3,3,3,3,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+41,21,14,10, 8, 6,5,5,4,4,3,3,3,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+42,21,14,10, 8, 7,6,5,4,4,3,3,3,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+43,22,15,11, 8, 7,6,5,4,4,3,3,3,3,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+44,22,15,11, 9, 7,6,5,4,4,3,3,3,3,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+45,23,15,11, 9, 7,6,5,5,4,4,3,3,3,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+46,23,16,11, 9, 7,6,5,5,4,4,3,3,3,3,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+47,24,16,12, 9, 7,6,5,5,4,4,3,3,3,3,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+48,24,16,12, 9, 8,6,6,5,4,4,3,3,3,3,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+49,25,17,12,10, 8,7,6,5,4,4,4,3,3,3,3,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+50,25,17,13,10, 8,7,6,5,5,4,4,3,3,3,3,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+51,26,18,13,10, 8,7,6,5,5,4,4,3,3,3,3,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
+52,26,18,13,10, 8,7,6,5,5,4,4,3,3,3,3,3,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,
+53,27,18,13,10, 9,7,6,5,5,4,4,4,3,3,3,3,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,
+54,27,19,14,11, 9,7,6,6,5,4,4,4,3,3,3,3,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,
+55,28,19,14,11, 9,7,6,6,5,5,4,4,3,3,3,3,3,2,2,2,2,2,2,2,2,2,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,
+56,28,19,14,11, 9,8,7,6,5,5,4,4,3,3,3,3,3,2,2,2,2,2,2,2,2,2,1,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,
+57,29,20,14,11, 9,8,7,6,5,5,4,4,4,3,3,3,3,2,2,2,2,2,2,2,2,2,2,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,
+58,29,20,15,11, 9,8,7,6,5,5,4,4,4,3,3,3,3,3,2,2,2,2,2,2,2,2,2,1,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,
+59,30,20,15,12,10,8,7,6,5,5,4,4,4,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,
+60,30,21,15,12,10,8,7,6,6,5,5,4,4,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,1,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,
+61,31,21,15,12,10,8,7,6,6,5,5,4,4,4,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,
+62,31,22,16,12,10,9,7,6,6,5,5,4,4,4,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,1,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
+63,32,22,16,13,10,9,7,7,6,5,5,4,4,4,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,
+64,32,22,16,13,10,9,8,7,6,5,5,4,4,4,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,1,
+ 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
+};
+#else
+#error No table for provided NBITS
+#endif
+
+static const unsigned char *tabp = tab - (1 << (NBITS - 1) << NBITS);
+
+static inline mp_double_limb_t
+div1 (mp_limb_t n0, mp_limb_t d0)
+{
+ int ncnt;
+ size_t nbi, dbi;
+ mp_limb_t q0;
+ mp_limb_t r0;
+ mp_limb_t mask;
+ mp_double_limb_t res;
+
+ ASSERT (n0 >= d0); /* Actually only msb position is critical. */
+
+ count_leading_zeros (ncnt, n0);
+ nbi = n0 << ncnt >> (GMP_LIMB_BITS - NBITS);
+ dbi = d0 << ncnt >> (GMP_LIMB_BITS - NBITS);
+
+ q0 = tabp[(nbi << NBITS) + dbi];
+ r0 = n0 - q0 * d0;
+ mask = -(mp_limb_t) (r0 >= d0);
+ q0 -= mask;
+ r0 -= d0 & mask;
+
+ if (UNLIKELY (r0 >= d0))
+ {
+ q0 = n0 / d0;
+ r0 = n0 - q0 * d0;
+ }
+
+ res.d1 = q0;
+ res.d0 = r0;
+ return res;
+}
+
+#elif HGCD2_DIV1_METHOD == 5
+
+/* Table inverses of divisors. We don't bother with suppressing the msb from
+ the tables. We index with the NBITS most significant divisor bits,
+ including the always-set highest bit, but use addressing trickery via tabp
+ to suppress it.
+
+ Possible improvements:
+
+ * Do first multiply using 32-bit operations on 64-bit computers. At least
+ on most Arm64 cores, that uses 3 times less resources. It also saves on
+ many x86-64 processors.
+*/
+
+#ifndef NBITS
+#define NBITS 7
+#endif
+
+#if NBITS == 5
+/* This needs full division about 1.63% of the time. */
+static const unsigned char tab[16] = {
+ 63, 59, 55, 52, 50, 47, 45, 43, 41, 39, 38, 36, 35, 34, 33, 32
+};
+static const unsigned char *tabp = tab - (1 << (NBITS - 1));
+#elif NBITS == 6
+/* This needs full division about 0.93% of the time. */
+static const unsigned char tab[32] = {
+127,123,119,116,112,109,106,104,101, 98, 96, 94, 92, 90, 88, 86,
+ 84, 82, 80, 79, 77, 76, 74, 73, 72, 70, 69, 68, 67, 66, 65, 64
+};
+static const unsigned char *tabp = tab - (1 << (NBITS - 1));
+#elif NBITS == 7
+/* This needs full division about 0.49% of the time. */
+static const unsigned char tab[64] = {
+255,251,247,243,239,236,233,229,226,223,220,217,214,211,209,206,
+203,201,198,196,194,191,189,187,185,183,181,179,177,175,173,171,
+169,167,166,164,162,161,159,158,156,155,153,152,150,149,147,146,
+145,143,142,141,140,139,137,136,135,134,133,132,131,130,129,128
+};
+static const unsigned char *tabp = tab - (1 << (NBITS - 1));
+#elif NBITS == 8
+/* This needs full division about 0.26% of the time. */
+static const unsigned short tab[128] = {
+511,507,503,499,495,491,488,484,480,477,473,470,467,463,460,457,
+454,450,447,444,441,438,435,433,430,427,424,421,419,416,413,411,
+408,406,403,401,398,396,393,391,389,386,384,382,380,377,375,373,
+371,369,367,365,363,361,359,357,355,353,351,349,347,345,343,342,
+340,338,336,335,333,331,329,328,326,325,323,321,320,318,317,315,
+314,312,311,309,308,306,305,303,302,301,299,298,296,295,294,292,
+291,290,288,287,286,285,283,282,281,280,279,277,276,275,274,273,
+272,270,269,268,267,266,265,264,263,262,261,260,259,258,257,256
+};
+static const unsigned short *tabp = tab - (1 << (NBITS - 1));
+#else
+#error No table for provided NBITS
+#endif
+
+static inline mp_double_limb_t
+div1 (mp_limb_t n0, mp_limb_t d0)
+{
+ int ncnt, dcnt;
+ size_t dbi;
+ mp_limb_t inv;
+ mp_limb_t q0;
+ mp_limb_t r0;
+ mp_limb_t mask;
+ mp_double_limb_t res;
+
+ count_leading_zeros (ncnt, n0);
+ count_leading_zeros (dcnt, d0);
+
+ dbi = d0 << dcnt >> (GMP_LIMB_BITS - NBITS);
+ inv = tabp[dbi];
+ q0 = ((n0 << ncnt) >> (NBITS + 1)) * inv >> (GMP_LIMB_BITS - 1 + ncnt - dcnt);
+ r0 = n0 - q0 * d0;
+ mask = -(mp_limb_t) (r0 >= d0);
+ q0 -= mask;
+ r0 -= d0 & mask;
+
+ if (UNLIKELY (r0 >= d0))
+ {
+ q0 = n0 / d0;
+ r0 = n0 - q0 * d0;
+ }
+
+ res.d1 = q0;
+ res.d0 = r0;
+ return res;
}
+#else
+#error Unknown HGCD2_DIV1_METHOD
+#endif
+
+#if HAVE_NATIVE_mpn_div_22
+
+#define div2 mpn_div_22
/* Two-limb division optimized for small quotients. */
-static inline mp_limb_t
+mp_limb_t div2 (mp_ptr, mp_limb_t, mp_limb_t, mp_limb_t, mp_limb_t);
+
+#elif HGCD2_DIV2_METHOD == 1
+
+static mp_limb_t
div2 (mp_ptr rp,
- mp_limb_t nh, mp_limb_t nl,
- mp_limb_t dh, mp_limb_t dl)
+ mp_limb_t n1, mp_limb_t n0,
+ mp_limb_t d1, mp_limb_t d0)
{
- mp_limb_t q = 0;
-
- if ((mp_limb_signed_t) nh < 0)
+ mp_double_limb_t rq = div1 (n1, d1);
+ if (UNLIKELY (rq.d1 > d1))
{
- int cnt;
- for (cnt = 1; (mp_limb_signed_t) dh >= 0; cnt++)
+ mp_limb_t n2, q, t1, t0;
+ int c;
+
+ /* Normalize */
+ count_leading_zeros (c, d1);
+ ASSERT (c > 0);
+
+ n2 = n1 >> (GMP_LIMB_BITS - c);
+ n1 = (n1 << c) | (n0 >> (GMP_LIMB_BITS - c));
+ n0 <<= c;
+ d1 = (d1 << c) | (d0 >> (GMP_LIMB_BITS - c));
+ d0 <<= c;
+
+ udiv_qrnnd (q, n1, n2, n1, d1);
+ umul_ppmm (t1, t0, q, d0);
+ if (t1 > n1 || (t1 == n1 && t0 > n0))
{
- dh = (dh << 1) | (dl >> (GMP_LIMB_BITS - 1));
- dl = dl << 1;
+ ASSERT (q > 0);
+ q--;
+ sub_ddmmss (t1, t0, t1, t0, d1, d0);
}
+ sub_ddmmss (n1, n0, n1, n0, t1, t0);
- while (cnt)
- {
- q <<= 1;
- if (nh > dh || (nh == dh && nl >= dl))
- {
- sub_ddmmss (nh, nl, nh, nl, dh, dl);
- q |= 1;
- }
- dl = (dh << (GMP_LIMB_BITS - 1)) | (dl >> 1);
- dh = dh >> 1;
- cnt--;
- }
+ /* Undo normalization */
+ rp[0] = (n0 >> c) | (n1 << (GMP_LIMB_BITS - c));
+ rp[1] = n1 >> c;
+
+ return q;
}
else
{
- int cnt;
- for (cnt = 0; nh > dh || (nh == dh && nl >= dl); cnt++)
+ mp_limb_t q, t1, t0;
+ n1 = rq.d0;
+ q = rq.d1;
+ umul_ppmm (t1, t0, q, d0);
+ if (UNLIKELY (t1 >= n1) && (t1 > n1 || t0 > n0))
{
- dh = (dh << 1) | (dl >> (GMP_LIMB_BITS - 1));
- dl = dl << 1;
- }
-
- while (cnt)
- {
- dl = (dh << (GMP_LIMB_BITS - 1)) | (dl >> 1);
- dh = dh >> 1;
- q <<= 1;
- if (nh > dh || (nh == dh && nl >= dl))
- {
- sub_ddmmss (nh, nl, nh, nl, dh, dl);
- q |= 1;
- }
- cnt--;
+ ASSERT (q > 0);
+ q--;
+ sub_ddmmss (t1, t0, t1, t0, d1, d0);
}
+ sub_ddmmss (rp[1], rp[0], n1, n0, t1, t0);
+ return q;
}
-
- rp[0] = nl;
- rp[1] = nh;
-
- return q;
}
-#if 0
-/* This div2 uses less branches, but it seems to nevertheless be
- slightly slower than the above code. */
-static inline mp_limb_t
+#elif HGCD2_DIV2_METHOD == 2
+
+/* Bit-wise div2. Relies on fast count_leading_zeros. */
+static mp_limb_t
div2 (mp_ptr rp,
- mp_limb_t nh, mp_limb_t nl,
- mp_limb_t dh, mp_limb_t dl)
+ mp_limb_t n1, mp_limb_t n0,
+ mp_limb_t d1, mp_limb_t d0)
{
mp_limb_t q = 0;
int ncnt;
int dcnt;
- count_leading_zeros (ncnt, nh);
- count_leading_zeros (dcnt, dh);
+ count_leading_zeros (ncnt, n1);
+ count_leading_zeros (dcnt, d1);
dcnt -= ncnt;
- dh = (dh << dcnt) + (-(dcnt > 0) & (dl >> (GMP_LIMB_BITS - dcnt)));
- dl <<= dcnt;
+ d1 = (d1 << dcnt) + (d0 >> 1 >> (GMP_LIMB_BITS - 1 - dcnt));
+ d0 <<= dcnt;
do
{
- mp_limb_t bit;
+ mp_limb_t mask;
q <<= 1;
- if (UNLIKELY (nh == dh))
- bit = (nl >= dl);
+ if (UNLIKELY (n1 == d1))
+ mask = -(n0 >= d0);
else
- bit = (nh > dh);
+ mask = -(n1 > d1);
- q |= bit;
+ q -= mask;
- sub_ddmmss (nh, nl, nh, nl, (-bit) & dh, (-bit) & dl);
+ sub_ddmmss (n1, n0, n1, n0, mask & d1, mask & d0);
- dl = (dh << (GMP_LIMB_BITS - 1)) | (dl >> 1);
- dh = dh >> 1;
+ d0 = (d1 << (GMP_LIMB_BITS - 1)) | (d0 >> 1);
+ d1 = d1 >> 1;
}
while (dcnt--);
- rp[0] = nl;
- rp[1] = nh;
+ rp[0] = n0;
+ rp[1] = n1;
return q;
}
+#else
+#error Unknown HGCD2_DIV2_METHOD
#endif
-#else /* GMP_NAIL_BITS != 0 */
-/* Check all functions for nail support. */
-/* hgcd2 should be defined to take inputs including nail bits, and
- produce a matrix with elements also including nail bits. This is
- necessary, for the matrix elements to be useful with mpn_mul_1,
- mpn_addmul_1 and friends. */
-#error Not implemented
-#endif /* GMP_NAIL_BITS != 0 */
-
/* Reduces a,b until |a-b| (almost) fits in one limb + 1 bit. Constructs
matrix M. Returns 1 if we make progress, i.e. can perform at least
one subtraction. Otherwise returns zero. */
@@ -341,9 +627,8 @@ mpn_hgcd2 (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl,
}
}
- /* NOTE: Since we discard the least significant half limb, we don't
- get a truly maximal M (corresponding to |a - b| <
- 2^{GMP_LIMB_BITS +1}). */
+ /* NOTE: Since we discard the least significant half limb, we don't get a
+ truly maximal M (corresponding to |a - b| < 2^{GMP_LIMB_BITS +1}). */
/* Single precision loop */
for (;;)
{
@@ -361,9 +646,10 @@ mpn_hgcd2 (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl,
}
else
{
- mp_limb_t r;
- mp_limb_t q = div1 (&r, ah, bh);
- ah = r;
+ mp_double_limb_t rq = div1 (ah, bh);
+ mp_limb_t q = rq.d1;
+ ah = rq.d0;
+
if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
{
/* A is too small, but q is correct. */
@@ -390,9 +676,10 @@ mpn_hgcd2 (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl,
}
else
{
- mp_limb_t r;
- mp_limb_t q = div1 (&r, bh, ah);
- bh = r;
+ mp_double_limb_t rq = div1 (bh, ah);
+ mp_limb_t q = rq.d1;
+ bh = rq.d0;
+
if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
{
/* B is too small, but q is correct. */