[GiNaC-list] lsolve problem
Richard B. Kreckel
kreckel at ginac.de
Mon Feb 18 23:52:20 CET 2008
Hi!
Alexei Sheplyakov wrote:
> On Mon, Feb 18, 2008 at 02:32:45PM +0100, Diego Conti wrote:
>
>> But then, wouldn't it be natural to reorganize the code so that
>> power::to_rational throws an exception when it is called as
>> a result of the call to to_rational in fraction_free_elimination?
>
> This approach sounds a bit wrong (to put it very mildly).
One of the reasons why this fix is bad is because it would only be
acceptable for the case solve_algo::automatic. But there are several
other ways of calling fraction_free_elimination more directly.
> First of all, power::to_rational does its work just fine, there's no
> need to change it. It's fraction_free_elimination which is buggy and
> needs to be fixed.
>
> Secondly, this would reject many valid expressions, i.e.
>
> #include <iostream>
> #include <ginac/ginac.h>
> using namespace std;
> using namespace GiNaC;
>
> int main(int argc, char** argv)
> {
> symbol a("A"),b("B"),c("C"),d("D"),x("X");
> lst syms;
> syms = a,b;
> lst eqns;
> eqns = sqrt(1-pow(x,2))*a + x*b == c,
> -x*a + sqrt(1-pow(x,2))*b == d;
> ex sol = lsolve(eqns, syms, solve_algo::bareiss);
> cout << "System: " << eqns << endl << "Solution: " << sol << endl;
> return 0;
> }
Right. This is one way of calling fraction_free_elimination more directly.
>> Then one could have matrix::solve catch the exception and switch
>> automatically to a different algorithm,
>
> The problem is not that simple (see the example above).
Clearly, calling .expand() on the numerator and denominator in the
elimination step is the culprit. It fails to simplify s^2-3 to zero.
On the other hand, I believe one must not try to be more clever in the
elimination step proper: it could lead to a situation where
m[k-1](k-1,k-1) fails to divide m[k+1](r,c) because the latter has been
simplified too much. And that division to work without having to embark
on a GCD computation is the crucial point about Bareiss' fraction free
elimination.
I am having a guts feeling that the problem is only that a zero pivot
element like s^2-3 with s==sqrt(3) is being chosen. This is bad in its
own right because a pivot should be non-zero by definition. If this is
correct, the attached patch would solve all these problems without
incurring much overhead. Maybe it can be written more elegantly?
It passes our testsuite without noticeable benchmark penalties and it
solves Diego's problem. But I haven't done much more thinking and
testing. Diego, can you please try if this works as expected for
whatever you are doing there?
Best wishes
-richy.
--
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: matrix.patch
Type: text/x-patch
Size: 1263 bytes
Desc: not available
Url : http://www.cebix.net/pipermail/ginac-list/attachments/20080218/fab77c3c/attachment.patch
More information about the GiNaC-list
mailing list