Karim BELABAS on Thu, 23 Jul 1998 15:18:21 +0200 (MET DST)


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

Re: bug report/techn. question: incorrect type in ...


[Maximilian Pitschi:]
> I'm trying to do some matix computations, e.g., matdet(...), with
> identical matrices M1 and M2, which depend on a parameter p.  See log
> below. To me, the matrices differ only in the way they were defined.
> 
> Why does matdet(M2) fail with "incorrect type in gexpo" while matdet(M1)
> works, although type(M2[...]) says the same type(M1[...])? --- Is this a
> bug or ...?
> [...]
> Log:
> parisize = 4000000, primelimit = 500000, buffersize = 30000
> ? M1[1,1] += 4999500/1000000/p;
> [...]
> ? M2[1,1] += 4.999500/p;

In all the standard linear algebra routines, `inexact real' matrices
(involving a real number somewhere in the structure) are processed using a
maximal pivot strategy, to ensure maximal numerical stability. Henceforth,
the routine more or less assume they can compare coefficients with each
other as if they were all scalars, and this fails badly when polynomials
are involved...

This is one of the cases where PARI's flexibility (user doesn't have to
specify precisely his objects, and routines rely on quick heuristics to
determine what exact routine to call) is actually harmful. It's easy to
change the heuristic, but probably impossible to make it foolproof.

[Change gauss_get_prec() in src/basemath/alglin1.c, inhibiting the maximal
pivot thing if a non-scalar occurs in the structure. I can do that.]

Karim.
--
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
--
PARI/GP Home Page: http://pari.home.ml.org