[GiNaC-list] lsolve problem

[Richard B. Kreckel]

Hi!

Diego Conti wrote:
> 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.
> 
> code:
> 
> #include "ginac/ginac.h"
> #include <iostream>
> using namespace std;
> using namespace GiNaC;
> 
> int main()
> {
> symbol 
> lambda1("lambda1"),lambda2("lambda2"),lambda3("lambda3"),lambda4("lambda4"),lambda5("lambda5"), 
> lambda6("lambda6"),lambda7("lambda7"),    
> lambda8("lambda8"),lambda9("lambda9"),lambda10("lambda10"),lambda11("lambda11"),lambda12("lambda12"),lambda13("lambda13"), 
> lambda14("lambda14"),    
> lambda15("lambda15"),lambda16("lambda16"),lambda17("lambda17"),lambda18("lambda18"),lambda19("lambda19"),lambda20("lambda20"),lambda21("lambda21"),    
> lambda22("lambda22"),lambda23("lambda23"),lambda24("lambda24"),lambda25("lambda25"),lambda26("lambda26"), 
> lambda27("lambda27"),lambda28("lambda28"),lambda29("lambda29");
> ex sqrt3=sqrt(ex(3));
> lst syms;
> syms=lambda1,lambda2,lambda3,lambda4,lambda5,lambda6,lambda7,lambda8,lambda9,lambda10,lambda11,
> lambda12,lambda13,lambda14,lambda15,lambda16,lambda17,lambda18,lambda19,lambda20,lambda21,lambda22,lambda23,lambda24,
> lambda25,lambda26,lambda27,lambda28,lambda29;   
> 
>     lst eqns;
>     eqns.append(-lambda20+lambda24-lambda2-lambda25*sqrt3==0);
>     eqns.append(-sqrt3*lambda23-lambda22==0);
>     eqns.append(lambda15-lambda9-sqrt3*lambda7-lambda6==0);
>     eqns.append(lambda23+lambda1-2*lambda21==0);
>     eqns.append(lambda15-lambda9-lambda6==0);
>     eqns.append(-2*lambda28==0);
>     eqns.append(sqrt3*lambda12-lambda13==0);
>     eqns.append(lambda26-sqrt3*lambda27+lambda19+lambda8==0);
>     eqns.append(-lambda28==0);
>     eqns.append(-sqrt3*lambda19+lambda27+sqrt3*lambda8==0);
>     eqns.append(sqrt3*lambda17+lambda18==0);
>     eqns.append(lambda4-2*lambda17-lambda11==0);
>     eqns.append(lambda20+2*lambda24+lambda2==0);
>     eqns.append(-lambda16+lambda3+2*lambda12==0);
>     eqns.append(-sqrt3*lambda22-lambda23-lambda1-lambda21==0);
>     eqns.append(lambda9*sqrt3-lambda7+lambda15*sqrt3==0);
>     eqns.append(-sqrt3*lambda3-lambda13==0);
>     eqns.append(-lambda15+lambda9-sqrt3*lambda7+lambda6==0);
>     eqns.append(lambda25-sqrt3*lambda24==0);
>     eqns.append(-lambda5-lambda10-lambda14==0);
>     eqns.append(2*lambda20+lambda24-lambda2==0);
>     cout<<GiNaC::lsolve(eqns,syms)<<endl;
> }
> 
> output:
> {lambda1==lambda23,lambda2==-1/3*sqrt(3)*lambda25,lambda3==-lambda16,lambda4==-2/3*sqrt(3)*lambda18+lambda11,lambda5==-lambda14-lambda10,lambda6==2*lambda15,lambda7==0,lambda8==-1/2*lambda26+1/3*sqrt(3)*lambda27,lambda9==-lambda15,lambda10==lambda10,lambda11==lambda11,lambda12==lambda16,lambda13==sqrt(3)*lambda16,lambda14==lambda14,lambda15==lambda15,lambda16==lambda16,lambda17==-1/3*sqrt(3)*lambda18,lambda18==lambda18,lambda19==-1/6*(sqrt(3)*lambda26-4*lambda27)*sqrt(3),lambda20==-1/3*sqrt(3)*lambda25,lambda21==lambda23,lambda22==-lambda23*sqrt(3),lambda23==lambda23,lambda24==1/3*sqrt(3)*lambda25,lambda25==lambda25,lambda26==lambda26,lambda27==lambda27,lambda28==lambda28,lambda29==lambda29}
> 
> I suspect the problem is the sqrt(3), since replacing it with an integer 
> the result appears to be correct. Or am I doing something wrong?

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.
-- 
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>