#	functional equation solver;
#	solve functional, difference and recurrence equation for function y(x) with independent variable x;
#	e.g. rsolve( difference(y(x))=x ); rsolve( y(x+1)=y(x)+2 );

#2;
rsolve(y(xx_)+b_,x_):= -replace(b,x,inverse(xx,x));
rsolve(a_*y(xx_)+b_,x_):= -replace(b/a,x,inverse(xx,x));
rsolve(c_*y(a_*x_)+b_,x_):= -replace(b/c,x,x/a);
rsolve(y(a_*x_)+b_,x_):= -replace(b,x,x/a);
rsolve(a_*y(x_+n_)+b_,x_):= -replace(b/a,x,x-n);
rsolve(y(x_+a_)+b_,x_):= -replace(b,x,x-a);
rsolve(y(x_)+b_,x_):= -b;
rsolve(a_*y(x_)+b_,x_):= -b/a;

rsolve(y(a_)+y(b_),x_):=if(a*b==1,C_1*log(x));
rsolve(y(a_)-y(b_),x_):=if(a*b==1,C_1*log(-x)+C_2,C_1);
rsolve(y(b_)-y(a_),x_):=if(a*b==1,C_1*log(-x)+C_2,C_1);

rsolve(y(x_^a_)+b_*y(x_),x_):=if(a==b and a<0,C_1*log(-x),if(a+b==0,C_1*log(x) ));
rsolve(y(a_*x_)+b_*y(x_),x_):=if(a+b==0,C_1*x);

rsolve(y(x_+1)+a_*y(x_),x_):= if(hasnot(a,x),C_1*(-a)^x );
rsolve(y(x_+1)-y(x_)*x_,x_):=C_1*Gamma(x);
rsolve(y(1-x_)+y(-x_)*x_,x_):=C_1*Gamma(x);
rsolve(y(x_+1)-y(x_)*(1+x_),x_):= C_1*x!;
rsolve((-x_)*y(-1+x_)+y(x_),x_):=C_1*x!;

rsolve(a_*y(x_)+c_*y(b_+x_),x_):=if(hasnot(a,c,x),(-a/c)^(x/b)*C_);
rsolve(y(x_)+a_*y(b_+x_),x_):=if(hasnot(a,x),(-a)^(-x/b)*C_1);
rsolve(y(x_+b_)+a_*y(x_),x_):=if(hasnot(a,x),(-a)^(x/b)*C_1);

rsolve(-y(a_+b_)+y(a_)*y(b_),x_):=exp(x);
rsolve(-y(a_-b_)+y(a_)/y(b_),x_):=exp(x);
rsolve(y(a_+b_)-y(a_)*y(b_),x_):=exp(x);
rsolve(y(a_-b_)-y(a_)/y(b_),x_):=exp(x);

rsolve(y(x_^n_)-y(x_)^n_,x_):=x^C_1;
rsolve(y(2x_)-2y(x_)^2+1,x_):=cos(x);
rsolve(y(2x_^2-1)-2y(x_),x_):=C_1*acos(x);
rsolve(y(x_)*y(1-x_)-pi/sin(pi*x_),x_):=gamma(x);

rsolve(y(x_+2pi)-y(x_),x_):=C_1*trig(x);
rsolve(y(x_)-y(x_+2pi),x_):=C_1*trig(x);
rsolve(y(x_+pi)-y(x_),x_):=C_1*trig(2x);
rsolve(y(x_)-y(x_+pi),x_):=C_1*trig(2x);

rsolve(y(-xx_)+y(xx_),x_):=C_1*odd(x);
rsolve(y(-xx_)-y(xx_),x_):=C_1*even(x);
rsolve(y(xx_)-y(-xx_),x_):=C_1*even(x);
rsolve(y(-xx_)-1/y(xx_),x_):=exp(x);
rsolve(1/y(xx_)-y(-xx_),x_):=exp(x);
rsolve(-1+y(-xx_)*y(xx_),x_):=exp(x);
rsolve(y(f(x_))-f(y(x_)),x_):=x;


#3;
rsolve(y(a_)+y(b_)+f_,x_):=if(a*b==1,C_1*log(x)-f/2);
rsolve(y(x_^a_)+b_*y(x_)+f_,x_):=if(a==b and a<0,C_1*log(-x),if(a+b==0,C_1*log(x)-b*f*x/2 ));
rsolve(y(x_+1)+a_*y(x_)+d_,x_):= if(hasnot(a,x), C_1*(-a)^x+rpsolution(a,d,x) );

rsolve(-y(a_+b_)+y(a_)+y(b_),x_):=C_1*x;
rsolve(y(a_+b_)-y(a_)-y(b_),x_):=C_1*x;
rsolve(-y(a_*b_)+y(a_)+y(b_),x_):=C_1*log(x);
rsolve(y(a_*b_)-y(a_)-y(b_),x_):=C_1*log(x);
rsolve(y(x_*b_)-y(x_)*y(b_),x_):=x^n;
rsolve(y(x_+b_)-y(x_)*y(b_)/y(x_+b_),x_):=C_1/x;
rsolve(y(x_+b_)+y(x_-b_)-2y(x_)*y(b_),x_):=cos(x);

#rsolve(d_+a_*y(x_)+c_*y(b_+x_),x_):=if(hasnot(a,c,d,x),(-a/c)^(x/b)*C_1-d/(a+c));
#rsolve(d_+y(x_)+c_*y(b_+x_),x_):=if(hasnot(c,d,x),(-1/c)^(x/b)*C_1-d/(1+c));
#rsolve(y(x_+b_)+a_*y(x_)+c_,x_):=if(hasnot(a,c,x),(-a)^(x/b)*C_1-c/(a+1));
#rsolve(y(x_+b_)+y(x_)+c_,x_):=if(hasnot(c,x),(-1)^(x/b)*C_1-c/2);
rsolve(y(a_*x_)+b_*y(x_)+c_,x_):=if(a+b==0,C_1*x-c);
rsolve(y(x_)-y(b_+x_)+c_,x_):= if(b<0,sum(c,x)/b+C_1,replace(sum(c,x),x,-x)/b+C_1 );
rsolve(y(x_+b_)-y(x_)+c_,x_):= if(b<0,sum(c,x)/b+C_1,replace(sum(c,x),x,-x)/b+C_1 );
rsolve(y(x_+b_)-y(x_)+f_*x_^n_,x_):= if(b<0,sum(c,x)/b+C_1,f*zeta(-n,x)-f*zeta(-n,0)+C_1 );

#rsolve(y(x_+a_)+d_*y(x_+b_)+f_*x_^n_,x_):= rsolve(y(x+b-a)+d*y(x)+f*x^n,x);
#rsolve(y(x_+a_)+y(x_+b_)+f_*x_^n_,x_):= rsolve(y(x+b-a)+y(x)+f*x^n,x);
rsolve(y(x_+a_)+y(x_+b_)+f_,x_):= replace(rsolve(y(x+b-a)+y(x)+f,x),x,x-a);
rsolve(y(x_+a_)+d_*y(x_+b_)+f_,x_):= replace(rsolve(y(x+b-a)+d*y(x)+f,x),x,x-a);
rsolve(c_*y(x_+a_)+y(x_+b_)+f_,x_):= replace(rsolve(y(x+b-a)+c*y(x)+f,x),x,x-a);
rsolve(c_*y(x_+a_)+d_*y(x_+b_)+f_,x_):= replace(rsolve(d*y(x+b-a)+c*y(x)+f,x),x,x-a);
rsolve(y(2+x_)+c_*y(x_)+f_,x_):= if(hasnot(c,f,x), rsolution(1,2,0,1,c)+ if(c==-1, 1/4f*((-1)^(2 x) - 2 x), -f/(1+c)) );
rsolve(y(2+x_)+y(x_)+f_,x_):= if(hasnot(f,x), rsolution(1,2,0,1,1)-f/2 );

rsolve(a_*y(x_)+b_*y(1/x_)+c_*x_,x_):=if(hasnot(a,b,c,x), b*c/(a^2-b^2)/x-a*c/(a^2-b^2)*x);
rsolve(y(2+x_)+y(x_)+f_*x_,x_):= rsolution(1,2,0,1,1)-f/2*(x-1);
rsolve(y(2+x_)-y(x_)+f_*x_^n_,x_):= c_2 *(-1)^x + c_1 + 1/2 f *(zeta(-n, x) + 2^n *(-1)^(2 x) *(zeta(-n, x/2) - zeta(-n, (x + 1)/2))) ;
rsolve(y(2+x_)+y(x_)+x_,x_):= rsolution(1,2,0,1,1)-x/2+1/2;
rsolve(y(2+x_)+y(x_),x_):= rsolution(1,2,0,1,1);
rsolve(y(2+x_)+c_*y(x_),x_):= rsolution(1,2,0,1,c);
rsolve(y(2+x_)+b_*y(1+x_),x_):=rsolution(1,2,b,1,0);
rsolve(y(2+x_)+y(1+x_)+y(x_),x_):= rsolution(1,2,1,1,1);
rsolve(y(2+x_)+y(1+x_)+c_*y(x_),x_):=rsolution(1,2,1,1,c);
rsolve(y(2+x_)+b_*y(1+x_)+y(x_),x_):=rsolution(1,2,b,1,1);
rsolve(y(2+x_)+b_*y(1+x_)+c_*y(x_),x_):=rsolution(1,2,b,1,c);
rsolve(a_*y(2+x_)+b_*y(1+x_)+c_*y(x_),x_):=rsolution(a,2,b,1,c);

rsolve(y(x_+1)-y(x_)*x_-y(x_),x_):= C_1*x!;
rsolve(y(x_+1)-y(x_)-y(-1+x_),x_):=C_1*fibonacci(x);


4;
rsolve(y(2+x_)+b_*y(1+x_)+y(x_)+f_,x_):= if(hasnot(f,x), rsolution(1,2,b,1,1) + if(b==-2,-f*x^2/2,-f/(2+b) ));
rsolve(y(2+x_)+b_*y(1+x_)+c_*y(x_)+f_,x_):= if(hasnot(f,x), rsolution(1,2,b,1,c)+ if(1+b+c==0, -f/(2+b)*x, -f/(1+b+c) ));

rsolve(y(2+x_)+b_*y(1+x_)+y(x_)+f_*x_,x_):= if(hasnot(f,x), rsolution(1,2,b,1,1)+ if(b==-2, -f/2*x^2, (f* (-x *(b + 2) + b + 2))/(b + 2)^2  ));
rsolve(y(2+x_)+b_*y(1+x_)+c_*y(x_)+f_*x_,x_):= if(hasnot(f,x), rsolution(1,2,b,1,c)+ if(1+b+c==0, -f/(2+b)*x, (f* (-x *(b + c + 1) + b + 2))/(b + c + 1)^2  ));
rsolve(y(2+x_)+y(1+x_)+y(x_)+f_*x_,x_):= rsolution(1,2,1,1,1) - f*(x - 1)/3 ;
rsolve(y(2+x_)+y(1+x_)+y(x_)+x_,x_):= rsolution(1,2,1,1,1) - x/3+1/3 ;


#0;
rsolve(a_=b_,x_):=rsolve(a-b,x);
rsolve(a_):=rsolve(a,x);