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