[GiNaC-list] lsolve problem

[Richard B. Kreckel]

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