Bill Allombert on Thu, 14 May 2009 18:48:49 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: Another problem with matrix inversion |
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.