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