Bill Allombert on Sat, 23 May 2009 17:53:06 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: Another problem with matrix inversion |
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.