[GiNaC-list] lsolve problem
diego.conti at unimib.it
Mon Feb 18 14:32:45 CET 2008
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.
> 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.
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.
More information about the GiNaC-list