gcdext, computing the biggest cofactor only [was: mpz reuse test takes too much time]
Marco Bodrato
bodrato at mail.dm.unipi.it
Sun Dec 18 17:00:42 UTC 2016
Ciao,
Il Sab, 3 Dicembre 2016 12:13 pm, Niels Möller ha scritto:
> In the case
>
> mpn_gcdext(G, S, NULL, A, B)
>
> with size(A) < size(B), IIRC we currently handle that by swapping
> arguments, calling mpn_gcdext to compute T (even though the caller
You remember correctly.
> That doesn't seem ideal. The alternative would be to do an up-front
> division, and call mpn_gcdext with inputs A, B mod A, and some different
Another alternative is to zero-pad the smallest operand. The _documented_
interface of mpn_gcdext allows this trick.
> One also needs to keep in mind that if |A| < |B|, then typically |T| <
> |S| too. Which could influence algorithm choice if the size difference
I'd rather delegate the algorithm choice to the mpn layer :-)
A possible patch follows, it swaps back A and B if only the biggest
cofactor is needed and zero pad the shortest operand. It works smoothly,
but I did not test its impact on speed...
diff -r 1cdb2f2d864b mpz/gcdext.c
--- a/mpz/gcdext.c Sun Dec 18 09:19:44 2016 +0100
+++ b/mpz/gcdext.c Sun Dec 18 18:58:15 2016 +0100
@@ -82,6 +82,25 @@
TMP_MARK;
+ if (s == NULL)
+ {
+ MPZ_SRCPTR_SWAP (a, b);
+ MP_SIZE_T_SWAP (asize, bsize);
+ MPZ_PTR_SWAP (s, t);
+ ASSERT (asize < bsize);
+
+ TMP_ALLOC_LIMBS_3 (tmp_gp, bsize, tmp_sp, bsize + 1, tmp_bp, bsize
* 2);
+ tmp_ap = tmp_bp + bsize;
+ MPN_COPY (tmp_ap, PTR (a), asize);
+ MPN_FILL (tmp_ap + asize, bsize - asize, 0);
+ MPN_COPY (tmp_bp, PTR (b), bsize);
+
+ gsize = mpn_gcdext (tmp_gp, tmp_sp, &tmp_ssize, tmp_ap, bsize,
tmp_bp, bsize);
+ ssize = ABS (tmp_ssize);
+ tmp_ssize = SIZ (a) >= 0 ? tmp_ssize : -tmp_ssize;
+ }
+ else
+ {
TMP_ALLOC_LIMBS_2 (tmp_gp, bsize, tmp_sp, asize + bsize + bsize + 1);
tmp_bp = tmp_sp + bsize + 1;
tmp_ap = tmp_bp + bsize;
@@ -112,8 +131,8 @@
mpz_sub (x, gtmp, x);
mpz_divexact (t, x, b);
}
+ }
- if (s != NULL)
{
mp_ptr sp;
Best regards,
m
--
http://bodrato.it/papers/
More information about the gmp-devel
mailing list