Thank you!
On Saturday, November 5, 2022 at 12:15:26 PM UTC+1 dim...@gmail.com wrote: > I've opened https://trac.sagemath.org/ticket/34724 to deal with this > > On Sat, Nov 5, 2022 at 11:08 AM Dima Pasechnik <dim...@gmail.com> wrote: > > > > Well, applying a naive exact linear algebra routine to inexact data, > > and that's what Sage is doing here, is prone to errors. > > sage: A=m.augment(identity_matrix(RR,2)) > > sage: A > > [-6.12323399573677e-17 -1.72508242466029 1.00000000000000 > > 0.000000000000000] > > [ 0.579682446302195 6.12323399573677e-17 0.000000000000000 > > 1.00000000000000] > > sage: A.echelonize(algorithm="classical");A # OOOPS! - that's the > default here > > [ 1.00000000000000 0.000000000000000 4.00000000000000 > > 1.72508242466029] > > [ 0.000000000000000 1.00000000000000 -0.579682446302195 > > -6.12323399573676e-17] > > sage: A=m.augment(identity_matrix(RR,2)) > > sage: A.echelonize(algorithm='scaled_partial_pivoting');A # that's how > > it should be > > [ 1.00000000000000 0.000000000000000 6.12323399573677e-17 > > 1.72508242466029] > > [ 0.000000000000000 1.00000000000000 -0.579682446302195 > > -6.12323399573677e-17] > > > > On Sat, Nov 5, 2022 at 10:21 AM Emmanuel Charpentier > > <emanuel.c...@gmail.com> wrote: > > > > > > something is definitely unhinged here : On 9.8.beta3 running on Debian > testing on core i7 + 16 GB RAM, after running : > > > > > > a = RR(-4967757600021511 / 2**106) > > > b = RR(-7769080564883485 / 2**52) > > > c = RR( 5221315298319565 / 2**53) > > > m = matrix([[a, b], [c, -a]]) > > > M = matrix([[var("p%d%d"%(u, v), latex_name="p_{%s,%d}"%(u, v)) > > > for v in range(2)] > > > for u in range(2)]) > > > S = dict(zip(M.list(), [a, b, c, -a])) > > > MN = M.apply_map(lambda u:u.subs(S)) > > > > > > one gets : > > > > > > sage: m.parent() > > > Full MatrixSpace of 2 by 2 dense matrices over Real Field with 53 bits > of precision > > > sage: m*~m > > > [ 1.00000000000000 -1.23259516440783e-32] > > > [ 2.31872978520878 1.00000000000000] > > > sage: (~m)*m > > > [ 1.00000000000000 -6.90032969864117] > > > [6.16297582203915e-33 1.00000000000000] > > > sage: MN.parent() > > > Full MatrixSpace of 2 by 2 dense matrices over Symbolic Ring > > > sage: MN*~MN > > > [ 1.00000000000000 -1.23259516440783e-32] > > > [ 2.31872978520878 1.00000000000000] > > > sage: (~MN)*MN > > > [ 1.00000000000000 -6.90032969864117] > > > [6.16297582203915e-33 1.00000000000000] > > > > > > all being wrong, wrong, wrong… > > > > > > However : > > > > > > sage: (M*~M).apply_map(lambda u:u.subs(S)) > > > [ 1.00000000000000 0] > > > [-3.54953126192945e-17 1.00000000000000] > > > sage: ((~M)*M).apply_map(lambda u:u.subs(S)) > > > [ 1.00000000000000 1.05630833481279e-16] > > > [ 0 1.00000000000000] > > > > > > both being acceptable. > > > > > > One also notes that the form of : > > > > > > sage: ~M > > > [1/p00 - p01*p10/(p00^2*(p01*p10/p00 - p11)) p01/(p00*(p01*p10/p00 - > p11))] > > > [ p10/(p00*(p01*p10/p00 - p11)) -1/(p01*p10/p00 - p11)] > > > sage: (~M).apply_map(simplify) > > > [1/p00 - p01*p10/(p00^2*(p01*p10/p00 - p11)) p01/(p00*(p01*p10/p00 - > p11))] > > > [ p10/(p00*(p01*p10/p00 - p11)) -1/(p01*p10/p00 - p11)] > > > > > > is somewhat unexpected ; one expects : > > > > > > sage: (~M).apply_map(lambda u:u.simplify_full()) > > > [-p11/(p01*p10 - p00*p11) p01/(p01*p10 - p00*p11)] > > > [ p10/(p01*p10 - p00*p11) -p00/(p01*p10 - p00*p11)] > > > > > > which is also the form returned by maxima : > > > > > > sage: maxima_calculus.invert(M).sage() > > > [-p11/(p01*p10 - p00*p11) p01/(p01*p10 - p00*p11)] > > > [ p10/(p01*p10 - p00*p11) -p00/(p01*p10 - p00*p11)] > > > > > > giac : > > > > > > sage: giac.inverse(giac(M)).sage() > > > [[-p11/(p01*p10 - p00*p11), p01/(p01*p10 - p00*p11)], > > > [p10/(p01*p10 - p00*p11), -p00/(p01*p10 - p00*p11)]] > > > > > > fricas : > > > > > > sage: fricas.inverse(M._fricas_()).sage() > > > [-p11/(p01*p10 - p00*p11) p01/(p01*p10 - p00*p11)] > > > [ p10/(p01*p10 - p00*p11) -p00/(p01*p10 - p00*p11)] > > > > > > mathematica : > > > > > > sage: mathematica.Inverse(M).sage() > > > [[-p11/(p01*p10 - p00*p11), p01/(p01*p10 - p00*p11)], > > > [p10/(p01*p10 - p00*p11), -p00/(p01*p10 - p00*p11)]] > > > > > > and (somewhat un-backconvertible) : > > > > > > sage: sympy.sympify(M)^-1._sage_() > > > Matrix([ > > > [ p11/(p00*p11 - p01*p10), -p01/(p00*p11 - p01*p10)], > > > [-p10/(p00*p11 - p01*p10), p00/(p00*p11 - p01*p10)]]) > > > > > > This is, IMNSHO, a critical bug. Could you open a tichet for this, and > mark it as such ? > > > > > > Le samedi 5 novembre 2022 à 07:59:27 UTC+1, Håkan Granath a écrit : > > >> > > >> Hi, there seems to be a problem with inverses of matrices with > elements in RR. It only occurs very sporadically for me, but here is an > example: > > >> > > >> a = RR(-4967757600021511 / 2**106) > > >> b = RR(-7769080564883485 / 2**52) > > >> c = RR( 5221315298319565 / 2**53) > > >> > > >> m = matrix([[a, b], [c, -a]]) > > >> > > >> print(m) > > >> print() > > >> print(~m) > > >> > > >> On my machines it produces the output > > >> > > >> [-6.12323399573677e-17 -1.72508242466029] > > >> [ 0.579682446302195 6.12323399573677e-17] > > >> > > >> [ 4.00000000000000 1.72508242466029] > > >> [ -0.579682446302195 -6.12323399573676e-17] > > >> > > >> Clearly the element 4 is wrong (the correct inverse is -m). Is this a > known bug? > > >> > > >> Some system information: > > >> > > >> SageMath version 9.7, using Python 3.10.5 > > >> OS: Ubuntu 20.04.5 LTS > > >> CPU: Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz > > >> > > >> Best regards, > > >> > > >> Håkan Granath > > > > > > -- > > > You received this message because you are subscribed to the Google > Groups "sage-devel" group. > > > To unsubscribe from this group and stop receiving emails from it, send > an email to sage-devel+...@googlegroups.com. > > > To view this discussion on the web visit > https://groups.google.com/d/msgid/sage-devel/bec2d18e-8034-440c-8076-9a5e9ec93d2cn%40googlegroups.com > . > -- You received this message because you are subscribed to the Google Groups "sage-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/b1098ba2-7883-434f-95ab-56221200ee94n%40googlegroups.com.
