Re: [Maxima-bugs] [ maxima-Bugs-988701 ] incorrect matrix inversion



On Tue, 2004-07-13 at 12:02, Camm Maguire wrote:
> Can we pin down the maxima and gcl versions a bit more?  If there is a
> bug, I want to fix it soon, as it appears we have to put out a 2.6.3
> errata release.

I've been looking at this a little. I must say that I am perplexed. In
all cases, I am using GCL 2.6.2 from the release tarball, compiled with
just --enable-ansi. I am using a Red Hat Enterprise Linux (-like)
system:

|> gcc --version
gcc (GCC) 3.2.3 20030502 (Red Hat Linux 3.2.3-24)

With Maxima 5.9.0, I do not see the problem:

---------------------------------------------------------------------
Maxima 5.9.0 http://maxima.sourceforge.net
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(C1) mf:matrix([1,0,0],[-.3,1,0],[0,-.3,1]);

                              [   1      0    0 ]
                              [                 ]
(D1)                          [ - 0.3    1    0 ]
                              [                 ]
                              [   0    - 0.3  1 ]
(C2) invert(mf);

                              [ 1.0    0    0  ]
                              [                ]
(D2)                          [ 0.3   1.0   0  ]
                              [                ]
                              [ 0.09  0.3  1.0 ]
(C3) :lisp (trace mul*)

(MUL*)
(C3) invert(mf);

  1> (MUL* 1 1)
  <1 (MUL* 1)
  1> (MUL* -1 0)
  <1 (MUL* 0)
  1> (MUL* 1 0)
  <1 (MUL* 0)
  1> (MUL* -1 -0.29999999999999999)
  <1 (MUL* 0.29999999999999999)
  1> (MUL* 1 1)
  <1 (MUL* 1)
  1> (MUL* -1 0)
  <1 (MUL* 0)
  1> (MUL* 1 0.089999999999999997)
  <1 (MUL* 0.089999999999999997)
  1> (MUL* -1 -0.29999999999999999)
  <1 (MUL* 0.29999999999999999)
  1> (MUL* 1 1)
  <1 (MUL* 1)
                              [ 1.0    0    0  ]
                              [                ]
(D3)                          [ 0.3   1.0   0  ]
                              [                ]
                              [ 0.09  0.3  1.0 ]
(C4) quit();
------------------------------------------------------------------

With Maxima 5.9.0.9beta1, I do see the problem. Note that mul* does not
appear when traced:

------------------------------------------------------------------
Maxima 5.9.0.9beta1 http://maxima.sourceforge.net
Using Lisp Kyoto Common Lisp GCL 2.6.2 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i1) mf:matrix([1,0,0],[-.3,1,0],[0,-.3,1]);

                              [   1      0    0 ]
                              [                 ]
(%o1)                         [ - 0.3    1    0 ]
                              [                 ]
                              [   0    - 0.3  1 ]
(%i2) invert(mf);

                        [      1            0       0 ]
                        [                             ]
(%o2)                   [ - 858993459       1       0 ]
                        [                             ]
                        [ 1889785610   - 858993459  1 ]
(%i3) :lisp (trace mul*)

(MUL*)
(%i3) invert(mf);

                        [      1            0       0 ]
                        [                             ]
(%o3)                   [ - 858993459       1       0 ]
                        [                             ]
                        [ 1889785610   - 858993459  1 ]
------------------------------------------------------------------


invert.lisp contains the code for invert and adjoint, which is called by
invert. The file contains little else. The bug seems to appear only
after $adjoint is compiled:

------------------------------------------------------------------
(%i4) load("invert.lisp");

(%o4)                             invert.lisp
(%i5) invert(mf);

  1> (MUL* 1 1)
  <1 (MUL* 1)
  1> (MUL* -1 0)
  <1 (MUL* 0)
  1> (MUL* 1 0)
  <1 (MUL* 0)
  1> (MUL* -1 -0.29999999999999999)
  <1 (MUL* 0.29999999999999999)
  1> (MUL* 1 1)
  <1 (MUL* 1)
  1> (MUL* -1 0)
  <1 (MUL* 0)
  1> (MUL* 1 0.089999999999999997)
  <1 (MUL* 0.089999999999999997)
  1> (MUL* -1 -0.29999999999999999)
  <1 (MUL* 0.29999999999999999)
  1> (MUL* 1 1)
  <1 (MUL* 1)
                              [ 1.0    0    0  ]
                              [                ]
(%o5)                         [ 0.3   1.0   0  ]
                              [                ]
                              [ 0.09  0.3  1.0 ]
(%i6) :lisp (compile '$adjoint)

Compiling gazonk0.lsp.
End of Pass 1.
End of Pass 2.
OPTIMIZE levels: Safety=2, Space=2, Speed=2
Finished compiling gazonk0.lsp.
#<compiled-function $ADJOINT>
(%i6) invert(mf);

                        [      1            0       0 ]
                        [                             ]
(%o6)                   [ - 858993459       1       0 ]
                        [                             ]
                        [ 1889785610   - 858993459  1 ]

------------------------------------------------------------------

The file invert.lisp has not changed since the 5.9.0 release. $adjoint
calls maset from transl.lisp and mul* from runtim.lisp. Both maset and
mul* are unchanged from 5.9.0. If I replace mul* with mul in
invert.lisp, the problem goes away: (invert.lisp has been modified in
the below example)

------------------------------------------------------------------
(%i7) load("invert.lisp");

(%o7)                             invert.lisp
(%i8) :lisp (untrace mul*)

(MUL*)
(%i8) :lisp (trace mul)

(MUL)
(%i8) invert(mf);

  1> (MUL 1 1)
  <1 (MUL 1)
  1> (MUL -1 0)
  <1 (MUL 0)
  1> (MUL 1 0)
  <1 (MUL 0)
  1> (MUL -1 -0.29999999999999999)
  <1 (MUL 0.29999999999999999)
  1> (MUL 1 1)
  <1 (MUL 1)
  1> (MUL -1 0)
  <1 (MUL 0)
  1> (MUL 1 0.089999999999999997)
  <1 (MUL 0.089999999999999997)
  1> (MUL -1 -0.29999999999999999)
  <1 (MUL 0.29999999999999999)
  1> (MUL 1 1)
  <1 (MUL 1)
                              [ 1.0    0    0  ]
                              [                ]
(%o8)                         [ 0.3   1.0   0  ]
                              [                ]
                              [ 0.09  0.3  1.0 ]
(%i9) :lisp (compile '$adjoint)

Compiling gazonk0.lsp.
End of Pass 1.
End of Pass 2.
OPTIMIZE levels: Safety=2, Space=2, Speed=2
Finished compiling gazonk0.lsp.
#<compiled-function $ADJOINT>
(%i9) invert(mf);

                              [ 1.0    0    0  ]
                              [                ]
(%o9)                         [ 0.3   1.0   0  ]
                              [                ]
                              [ 0.09  0.3  1.0 ]
(%i10)
------------------------------------------------------------------

I see that the trace on mul doesn't work after compiling. Since we now
get the right answer, I'll guess that the earlier failure of trace on
mul* was a red herring.

--Jim