# 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(a_=b_,c_=d_, x_,y_) := solve(a-b,c-d, x,y);
solve(a_=b_, x_,y_) := solve(a-b, x,y);

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

solve(c_*g_=b_,x_):= If(isfree(b,c,x), replace(inverse(g,x),x,b/c) );
solve(e^a_=b_,x_):= If(isfree(b,x), solve(a=log(b),x) );
solve(e^x_=b_,x_):= If(isfree(b,x), log(b) );

solve(a_*b_,x_) := solve(a,x) and solve(b,x);

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,(-c/a)^(1/n) )));
solve(c_+x_^n_,x_) := If(isfree(c,x),  If(n==2, sqrt(-c) and -sqrt(-c), if(n>4,(-c)^(1/n) )));

solve(polys(a_,b_,x_),x_):= -b/a;
solve(polys(a_,b_,c_,x_),x_):= (-b+sqrt(b^2-4*a*c))/2/a and (-b-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(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(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(x_^3+x_, x_) := 0 and i and -i;
solve(x_^3-x_, x_) := 0 and 1 and -1;
solve(x_^3+c_*x_, x_) := 0 and sqrt(-c) and - sqrt(-c);
solve(x_^3+x_^2+x_+1, x_) := -1 and i and -i;

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

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(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(mod(p_,b_),x_):=mod(solve(p-b,x),b)+b*n;
solve(a_+mod(p_,b_),x_):=mod(solve(a+p-b,x),b)+b*n;
solve(mod(p_,b_)+z_,x_):=mod(solve(z+p-b,x),b)+b*n;
solve(a_+mod(p_,b_)+z_,x_):=mod(solve(a+p-b+z,x),b)+b*n;
solve(a_+c_*mod(p_,b_),x_):= mod(solve(a/c+p-b,x),b)+b*n;
solve(c_*mod(p_,b_)+z_,x_):= mod(solve(p-b+z/c,x),b)+b*n;
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_):= left(solve(c+b*yy+z,x,yy));
solve(a_+c_*mod(b_),x_):= left(solve(a+c+b*yy,x,yy));
solve(a_+c_*mod(b_)+z_,x_):= left(solve(a+c+b*yy+z,x,yy));

#solve(y(p_,x_)+b_,x_):=dsolve(y(p,x)+b,y);
#solve(a_*y(p_,x_)+b_,x_):=dsolve(a*y(p,x)+b,y);
#solve(y(p_,x_)^n_+b_,x_):= dsolve(solve(y(p,x)^n+b,y(p,x)),y,x,p);
#solve(a_*y(p_,x_)^n_+b_,x_):=dsolve(solve(y(p,x)^n+expand(b/a),y(p,x)),y,x,p);

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(l_=(r1_ and r2_),x_) := solve(l-r1,x) and solve(l-r2,x);
solve(left_=r_, x_) := solve(left-r, x);

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


#solve(f_) := solve(replace(toyn(f),y(x),y),x);
solve(f_) := solve(f,x);