#5655 linsolve returns wrong solutions for some systems with %i or radicals

[David Scherfgen]

By Claude:

linsolve returns a wrong solution for some systems whose coefficients contain %i or radicals together with a symbolic parameter. Substituting the solution back should give [0, 0, 0], here at c = 2:

(%i1) display2d:false
(%i2) eqs:[x+c*y+%i*z = 1+%i,c*x+sqrt(2)*y+z = 2,-x+y+sqrt(c)*z = %i]
(%i3) sol:linsolve(eqs,[x,y,z])
(%i4) float(rectform(subst(c = 2,subst(sol,map(lhs,eqs)-map(rhs,eqs)))))
(%o4) [3.3502382219814892*%i+3.7625658698293707,
       1.3322676295501878e-15*%i-6.38378239159465e-16,
       2.249515306418351-8.327112270558679*%i]

This is not rare: in random 3x3 systems with coefficients like these, 3 to 5 % of the solutions were wrong.

The Risch integrator solves linear systems the same way. Since the fix for bug [#5241], the following integral takes about 90 seconds (see bug [#5651]), and the result is wrong where cos(x) < 0, where its derivative is the negated integrand (relative error 2):

(%i5) g:((a+%i*a*tan(x))^(7/2)*(A+B*tan(x)))/(c-%i*c*tan(x))^(13/2)
(%i6) showtime:true
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%i7) G:integrate(g,x)
Evaluation took 93.1241 seconds (93.1640 elapsed) using 17003.017 MB.
(%i8) showtime:false
(%i9) vals:[a = 2,c = 3,A = 1,B = -1]
(%i10) relerr(p):=float(
              cabs(subst(x = p,subst(vals,diff(G,x)-g)))
               /cabs(subst(x = p,subst(vals,g))))
(%i11) map(relerr,[1/2,2,4])
(%o11) [4.403230523413117e-14,1.999999999999986,2.0000000000001528]

[David Scherfgen]

Fixed by commit [386344].