Re: Another problem with matrix inversion

[Bill Allombert <Bill.Allombert@math.u-bordeaux1.fr>]

On Thu, May 21, 2009 at 04:43:46AM +0200, Lorenz Minder wrote:
> Hi,
> 
> BA:
> > Sorry for the rushed our reply, but maybe you would prefer it to 
> > a long delayed one.
> 
> I absolutely do, yes.  Thanks for taking the time.

But I am afraid this lead to misunderstanding. Sorry about that.

> I'll fix the trivial stuff you mentionned (and also a bug I spotted),
> but there are more fundamental issues that I think need discussion.
> 
> > You also need to add a test-suite, and the documentation of the new public
> > function.
> 
> I'll look in the test suite, not sure how this works.  (Documentation
> is noted, let's talk about this later.)

Test-suite is merely a matter of adding a file src/test/in/foo or
expanding an existing file that is not used by 'make bench'.

> > I am a bit concerned that your patch is large, fixing the interface issue
> > would make it larger yet,

First, when I said the patch was large, I really meant that it added a lot of
code to PARI, not so much that it made a lot of change.

> I strongly suggest to go back to the way I did it originally, namely
> simply replace FpM_gauss() and Flm_gauss().  Let me quote from the

I have an alternative:
1) We add function FpM_invsafe/Flm_invsafe that would do that same for matrix
inversion that Fp_invsafe do for integers. Maybe even FpM_gausssafe if needed.
This is a minimal change, that use a well-defined interface.
I can do it if you want.

2) We add a function that take a factored modulus and use chinese remainder and
Hensel lifting to use FpM_invsafe to inverse a matrix. Hensel lifting alleviate
the need to change FpM_gauss to look for an invertible pivot. This is similar
to how gen_Shanks_sqrtn works.

3) We add a function that take an unfactored modulus, call the function above
as if the modulus is prime, and if not, refine the factorisation and iterate.

> libpari manual, the part that talks about Fp* and Fq* functions:
> 
>   7.2 Modular arithmetic. 
>   
>   These routines implement univariate polynomial arithmetic and linear
>   algebra over finite fields, in fact over finite rings of the form
>   (Z/pZ)[X]/(T), where p is not necessarily prime and T \in (Z/pZ)[X]
>   is possibly reducible; and finite extensions thereof. All this can be
>   emulated with t_INTMOD and t_POLMOD coefficients and using generic
>   routines, at a considerable loss of efficiency.
> 
> As far as I can see, this states unambiguously that p does not have to
> be prime.  Don't you agree?

I do not. I designed this interface and this says that the Fp operations 
perform as if the ring is a finite field. This means they do not go out 
of their way to handle non-fields.  The whole point of the Fp operations being
to get rid of genericity, which is quite costly.
In practice inside PARI this is used for non-field only for prime power
to handle operation over Z/p^mZ (e.g. Hensel lifting).

> I suggest we start thinking of this patch as what it is, a bug fix.

it is not a "bug fix": the PARI documentation is quite clear on that:
==============================================================================
Vectors, matrices, linear algebra and sets:

   Since PARI does not have a strong typing system, scalars live in unspecified
commutative  base  rings.   It is very difficult to write robust linear algebra
routines  in  such  a general setting.   We thus assume that the base ring is a
domain and work over its field of fractions.  If the base ring is not a domain,
one  gets  an error as soon as a non-zero pivot turns out to be non-invertible.
Some  functions,   e.g. mathnf or mathnfmod,  specifically assume that the base
ring is Z.
==============================================================================
This describes exactly the situation here, so it is not a bug.

Now I agree that, as a matter of "quality of implementation", it would be
nice if PARI was able to inverse all invertible matrices. But this is not a
bug.

> >and does not handle the general matrix inversion
> > problem. 
> 
> That was never the intention, and there are plenty of other patches
> that don't do that either.  Of course, if would be great if PARI would
> more generically invert matrices, but that's an orthogonal issue,
> really.  (Unless the generic algorithm would be as fast as what is
> there now, but frankly I don't buy that without having seen it.  For
> the record, the matadjoint(,1) thing is so slow it's useless for me
> currently.)

OK, I have made some test and I agree that matadjoint is very slow,
even compared to plain RgM_inv without the is_FpM otpimisation.

> Is the patch large? Possibly, but unless/until someone comes up with a
> significantly shorter algorithm that performs equally well, why does
> it matter?

There is always a cost (in storage/ memory/ testing/ maintainance/
consistentcy/ etc.) associate to extra code in PARI. This has to be
balanced with use case.

> > It would be nice if the whole strategy was abstracted to be reused in
> > similar
> > situation.
> 
> FqM_gauss() has the analogous bug of course, and it may be possible to
> reuse the same kind of logic, and possibly share some code if p is a
> prime but the polynomial reducible.  But if p=9, we get things like
> (3x + 1)^3 = 1, and I'm pretty sure that will screw us, although I
> admit not having thought about this carefully.

Agreed. However, factoring polynomials over a field is much easier than
factoring integers, so this trink is not very interesting for polynomials.

> My take on it is that it is not hard to make it more abstract later
> if it proves useful, but I think there is really not much of a point
> doing it just now.
> 
> > Do you think there is really a need for 'ZlM' ?  In that range, factoring
> > the modulus is trivial, and that would reduce the size of the patch 
> > somewhat.
> 
> For Flm_gauss() the problem could indeed be solved differently, but
> I'm not sure it's so much easier, especially since the case of p having
> square factors needs to be considered, and the elimination routines
> fixed accordingly in any case.

This is why I suggest to use Hensel lifting.

Cheers,
Bill.

PS:
> +    {
> +      GEN piv = remii(gcoeff(a,k,i), p);
> +      if (signe(piv)) {
> +        GEN res;
> +     if(!invmod(piv, p, &res)) {
> +       if(mod_is_prime) {
> +         /* Trigger error */
> +         Fp_inv(piv, p);
> +       }

Please use Fp_invsafe() (and implement a function Fl_invsafe for the Fl case).

Could you avoid the code duplication between Z_find_coprime_factors and
find_coprime_factors ? I do not think this is performance critical enough
to warrant it.