gcdext, computing the biggest cofactor only [was: mpz reuse test takes too much time]

[Marco Bodrato]

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/