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/>