Fix determinant for RR

[ppurka]

Refer to ​this sage-devel post

A = matrix([[              1.0,    -1.50614628068,     2.26847661882,    -3.41665762226,     5.14598617013,     -7.7506079306,     11.6735493077,    -17.5820728722],
[              1.0,   -0.936702701875,    0.877411951699,   -0.821874145813,    0.769851732984,   -0.721122198329,    0.675477111557,   -0.632721235449],
[              1.0,   -0.443181140009,     0.19640952286,  -0.0870449962496,     0.03857670067,  -0.0170964661807,  0.00757683137208, -0.00335790876514],
[              1.0,    0.352786603689,    0.124458387743,   0.0439072519123,   0.0154898902795,  0.00546462578321,  0.00192784677049, 0.000680118514595],
[              1.0,    0.647213396311,    0.418885180364,    0.271108100248,    0.175464794329,     0.11356316547,     0.07349960202,   0.0475699270508],
[              1.0,     1.44318114001,     2.08277180288,     3.00581698486,     4.33793838286,     6.26043086067,     9.03493574645,     13.0390488705],
[              1.0,     1.93670270187,     3.75081735545,     7.26421810653,     14.0686308339,     27.2467553477,     52.7688646993,     102.197602838],
[              1.0,     2.50614628068,     6.28076918019,     15.7405263208,     39.4480614948,     98.8626125954,     247.764168855,     620.933250262]])

B = A.change_ring(RDF)
print "det(A) = {}, det(B) = {}".format(A.determinant(), B.determinant())

det(A) = -4.19430400000000e6, det(B) = 16801.7979988

According to Peter Bruin, a possible fix is to use pari:

Well, it should also be fixed for a RealField? of higher precision. An easy solution for that is to use PARI, which uses a numerically more stable algorithm (Gaussian elimination, choosing pivots of maximal absolute value; I don't know about proven error bounds). Example:

sage: A._pari_().matdet()
16801.7979988279  # same as when doing the computation over QQ

Sage's determinant() already uses PARI over Z/nZ for n less than the machine word size; it would be trivial to adapt it to work also over the reals.

According to Dima, the fix should be:

Is the default choice of the algorithm the right one? One can see that

sage: A.determinant(algorithm="hessenberg")
16801.7979988558

is quite good...

Sage computes det() by computing charpoly(0), too... The division-free algorithm is obviously meant for more general setting of rings, not fields. I don't know why it was made default here, perhaps just an oversight.

[jakobkroeker]

now (sage 7.6 beta 4)

A.determinant(algorithm="hessenberg")

also returns the faulty -4.19430400000000e6