[patch] Fix a failure in matrix::division_free_elimination()
Hi, folks. I've run into another case of matrix::solve() failing. This time it seems like the division-free elimination doesn't handle un-normal zeros properly: currently it only does expand() (instead of normal()) on the results of it's calculations, but then uses matrix::pivot() to find the pivot, and that routine also uses expand().is_zero() for zero testing. This breaks when fractions are present in the initial matrix. Not every time, but at least in one example. Note that this problem also affects matrix::solve() with 'algo' set to solve_algo::automatic, which is how I encountered it. I'm attaching the example added to the tests, and also the fix that makes matrix::division_free_elimination() use normal() in stead of expand(). An alternative would be to use normal() in the matrix::pivot() routine, but I decided against it because: 1) that way the results of the normal() call would be discarded; 2) that could slow down e.g. matrix::gauss_elimination() which also uses matrix::pivot(), but normalizes expressions itself. Of course, using normal() makes division_free_elimination only truly "division-free" if the input matrix didn't have fractions. In other cases, GCDs will be computed.
It appears that I've messed up the formatting and left the useless "matrix cmp" variable in the patch. Sorry about that. Here are the fixed versions.
On 09.10.2018 14:59, Vitaly Magerya wrote:
It appears that I've messed up the formatting and left the useless "matrix cmp" variable in the patch. Sorry about that. Here are the fixed versions.
Thank you, Vitaly! Maybe it is just impossible to avoid unnormal numbers without calling .normal() in these matrix algorithms. Your recent patches make me wonder if the heuristics in matrix::determinant() surrounding variable normal_flag are really a good idea at all. In any case, Laplace elimination does return an unnormal 1 on matrix{{a/(a-b), 1}, {b/(a-b), 1}}, and this is only remedied by the very strange line 795 in matrix.cpp. What do you think? All my best, -richy.
participants (2)
-
Richard B. Kreckel
-
Vitaly Magerya