Re: Another problem with matrix inversion

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

On Wed, May 13, 2009 at 08:50:53AM +0200, Lorenz Minder wrote:
> Hi,
> 
> > PARI only know perform polynomial and linear algebra over a field, so
> > it assumes that moduli are prime.
> 
> Yes, of course.  My point is that it would be better if it worked also
> if the base-ring is not a field. Since even in that case, this is a
> meaningful question, a useful answer would be preferable, in my opinion.

There is the major issue that module over a ring do not generally have a
basis.

> I did a sketchy implementation of what I outlined in the other mail, and
> a patch is attached.  This patch only works for small integers right now,
> as I have only modified Flm_gauss(), but not FpM_gauss().  It's not
> yet production code, only a prototype for discussion.

Please do not change Flm and FpM function: they are meant to be restricted
to prime modulus. You will need to create new functions ZnM_xxx
to handle matrices over Z/nZ.

> It seems to work fine for what I need it, and so I'd like to
> ask if a more complete and tested patch among the same lines would be
> considered for inclusion in PARI.

I trust that you can write a PARI function ZnM_gauss that would work correctly.

The issue is how do you plan to integrate that with GP and keep some 
consistency ? i.e. it is always annoying when a feature (support for rings)
is only supported for a particular ring (Z/nZ) and for a particular function.

And furthermore what application do you have in mind ?

> Baring any mistakes I may have made, the code should work for any
> modulus, and claim that a matrix is singular only if it is indeed
> singular.
> 
> As is, the code in the patch reduces to the original algorithm if p is
> a prime. For composite moduli, O(log(m)) copies of the matrix are kept
> on the stack in the worst case.
> 
> Incidentially, this patch seems to fix a nasty bug with Flm_gauss() when
> the right hand side is a t_COL instead of a t_MAT; but I could not find
> any caller using it that way in any case. (And I _might_ be
> misunderstanding the old code, but I think that is unlikely.)

A function named FpM_gauss is not supposed to handle arguments that are not
t_MAT.

> I have now also noticed that there is matsolvemod(), which seems to work
> fine, too. So alternatively it would be possible to implement an
> algorithm on top of modsolvemod(); would that be better or does it have
> other drawbacks such as slower performance?  (Is it possible to give a
> bunch of right hand sides to matsolvemod() at once?)

I would not know.

Cheers,
Bill.