Re: Another problem with matrix inversion

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

On Thu, May 14, 2009 at 04:03:14AM +0200, Lorenz Minder wrote:
> Hi,
> 
> BA:
> > On Wed, May 13, 2009 at 08:50:53AM +0200, Lorenz Minder wrote:
> > > > 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.
> 
> Right.  But I can't see what you're getting at here.

It would be a problem for matker and other linear algebra operation.
 
> > > 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.
> 
> Ok. (What would you do with, e.g., FpM_gauss()?  Leave it there
> or have it simply call ZnM_gauss()?)

I would prefer if ZnM_gauss() would call FpM_gauss() that would not be
changed too much (maybe only make it return NULL for division by zero 
error).

> > > 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.
> 
> Right now PARI knows how to add, subtract, multiply and divide in the
> ring A, where A is any of Z, Z[i], F[X], M_n(R) and others. (In
> particular, for any of these, 1/x gives me the A-inverse of x if it
> exists.) It also knows to add, subtract and multiply in M_n(Z/mZ).

Because they are domain and the computation is done in the fraction field.

> So as far as I can see, the inconsistency here is that it does not
> know in general how to divide in M_n(Z/mZ).

or matrices over any other non-domain ring, e.g. Q[X]/(X^2-1),
Z/5Z[I] etc.

In any case if we restrict the problem to matrix inversion, we could implement
a division free algorithm, like 
matinv(A)=matadjoint(A,1)/matdet(A)

At this point the only trouble is to pick the fastest algorithm depending 
on the matrices entries.

> It would of course be integrated into GP the same way as division over
> the other rings is implemented, by overloading the "/" operator, and by
> having matsolve() do something sane.  So it would behave as it does
> with the current patch (except for the fact that the current patch
> does not handle large moduli).
> 
> > And furthermore what application do you have in mind ?
> 
> In my case, I need to solve a set of multivariate equations over GF(q).
> With the right change of variables I can (in some cases) reduce this to
> equations with only products, e.g. x*y*z = const_1, y*z*v = const_2,
> etc.  Then I take the log (to some fixed basis) which gives me a set
> of linear equations in Z/(q - 1)Z.
> 
> Of course, q - 1 is not prime.

So I expect you are using Pohlig-Hellman algorithm along the way,
so you know how to factor q-1, and it might be more efficient to 
use full chinese remainder than this trick.

> > > 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.)

Ah sorry, I looked at FpM_gauss(). Flm_gauss cannot possibly take a t_COL as
second argument, only a t_VECSMALL. There might be a fix in your code for some
bug still but I could not find it.

> > A function named FpM_gauss is not supposed to handle arguments that are
> > not
> > t_MAT.
> 
> Then I think the non-functional support that someone has coded into these
> functions should be removed.

Agreed. We should really provide FpM_FpC_gauss (etc.) functions instead.

Cheers,
Bill.