[Gmp-commit] /var/hg/gmp: hgcd simplifications. New file hgcd_matrix.c
mercurial at gmplib.org
mercurial at gmplib.org
Thu May 19 21:18:53 CEST 2011
details: /var/hg/gmp/rev/bd597269199a
changeset: 14188:bd597269199a
user: Niels M?ller <nisse at lysator.liu.se>
date: Thu May 19 21:09:39 2011 +0200
description:
hgcd simplifications. New file hgcd_matrix.c
diffstat:
ChangeLog | 14 +
configure.in | 2 +-
gmp-impl.h | 6 +
mpn/generic/hgcd.c | 396 ++-------------------------------------------
mpn/generic/hgcd_matrix.c | 253 +++++++++++++++++++++++++++++
5 files changed, 293 insertions(+), 378 deletions(-)
diffs (truncated from 747 to 300 lines):
diff -r f901b9fcfdd3 -r bd597269199a ChangeLog
--- a/ChangeLog Thu May 19 20:55:23 2011 +0200
+++ b/ChangeLog Thu May 19 21:09:39 2011 +0200
@@ -1,5 +1,19 @@
2011-05-19 Niels Möller <nisse at lysator.liu.se>
+ * configure.in (gmp_mpn_functions): Added hgcd_matrix.
+
+ * mpn/generic/hgcd.c (hgcd_matrix_update_1): Deleted. Several other
+ helper functions moved to hgcd_matrix.c, see below.
+ (hgcd_hook): New function.
+ (hgcd_step): Simplified, using mpn_gcd_subdiv_step and hgcd_hook.
+
+ * mpn/generic/hgcd_matrix.c: New file.
+ (mpn_hgcd_matrix_init): Moved here, from hgcd.c.
+ (mpn_hgcd_matrix_update_q): Likewise.
+ (mpn_hgcd_matrix_mul_1): Likewise.
+ (mpn_hgcd_matrix_mul): Likewise.
+ (mpn_hgcd_matrix_adjust): Likewise.
+
* mpn/generic/gcd_subdiv_step.c (mpn_gcd_subdiv_step): New
argument s, for use by hgcd.
* gmp-impl.h (mpn_gcd_subdiv_step): Update declaration.
diff -r f901b9fcfdd3 -r bd597269199a configure.in
--- a/configure.in Thu May 19 20:55:23 2011 +0200
+++ b/configure.in Thu May 19 21:09:39 2011 +0200
@@ -2540,7 +2540,7 @@
gcdext_lehmer \
div_q tdiv_qr jacbase jacobi_lehmer get_d \
matrix22_mul matrix22_mul1_inverse_vector \
- hgcd2 hgcd mullo_n mullo_basecase \
+ hgcd_matrix hgcd2 hgcd mullo_n mullo_basecase \
toom22_mul toom32_mul toom42_mul toom52_mul toom62_mul \
toom33_mul toom43_mul toom53_mul toom63_mul \
toom44_mul \
diff -r f901b9fcfdd3 -r bd597269199a gmp-impl.h
--- a/gmp-impl.h Thu May 19 20:55:23 2011 +0200
+++ b/gmp-impl.h Thu May 19 21:09:39 2011 +0200
@@ -3926,6 +3926,12 @@
#define mpn_hgcd_matrix_init __MPN (hgcd_matrix_init)
__GMP_DECLSPEC void mpn_hgcd_matrix_init __GMP_PROTO ((struct hgcd_matrix *, mp_size_t, mp_ptr));
+#define mpn_hgcd_matrix_update_q __MPN (hgcd_matrix_update_q)
+__GMP_DECLSPEC void mpn_hgcd_matrix_update_q __GMP_PROTO ((struct hgcd_matrix *, mp_srcptr, mp_size_t, unsigned, mp_ptr));
+
+#define mpn_hgcd_matrix_mul_1 __MPN (hgcd_matrix_mul_1)
+__GMP_DECLSPEC void mpn_hgcd_matrix_mul_1 __GMP_PROTO ((struct hgcd_matrix *, const struct hgcd_matrix1 *, mp_ptr));
+
#define mpn_hgcd_matrix_mul __MPN (hgcd_matrix_mul)
__GMP_DECLSPEC void mpn_hgcd_matrix_mul __GMP_PROTO ((struct hgcd_matrix *, const struct hgcd_matrix *, mp_ptr));
diff -r f901b9fcfdd3 -r bd597269199a mpn/generic/hgcd.c
--- a/mpn/generic/hgcd.c Thu May 19 20:55:23 2011 +0200
+++ b/mpn/generic/hgcd.c Thu May 19 21:09:39 2011 +0200
@@ -4,7 +4,7 @@
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 2003, 2004, 2005, 2008 Free Software Foundation, Inc.
+Copyright 2003, 2004, 2005, 2008, 2011 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
@@ -25,129 +25,25 @@
#include "gmp-impl.h"
#include "longlong.h"
-/* For input of size n, matrix elements are of size at most ceil(n/2)
- - 1, but we need two limbs extra. */
-void
-mpn_hgcd_matrix_init (struct hgcd_matrix *M, mp_size_t n, mp_ptr p)
+
+static void
+hgcd_hook (void *p, mp_srcptr gp, mp_size_t gn,
+ mp_srcptr qp, mp_size_t qn, int d)
{
- mp_size_t s = (n+1)/2 + 1;
- M->alloc = s;
- M->n = 1;
- MPN_ZERO (p, 4 * s);
- M->p[0][0] = p;
- M->p[0][1] = p + s;
- M->p[1][0] = p + 2 * s;
- M->p[1][1] = p + 3 * s;
+ ASSERT (gp == NULL);
+ ASSERT (d >= 0);
+ ASSERT (d <= 1);
- M->p[0][0][0] = M->p[1][1][0] = 1;
-}
-
-/* Updated column COL, adding in column (1-COL). */
-static void
-hgcd_matrix_update_1 (struct hgcd_matrix *M, unsigned col)
-{
- mp_limb_t c0, c1;
- ASSERT (col < 2);
-
- c0 = mpn_add_n (M->p[0][col], M->p[0][0], M->p[0][1], M->n);
- c1 = mpn_add_n (M->p[1][col], M->p[1][0], M->p[1][1], M->n);
-
- M->p[0][col][M->n] = c0;
- M->p[1][col][M->n] = c1;
-
- M->n += (c0 | c1) != 0;
- ASSERT (M->n < M->alloc);
-}
-
-/* Updated column COL, adding in column Q * (1-COL). Temporary
- * storage: qn + n <= M->alloc, where n is the size of the largest
- * element in column 1 - COL. */
-static void
-hgcd_matrix_update_q (struct hgcd_matrix *M, mp_srcptr qp, mp_size_t qn,
- unsigned col, mp_ptr tp)
-{
- ASSERT (col < 2);
-
- if (qn == 1)
+ MPN_NORMALIZE (qp, qn);
+ if (qn > 0)
{
- mp_limb_t q = qp[0];
- mp_limb_t c0, c1;
-
- c0 = mpn_addmul_1 (M->p[0][col], M->p[0][1-col], M->n, q);
- c1 = mpn_addmul_1 (M->p[1][col], M->p[1][1-col], M->n, q);
-
- M->p[0][col][M->n] = c0;
- M->p[1][col][M->n] = c1;
-
- M->n += (c0 | c1) != 0;
+ struct hgcd_matrix *M = (struct hgcd_matrix *) p;
+ /* NOTES: This is a bit ugly. A tp area is passed to
+ gcd_subdiv_step, which stores q at the start of that area. We
+ now use the rest. */
+ mp_ptr tp = (mp_ptr) qp + qn;
+ mpn_hgcd_matrix_update_q (M, qp, qn, d, tp);
}
- else
- {
- unsigned row;
-
- /* Carries for the unlikely case that we get both high words
- from the multiplication and carries from the addition. */
- mp_limb_t c[2];
- mp_size_t n;
-
- /* The matrix will not necessarily grow in size by qn, so we
- need normalization in order not to overflow M. */
-
- for (n = M->n; n + qn > M->n; n--)
- {
- ASSERT (n > 0);
- if (M->p[0][1-col][n-1] > 0 || M->p[1][1-col][n-1] > 0)
- break;
- }
-
- ASSERT (qn + n <= M->alloc);
-
- for (row = 0; row < 2; row++)
- {
- if (qn <= n)
- mpn_mul (tp, M->p[row][1-col], n, qp, qn);
- else
- mpn_mul (tp, qp, qn, M->p[row][1-col], n);
-
- ASSERT (n + qn >= M->n);
- c[row] = mpn_add (M->p[row][col], tp, n + qn, M->p[row][col], M->n);
- }
- if (c[0] | c[1])
- {
- M->n = n + qn + 1;
- M->p[0][col][n-1] = c[0];
- M->p[1][col][n-1] = c[1];
- }
- else
- {
- n += qn;
- n -= (M->p[0][col][n-1] | M->p[1][col][n-1]) == 0;
- if (n > M->n)
- M->n = n;
- }
- }
-
- ASSERT (M->n < M->alloc);
-}
-
-/* Multiply M by M1 from the right. Since the M1 elements fit in
- GMP_NUMB_BITS - 1 bits, M grows by at most one limb. Needs
- temporary space M->n */
-static void
-hgcd_matrix_mul_1 (struct hgcd_matrix *M, const struct hgcd_matrix1 *M1,
- mp_ptr tp)
-{
- mp_size_t n0, n1;
-
- /* Could avoid copy by some swapping of pointers. */
- MPN_COPY (tp, M->p[0][0], M->n);
- n0 = mpn_hgcd_mul_matrix1_vector (M1, M->p[0][0], tp, M->p[0][1], M->n);
- MPN_COPY (tp, M->p[1][0], M->n);
- n1 = mpn_hgcd_mul_matrix1_vector (M1, M->p[1][0], tp, M->p[1][1], M->n);
-
- /* Depends on zero initialization */
- M->n = MAX(n0, n1);
- ASSERT (M->n < M->alloc);
}
/* Perform a few steps, using some of mpn_hgcd2, subtraction and
@@ -159,7 +55,7 @@
limbs, and hgcd_mul_matrix1_inverse_vector, n limbs. If hgcd2
fails, needs space for the quotient, qn <= n - s + 1 limbs, for and
hgcd_matrix_update_q, qn + (size of the appropriate column of M) <=
- resulting size of $.
+ resulting size of M.
If N is the input size to the calling hgcd, then s = floor(N/2) +
1, M->n < N, qn + matrix size <= n - s + 1 + n - s = 2 (n - s) + 1
@@ -209,7 +105,7 @@
if (mpn_hgcd2 (ah, al, bh, bl, &M1))
{
/* Multiply M <- M * M1 */
- hgcd_matrix_mul_1 (M, &M1, tp);
+ mpn_hgcd_matrix_mul_1 (M, &M1, tp);
/* Can't swap inputs, so we need to copy. */
MPN_COPY (tp, ap, n);
@@ -218,262 +114,8 @@
}
subtract:
- /* There are two ways in which mpn_hgcd2 can fail. Either one of ah and
- bh was too small, or ah, bh were (almost) equal. Perform one
- subtraction step (for possible cancellation of high limbs),
- followed by one division. */
- /* Since we must ensure that #(a-b) > s, we handle cancellation of
- high limbs explicitly up front. (FIXME: Or is it better to just
- subtract, normalize, and use an addition to undo if it turns out
- the the difference is too small?) */
- for (an = n; an > s; an--)
- if (ap[an-1] != bp[an-1])
- break;
-
- if (an == s)
- return 0;
-
- /* Maintain a > b. When needed, swap a and b, and let col keep track
- of how to update M. */
- if (ap[an-1] > bp[an-1])
- {
- /* a is largest. In the subtraction step, we need to update
- column 1 of M */
- col = 1;
- }
- else
- {
- MP_PTR_SWAP (ap, bp);
- col = 0;
- }
-
- bn = n;
- MPN_NORMALIZE (bp, bn);
- if (bn <= s)
- return 0;
-
- /* We have #a, #b > s. When is it possible that #(a-b) < s? For
- cancellation to happen, the numbers must be of the form
-
- a = x + 1, 0, ..., 0, al
- b = x , GMP_NUMB_MAX, ..., GMP_NUMB_MAX, bl
-
- where al, bl denotes the least significant k limbs. If al < bl,
- then #(a-b) < k, and if also high(al) != 0, high(bl) != GMP_NUMB_MAX,
- then #(a-b) = k. If al >= bl, then #(a-b) = k + 1. */
-
- if (ap[an-1] == bp[an-1] + 1)
- {
- mp_size_t k;
- int c;
- for (k = an-1; k > s; k--)
- if (ap[k-1] != 0 || bp[k-1] != GMP_NUMB_MAX)
- break;
-
- MPN_CMP (c, ap, bp, k);
- if (c < 0)
- {
- mp_limb_t cy;
-
- /* The limbs from k and up are cancelled. */
- if (k == s)
- return 0;
- cy = mpn_sub_n (ap, ap, bp, k);
- ASSERT (cy == 1);
- an = k;
- }
- else
- {
- ASSERT_NOCARRY (mpn_sub_n (ap, ap, bp, k));
- ap[k] = 1;
- an = k + 1;
More information about the gmp-commit
mailing list