forget()
def deg(p,n):
if p.degree(n) == -1:
return -Infinity
else:
return p.degree(n)
def const_coeffic(d, D, ff):
var ('i n L Lc LL CC c')
C=[var('C%s'%i) for i in range(d+1)]
DD=range(d+1)
DD[0] = expand(sum(C[i]*n^(d-i) for i in range(d+1)))
for j in range(1,d+1):
DD[j]=0
for i in range(d+1):
f(n)=DD[j-1].coefficient(C[i])
DD[j]+=C[i]*expand(f(n+1)-f(n))
L=sum(DD[i]*D[i] for i in range(d+1))-ff
Lc=L.coefficients(n)
LL=[]
for i in range(len(Lc)):
LL+=[Lc[i][0] == 0]
c=list(L.variables())
if n in c: c.remove(n)
if len(c)>0:
if len(c)==1:
c=c[0]
print LL,c
if len(solve(LL,c))>0:
solutions=solve(LL,c,solution_dict=True)
else:
return []
else:
return []
LL=range(len(solutions))
for j in range(len(solutions)):
LL[j]=0
for i in range(d+1):
if C[i] in solutions[j]:
LL[j]+=solutions[j][C[i]]*n^(d-i)
else:
LL[j]+=C[i]*n^(d-i)
return LL
def poly_solve(data, f):
var ('n i r f d dd c F L E C D DD N Nf B b P')
r = data.nrows()
c = data.ncols()
rf = len(f)
N = []
Nf = []
B = []
DD = []
R.<n,D,E> = QQ[]
S.<d> = QQ[]
for a in range(c):
N = [n^a]+N
for a in range(rf):
Nf = [n^a]+Nf
ee = data*vector(N)
f = sum(Nf[i]*f[i] for i in range(rf))
E = D+1
L= sum(ee[i]*E^(r-i-1) for i in range(r))
for j in range(r):
DD+=[L.coefficient({D:j})]
B=[deg(L.coefficient({D:j}),n)-j]+B
b = max(B)
P=0
for j in range(r):
if b==B[j]:
P=P+L.coefficient({D:r-j-1}).lc()
if j!=r-1:
P=P*(d-r+j+2)
a=P.roots(multiplicities=False,ring=ZZ)
if len(a) == 0:
a=[-1]
dd=max([deg(f,n)-b]+a)
if dd < 0 :
if f != 0 :
return []
else :
return [0]
else:
return const_coeffic(dd, DD, f)
#Testni primeri:
data = Matrix([[0,3],[-1,0],[1,-1]]);f=[0]
print data
print f
y = poly_solve(data, f)
print ' ...'
print y
print '----------'
data = Matrix([[2,1],[-4,-4],[2,3]]);f=[0]
print data
print f
y = poly_solve(data, f)
print ' ...'
print y
print '----------'
data = Matrix([[1,2],[-8,-10],[16,8]]);f=[9,0,-6]
print data
print f
y = poly_solve(data, f)
print ' ...'
print y
print '----------'
data = Matrix([[1,0],[0,-1],[2,1]]);f=[6,25,-2]
print data
print f
y = poly_solve(data, f)
print ' ...'
print y
print '----------'
data = Matrix([[2,2],[2,-1]]);f=[1]
print data
print f
y = poly_solve(data, f)
print ' ...'
print y
print '----------'
data = Matrix([[4,2],[1,0]]);f=[3,1]
print data
print f
y = poly_solve(data, f)
print ' ...'
print y
print '----------'
#(n-1) n^2 y(n+3) + (3 n^3 + 6 n^2 - 3 n - 2)(n + 1) y(n+1) - (n - 1)(n^3 + 6 n^2 + 4 n + 1) y(n+2) - 2 n (n + 1)^3 y(n) = n^5 + n^4 - n^2 - n
# + n^3*y(n + 3)
# - n^2*y(n + 3)
#
# - n^4*y(n + 2)
# - 5*n^3*y(n + 2)
# + 2*n^2*y(n + 2)
# + 3*n *y(n + 2)
# + y(n + 2)
#
# + 3*n^4*y(n + 1)
# + 9*n^3*y(n + 1)
# + 3*n^2*y(n + 1)
# - 5*n *y(n + 1)
# - 2 *y(n + 1)
#
# - 2*n^4*y(n)
# - 6*n^3*y(n)
# - 6*n^2*y(n)
# - 2*n *y(n)
data = Matrix([[0,1,-1,0,0],[-1,-5,2,3,1],[3,9,3,-5,-2],[-2,-6,-6,-2,0]]);f=[1,1,0,-1,-1,0]
print data
print f
y = poly_solve(data, f)
print ' ...'
print y
print '----------'
|
|
[ 0 3]
[-1 0]
[ 1 -1]
[0]
[12*C0 + 6*C1 + 2*C2 == 0, 11*C0 + C1 == 0] [C0, C1, C2]
...
[1/27*n^2*r2 - 11/27*n*r2 + r2]
----------
[ 2 1]
[-4 -4]
[ 2 3]
[0]
[-2*C1 == 0] C1
...
[C0*n^2 + C2]
----------
[ 1 2]
[ -8 -10]
[ 16 8]
[9, 0, -6]
[-6*C0 + 6 == 0, -6*C0 + 9*C1 == 0, 9*C0 - 9 == 0] [C0, C1]
...
[n + 2/3]
----------
[ 1 0]
[ 0 -1]
[ 2 1]
[6, 25, -2]
[-C0 + 2 == 0, 2*C0 + 3*C1 - 25 == 0, 3*C0 - 6 == 0] [C0, C1]
...
[2*n + 7]
----------
[ 2 2]
[ 2 -1]
[1]
...
[]
----------
[4 2]
[1 0]
[3, 1]
[2*C0 - 1 == 0, 5*C0 - 3 == 0] C0
...
[]
----------
[ 0 1 -1 0 0]
[-1 -5 2 3 1]
[ 3 9 3 -5 -2]
[-2 -6 -6 -2 0]
[1, 1, 0, -1, -1, 0]
[2*C0 - C2 == 0, 7*C0 - 4*C2 + 1 == 0, 3*C0 - 2*C2 + 1 == 0, 2*C0 - C2
== 0, C0 - 1 == 0, C0 - 1 == 0] [C0, C1, C2]
...
[n^2 + n*r4 + 2]
----------
[ 0 3]
[-1 0]
[ 1 -1]
[0]
[12*C0 + 6*C1 + 2*C2 == 0, 11*C0 + C1 == 0] [C0, C1, C2]
...
[1/27*n^2*r2 - 11/27*n*r2 + r2]
----------
[ 2 1]
[-4 -4]
[ 2 3]
[0]
[-2*C1 == 0] C1
...
[C0*n^2 + C2]
----------
[ 1 2]
[ -8 -10]
[ 16 8]
[9, 0, -6]
[-6*C0 + 6 == 0, -6*C0 + 9*C1 == 0, 9*C0 - 9 == 0] [C0, C1]
...
[n + 2/3]
----------
[ 1 0]
[ 0 -1]
[ 2 1]
[6, 25, -2]
[-C0 + 2 == 0, 2*C0 + 3*C1 - 25 == 0, 3*C0 - 6 == 0] [C0, C1]
...
[2*n + 7]
----------
[ 2 2]
[ 2 -1]
[1]
...
[]
----------
[4 2]
[1 0]
[3, 1]
[2*C0 - 1 == 0, 5*C0 - 3 == 0] C0
...
[]
----------
[ 0 1 -1 0 0]
[-1 -5 2 3 1]
[ 3 9 3 -5 -2]
[-2 -6 -6 -2 0]
[1, 1, 0, -1, -1, 0]
[2*C0 - C2 == 0, 7*C0 - 4*C2 + 1 == 0, 3*C0 - 2*C2 + 1 == 0, 2*C0 - C2 == 0, C0 - 1 == 0, C0 - 1 == 0] [C0, C1, C2]
...
[n^2 + n*r4 + 2]
----------
|