[GiNaC-list] lsolve problem

[Diego Conti]

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