Lorenz Minder on Thu, 21 May 2009 04:58:52 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: Another problem with matrix inversion


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.

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

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

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
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'm quite relieved, this makes much more sense to me anyways.  And if
p had to be prime, then other things such as is_FpM() would be broken.)

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

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

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

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

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.

So in summary, is it OK if I move back to simply changing Flm_gauss()
and FpM_gauss() instead?  And can we finally agree that this is a bug
fix, not some new functionality thing?

Best regards,
--Lorenz
-- 
Neu: GMX FreeDSL Komplettanschluss mit DSL 6.000 Flatrate + Telefonanschluss für nur 17,95 Euro/mtl.!* http://dslspecial.gmx.de/freedsl-aktionspreis/?ac=OM.AD.PD003K11308T4569a