fichier programme bezout.py

1
#######################" fonctions python numpy scipy #########################
2

3
def _GH():
4
	Gx = []; Hx = []; Gy = []; Hy = []
5
	for j in range(n):
6
		dft = ft.fft(np.eye(2*fn*deg[j]))
7
		gx = dft[0::2*fn/(j+1), :deg[j]+1]
8
		gy = dft[1::2*fn/(n-j), :deg[j]+1]
9
		hx = dft[0::2*fn/(j+1), :dx[j]]
10
		hy = dft[1::2*fn/(n-j), :dy[j]]
11
		Gx.append(gx); Hx.append(hx); Gy.append(gy); Hy.append(hy)
12
	return Gx, Gy, Hx, Hy
13

14
def _HK():
15
	H = 1; K = 1
16
	for k in range(n):
17
		H = np.kron(Hx[k], H); K = np.kron(Hy[k], K)
18
	return H, K
19

20
def evaluate(i, j):
21
	f = F[i].transpose(range(n-1, -1, -1)).reshape(prod(fshape[j:]), prod(fshape[:j]))
22
	L = 1; R = 1
23
	for k in range(j):
24
		o = np.ones((dx[k], 1))
25
		L = np.kron(o, L); R = np.kron(Gy[k], R)
26
	for k in range(j, n):
27
		o = np.ones((dy[k], 1))
28
		R = np.kron(o, R); L = np.kron(Gx[k], L)
29
	return L.dot(f).dot(R.T)
30

31
def _J():
32
	J = np.empty((Dx, Dy, 2*n+1, n+1), dtype=complex)
33
	for i in range(2*n+1):
34
		for j in range(n+1):
35
			J[:, :, i, j] = evaluate(i, j)
36
	return J
37

38
def _C():
39
	C = np.empty((n+1, Dx, Dy), dtype=complex)
40
	indices_vol = range(n, 2*n+1)
41
	for i in range(n+1):
42
		indices = range(n) + [n+i]
43
		vol = np.linalg.det(J[:, :, indices_vol, :])
44
		num = np.linalg.det(J[:, :, indices, :])
45
		C[i] = num/vol
46
	return C
47

48
def _B():
49
	B = np.empty((n+1, Dx, Dy), dtype=complex)
50
	for i in range(n+1):
51
		HC  = conjugate(H.T).dot(C[i])/Dx
52
		KHC = conjugate(K.T).dot(HC.T)/Dy
53
		B[i] = KHC.T
54
	return np.around(B).real.astype(int)
55

56
def build_ixy(k):
57
	ix = ()
58
	iy = ()
59
	for j in range(n):
60
		if j < k:
61
			ix = ix + (slice(None),)
62
			iy = (slice(None),) + iy
63
		elif j == k:
64
			ix = ix + (slice(-deg[k], None),)
65
			iy = (slice(-deg[n-1-k], None),) + iy
66
		else:
67
			ix = ix + (slice(None, -deg[j]),)
68
			iy = (slice(None, -deg[n-1-j]),) + iy
69
	return ix, iy
70

71
def permut():
72
	aax = np.arange(Dx, dtype=int).reshape(dx[::-1]).transpose(range(n-1,-1,-1))
73
	aay = np.arange(Dy, dtype=int).reshape(dy[::-1]).transpose(range(n-1,-1,-1))
74
	ax = np.zeros(0, dtype=int)
75
	ay = np.zeros(0, dtype=int)
76
	for k in range(n):
77
		ix, iy = build_ixy(k)
78
		tx = (np.arange(n)*(n-1) + k) % n
79
		ty = (n-1 - tx) % n
80
		ay = np.concatenate(( ay, np.transpose(aay[iy], ty).reshape(Dy/n) ))
81
		ax = np.concatenate(( ax, np.transpose(aax[ix], tx).reshape(Dx/n) ))
82
	return ax, ay
83
	
84
def block_triang():
85
	ax, ay = permut()
86
	for k in range(n+1):
87
		B[k] = np.fliplr(B[k][np.ix_(ax,ay)])
88
	return B
89

90
def block_size():
91
	bls = np.zeros(0, dtype=int)
92
	for i in range(n):
93
		bls = np.concatenate(( bls, np.ones(deg[i], dtype=int)*Dx/(n*deg[i]) ))
94
	return bls
95

96
############################################# fcts sage #######################################
97

98
def rand_poly(j, m):
99
	# construit un polynome random en j+1 variables 0..j, de degres max(deg, m)
100
	p = 0
101
	for k in range(min(deg[j], m) + 1):
102
		if j > 0:
103
			p += rand_poly(j-1, m-k)*x[j]**k
104
		else:
105
			coeff = ZZ.random_element(0.45, distribution='gaussian')
106
			p += coeff*x[j]**k		
107
	return p
108

109
def poly2prism(p):
110
	z = np.zeros(fshape)
111
	mns = p.monomials(); cfs = p.coefficients()
112
	for k in range(len(mns)):
113
		z[mns[k].degrees()] = cfs[k]
114
	return z
115

116
def poly2sparse(P):
117
	degrees = np.zeros((0, n))
118
	coeffs = np.zeros((0))
119
	nb_monomials = np.zeros(n, dtype=int)
120
	for j in range(n):
121
		p = P[j]
122
		mns = p.monomials(); 
123
		degs = np.zeros((len(mns), n))
124
		for k in range(len(mns)):
125
			degs[k, :] = np.asarray(mns[k].degrees())
126
		cfs = np.asarray(p.coefficients())
127
		degrees = np.concatenate((degrees, degs))
128
		nb_monomials[j] = len(mns)
129
		coeffs = np.concatenate((coeffs, cfs))
130
	return degrees, coeffs, nb_monomials
131
		
132
def sparse2prism(sp):
133
	degrees = sp[0]
134
	l = degrees.shape[0]
135
	coeffs = sp[1]
136
	f = np.zeros(np.asarray(deg)+1, dtype=float)
137
	for k in range(l):
138
		f[tuple(degrees[k])] = coeffs[k]
139
	return f
140

141
def shrink(axe):
142
	iz = np.all(B == 0, axis=axe)
143
	pix = np.prod(iz, axis=0)
144
	ix = np.nonzero(1 - pix)[0]
145
	return list(ix)
146

147
def reduct():
148
	nr, nc = BB[0].nrows(), BB[0].ncols()
149
	ker_basis = BB[0].left_kernel().basis_matrix()
150
	print 'ker_size = ', ker_basis.nrows()
151
	relations = matrix(QQ, 0 , nc)
152
	for k in range(1, n+1):
153
		new_relat = ker_basis*BB[k]
154
		relations = block_matrix(2, 1, [relations, new_relat])
155
	r = rank(relations)
156
	print 'r = ', r
157
	rkbt = relations.right_kernel().basis_matrix().transpose()
158
	BB_out = []
159
	for k in range(n+1):
160
		BB_out.append(BB[k]*rkbt)
161
	den = BB_out[0].denominator()
162
	print 'den = ', den
163
	return r, BB_out
164

165
def _Bred():
166
	nr, nc = BB[0].nrows(), BB[0].ncols()
167
	d = rank(BB[0])
168
	print 'd =', d
169
	piv_cols, piv_rows = BB[0].pivots(), BB[0].transpose().pivots()
170
	Bred = np.empty((n+1, d, d), dtype=float)
171
	for j in range(n+1):
172
		Bred[j] = BB[j][piv_rows, piv_cols] + matrix(RR, d, d)
173
	return Bred
174

175
def _XX_chow():
176
	nr, nc = BB[0].nrows(), BB[0].ncols()
177
	d = rank(BB[0])
178
	print 'd =', d
179
	piv_cols, piv_rows = BB[0].pivots(), BB[0].transpose().pivots()
180
	nBB = np.empty((n+1, d, d), dtype=float)
181
	XX = np.empty((n, d, d), dtype=float)
182
	for k in range(n+1):
183
		nBB[k, :, :] = np.array(BB[k][piv_rows, piv_cols])
184
	for k in range(n):
185
		XX[k, :, :] = la.solve(nBB[0], nBB[k+1])
186
	chow_mat = np.zeros((d, d))
187
	for k in range(n):
188
		chow_mat += np.random.random()*nBB[k+1]
189
	Ex = la.eig(chow_mat, nBB[0])
190
	Wx = Ex[1]
191
	X = np.empty((n, d), dtype=complex)
192
	for k in range(n):
193
		X[k] = la.solve(Wx, XX[k].dot(Wx)).diagonal(offset=0)
194
	return XX, X
195

196

197
def _jPZ():
198
	d = X.shape[1]
199
	jz = np.empty((n, n), dtype=complex)
200
	Pz = np.empty((n, 1), dtype=complex)
201
	jPZ = np.empty((n, d), dtype=complex)
202
	for k in range(d):
203
		z = tuple(X[:, k])
204
		for i in range(n):
205
			Pz[i] = P[i](z)
206
#			for j in range(n):
207
#				jz[i, j] = jac[i, j](z)
208
#		if abs(la.det(jz)) > 1e-4:
209
#			jPZ[:, k:k+1] = la.solve(jz, Pz)
210
#		else:
211
		jPZ[:, k:k+1] = Pz
212
	return jPZ
213

214
##################" debut programme python / sage ######################
215

216
import matplotlib.pyplot as plt
217
import numpy as np
218
import scipy.linalg as la
219
import scipy.fftpack as ft
220
import time
221
import sys
222
import scipy.io as sio
223

224
deg = [5, 4, 4]
225
m = 20
226
n = len(deg)
227
R = PolynomialRing(QQ, 'x', n)
228
x = R.gens()
229
xx = [x[0]**0] + list(x)
230

231
fshape = [d+1 for d in deg]
232
dx = [(i+1)*deg[i] for i in range(n)]
233
dy = [(n-i)*deg[i] for i in range(n)]
234
fn, Dx, Dy = factorial(n), prod(dx), prod(dy)
235

236
P = [rand_poly(n-1, m) for i in range(n)] + xx
237
#~ a = [ZZ.random_element() for _ in range(7)]
238
#~ P = [a[0]*x[0]**3*x[1]**2 + a[1]*x[0] + a[2]*x[1]**2 + a[3], a[4]*x[0]*x[1]**4 + a[5]*x[0]**3 + a[6]*x[1]] + xx
239
degrees, coeffs, nb_monomials = poly2sparse(P)
240

241
jac = jacobian(P[:n], x)
242
F = [poly2prism(p) for p in P]
243

244
Gx, Gy, Hx, Hy = _GH()
245
H, K = _HK()
246
J = _J()
247
C = _C()
248
B = _B()
249

250
B = block_triang()
251

252

253
print 'debut sage'
254
BB = []
255
for k in range(n+1):
256
	Bk = matrix(QQ, B[k])
257
	BB.append(Bk[:, :])
258
r = 1
259
while r > 0:
260
	r, BB = reduct()
261

262
Bred = _Bred()
263

264
bls = block_size()
265

266
sio.savemat('np_B.mat', {'Bred':np.transpose(Bred.astype(float), (1, 2, 0)), 'B':np.transpose(B.astype(float), (1, 2, 0)), 'deg':deg, 'bls':bls, 'degrees':degrees, 'coeffs':coeffs, 'nb_monomials':nb_monomials})
267

268
XX, X = _XX_chow()
269

270
jPZ = _jPZ()
271

272
plt.close()
273
#plt.plot(np.log10(2**-52 + abs(np.transpose(jPZ[:,:]))))
274
for i in range(n):
275
	hh = plt.hist(np.log10(2**-52 + abs(jPZ[i])), 50)
276
plt.grid();plt.savefig('ref.png'); plt.show()
277

278

279
I = R.ideal(P[:n])
280
dim = I.vector_space_dimension()
281
print 'dim =', dim
282

283

284