bug#41223: 28.0.50; Calc gives wrong determinant for some matrices with symbolic variables

bug#41223: 28.0.50; Calc gives wrong determinant for some matrices with symbolic variables

[Mauro Aranda]

[-- Attachment #1: Type: text/plain, Size: 4152 bytes --]

I found this problem with matrices larger than 3x3.

Consider a matrix like this:
1 0 0 0
0 1 0 0
0 0 0 1
0 0 a 0

Its determinant is -a, but Calc says its 0.

Steps to reproduce:

1) emacs -Q
2) Start Calc: C-x * c
3) Enter algebraic mode by typing '
4) Enter the above matrix, like this:
[[1 0 0 0][0 1 0 0][0 0 0 1][0 0 a 0]]
5) Calculate the determinant: V D

Calc says that the determinant is 0, which is wrong.

If, however, a is replaced with a value, say 5, Calc correctly solves it
and reports the determinant as -5.

With a 3x3 matrix, the result is correct.  That is:
V D on [[1 0 0][0 0 1][0 a 0]] gives -a as the result.

I think the algorithm for choosing the pivots should be improved to be
able to permutate the rows and choose a symbolic variable as the pivot.

In GNU Emacs 28.0.50 (build 7, x86_64-pc-linux-gnu, GTK+ Version 3.22.30,
cairo version 1.15.10)
 of 2020-05-12 built on tbb-desktop
Repository revision: 703115829b35de6a90d7bafb7931f905e79d0d35
Repository branch: master
Windowing system distributor 'The X.Org Foundation', version 11.0.12005000
System Description: Ubuntu 18.04.4 LTS

Recent messages:
For information about GNU Emacs and the GNU system, type C-h C-a.
C-x *  (Type ? for a list of Calc options)
Welcome to the GNU Emacs Calculator!  Press ‘?’ or ‘h’ for help, ‘q’ to quit
Working... LUD step = 3/4
C-x *  (Type ? for a list of Calc options)

Configured features:
XPM JPEG TIFF GIF PNG CAIRO SOUND DBUS GSETTINGS GLIB NOTIFY INOTIFY
GNUTLS FREETYPE HARFBUZZ ZLIB TOOLKIT_SCROLL_BARS GTK3 X11 XDBE XIM
MODULES THREADS PDUMPER GMP

Important settings:
  value of $LC_MONETARY: es_AR.UTF-8
  value of $LC_NUMERIC: es_AR.UTF-8
  value of $LC_TIME: es_AR.UTF-8
  value of $LANG: en_US.UTF-8
  locale-coding-system: utf-8-unix

Major mode: Lisp Interaction

Minor modes in effect:
  tooltip-mode: t
  global-eldoc-mode: t
  eldoc-mode: t
  electric-indent-mode: t
  mouse-wheel-mode: t
  tool-bar-mode: t
  menu-bar-mode: t
  file-name-shadow-mode: t
  global-font-lock-mode: t
  font-lock-mode: t
  blink-cursor-mode: t
  auto-composition-mode: t
  auto-encryption-mode: t
  auto-compression-mode: t
  line-number-mode: t
  transient-mark-mode: t

Load-path shadows:
None found.

Features:
(shadow sort mail-extr emacsbug message rmc puny dired dired-loaddefs
format-spec rfc822 mml mml-sec password-cache epa derived epg epg-config
gnus-util rmail rmail-loaddefs text-property-search time-date subr-x seq
byte-opt gv bytecomp byte-compile cconv mm-decode mm-bodies mm-encode
mail-parse rfc2231 mailabbrev gmm-utils mailheader sendmail rfc2047
rfc2045 ietf-drums mm-util mail-prsvr mail-utils calc-arith calc-mtx
calccomp calc-alg calc-vec calc-aent calc-menu easymenu calc-misc
calc-prog calc-ext cl-loaddefs cl-lib calc calc-loaddefs calc-macs
tooltip eldoc electric uniquify ediff-hook vc-hooks lisp-float-type
mwheel term/x-win x-win term/common-win x-dnd tool-bar dnd fontset image
regexp-opt fringe tabulated-list replace newcomment text-mode elisp-mode
lisp-mode prog-mode register page tab-bar menu-bar rfn-eshadow isearch
timer select scroll-bar mouse jit-lock font-lock syntax facemenu
font-core term/tty-colors frame minibuffer cl-generic cham georgian
utf-8-lang misc-lang vietnamese tibetan thai tai-viet lao korean
japanese eucjp-ms cp51932 hebrew greek romanian slovak czech european
ethiopic indian cyrillic chinese composite charscript charprop
case-table epa-hook jka-cmpr-hook help simple abbrev obarray
cl-preloaded nadvice loaddefs button faces cus-face macroexp files
text-properties overlay sha1 md5 base64 format env code-pages mule
custom widget hashtable-print-readable backquote threads dbusbind
inotify dynamic-setting system-font-setting font-render-setting cairo
move-toolbar gtk x-toolkit x multi-tty make-network-process emacs)

Memory information:
((conses 16 76586 9357)
 (symbols 48 9259 1)
 (strings 32 22656 1028)
 (string-bytes 1 685264)
 (vectors 16 11989)
 (vector-slots 8 156067 9126)
 (floats 8 28 37)
 (intervals 56 207 0)
 (buffers 992 12))

[-- Attachment #2: Type: text/html, Size: 4469 bytes --]

bug#41223: 28.0.50; Calc gives wrong determinant for some matrices with symbolic variables

[Mattias Engdegård]

Thank you for reporting this. If Calc cannot compute something, at least it shouldn't lie.

In this case I'm tempted to apply the expedient below. Would that suffice for you?

--- a/lisp/calc/calc-mtx.el
+++ b/lisp/calc/calc-mtx.el
@@ -275,7 +275,7 @@ math-do-matrix-lud
                k (1+ k)))
        (setcar (nthcdr j (nth i lu)) sum)
        (let ((dum (math-lud-pivot-check sum)))
-         (if (Math-lessp big dum)
+         (if (or (math-zerop big) (Math-lessp big dum))
              (setq big dum
                    imax i)))
        (setq i (1+ i)))

bug#41223: 28.0.50; Calc gives wrong determinant for some matrices with symbolic variables

[Mauro Aranda]

[-- Attachment #1: Type: text/plain, Size: 800 bytes --]

Mattias Engdegård <mattiase@acm.org> writes:

Hello Mattias.

> Thank you for reporting this. If Calc cannot compute something, at least
it shouldn't lie.

Agreed.  And thank you for taking care of this so quickly.

> In this case I'm tempted to apply the expedient below. Would that suffice
for you?
>
> --- a/lisp/calc/calc-mtx.el
> +++ b/lisp/calc/calc-mtx.el
> @@ -275,7 +275,7 @@ math-do-matrix-lud
>                 k (1+ k)))
>         (setcar (nthcdr j (nth i lu)) sum)
>         (let ((dum (math-lud-pivot-check sum)))
> -         (if (Math-lessp big dum)
> +         (if (or (math-zerop big) (Math-lessp big dum))
>               (setq big dum
>                     imax i)))
>         (setq i (1+ i)))

It is good enough for me.  Thanks.

Best regards,
Mauro.

[-- Attachment #2: Type: text/html, Size: 1010 bytes --]

bug#41223: 28.0.50; Calc gives wrong determinant for some matrices with symbolic variables

[Mattias Engdegård]

14 maj 2020 kl. 00.15 skrev Mauro Aranda <maurooaranda@gmail.com>:

> It is good enough for me.  Thanks.

The pleasure was all mine! It has now been pushed to master.

No doubt there are better ways of computing the determinant for non-numeric matrices. In particular, division shouldn't be necessary at all, and the current algorithm results in some untidy expressions. For example, the determinant of

[a 0 0 1]
[0 b 0 0]
[0 0 c 0]
[1 0 0 d]

is given as abc*(d-1/a), falsely giving the impression of being undefined for a=0, rather than bc(ad-1).

bug#41223: 28.0.50; Calc gives wrong determinant for some matrices with symbolic variables

[Mattias Engdegård]

By the way, the change also fixes inversion of matrices such as the one in your example. Previously, Calc would (erroneously) complain that the matrix was singular.

bug#41223: 28.0.50; Calc gives wrong determinant for some matrices with symbolic variables

[Mauro Aranda]

[-- Attachment #1: Type: text/plain, Size: 312 bytes --]

Mattias Engdegård <mattiase@acm.org> writes:

> By the way, the change also fixes inversion of matrices such as the
> one in your example. Previously, Calc would (erroneously) complain
> that the matrix was singular.

Ah, yes.  I tested that, but forgot to include it in the bug report.

Thanks again!

[-- Attachment #2: Type: text/html, Size: 436 bytes --]

bug#41223: 28.0.50; Calc gives wrong determinant for some matrices with symbolic variables

[Mauro Aranda]

[-- Attachment #1: Type: text/plain, Size: 600 bytes --]

Mattias Engdegård <mattiase@acm.org> writes:

> No doubt there are better ways of computing the determinant for
> non-numeric matrices. In particular, division shouldn't be necessary
> at all, and the current algorithm results in some untidy
> expressions. For example, the determinant of
>
> [a 0 0 1]
> [0 b 0 0]
> [0 0 c 0]
> [1 0 0 d]
>
> is given as abc*(d-1/a), falsely giving the impression of being
> undefined for a=0, rather than bc(ad-1).

I see.  Yes, using a fraction-free algorithm would be good.  But for the
moment, getting bc(ad-1) as the result is one `j M' away.

[-- Attachment #2: Type: text/html, Size: 743 bytes --]