restart;


po:=(a,b)->pochhammer(a,b):


parterrdenom1:=proc()
global u1, denrrp1, derrp1:
local vv, numfacdif, ii, hh, facdenon, facden, denomu1: 
u1:=(n,k)->simplify(bb(n+1,k)/bb(n,k)):
denomu1:=(n,k)->denom(u1(n,k)); facdenon:=factors(denomu1(n,k))[2]; vv:=1:
numfacdif:=nops(facdenon(n,k));
for ii from 1 to numfacdif do; facden:=coeff(facdenon[ii][1],k,0):
if subs({n=-1,k=0},facden)=0 then; hh:=1; else: 
hh:=(facdenon[ii][1])^(facdenon[ii][2]); fi:
vv:=vv*hh; denrrp1:=(n,k)->vv; 
od; derrp1:=(o,i)->subs(n=o,k=i,denrrp1(n,k)):
end:


parterrdenom2:=proc()
global u3, denrrp2, derrp2:
local vv, numfacdif, ii, hh, facnumer, facnum, denomu2:
u3:=(n,k)->simplify(bb(n,k)/bb(n-1,k)):
denomu2:=(n,k)->numer(u3(n,k)); facnumer:=factors(denomu2(n,k))[2]; vv:=1:
numfacdif:=nops(facnumer(n,k)):
for ii from 1 to numfacdif do; facnum:=coeff(facnumer[ii][1],k,0):
if subs({n=0,k=0},facnum)=0 then:  
hh:=(facnumer[ii][1])^(facnumer[ii][2]); else; hh:=1; fi:
vv:=vv*hh; denrrp2:=(n,k)->vv: 
od; derrp2:=(o,i)->subs(n=o,k=i,denrrp2(n,k)):
end:


parterrdenom:=proc()
global derr:
parterrdenom1(); parterrdenom2(): 
derr:=(n,k)->derrp1(n,k)*derrp2(n,k):
end:


parterrnumer:=proc(aaa,bbb,ccc,ddd)
global numrr1, rr1, rrr1, inderr1; local gradodenrr:
gradodenrr:=grado+degree(derr(n,k)):
numrr1:=(n,k)->k*sum(sum(dd(i,j)*n^i*k^(j-1),i=0..gradodenrr-j),j=1..gradodenrr)+(aaa+bbb*n+ccc*n^2+ddd*n^3)*derr(n,0);
rr1:=(n,k)->numrr1(n,k)/derr(n,k); rrr1:=(o,x)->subs({n=o,k=x},rr1(n,k)):
inderr1:={seq(seq(dd(ii,jj),ii=0..gradodenrr-jj),jj=1..gradodenrr)}:
return(rrr1(n,k)):
end:


parterr:=proc(aaa,bbb,ccc,ddd)
parterrdenom(); parterrnumer(aaa,bbb,ccc,ddd):
end:


partessdenom1:=proc()
global u2, indess1, denssp1, dessp1:
local vv, numfacdif, ii, hh, facdenon, facden, denomu2: 
u2:=(n,k)->simplify(bb(n,k+1)/bb(n,k)); 
denomu2:=(n,k)->denom(u2(n,k)); facdenon:=factors(denomu2(n,k))[2]; vv:=1;
numfacdif:=nops(facdenon(n,k)):
for ii from 1 to numfacdif do; facden:=coeff(facdenon[ii][1],n,0);
if subs({n=0,k=-1},facden)=0 then; hh:=1; else; hh:=(facdenon[ii][1])^(facdenon[ii][2]); fi; 
vv:=vv*hh; denssp1:=(n,k)->vv:
od; dessp1:=(o,i)->subs(n=o,k=i,denssp1(n,k)):
end:


partessdenom2:=proc()
global u4, denssp2, dessp2:
local vv, numfacdif, ii, hh, facnumer, facnum, denomu2:
u4:=(n,k)->simplify(bb(n,k)/bb(n,k-1)); 
denomu2:=(n,k)->numer(u4(n,k)); facnumer:=factors(denomu2(n,k))[2]; vv:=1;
numfacdif:=nops(facnumer(n,k)):
for ii from 1 to numfacdif do; facnum:=coeff(facnumer[ii][1],n,0):
if subs({n=0,k=0},facnum)=0 then; hh:=(facnumer[ii][1])^(facnumer[ii][2]); else; hh:=1; fi: 
vv:=vv*hh; denssp2:=(n,k)->vv:
od; dessp2:=(o,i)->subs(n=o,k=i,denssp2(n,k)):
end:


partessdenom:=proc()
global dess:
partessdenom1(); partessdenom2():
dess:=(n,k)->dessp1(n,k)*dessp2(n,k):
end:


partessnumer:=proc()
global numss1, ss1, sss1, indess1: 
local gradodenss:
gradodenss:=grado-1+degree(dess(n,k)):
numss1:=(n,k)->sum(sum(pp(i,j)*n^(i+1)*k^j,i=0..gradodenss-j),j=0..gradodenss):
ss1:=(n,k)->numss1(n,k)/dess(n,k); sss1:=(o,x)->subs({n=o,k=x},ss1(n,k)):
indess1:={seq(seq(pp(ii,jj),ii=0..gradodenss-jj),jj=0..gradodenss)}:
return(sss1(n,k)):
end:


partess:=proc()
partessdenom(); partessnumer():
end:


resolver:=proc(zz,aaa,bbb,ccc,ddd)
global lafor, ecus, indeter, numincog, sols, yyy, zzz:
parterr(aaa,bbb,ccc,ddd); partess():
lafor:=(n,k)->u1(n,k)*zz*sss1(n+1,k)-sss1(n,k)-u2(n,k)*yy*rrr1(n,k+1)+rrr1(n,k):
indeter:=inderr1 union indess1 union {yy}; numincog:=nops(indeter):
ecus:={coeffs(collect(numer(lafor(n,k)),{n,k},distributed),{n,k})}:
sols:=solve(ecus,indeter); yyy:=subs(sols,yy); zzz:=subs(sols,zz):
if type(yyy,numeric) then print(subs(sols,{y=yy,z=zz})); print(R(n,k)=factor(simplify(subs(sols,rrr1(n,k))))); 
print(S(n,k)=factor(simplify(subs(sols,sss1(n,k))))); else; print("no tiene solucion"); fi;
end: