# solve(eq) equation for default unknown x;
# solve(eq,x) equation for unknown x;
# solve(eq,x,y) equation for unknown x and y;
# solve(eq,eq2,x,y) solve system of equations eq and eq2 for unknown x and y;
#solve(c_+x_^2+x_^4, x_) := If(isfree(c,x), (sqrt(-1/2+sqrt(1-4*c)/2) and sqrt(-1/2-sqrt(1-4*c)/2) and -sqrt(-1/2+sqrt(1-4*c)/2) and -sqrt(-1/2-sqrt(1-4*c)/2)));
#solve(c_+b_*x_^2+x_^4, x_) := If(isfree(b,c,x), (sqrt((-b+sqrt(b^2-4*c))/2) and sqrt((-b-sqrt(b^2-4*c))/2) and -sqrt((-b+sqrt(b^2-4*c))/2) and -sqrt((-b-sqrt(b^2-4*c))/2)));
#solve(c_+b_*x_^2+a_*x_^4, x_) := If(isfree(a,b,c,x), (sqrt((-b+sqrt(b^2-4*c*a))/(2a)) and sqrt((-b-sqrt(b^2-4*c*a))/(2a)) and -sqrt((-b+sqrt(b^2-4*c*a))/(2a)) and -sqrt((-b-sqrt(b^2-4*c*a))/(2a))));


solve(a_=b_, x_,y_,z_) := solve(a-b, x,y,z);
solve(x_^2+y_^2-z_^2, x_,y_,z_):= x=3 and y=4 and z=5;

solve(a_=b_,c_=d_, x_,y_,z_) := solve(a-b,c-d, x,y,z);
solve(x_^2-y_^2+b_,y^2-z^2+(-9), x_,y_,z_):= x=sqrt(25-b) and y=5 and z=4;
solve(x_^2-y_^2+b_,y^2-z^2+(-16), x_,y_,z_):= x=sqrt(25-b) and y=5 and z=3;
solve(x_^2-y_^2+(-16),y^2-z^2+c_, x_,y_,z_):= x=5 and y=3 and z=sqrt(9+c);


solve(a_=b_,c_=d_, x_,y_) := solve(a-b,c-d, x,y);
solve(x_+y_+b_,x_*y_+c_, x_,y_) := x=-b/2+sqrt(b^2+4c)/2 and y=-b/2-sqrt(b^2+4c)/2;
solve(x_^2-y_^2+b_,x_*y_+c_, x_,y_) := x=-b/2+sqrt(b^2+4c)/2 and y=-b/2-sqrt(b^2+4c)/2;
solve(x_^2+2x_*y_+y_^2+b_,x_*y_+c_, x_,y_) := x=-sqrt(-b)/2+sqrt(-b+4c)/2 and y=-sqrt(-b)/2-sqrt(-b+4c)/2;

solve(a_=b_, x_,y_) := solve(a-b, x,y);
solve(a_+y_^2, x_,y_) := (x=solve(replace(a,y,0),x) and y=0, x=0 and y=solve(replace(a,x,0)+y^2,y));
solve(a_+c_*y_^2, x_,y_) := (x=solve(replace(a,y,0),x) and y=0, x=0 and y=solve(replace(a,x,0)+c*y^2,y));
solve(a_+c_*y_, x_,y_) :=  block(xx:=mod(solve(a-c,x),c), if(isinteger(left(xx)), (x=left(xx)+c*n and y=expand(-replace(a/c,x,left(xx)+c*n))) ));
solve(a_+y_, x_,y_) := block(xx:=solve(a-1,x), if(isinteger(left(xx)), x=left(xx)+n and y=expand(-replace(a,x,left(xx)+n)) ));
solve(x_+y_+f_, x_,y_) := x=1+n and y=-f-1-n;
solve(c_+b_*x_+a_*y_+x_*y_, x_,y_) := If(isconstant(a,b,c), x=sqrt(a*b-c)-a and y=sqrt(a*b-c)-b);
solve(c_+b_*x_+a_*y_+d_*x_*y_, x_,y_) := If(isconstant(a,b,c,d), x=sqrt(a*b-c*d)/d-a/d and y=sqrt(a*b-c*d)/d-b/d);
solve(x_^2+b_*x_*y_+c_*y_^2+f_, x_,y_):=if(b==2*(c)^0.5,solve(x+sqrt(c)*y+sqrt(-f),x,y));
solve(a_*x_^2+b_*x_*y_+y_^2+f_, x_,y_):=if(b==2*(a)^0.5,solve(sqrt(a)*x+y+sqrt(-f),x,y));
solve(a_*x_^2+b_*x_*y_+c_*y_^2+f_, x_,y_):=if(b==2*(a*c)^0.5,solve(sqrt(a)*x+sqrt(c)*y+sqrt(-f),x,y));
solve((a_+x_)*(b_+y_)+c_, x_,y_) := If(isconstant(a,b,c), x=sqrt(-c)-a and y=sqrt(-c)-b);
solve(a_*x_+b_*y_+c_*x_*y_,x_,y_):= (x=0 and y=0, x=-(a+b)/c and y=-(a+b)/c);
solve(x_^2+b_*y_^2+f_,x_,y_):= x=sqrt(-f-b*n^2) and y=n;
solve(x_^2+2x_*y_+y_^2+f_,x_,y_):= solve(x+y-sqrt(-f),x,y);
solve(x_^2+(-2)*x_*y_+y_^2+f_,x_,y_):= solve(x-y-sqrt(-f),x,y);
solve(x_^2+y_^2+(-25),x_,y_):= x=3 and y=4;
solve(x_^2-y_^2+(-9),x_,y_):= x=5 and y=4;
solve(x_^2-y_^2+(-16),x_,y_):= x=5 and y=3;

solve(a_=b_,c_=d_,x_) := solve(a-b,c-d, x);

solve(r1_ and r2_,x_) := (solve(r1,x) and solve(r2,x));

solve(b_+a_*x_, x_) := If(isfree(a,b,x), -b/a);
solve(c_+a_*x_^n_, x_) := If(isfree(c,a,n,x), If(n==2, sqrt(-c/a) and -sqrt(-c/a), 
	if(n==3, -cbrt(c/a) and (-1)^(1/3)*cbrt(c/a) and -(-1)^(2/3)*cbrt(c/a), 
	If(n==4, sqrt(-c/a) and -sqrt(-c/a) and -i*sqrt(-c/a) and i*sqrt(-c/a), 
	(-c/a)^(1/n) ))));
solve(c_+x_^n_,x_) := If(isfree(c,x), If(n==2, sqrt(-c) and -sqrt(-c), 
	if(n==3, -cbrt(c) and (-1)^(1/3)*cbrt(c) and -(-1)^(2/3)*cbrt(c), 
	If(n==4, (-c)^(1/4) and -(-c)^(1/4) and -i*(-c)^(1/4) and i*(-c)^(1/4), 
	(-c)^(1/n) ))));
solve(c_*x_+a_*x_^n_,x_) := If(isfree(a,c,x), (-c/a)^(1/(n-1)) );
solve(c_*x_^m_+a_*x_^n_,x_) := If(isfree(a,c,x), (-c/a)^(1/(n-m)) );
solve(c_*x_+x_^n_,x_) := If(isfree(c,x), (-c)^(1/(n-1)) );
solve(c_*x_^m_+x_^n_,x_) := If(isfree(c,x), (-c)^(1/(n-m)) );
solve(x_+a_*x_^n_,x_) := If(isfree(a,x), (-1/a)^(1/(n-1)) );

solve(polys(a_,b_,x_),x_):= -b/a;
solve(polys(a_,b_,c_,x_),x_):= -b/2/a+sqrt(b^2-4*a*c)/2/a and -b/2/a-sqrt(b^2-4*a*c)/2/a;
solve(polys(a_,b_,c_,d_,x_),x_):=   block( 
  r1:=(sqrt((-27 a^2 *d + 9 a *b *c - 2 b^3)^2 + 4 *(3 a *c - b^2)^3) - 27 a^2 *d + 9 a *b* c - 2 b^3)^(1/3)/(3* 2^(1/3) *a), 
  r2:=(sqrt((-27 a^2* d + 9 a* b* c - 2 b^3)^2 + 4* (3 a* c - b^2)^3) + 27 a^2* d - 9 a* b* c + 2 b^3)^(1/3)/(3* 2^(1/3) *a), 
  r1-r2 - b/(3 a) and r1*(-1 + I *sqrt(3))/2-r2*(-1 - I *sqrt(3))/2 - b/(3 a) and r1*(-1 - I *sqrt(3))/2-r2*(-1 + I *sqrt(3))/2 - b/(3 a));
solve(polys(a_,0,b_,0,c_,x_), x_) := If(isfree(a,b,c,x), (sqrt((-b+sqrt(b^2-4*c*a))/(2a)) and sqrt((-b-sqrt(b^2-4*c*a))/(2a)) and -sqrt((-b+sqrt(b^2-4*c*a))/(2a)) and -sqrt((-b-sqrt(b^2-4*c*a))/(2a))));
solve(polys(1,0,0,0,c_,x_), x_) :=  (-c)^(1/4) and -(-c)^(1/4)  and (-c)^(1/4)*i  and -(-c)^(1/4)*i;
solve(polys(1,0,1,0), x_) := (0 and i and -i);
solve(polys(1,0,c_,0), x_) := (0 and sqrt(-c) and -sqrt(-c));
solve(polys(1,1,1,1), x_) := (-1 and i and -i);

solve(y_+c_*y_+a_,y_):= if(hasnot(c,a,y), -a/(1+c) );
solve(b_*y_+1/y_+a_,y_):= if(hasnot(b,a,y), (sqrt(a^2-4b)-a)/2/b and (-sqrt(a^2-4b)-a)/2/b );
solve(b_*y_+c_*y_+a_,y_):= if(hasnot(b,c,a,y), -a/(b+c) );
solve(b_*y_+c_/y_+a_,y_):= if(hasnot(b,c,a,y), (sqrt(a^2-4b*c)-a)/2/b and (-sqrt(a^2-4b*c)-a)/2/b );

solve(x_^3+c_*x_, x_) := (0 and sqrt(-c) and -sqrt(-c));
solve(x_^3+x_, x_) := (0 and i and -i);
solve(x_^3+x_^2+x_+1, x_) :=  (-1 and i and -i);

solve(a_*x_, x_) := 0;
solve(x_*y_, x_) := 0;
solve(b*(c_+x_), x_) := If(isfree(b,c,x), -c);


#solve(a_+abs(f_),x_) := If(a<0, solve(a-f,x) and solve(a+f,x));
#solve(a_+abs(x_),x_) := If(a<0, a and (-a));
#solve(a_+c_*abs(f_),x_) := If(c*a<0, solve(a/c-f,x) and solve(f+a/c,x));
#solve(a_+c_*abs(x_),x_) := If(c*a<0, a/c and -a/c);
solve(a_+abs(f_)+abs(g_),x_) := If(a<0, solve(a-f-g,x) and solve(a+f+g,x) );
solve(a_+cos(x_)+sin(x_),x_) :=asin(-sqrt(2)/2*a)-pi/4;

solve(vector(a_,b_)*x_+c_,x_):= -c/vector(a,b);
solve(vector(a_,b_)*x_+vector(c_,d_)*y_+vector(f_,g_),x_,y_):= solve(a*x+c*y+f,b*x+d*y+g,x,y);
#solve(vector(a_,b_),x_):= if(hasnot(a,b,y),if(solve(a,x)==solve(b,x),solve(a,x)),solve(a,b,x,y) );

solve(c_*mod(p_,b_)+z_,x_):= mod(solve(p-b+z/c,x),b);
solve(a_+c_*mod(p_,b_)+z_,x_):= mod(solve(a/c+p-b+z/c,x),b)+b*n;
solve(c_*mod(b_)+z_,x_):= solve(z+c-b,x)+b*n;
solve(mod(p_,b_)+z_,x_):=mod(solve(z+p-b,x),b);
solve(a_+mod(p_,b_)+z_,x_):=mod(solve(a+p-b+z,x),b)+b*n;
solve(a_+c_*mod(b_)+z_,x_):= solve(a+c+z-b,x)+b*n;
solve(a_+c_*mod(b_)+y_+z_,x_):= solve(a+c+y+z-b,x)+b*n;
solve(mod(p_,b_),x_):= if(has(p,x),mod(solve(p-b,x),b)+b*n, if(has(b,x),solve(b-p,x) ));
solve(a_+mod(p_,b_),x_):=mod(solve(p+a-b,x),b)+b*n;
solve(a_+mod(b_+c_*x_,m_),x_):= if(mod(a+b,c)==0, -(a+b)/c+m*n, if(a+b+1==0,inversemod(c,m)+m*n, mod((m-a-b)/c,m)+m*n ));
solve(a_+mod(c_*x_,m_),x_):= if(mod(a,c)==0, -a/c+m*n, if(a+1==0,inversemod(c,m)+m*n, mod((m-a)/c,m)+m*n ));
solve(a_+mod(x_,b_),x_):= -a+b*n;
solve(a_+c_*mod(p_,b_),x_):= mod(solve(a/c+p-b,x),b)+b*n;
solve(a_+c_*mod(b_),x_):= solve(a+c-b,x)+b*n;
solve(a_*x_+c_*mod(b_),x_):=mod(-c*inversemod(a,b),b)+b*n;
solve(a_*x_*mod(b_)+c_,x_):=mod(-c*inversemod(a,b),b)+b*n;
solve(mod(x_,b_)+mod(x_,c_),x_):= b*c*n;
solve(a_+mod(x_,b_)+mod(x_,c_),x_):= if(iseven(a) and b> -a/2,-a/2+b*c*n, 
	if(a+b==0,b+b*c*n, 
	if((-a+b)/2-c==0,b*c+a+b+b*c*n 
	)));

solve(a_+mod(x_,b_),a2_+mod(x_,c_),x_):= if(a+b==a2+c,b*c-a-b+b*c*n);
solve(a_+mod(x_,b_),a_+mod(x_,c_),x_):= -a+b*c*n;

solve(a_+b_*exp(c_+x_)+z_,x_) := -c-log((-a-z)/b);

#solve(log(f_)+z_,y_) :=solve(f-exp(-z),y);
#solve(a_+log(f_),y_) :=solve(f-exp(-a),y);
#solve(a_+log(f_)+z_,y_) :=solve(f-exp(-a-z),y);
#solve(a_+b_*log(f_),y_) :=solve(f-exp(-a/b),y);
solve(c_+b_*log(y_)+y_,y_) := b*W(exp(-c/b)/b);
solve(c_+log(y_)+y_,y_) := W(exp(-c));
solve(c_+log(y_)+a_*y_,y_) := W(a*exp(-c))/a;
solve(c_+log(y_)+x_+y_,y_) := W(exp(-c-x));
solve(c_+b_*log(y_)+x_+y_,y_) := b*W(exp(-c/b-x/b)/b);


solve(left_>r_, x_) := block(s:=solve(left-r, x),  if(max(s)<>min(s),x< min(s) and x>max(s), x>s ));
solve(left_>=r_, x_) := block(s:=solve(left-r, x), if(max(s)<>min(s),x<=min(s) and x>=max(s), x>=s ));
solve(left_<r_, x_) := block(s:=solve(left-r, x),  if(max(s)<>min(s),x> min(s) and x<max(s), x<s ));
solve(left_<=r_, x_) := block(s:=solve(left-r, x), if(max(s)<>min(s),x>=min(s) and x<=max(s), x<=s ));

solve(eq2_=eq1_, y_) := if(eq1<>0 and hasnot(eq1,y) and hasnot(eq2,abs(_)) and hasnot(eq2-eq1,mod(_)), replace(inverse(eq2,y),y,eq1), solve(eq2-eq1,y) );
solve(a_+b_=c_,y_):=solve(a+b-c,y);
solve(a_*x_+b_*x_=c_,y_):=solve((a+b)*x=c,y);
#solve(left_=r_, x_) :=  solve(left-r, x);
solve(c_*g_=b_,x_):= if(has(b,mod()),solve(c*g-b),
	if(isfree(b,c,x), solve(g=b/c,x), 
	If(isfree(b,g,x), solve(c=b/g,x), 
	If(isfree(b,x),  replace(inverse(c*g,x),x,b),
	solve(c*g-b) ))));
solve(c_*log(b_+exp(y_))+d_*y_=x_,y_) :=if(c+d==0, -log(-exp(-b*x)/b-1/b) );
solve(c_*log(b_+a_*exp(y_))+d_*y_=x_,y_) :=if(c+d==0,-log(-exp(-b*x)/b-a/b) );
solve(-log(1+a_*exp(y_))+y_=x_,y_) := -log(-exp(-x)-a);
solve(-log(1+exp(y_))+y_=x_,y_) := -log(-exp(-x)-1);
solve(log(-1-exp(y_))-y_=x_,y_) := -log(exp(x)-1);
solve(log(1+exp(y_))-y_=x_,y_) := -log(exp(x)-1);
solve(log(-1+exp(y_))-y_=x_,y_) := -log(exp(x)+1);
solve(log(y_)+b_*y_=x_,y_):= 1/b*W(b*exp(x));
solve(a_*log(y_)+y_=x_,y_):= a*W(1/a*exp(x/a));
solve(a_*log(y_)+b_*y_=x_,y_):= a/b*W(b/a*exp(x/a));
solve(log(y_)*y_-y_=x_,y_):= exp(W(x/exp(1))+1);
solve(log(y_)*y_+y_=x_,y_):= exp(W(x*exp(1))-1);
solve(exp(b_*y_)+a_*y_=x_,y_):= if(a+b==0, x/a+W(-exp(-x))/a );
solve(exp(a_*y_)+a_*y_=x_,y_):= x/a-W(exp(x))/a;
solve(exp(y_)+y_=x_,y_):= x-W(exp(x));
solve(exp(-y_)+y_=x_,y_):= x+W(-exp(-x));

solve(l_=(r1_ and r2_),x_) := solve(l-r1,x) and solve(l-r2,x);

solve(a_=b_,c_=d_) := if(has(a=b,y), solve(a-b,c-d, x,y), solve(a-b,c-d, x) );
solve(f_) := solve(f,x);