2-value divmod

[hvds]

Given eg x = 2y and x == 6 (mod 8), we do not know the value of y (mod 8), but we can say y == 3 (mod 4). Do we already have some variant of divmod, say div2mod, that would return (3, 4) for div2mod(6, 2, 8)? As far as I can see that just involves dividing all arguments by the common gcd before invoking divmod, so if we don't already have something like that it might be a useful convenience function to add.

[danaj]

Interesting. I see in util.h:

extern UV modinverse(UV a, UV p);              /* Returns 1/a mod p */
extern UV divmod(UV a, UV b, UV n);            /* Returns a/b mod n */
extern UV gcddivmod(UV a, UV b, UV n);         /* divmod(a/gcd,b/gcd,n) */

which simplified things in the rootmod code. So clearly there are uses for it.

Here it's just doing a gcd of the two values A,B in A/B mod N, leaving the modulo alone. So a single value is returned. Is there a good use for the two value return, and is there a good reason for or against including the modulo in the gcd?

[hvds]

Is there a good use for the two value return, and is there a good reason for or against including the modulo in the gcd?

I think the two questions are coupled: in the above example given a relationship between two variables (x = 2y) and a known constraint on one (x == 6 mod 8) we want to find the corresponding constraint on the other. I'm not sure when it is useful to find gcddivmod(6, 2, 8) == 3 without also knowing that the relevant gcd was 2.

For generality it's useful to consider separately each of the cases of pairwise gcd for this purported div2mod. Assuming (p, q, r) pairwise coprime, a > 1, I think I would expect:

• f(ap, aq, r) => (p * invmod(q, r), r) == (ap * invmod(aq, r), r);
• f(ap, q, ar) => (ap * invmod(q, ar), ar);
• f(p, aq, ar) => (error value);
• f(ap, aq, ar) => (p * invmod(q, r), r).

The third case reflects the fact that the relationship x = 2y and the constraint x == 5 (mod 8) are incompatible. But in overview, that means nothing interesting happens unless gcd(b, n) > 1, and I think an unoptimized implementation would be something like:

sub div2mod {
  my($a, $b, $n) = @_;
  my $g = gcd($b, $n);
  return error if $a % $g;
  my($a2, $b2, $n2) = map $_ / $g, ($a, $b, $n);
  return (divmod($a2, $b2, $n2), $n2);
}

I use this sort of transformation in my search for APs of integers with various properties (eg OEIS A165497, A165500, A292580): I'm building up a set of modular constraints (such as the x == 6 (mod 8) above); in some cases I then need to do arithmetic on them, eg in the above case to turn it into a constraint on y.