Hi, thanks for your reply.
I think I have found a bug in lsolve. I have tested the following program under ginac 1.3.9 and 1.4.1, and I obtain an incorrect result in both cases: the linear system contains the equation lambda28=0, but I obtain lambda28==lambda28 in the solution.
<snip> An interesting problem. It is due to the division free elimination algorithm replacing sqrt(3) with a temporary symbol s (see call of .to_rational(srl).numer_denom() in function matrix::fraction_free_elimination(bool), file matrix.cpp). Later on, that symbol appears in expanded polynomials like s^2-3. It's not immediately apparent for the program that this is zero. The program does this replacement because the divide function expects polynomials over the rationals as numerator and denominator.
Calling lsolve(eqns,syms,solve_algo::divfree) or lsolve(eqns,syms,solve_algo::gauss) will switch to alternative algorithms which are by construction not susceptible to this kind of problem. They may, however, turn out to be too slow for more densely populated systems of equations. You should give them a try.
Cheers -richy.
OK, now I understand. The alternative algorithms work fine. 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? Then one could have matrix::solve catch the exception and switch automatically to a different algorithm, ensuring the correct result is returned. Thanks again Diego