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



Greetings!  Just a quick note to see if you saw my recent post on
this.  I think the problem is that the logic

#+cl ((and (numberp x) (numberp y)) (f* x y))

in mul2* of opers.lisp is not right, given that (f* ...) expands into
(the fixnum ...)

I'm guessing someone is pushing cl to *features* in the latest maxima
cvs build, but not in 5.9.0.

Also would appreciate your thoughts on the function proclamation issue
I raised if/when you get a chance.

Take care,

James Amundson <amundson@users.sourceforge.net> writes:

> 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
> 
> 
> 
> 
> _______________________________________________
> Maxima mailing list
> Maxima@www.math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
> 
> 
> 

-- 
Camm Maguire			     			camm@enhanced.com
==========================================================================
"The earth is but one country, and mankind its citizens."  --  Baha'u'llah