feature: allsqrtmod(a, n)

[hvds]

I found a need for allsqrtmod, and wrote (perl code for) the functions documented below. Would you be interested in incorporating some or all of these in MPU? If you would, I'd be happy to add some tests (and maybe XS variants) and put together a PR.

allsqrtmod is implemented in terms of allsqrtmodfact which in turn uses allsqrtmodpk; the latter two could also be private functions.

=head2 allsqrtmod( $a, $n )

Given integer C<a> and positive integer C<n>, return a list of the square
roots C<x> of C<a mod n> having C<< 0 <= x < n >>. If no square root
exists, an empty list is returned.

=head2 allsqrtmodfact( $a, $n, $f )

Given integer C<a> and positive integer C<n>, and the factorization C<f>
of C<n> as returned by L<factor_exp>, return a list of the square roots C<x>
of C<a mod n> having C<< 0 <= x < n >>. If no square root exists, an empty
list is returned.

=head2 allsqrtmodpk( $a, $p, $k )

Given integer C<a>, prime C<p> and positive integer C<k>, return a list of
the square roots C<x> of C<a mod p^k> having C<< 0 <= x < p^k >>. If no 
square root exists, an empty list is returned.

=cut

[danaj]

I'm definitely interested. I looked at it a little a few months ago, and added "rootmod: a way to get all the roots" to the TODO. It should be easy for primes, reasonable for sqrt, not sure about general roots.

I really need to get things added to the GMP back end and release something before I add too much more. But allsqrtmod (and maybe allrootmod) would be great.

I did think a bit about rootmodp, but it only saves the primality test. The prime power case saves a generic factoring. In both cases one worry is what to do if they lied. I realize that can be undefined, but it's bothersome. Other than that, I don't think there's a big savings in XS, but in pure Perl there is.

[hvds]

Ok, I've created PR #55 with allsqrtmod support. I haven't thought too much about general roots, not sure what's involved with that.

There's still a bit of duplication of effort, eg having to recalculate p^k, but I think it should be reasonably efficient. I'm mainly unsure how applying chinese() recursively on pairs would compare to applying it in a single pass over sets of k elements (for n divisible by k distinct primes) - I'm guessing it's more efficient this way, but in any case doing it as a single pass would need fairly complex logic (or a dependency on Algorithm::Loops).

[danaj]

@hvds I added XS for the general root case. I'll work on the PP code for it.

Currently it shortcuts to the allsqrtmod code when given k=2. I'm thinking I might keep the two separate. I haven't found any differences in the results (of course) and there doesn't seem to be much performance difference.

[hvds]

Keeping the two separate makes sense to me; it's always possible that particular patterns of use might find a more noticeable difference, so this would leave users the opportunity to choose what works best for them.