1 | /////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="$Id: poly.lib,v 1.36 2005-05-06 14:38:55 hannes Exp $"; |
---|
3 | category="General purpose"; |
---|
4 | info=" |
---|
5 | LIBRARY: poly.lib Procedures for Manipulating Polys, Ideals, Modules |
---|
6 | AUTHORS: O. Bachmann, G.-M: Greuel, A. Fruehbis |
---|
7 | |
---|
8 | PROCEDURES: |
---|
9 | cyclic(int); ideal of cyclic n-roots |
---|
10 | katsura([i]); katsura [i] ideal |
---|
11 | freerank(poly/...) rank of coker(input) if coker is free else -1 |
---|
12 | is_zero(poly/...); int, =1 resp. =0 if coker(input) is 0 resp. not |
---|
13 | lcm(ideal); lcm of given generators of ideal |
---|
14 | maxcoef(poly/...); maximal length of coefficient occurring in poly/... |
---|
15 | maxdeg(poly/...); int/intmat = degree/s of terms of maximal order |
---|
16 | maxdeg1(poly/...); int = [weighted] maximal degree of input |
---|
17 | mindeg(poly/...); int/intmat = degree/s of terms of minimal order |
---|
18 | mindeg1(poly/...); int = [weighted] minimal degree of input |
---|
19 | normalize(poly/...); normalize poly/... such that leading coefficient is 1 |
---|
20 | rad_con(p,I); check radical containment of poly p in ideal I |
---|
21 | content(f); content of polynomial/vector f |
---|
22 | numerator(n); numerator of number n |
---|
23 | denominator(n) denominator of number n |
---|
24 | mod2id(M,iv); conversion of a module M to an ideal |
---|
25 | id2mod(i,iv); conversion inverse to mod2id |
---|
26 | substitute(I,...) substitute in I variables by polynomials |
---|
27 | subrInterred(i1,i2,iv);interred w.r.t. a subset of variables |
---|
28 | hilbPoly( I) Hilbert polynomial of basering/I |
---|
29 | (parameters in square brackets [] are optional) |
---|
30 | "; |
---|
31 | |
---|
32 | LIB "general.lib"; |
---|
33 | LIB "ring.lib"; |
---|
34 | /////////////////////////////////////////////////////////////////////////////// |
---|
35 | static proc bino(int a, int b) |
---|
36 | { |
---|
37 | //computes binomial var(1)+a over b |
---|
38 | int i; |
---|
39 | if(b==0){return(1);} |
---|
40 | poly p=(var(1)+a)/b; |
---|
41 | for(i=1;i<=b-1;i++) |
---|
42 | { |
---|
43 | p=p*(var(1)+a-i)/i; |
---|
44 | } |
---|
45 | return(p); |
---|
46 | } |
---|
47 | |
---|
48 | proc hilbPoly(ideal I) |
---|
49 | "USAGE: hilbPoly(I) I a homogeneous ideal |
---|
50 | RETURN: the Hilbert polynomial of basering/I as an intvec v=v_0,...,v_r |
---|
51 | such that the Hilbert polynomial is (v_0+v_1*t+...v_r*t^r)/r! |
---|
52 | EXAMPLE: example hilbPoly; shows an example |
---|
53 | " |
---|
54 | { |
---|
55 | def R=basering; |
---|
56 | if(!attrib(I,"isSB")){I=std(I);} |
---|
57 | intvec v=hilb(I,2); |
---|
58 | int s=dim(I); |
---|
59 | intvec hp; |
---|
60 | if(s==0){return(hp);} |
---|
61 | int d=size(v)-2; |
---|
62 | ring S=0,t,dp; |
---|
63 | poly p=v[1+d]*bino(s-1-d,s-1); |
---|
64 | int i; |
---|
65 | for(i=1;i<=d;i++) |
---|
66 | { |
---|
67 | p=p+v[d-i+1]*bino(s-1-d+i,s-1); |
---|
68 | } |
---|
69 | int n=1; |
---|
70 | for(i=2;i<=s-1;i++){n=n*i;} |
---|
71 | p=n*p; |
---|
72 | for(i=1;i<=s;i++) |
---|
73 | { |
---|
74 | hp[i]=int(leadcoef(p[s-i+1])); |
---|
75 | } |
---|
76 | setring R; |
---|
77 | return(hp); |
---|
78 | } |
---|
79 | example |
---|
80 | { "EXAMPLE:"; echo = 2; |
---|
81 | ring r = 0,(b,c,t,h),dp; |
---|
82 | ideal I= |
---|
83 | bct-t2h+2th2+h3, |
---|
84 | bt3-ct3-t4+b2th+c2th-2bt2h+2ct2h+2t3h-bch2-2bth2+2cth2+2th3, |
---|
85 | b2c2+bt2h-ct2h-t3h+b2h2+2bch2+c2h2-2bth2+2cth2+t2h2-2bh3+2ch3+2th3+3h4, |
---|
86 | c2t3+ct4-c3th-2c2t2h-2ct3h-t4h+bc2h2-2c2th2-bt2h2+4t3h2+2bth3-2cth3-t2h3 |
---|
87 | +bh4-6th4-2h5; |
---|
88 | hilbPoly(I); |
---|
89 | } |
---|
90 | |
---|
91 | /////////////////////////////////////////////////////////////////////////////// |
---|
92 | proc substitute (I,list #) |
---|
93 | "USAGE: - case 1: typeof(#[1])==poly: |
---|
94 | substitute (I,v,f[,v1,f1,v2,f2,...]); I object of basering which |
---|
95 | can be mapped, v,v1,v2,.. ring variables, f,f1,f2,... poly |
---|
96 | @* - case 2: typeof(#[1])==ideal: |
---|
97 | substitute1 (I,v,f); I object of basering which can be mapped, |
---|
98 | v ideal of ring variables, f ideal |
---|
99 | RETURN: object of same type as I, |
---|
100 | @* - case 1: ring variable v,v1,v2,... substituted by polynomials |
---|
101 | f,f1,f2,..., in this order |
---|
102 | @* - case 2: ring variables in v substituted by polynomials in f: |
---|
103 | v[i] is substituted by f[i], i=1,...,i=min(size(v),ncols(f)) |
---|
104 | NOTE: this procedure extends the built-in command subst which substitutes |
---|
105 | ring variables only by monomials |
---|
106 | EXAMPLE: example substitute; shows an example |
---|
107 | " |
---|
108 | |
---|
109 | { |
---|
110 | def bas = basering; |
---|
111 | ideal m = maxideal(1); |
---|
112 | int i,ii; |
---|
113 | if(typeof(#[1])=="poly") |
---|
114 | { |
---|
115 | poly v = #[1]; |
---|
116 | poly f = #[2]; |
---|
117 | map phi = bas,m; |
---|
118 | def J = I; |
---|
119 | for (ii=1; ii<=size(#) - 1; ii=ii+2) |
---|
120 | { |
---|
121 | m = maxideal(1); |
---|
122 | i=rvar(#[ii]); |
---|
123 | m[i] = #[ii+1]; |
---|
124 | phi = bas,m; |
---|
125 | J = phi(J); |
---|
126 | } |
---|
127 | return(J); |
---|
128 | } |
---|
129 | if(typeof(#[1])=="ideal") |
---|
130 | { |
---|
131 | ideal v = #[1]; |
---|
132 | ideal f = #[2]; |
---|
133 | int mi = size(v); |
---|
134 | if(ncols(f)<mi) |
---|
135 | { |
---|
136 | mi = ncols(f); |
---|
137 | } |
---|
138 | m[rvar(v[1])]=f[1]; |
---|
139 | map phi = bas,m; |
---|
140 | def J = phi(I); |
---|
141 | for (ii=2; ii<=mi; ii++) |
---|
142 | { |
---|
143 | m = maxideal(1); |
---|
144 | m[rvar(v[ii])]=f[ii]; |
---|
145 | phi = bas,m; |
---|
146 | J = phi(J); |
---|
147 | } |
---|
148 | return(J); |
---|
149 | } |
---|
150 | } |
---|
151 | example |
---|
152 | { "EXAMPLE:"; echo = 2; |
---|
153 | ring r = 0,(b,c,t),dp; |
---|
154 | ideal I = -bc+4b2c2t,bc2t-5b2c; |
---|
155 | substitute(I,c,b+c,t,0,b,b-1); |
---|
156 | ideal v = c,t,b; |
---|
157 | ideal f = b+c,0,b-1; |
---|
158 | substitute(I,v,f); |
---|
159 | } |
---|
160 | /////////////////////////////////////////////////////////////////////////////// |
---|
161 | proc cyclic (int n) |
---|
162 | "USAGE: cyclic(n); n integer |
---|
163 | RETURN: ideal of cyclic n-roots from 1-st n variables of basering |
---|
164 | EXAMPLE: example cyclic; shows examples |
---|
165 | " |
---|
166 | { |
---|
167 | //----------------------------- procedure body -------------------------------- |
---|
168 | ideal m = maxideal(1); |
---|
169 | m = m[1..n],m[1..n]; |
---|
170 | int i,j; |
---|
171 | ideal s; poly t; |
---|
172 | for ( j=0; j<=n-2; j=j+1 ) |
---|
173 | { |
---|
174 | t=0; |
---|
175 | for( i=1;i<=n;i=i+1 ) { t=t+product(m,i..i+j); } |
---|
176 | s=s+t; |
---|
177 | } |
---|
178 | s=s,product(m,1..n)-1; |
---|
179 | return (s); |
---|
180 | } |
---|
181 | //-------------------------------- examples ----------------------------------- |
---|
182 | example |
---|
183 | { "EXAMPLE:"; echo = 2; |
---|
184 | ring r=0,(u,v,w,x,y,z),lp; |
---|
185 | cyclic(nvars(basering)); |
---|
186 | homog(cyclic(5),z); |
---|
187 | } |
---|
188 | /////////////////////////////////////////////////////////////////////////////// |
---|
189 | |
---|
190 | proc katsura |
---|
191 | "USAGE: katsura([n]): n integer |
---|
192 | RETURN: katsura(n) : n-th katsura ideal of |
---|
193 | (1) newly created and set ring (32003, x(0..n), dp), if |
---|
194 | nvars(basering) < n |
---|
195 | (2) basering, if nvars(basering) >= n |
---|
196 | katsura() : katsura ideal of basering |
---|
197 | EXAMPLE: example katsura; shows examples |
---|
198 | " |
---|
199 | { |
---|
200 | int n; |
---|
201 | if ( size(#) == 1 && typeof(#[1]) == "int") |
---|
202 | { |
---|
203 | n = #[1] - 1; |
---|
204 | while (1) |
---|
205 | { |
---|
206 | if (defined(basering)) |
---|
207 | { |
---|
208 | if (nvars(basering) >= #[1]) {break;} |
---|
209 | } |
---|
210 | ring katsura_ring = 32003, x(0..#[1]), dp; |
---|
211 | keepring katsura_ring; |
---|
212 | break; |
---|
213 | } |
---|
214 | } |
---|
215 | else |
---|
216 | { |
---|
217 | n = nvars(basering) -1; |
---|
218 | } |
---|
219 | |
---|
220 | ideal s; |
---|
221 | int i, j; |
---|
222 | poly p; |
---|
223 | |
---|
224 | p = -1; |
---|
225 | for (i = -n; i <= n; i++) |
---|
226 | { |
---|
227 | p = p + kat_var(i, n); |
---|
228 | } |
---|
229 | s[1] = p; |
---|
230 | |
---|
231 | for (i = 0; i < n; i++) |
---|
232 | { |
---|
233 | p = -1 * kat_var(i,n); |
---|
234 | for (j = -n; j <= n; j++) |
---|
235 | { |
---|
236 | p = p + kat_var(j,n) * kat_var(i-j, n); |
---|
237 | } |
---|
238 | s = s,p; |
---|
239 | } |
---|
240 | return (s); |
---|
241 | } |
---|
242 | //-------------------------------- examples ----------------------------------- |
---|
243 | example |
---|
244 | { |
---|
245 | "EXAMPLE:"; echo = 2; |
---|
246 | ring r; basering; |
---|
247 | katsura(); |
---|
248 | katsura(4); basering; |
---|
249 | } |
---|
250 | |
---|
251 | proc kat_var(int i, int n) |
---|
252 | { |
---|
253 | poly p; |
---|
254 | if (i < 0) { i = -i;} |
---|
255 | if (i <= n) { p = var(i+1); } |
---|
256 | return (p); |
---|
257 | } |
---|
258 | /////////////////////////////////////////////////////////////////////////////// |
---|
259 | |
---|
260 | proc freerank |
---|
261 | "USAGE: freerank(M[,any]); M=poly/ideal/vector/module/matrix |
---|
262 | COMPUTE: rank of module presented by M in case it is free. |
---|
263 | By definition this is vdim(coker(M)/m*coker(M)) if coker(M) |
---|
264 | is free, where m = maximal ideal of the variables of the |
---|
265 | basering and M is considered as matrix. |
---|
266 | (the 0-module is free of rank 0) |
---|
267 | RETURN: rank of coker(M) if coker(M) is free and -1 else; |
---|
268 | in case of a second argument return a list: |
---|
269 | L[1] = rank of coker(M) or -1 |
---|
270 | L[2] = minbase(M) |
---|
271 | NOTE: freerank(syz(M)); computes the rank of M if M is free (and -1 else) |
---|
272 | EXAMPLE: example freerank; shows examples |
---|
273 | " |
---|
274 | { |
---|
275 | int rk; |
---|
276 | def M = simplify(#[1],10); |
---|
277 | resolution mre = res(M,2); |
---|
278 | intmat B = betti(mre); |
---|
279 | if ( ncols(B)>1 ) { rk = -1; } |
---|
280 | else { rk = sum(B[1..nrows(B),1]); } |
---|
281 | if (size(#) == 2) { list L=rk,mre[1]; return(L);} |
---|
282 | return(rk); |
---|
283 | } |
---|
284 | example |
---|
285 | {"EXAMPLE"; echo=2; |
---|
286 | ring r; |
---|
287 | ideal i=x; |
---|
288 | module M=[x,0,1],[-x,0,-1]; |
---|
289 | freerank(M); // should be 2, coker(M) is not free |
---|
290 | freerank(syz (M),""); |
---|
291 | // [1] should be 1, coker(syz(M))=M is free of rank 1 |
---|
292 | // [2] should be gen(2)+gen(1) (minimal relation of M) |
---|
293 | freerank(i); |
---|
294 | freerank(syz(i)); // should be 1, coker(syz(i))=i is free of rank 1 |
---|
295 | } |
---|
296 | /////////////////////////////////////////////////////////////////////////////// |
---|
297 | |
---|
298 | proc is_zero |
---|
299 | "USAGE: is_zero(M[,any]); M=poly/ideal/vector/module/matrix |
---|
300 | RETURN: integer, 1 if coker(M)=0 resp. 0 if coker(M)!=0, where M is |
---|
301 | considered as matrix. |
---|
302 | If a second argument is given, return a list: |
---|
303 | L[1] = 1 if coker(M)=0 resp. 0 if coker(M)!=0 |
---|
304 | L[2] = dim(M) |
---|
305 | EXAMPLE: example is_zero; shows examples |
---|
306 | " |
---|
307 | { |
---|
308 | int d=dim(std(#[1])); |
---|
309 | int a = ( d==-1 ); |
---|
310 | if( size(#) >1 ) { list L=a,d; return(L); } |
---|
311 | return(a); |
---|
312 | } |
---|
313 | example |
---|
314 | { "EXAMPLE:"; echo=2; |
---|
315 | ring r; |
---|
316 | module m = [x],[y],[1,z]; |
---|
317 | is_zero(m,1); |
---|
318 | qring q = std(ideal(x2+y3+z2)); |
---|
319 | ideal j = x2+y3+z2-37; |
---|
320 | is_zero(j); |
---|
321 | } |
---|
322 | /////////////////////////////////////////////////////////////////////////////// |
---|
323 | |
---|
324 | proc maxcoef (f) |
---|
325 | "USAGE: maxcoef(f); f poly/ideal/vector/module/matrix |
---|
326 | RETURN: maximal length of coefficient of f of type int (by counting the |
---|
327 | length of the string of each coefficient) |
---|
328 | EXAMPLE: example maxcoef; shows examples |
---|
329 | " |
---|
330 | { |
---|
331 | //----------------------------- procedure body -------------------------------- |
---|
332 | int max,s,ii,jj; string t; |
---|
333 | ideal i = ideal(matrix(f)); |
---|
334 | i = simplify(i,6); // delete 0's and keep first of equal elements |
---|
335 | poly m = var(1); matrix C; |
---|
336 | for (ii=2;ii<=nvars(basering);ii=ii+1) { m = m*var(ii); } |
---|
337 | for (ii=1; ii<=size(i); ii=ii+1) |
---|
338 | { |
---|
339 | C = coef(i[ii],m); |
---|
340 | for (jj=1; jj<=ncols(C); jj=jj+1) |
---|
341 | { |
---|
342 | t = string(C[2,jj]); s = size(t); |
---|
343 | if ( t[1] == "-" ) { s = s - 1; } |
---|
344 | if ( s > max ) { max = s; } |
---|
345 | } |
---|
346 | } |
---|
347 | return(max); |
---|
348 | } |
---|
349 | //-------------------------------- examples ----------------------------------- |
---|
350 | example |
---|
351 | { "EXAMPLE:"; echo = 2; |
---|
352 | ring r= 0,(x,y,z),ds; |
---|
353 | poly g = 345x2-1234567890y+7/4z; |
---|
354 | maxcoef(g); |
---|
355 | ideal i = g,10/1234567890; |
---|
356 | maxcoef(i); |
---|
357 | // since i[2]=1/123456789 |
---|
358 | } |
---|
359 | /////////////////////////////////////////////////////////////////////////////// |
---|
360 | |
---|
361 | proc maxdeg (id) |
---|
362 | "USAGE: maxdeg(id); id poly/ideal/vector/module/matrix |
---|
363 | RETURN: int/intmat, each component equals maximal degree of monomials in the |
---|
364 | corresponding component of id, independent of ring ordering |
---|
365 | (maxdeg of each var is 1). |
---|
366 | Of type int if id is of type poly, of type intmat else |
---|
367 | NOTE: proc maxdeg1 returns 1 integer, the absolute maximum; moreover, it has |
---|
368 | an option for computing weighted degrees |
---|
369 | EXAMPLE: example maxdeg; shows examples |
---|
370 | " |
---|
371 | { |
---|
372 | //-------- subprocedure to find maximal degree of given component ---------- |
---|
373 | proc findmaxdeg |
---|
374 | { |
---|
375 | poly c = #[1]; |
---|
376 | if (c==0) { return(-1); } |
---|
377 | //--- guess upper 'o' and lower 'u' bound, in case of negative weights ----- |
---|
378 | int d = (deg(c)>=0)*deg(c)-(deg(c)<0)*deg(c); |
---|
379 | int i = d; |
---|
380 | while ( c-jet(c,i) != 0 ) { i = 2*(i+1); } |
---|
381 | int o = i-1; |
---|
382 | int u = (d != i)*((i / 2)-1); |
---|
383 | //----------------------- "quick search" for maxdeg ------------------------ |
---|
384 | while ( (c-jet(c,i)==0)*(c-jet(c,i-1)!=0) == 0) |
---|
385 | { |
---|
386 | i = (o+1+u) / 2; |
---|
387 | if (c-jet(c,i)!=0) { u = i+1; } |
---|
388 | else { o = i-1; } |
---|
389 | } |
---|
390 | return(i); |
---|
391 | } |
---|
392 | //------------------------------ main program --------------------------------- |
---|
393 | matrix M = matrix(id); |
---|
394 | int r,c = nrows(M), ncols(M); int i,j; |
---|
395 | intmat m[r][c]; |
---|
396 | for (i=r; i>0; i=i-1) |
---|
397 | { |
---|
398 | for (j=c; j>0; j=j-1) { m[i,j] = findmaxdeg(M[i,j]); } |
---|
399 | } |
---|
400 | if (typeof(id)=="poly") { return(m[1,1]); } |
---|
401 | return(m); |
---|
402 | } |
---|
403 | //-------------------------------- examples ----------------------------------- |
---|
404 | example |
---|
405 | { "EXAMPLE:"; echo = 2; |
---|
406 | ring r = 0,(x,y,z),wp(1,2,3); |
---|
407 | poly f = x+y2+z3; |
---|
408 | deg(f); //deg; returns weighted degree (in case of 1 block)! |
---|
409 | maxdeg(f); |
---|
410 | matrix m[2][2]=f+x10,1,0,f^2; |
---|
411 | maxdeg(m); |
---|
412 | } |
---|
413 | /////////////////////////////////////////////////////////////////////////////// |
---|
414 | |
---|
415 | proc maxdeg1 (id,list #) |
---|
416 | "USAGE: maxdeg1(id[,v]); id=poly/ideal/vector/module/matrix, v=intvec |
---|
417 | RETURN: integer, maximal [weighted] degree of monomials of id independent of |
---|
418 | ring ordering, maxdeg1 of i-th variable is v[i] (default: v=1..1). |
---|
419 | NOTE: This proc returns one integer while maxdeg returns, in general, |
---|
420 | a matrix of integers. For one polynomial and if no intvec v is given |
---|
421 | maxdeg is faster |
---|
422 | EXAMPLE: example maxdeg1; shows examples |
---|
423 | " |
---|
424 | { |
---|
425 | //-------- subprocedure to find maximal degree of given component ---------- |
---|
426 | proc findmaxdeg |
---|
427 | { |
---|
428 | poly c = #[1]; |
---|
429 | if (c==0) { return(-1); } |
---|
430 | intvec v = #[2]; |
---|
431 | //--- guess upper 'o' and lower 'u' bound, in case of negative weights ----- |
---|
432 | int d = (deg(c)>=0)*deg(c)-(deg(c)<0)*deg(c); |
---|
433 | int i = d; |
---|
434 | if ( c == jet(c,-1,v)) //case: maxdeg is negative |
---|
435 | { |
---|
436 | i = -d; |
---|
437 | while ( c == jet(c,i,v) ) { i = 2*(i-1); } |
---|
438 | int o = (d != -i)*((i / 2)+2) - 1; |
---|
439 | int u = i+1; |
---|
440 | int e = -1; |
---|
441 | } |
---|
442 | else //case: maxdeg is nonnegative |
---|
443 | { |
---|
444 | while ( c != jet(c,i,v) ) { i = 2*(i+1); } |
---|
445 | int o = i-1; |
---|
446 | int u = (d != i)*((i / 2)-1); |
---|
447 | int e = 1; |
---|
448 | } |
---|
449 | //----------------------- "quick search" for maxdeg ------------------------ |
---|
450 | while ( ( c==jet(c,i,v) )*( c!=jet(c,i-1,v) ) == 0 ) |
---|
451 | { |
---|
452 | i = (o+e+u) / 2; |
---|
453 | if ( c!=jet(c,i,v) ) { u = i+1; } |
---|
454 | else { o = i-1; } |
---|
455 | } |
---|
456 | return(i); |
---|
457 | } |
---|
458 | //------------------------------ main program --------------------------------- |
---|
459 | ideal M = simplify(ideal(matrix(id)),8); //delete scalar multiples from id |
---|
460 | int c = ncols(M); |
---|
461 | int i,n; |
---|
462 | if( size(#)==0 ) |
---|
463 | { |
---|
464 | int m = maxdeg(M[c]); |
---|
465 | for (i=c-1; i>0; i=i-1) |
---|
466 | { |
---|
467 | n = maxdeg(M[i]); |
---|
468 | m = (m>=n)*m + (m<n)*n; //let m be the maximum of m and n |
---|
469 | } |
---|
470 | } |
---|
471 | else |
---|
472 | { |
---|
473 | intvec v=#[1]; //weight vector for the variables |
---|
474 | int m = findmaxdeg(M[c],v); |
---|
475 | for (i=c-1; i>0; i--) |
---|
476 | { |
---|
477 | n = findmaxdeg(M[i],v); |
---|
478 | if( n>m ) { m=n; } |
---|
479 | } |
---|
480 | } |
---|
481 | return(m); |
---|
482 | } |
---|
483 | //-------------------------------- examples ----------------------------------- |
---|
484 | example |
---|
485 | { "EXAMPLE:"; echo = 2; |
---|
486 | ring r = 0,(x,y,z),wp(1,2,3); |
---|
487 | poly f = x+y2+z3; |
---|
488 | deg(f); //deg returns weighted degree (in case of 1 block)! |
---|
489 | maxdeg1(f); |
---|
490 | intvec v = ringweights(r); |
---|
491 | maxdeg1(f,v); //weighted maximal degree |
---|
492 | matrix m[2][2]=f+x10,1,0,f^2; |
---|
493 | maxdeg1(m,v); //absolute weighted maximal degree |
---|
494 | } |
---|
495 | /////////////////////////////////////////////////////////////////////////////// |
---|
496 | |
---|
497 | proc mindeg (id) |
---|
498 | "USAGE: mindeg(id); id poly/ideal/vector/module/matrix |
---|
499 | RETURN: minimal degree/s of monomials of id, independent of ring ordering |
---|
500 | (mindeg of each variable is 1) of type int if id of type poly, else |
---|
501 | of type intmat. |
---|
502 | NOTE: proc mindeg1 returns one integer, the absolute minimum; moreover it |
---|
503 | has an option for computing weighted degrees. |
---|
504 | EXAMPLE: example mindeg; shows examples |
---|
505 | " |
---|
506 | { |
---|
507 | //--------- subprocedure to find minimal degree of given component --------- |
---|
508 | proc findmindeg |
---|
509 | { |
---|
510 | poly c = #[1]; |
---|
511 | if (c==0) { return(-1); } |
---|
512 | //--- guess upper 'o' and lower 'u' bound, in case of negative weights ----- |
---|
513 | int d = (ord(c)>=0)*ord(c)-(ord(c)<0)*ord(c); |
---|
514 | int i = d; |
---|
515 | while ( jet(c,i) == 0 ) { i = 2*(i+1); } |
---|
516 | int o = i-1; |
---|
517 | int u = (d != i)*((i / 2)-1); |
---|
518 | //----------------------- "quick search" for mindeg ------------------------ |
---|
519 | while ( (jet(c,u)==0)*(jet(c,o)!=0) ) |
---|
520 | { |
---|
521 | i = (o+u) / 2; |
---|
522 | if (jet(c,i)==0) { u = i+1; } |
---|
523 | else { o = i-1; } |
---|
524 | } |
---|
525 | if (jet(c,u)!=0) { return(u); } |
---|
526 | else { return(o+1); } |
---|
527 | } |
---|
528 | //------------------------------ main program --------------------------------- |
---|
529 | matrix M = matrix(id); |
---|
530 | int r,c = nrows(M), ncols(M); int i,j; |
---|
531 | intmat m[r][c]; |
---|
532 | for (i=r; i>0; i=i-1) |
---|
533 | { |
---|
534 | for (j=c; j>0; j=j-1) { m[i,j] = findmindeg(M[i,j]); } |
---|
535 | } |
---|
536 | if (typeof(id)=="poly") { return(m[1,1]); } |
---|
537 | return(m); |
---|
538 | } |
---|
539 | //-------------------------------- examples ----------------------------------- |
---|
540 | example |
---|
541 | { "EXAMPLE:"; echo = 2; |
---|
542 | ring r = 0,(x,y,z),ls; |
---|
543 | poly f = x5+y2+z3; |
---|
544 | ord(f); // ord returns weighted order of leading term! |
---|
545 | mindeg(f); // computes minimal degree |
---|
546 | matrix m[2][2]=x10,1,0,f^2; |
---|
547 | mindeg(m); // computes matrix of minimum degrees |
---|
548 | } |
---|
549 | /////////////////////////////////////////////////////////////////////////////// |
---|
550 | |
---|
551 | proc mindeg1 (id, list #) |
---|
552 | "USAGE: mindeg1(id[,v]); id=poly/ideal/vector/module/matrix, v=intvec |
---|
553 | RETURN: integer, minimal [weighted] degree of monomials of id independent of |
---|
554 | ring ordering, mindeg1 of i-th variable is v[i] (default v=1..1). |
---|
555 | NOTE: This proc returns one integer while mindeg returns, in general, |
---|
556 | a matrix of integers. For one polynomial and if no intvec v is given |
---|
557 | mindeg is faster. |
---|
558 | EXAMPLE: example mindeg1; shows examples |
---|
559 | " |
---|
560 | { |
---|
561 | //--------- subprocedure to find minimal degree of given component --------- |
---|
562 | proc findmindeg |
---|
563 | { |
---|
564 | poly c = #[1]; |
---|
565 | intvec v = #[2]; |
---|
566 | if (c==0) { return(-1); } |
---|
567 | //--- guess upper 'o' and lower 'u' bound, in case of negative weights ----- |
---|
568 | int d = (ord(c)>=0)*ord(c)-(ord(c)<0)*ord(c); |
---|
569 | int i = d; |
---|
570 | if ( jet(c,-1,v) !=0 ) //case: mindeg is negative |
---|
571 | { |
---|
572 | i = -d; |
---|
573 | while ( jet(c,i,v) != 0 ) { i = 2*(i-1); } |
---|
574 | int o = (d != -i)*((i / 2)+2) - 1; |
---|
575 | int u = i+1; |
---|
576 | int e = -1; i=u; |
---|
577 | } |
---|
578 | else //case: inded is nonnegative |
---|
579 | { |
---|
580 | while ( jet(c,i,v) == 0 ) { i = 2*(i+1); } |
---|
581 | int o = i-1; |
---|
582 | int u = (d != i)*((i / 2)-1); |
---|
583 | int e = 1; i=u; |
---|
584 | } |
---|
585 | //----------------------- "quick search" for mindeg ------------------------ |
---|
586 | while ( (jet(c,i-1,v)==0)*(jet(c,i,v)!=0) == 0 ) |
---|
587 | { |
---|
588 | i = (o+e+u) / 2; |
---|
589 | if (jet(c,i,v)==0) { u = i+1; } |
---|
590 | else { o = i-1; } |
---|
591 | } |
---|
592 | return(i); |
---|
593 | } |
---|
594 | //------------------------------ main program --------------------------------- |
---|
595 | ideal M = simplify(ideal(matrix(id)),8); //delete scalar multiples from id |
---|
596 | int c = ncols(M); |
---|
597 | int i,n; |
---|
598 | if( size(#)==0 ) |
---|
599 | { |
---|
600 | int m = mindeg(M[c]); |
---|
601 | for (i=c-1; i>0; i=i-1) |
---|
602 | { |
---|
603 | n = mindeg(M[i]); |
---|
604 | m = (m<=n)*m + (m>n)*n; //let m be the maximum of m and n |
---|
605 | } |
---|
606 | } |
---|
607 | else |
---|
608 | { |
---|
609 | intvec v=#[1]; //weight vector for the variables |
---|
610 | int m = findmindeg(M[c],v); |
---|
611 | for (i=c-1; i>0; i=i-1) |
---|
612 | { |
---|
613 | n = findmindeg(M[i],v); |
---|
614 | m = (m<=n)*m + (m>n)*n; //let m be the maximum of m and n |
---|
615 | } |
---|
616 | } |
---|
617 | return(m); |
---|
618 | } |
---|
619 | //-------------------------------- examples ----------------------------------- |
---|
620 | example |
---|
621 | { "EXAMPLE:"; echo = 2; |
---|
622 | ring r = 0,(x,y,z),ls; |
---|
623 | poly f = x5+y2+z3; |
---|
624 | ord(f); // ord returns weighted order of leading term! |
---|
625 | intvec v = 1,-3,2; |
---|
626 | mindeg1(f,v); // computes minimal weighted degree |
---|
627 | matrix m[2][2]=x10,1,0,f^2; |
---|
628 | mindeg1(m,1..3); // computes absolute minimum of weighted degrees |
---|
629 | } |
---|
630 | /////////////////////////////////////////////////////////////////////////////// |
---|
631 | |
---|
632 | proc normalize (id) |
---|
633 | "USAGE: normalize(id); id=poly/vector/ideal/module |
---|
634 | RETURN: object of same type with leading coefficient equal to 1 |
---|
635 | EXAMPLE: example normalize; shows an example |
---|
636 | " |
---|
637 | { |
---|
638 | return(simplify(id,1)); |
---|
639 | } |
---|
640 | //-------------------------------- examples ----------------------------------- |
---|
641 | example |
---|
642 | { "EXAMPLE:"; echo = 2; |
---|
643 | ring r = 0,(x,y,z),ls; |
---|
644 | poly f = 2x5+3y2+4z3; |
---|
645 | normalize(f); |
---|
646 | module m=[9xy,0,3z3],[4z,6y,2x]; |
---|
647 | normalize(m); |
---|
648 | ring s = 0,(x,y,z),(c,ls); |
---|
649 | module m=[9xy,0,3z3],[4z,6y,2x]; |
---|
650 | normalize(m); |
---|
651 | } |
---|
652 | /////////////////////////////////////////////////////////////////////////////// |
---|
653 | |
---|
654 | /////////////////////////////////////////////////////////////////////////////// |
---|
655 | // Input: <ideal>=<f1,f2,...,fm> and <polynomial> g |
---|
656 | // Question: Does g lie in the radical of <ideal>? |
---|
657 | // Solution: Compute a standard basis G for <f1,f2,...,fm,gz-1> where z is a |
---|
658 | // new variable. Then g is contained in the radical of <ideal> <=> |
---|
659 | // 1 is generator in G. |
---|
660 | /////////////////////////////////////////////////////////////////////////////// |
---|
661 | proc rad_con (poly g,ideal I) |
---|
662 | "USAGE: rad_con(g,I); g polynomial, I ideal |
---|
663 | RETURN: 1 (TRUE) (type int) if g is contained in the radical of I |
---|
664 | @* 0 (FALSE) (type int) otherwise |
---|
665 | EXAMPLE: example rad_con; shows an example |
---|
666 | " |
---|
667 | { def br=basering; |
---|
668 | int n=nvars(br); |
---|
669 | int dB=degBound; |
---|
670 | degBound=0; |
---|
671 | string mp=string(minpoly); |
---|
672 | execute("ring R=("+charstr(br)+"),(x(1..n),z),dp;"); |
---|
673 | execute("minpoly=number("+mp+");"); |
---|
674 | ideal irrel=x(1..n); |
---|
675 | map f=br,irrel; |
---|
676 | poly p=f(g); |
---|
677 | ideal J=f(I)+ideal(p*z-1); |
---|
678 | J=std(J); |
---|
679 | degBound=dB; |
---|
680 | if (J[1]==1) |
---|
681 | { return(1); |
---|
682 | } |
---|
683 | else |
---|
684 | { return(0); |
---|
685 | } |
---|
686 | } |
---|
687 | example |
---|
688 | { "EXAMPLE:"; echo=2; |
---|
689 | ring R=0,(x,y,z),dp; |
---|
690 | ideal I=x2+y2,z2; |
---|
691 | poly f=x4+y4; |
---|
692 | rad_con(f,I); |
---|
693 | ideal J=x2+y2,z2,x4+y4; |
---|
694 | poly g=z; |
---|
695 | rad_con(g,I); |
---|
696 | } |
---|
697 | /////////////////////////////////////////////////////////////////////////////// |
---|
698 | |
---|
699 | proc lcm (id, list #) |
---|
700 | "USAGE: lcm(p[,q]); p int/intvec q a list of integers or |
---|
701 | p poly/ideal q a list of polynomials |
---|
702 | RETURN: the least common multiple of the common entries of p and q: |
---|
703 | @* - of type int if p is an int/intvec |
---|
704 | @* - of type poly if p is a poly/ideal |
---|
705 | EXAMPLE: example lcm; shows an example |
---|
706 | " |
---|
707 | { |
---|
708 | int k,j; |
---|
709 | //------------------------------ integer case -------------------------------- |
---|
710 | if( typeof(id) == "int" or typeof(id) == "intvec" ) |
---|
711 | { |
---|
712 | intvec i = id; |
---|
713 | if (size(#)!=0) |
---|
714 | { |
---|
715 | for (j = 1; j<=size(#); j++) |
---|
716 | { |
---|
717 | if (typeof(#[j]) !="int" and typeof(#[j]) !="intvec") |
---|
718 | { ERROR("// ** list element must be an integer");} |
---|
719 | else |
---|
720 | { i = i,#[j]; } |
---|
721 | } |
---|
722 | } |
---|
723 | int p,q; |
---|
724 | if( i == 0 ) |
---|
725 | { |
---|
726 | return(0); |
---|
727 | } |
---|
728 | for(j=1;j<=size(i);j++) |
---|
729 | { |
---|
730 | if( i[j] != 0 ) |
---|
731 | { |
---|
732 | p=i[j]; |
---|
733 | break; |
---|
734 | } |
---|
735 | } |
---|
736 | for (k=j+1;k<=size(i);k++) |
---|
737 | { |
---|
738 | if( i[k] !=0) |
---|
739 | { |
---|
740 | q=gcd(p,i[k]); |
---|
741 | p=p/q; |
---|
742 | p=p*i[k]; |
---|
743 | } |
---|
744 | } |
---|
745 | if(p <0 ) |
---|
746 | {return(-p);} |
---|
747 | return(p); |
---|
748 | } |
---|
749 | |
---|
750 | //----------------------------- polynomial case ------------------------------ |
---|
751 | if( typeof(id) == "poly" or typeof(id) == "ideal" ) |
---|
752 | { |
---|
753 | ideal i = id; |
---|
754 | if (size(#)!=0) |
---|
755 | { |
---|
756 | for (j = 1; j<=size(#); j++) |
---|
757 | { |
---|
758 | if (typeof(#[j]) !="poly" and typeof(#[j]) !="ideal" |
---|
759 | and typeof(#[j]) !="int" and typeof(#[j]) !="intvec") |
---|
760 | { ERROR("// ** list element must be a polynomial");} |
---|
761 | else |
---|
762 | { i = i,#[j]; } |
---|
763 | } |
---|
764 | } |
---|
765 | poly p,q; |
---|
766 | i=simplify(i,10); |
---|
767 | if(size(i) == 0) |
---|
768 | { |
---|
769 | return(0); |
---|
770 | } |
---|
771 | for(j=1;j<=size(i);j++) |
---|
772 | { |
---|
773 | if(maxdeg(i[j])!= 0) |
---|
774 | { |
---|
775 | p=i[j]; |
---|
776 | break; |
---|
777 | } |
---|
778 | } |
---|
779 | if(deg(p)==-1) |
---|
780 | { |
---|
781 | return(1); |
---|
782 | } |
---|
783 | for (k=j+1;k<=size(i);k++) |
---|
784 | { |
---|
785 | if(maxdeg(i[k])!=0) |
---|
786 | { |
---|
787 | q=gcd(p,i[k]); |
---|
788 | if(maxdeg(q)==0) |
---|
789 | { |
---|
790 | p=p*i[k]; |
---|
791 | } |
---|
792 | else |
---|
793 | { |
---|
794 | p=p/q; |
---|
795 | p=p*i[k]; |
---|
796 | } |
---|
797 | } |
---|
798 | } |
---|
799 | return(p); |
---|
800 | } |
---|
801 | } |
---|
802 | example |
---|
803 | { "EXAMPLE:"; echo = 2; |
---|
804 | ring r = 0,(x,y,z),lp; |
---|
805 | poly p = (x+y)*(y+z); |
---|
806 | poly q = (z4+2)*(y+z); |
---|
807 | lcm(p,q); |
---|
808 | ideal i=p,q,y+z; |
---|
809 | lcm(i,p); |
---|
810 | lcm(2,3,6); |
---|
811 | lcm(2..6); |
---|
812 | } |
---|
813 | |
---|
814 | /////////////////////////////////////////////////////////////////////////////// |
---|
815 | |
---|
816 | proc content(f) |
---|
817 | "USAGE: content(f); f polynomial/vector |
---|
818 | RETURN: number, the content (greatest common factor of coefficients) |
---|
819 | of the polynomial/vector f |
---|
820 | EXAMPLE: example content; shows an example |
---|
821 | " |
---|
822 | { |
---|
823 | return(leadcoef(f)/leadcoef(cleardenom(f))); |
---|
824 | } |
---|
825 | example |
---|
826 | { "EXAMPLE:"; echo = 2; |
---|
827 | ring r=0,(x,y,z),(c,lp); |
---|
828 | content(3x2+18xy-27xyz); |
---|
829 | vector v=[3x2+18xy-27xyz,15x2+12y4,3]; |
---|
830 | content(v); |
---|
831 | } |
---|
832 | /////////////////////////////////////////////////////////////////////////////// |
---|
833 | |
---|
834 | proc numerator(number n) |
---|
835 | "USAGE: numerator(n); n number |
---|
836 | RETURN: number, the numerator of n |
---|
837 | SEE ALSO: denominator, content, cleardenom |
---|
838 | EXAMPLE: example numerator; shows an example |
---|
839 | " |
---|
840 | { |
---|
841 | poly p = cleardenom(n+var(1)); |
---|
842 | return (number(coeffs(p,var(1))[1,1])); |
---|
843 | } |
---|
844 | example |
---|
845 | { |
---|
846 | "EXAMPLE:"; echo = 2; |
---|
847 | ring r = 0,x, dp; |
---|
848 | number n = 3/2; |
---|
849 | numerator(n); |
---|
850 | } |
---|
851 | /////////////////////////////////////////////////////////////////////////////// |
---|
852 | |
---|
853 | proc denominator(number n) |
---|
854 | "USAGE: denominator(n); n number |
---|
855 | RETURN: number, the denominator of n |
---|
856 | SEE ALSO: denominator, content, cleardenom |
---|
857 | EXAMPLE: example denominator; shows an example |
---|
858 | " |
---|
859 | { |
---|
860 | poly p = cleardenom(n+var(1)); |
---|
861 | return (number(coeffs(p,var(1))[2,1])); |
---|
862 | } |
---|
863 | example |
---|
864 | { |
---|
865 | "EXAMPLE:"; echo = 2; |
---|
866 | ring r = 0,x, dp; |
---|
867 | number n = 3/2; |
---|
868 | denominator(n); |
---|
869 | } |
---|
870 | //////////////////////////////////////////////////////////////////////// |
---|
871 | |
---|
872 | //////////////////////////////////////////////////////////////////////// |
---|
873 | // The idea of the procedures mod2id, id2mod and subrInterred is, to |
---|
874 | // perform standard basis computation or interreduction of a submodule |
---|
875 | // of a free module with generators gen(1),...,gen(n) over a ring R |
---|
876 | // in a ring R[t1,...,tn]/<ti*tj> with gen(i) maped to ti |
---|
877 | //////////////////////////////////////////////////////////////////////// |
---|
878 | |
---|
879 | proc mod2id(matrix M,intvec vpos) |
---|
880 | "USAGE: mod2id(M,vpos); M matrix, vpos intvec |
---|
881 | ASSUME: vpos is an integer vector such that gen(i) corresponds |
---|
882 | to var(vpos[i]). |
---|
883 | The basering contains variables var(vpos[i]) which do not occur |
---|
884 | in M. |
---|
885 | RETURN: ideal I in which each gen(i) from the module is replaced by |
---|
886 | var(vpos[i]) and all monomials var(vpos[i])*var(vpos[j]) have |
---|
887 | been added to the generating set of I. |
---|
888 | NOTE: This procedure should be used in the following situation: |
---|
889 | one wants to pass to a ring with new variables, say e(1),..,e(s), |
---|
890 | which correspond to the components gen(1),..,gen(s) of the |
---|
891 | module M such that e(i)*e(j)=0 for all i,j. |
---|
892 | The new ring should already exist and be the current ring |
---|
893 | EXAMPLE: example mod2id; shows an example" |
---|
894 | { |
---|
895 | int p = printlevel-voice+3; // p=printlevel+1 (default: p=1) |
---|
896 | //---------------------------------------------------------------------- |
---|
897 | // define the ideal generated by the var(vpos[i]) and set up the matrix |
---|
898 | // for the multiplication |
---|
899 | //---------------------------------------------------------------------- |
---|
900 | ideal vars=var(vpos[1]); |
---|
901 | for(int i=2;i<=size(vpos);i++) |
---|
902 | { |
---|
903 | vars=vars,var(vpos[i]); |
---|
904 | } |
---|
905 | matrix varm[1][size(vpos)]=vars; |
---|
906 | if (size(vpos) > nrows(M)) |
---|
907 | { |
---|
908 | matrix Mt[size(vpos)][ncols(M)]; |
---|
909 | Mt[1..nrows(M),1..ncols(M)]=M; |
---|
910 | kill M; |
---|
911 | matrix M=Mt; |
---|
912 | } |
---|
913 | //---------------------------------------------------------------------- |
---|
914 | // define the desired ideal |
---|
915 | //---------------------------------------------------------------------- |
---|
916 | ideal ret=vars^2,varm*M; |
---|
917 | return(ret); |
---|
918 | } |
---|
919 | example |
---|
920 | { "EXAMPLE:"; echo=2; |
---|
921 | ring r=0,(e(1),e(2),x,y,z),(dp(2),ds(3)); |
---|
922 | module mo=x*gen(1)+y*gen(2); |
---|
923 | intvec iv=2,1; |
---|
924 | mod2id(mo,iv); |
---|
925 | } |
---|
926 | //////////////////////////////////////////////////////////////////////// |
---|
927 | |
---|
928 | proc id2mod(ideal i,intvec vpos) |
---|
929 | "USAGE: id2mod(I,vpos); I ideal, vpos intvec |
---|
930 | RETURN: module corresponding to the ideal by replacing var(vpos[i]) by |
---|
931 | gen(i) and omitting all generators var(vpos[i])*var(vpos[j]) |
---|
932 | NOTE: * This procedure only makes sense if the ideal contains |
---|
933 | all var(vpos[i])*var(vpos[j]) as monomial generators and |
---|
934 | all other generators of I are linear combinations of the |
---|
935 | var(vpos[i]) over the ring in the other variables. |
---|
936 | * This is the inverse procedure to mod2id and should be applied |
---|
937 | only to ideals created by mod2id using the same intvec vpos |
---|
938 | (possibly after a standard basis computation) |
---|
939 | EXAMPLE: example id2mod; shows an example" |
---|
940 | { |
---|
941 | int p = printlevel-voice+3; // p=printlevel+1 (default: p=1) |
---|
942 | //--------------------------------------------------------------------- |
---|
943 | // Initialization |
---|
944 | //--------------------------------------------------------------------- |
---|
945 | int n=size(i); |
---|
946 | int v=size(vpos); |
---|
947 | matrix tempmat; |
---|
948 | matrix mm[v][n]; |
---|
949 | //--------------------------------------------------------------------- |
---|
950 | // Conversion |
---|
951 | //--------------------------------------------------------------------- |
---|
952 | for(int j=1;j<=v;j++) |
---|
953 | { |
---|
954 | tempmat=coeffs(i,var(vpos[j])); |
---|
955 | mm[j,1..n]=tempmat[2,1..n]; |
---|
956 | } |
---|
957 | for(j=1;j<=v;j++) |
---|
958 | { |
---|
959 | mm=subst(mm,var(vpos[j]),0); |
---|
960 | } |
---|
961 | module ret=simplify(mm,10); |
---|
962 | return(ret); |
---|
963 | } |
---|
964 | example |
---|
965 | { "EXAMPLE:"; echo=2; |
---|
966 | ring r=0,(e(1),e(2),x,y,z),(dp(2),ds(3)); |
---|
967 | ideal i=e(2)^2,e(1)*e(2),e(1)^2,e(1)*y+e(2)*x; |
---|
968 | intvec iv=2,1; |
---|
969 | id2mod(i,iv); |
---|
970 | } |
---|
971 | /////////////////////////////////////////////////////////////////////// |
---|
972 | |
---|
973 | proc subrInterred(ideal mon, ideal sm, intvec iv) |
---|
974 | "USAGE: subrInterred(mon,sm,iv); |
---|
975 | sm: ideal in a ring r with n + s variables, |
---|
976 | e.g. x_1,..,x_n and t_1,..,t_s |
---|
977 | mon: ideal with monomial generators (not divisible by |
---|
978 | any of the t_i) such that sm is contained in the module |
---|
979 | k[t_1,..,t_s]*mon[1]+..+k[t_1,..,t_s]*mon[size(mon)] |
---|
980 | iv: intvec listing the variables which are supposed to be used |
---|
981 | as x_i |
---|
982 | RETURN: list l: |
---|
983 | l[1]=the monomials from mon in the order used |
---|
984 | l[2]=their coefficients after interreduction |
---|
985 | l[3]=l[1]*l[2] |
---|
986 | PURPOSE: Do interred only w.r.t. a subset of variables. |
---|
987 | The procedure returns an interreduced system of generators of |
---|
988 | sm considered as a k[t_1,..,t_s]-submodule of the free module |
---|
989 | k[t_1,..,t_s]*mon[1]+..+k[t_1,..,t_s]*mon[size(mon)]). |
---|
990 | EXAMPLE: example subrInterred; shows an example |
---|
991 | " |
---|
992 | { |
---|
993 | int p = printlevel-voice+3; // p=printlevel+1 (default: p=1) |
---|
994 | //----------------------------------------------------------------------- |
---|
995 | // check that mon is really generated by monomials |
---|
996 | // and sort its generators with respect to the monomial ordering |
---|
997 | //----------------------------------------------------------------------- |
---|
998 | int err; |
---|
999 | for(int i=1;i<=ncols(mon);i++) |
---|
1000 | { |
---|
1001 | if ( size(mon[i]) > 1 ) |
---|
1002 | { |
---|
1003 | err=1; |
---|
1004 | } |
---|
1005 | } |
---|
1006 | if (err==1) |
---|
1007 | { |
---|
1008 | ERROR("mon has to be generated by monomials"); |
---|
1009 | } |
---|
1010 | intvec sv=sortvec(mon); |
---|
1011 | ideal mons; |
---|
1012 | for(i=1;i<=size(sv);i++) |
---|
1013 | { |
---|
1014 | mons[i]=mon[sv[i]]; |
---|
1015 | } |
---|
1016 | ideal itemp=maxideal(1); |
---|
1017 | for(i=1;i<=size(iv);i++) |
---|
1018 | { |
---|
1019 | itemp=subst(itemp,var(iv[i]),0); |
---|
1020 | } |
---|
1021 | itemp=simplify(itemp,10); |
---|
1022 | def r=basering; |
---|
1023 | string tempstr="ring rtemp=" + charstr(basering) + ",(" + string(itemp) |
---|
1024 | + "),(C,dp);"; |
---|
1025 | //----------------------------------------------------------------------- |
---|
1026 | // express m in terms of the generators of mon and check whether m |
---|
1027 | // can be considered as a submodule of k[t_1,..,t_n]*mon |
---|
1028 | //----------------------------------------------------------------------- |
---|
1029 | module motemp; |
---|
1030 | motemp[ncols(sm)]=0; |
---|
1031 | poly ptemp; |
---|
1032 | int j; |
---|
1033 | for(i=1;i<=ncols(mons);i++) |
---|
1034 | { |
---|
1035 | for(j=1;j<=ncols(sm);j++) |
---|
1036 | { |
---|
1037 | ptemp=sm[j]/mons[i]; |
---|
1038 | motemp[j]=motemp[j]+ptemp*gen(i); |
---|
1039 | } |
---|
1040 | } |
---|
1041 | for(i=1;i<=size(iv);i++) |
---|
1042 | { |
---|
1043 | motemp=subst(motemp,var(iv[i]),0); |
---|
1044 | } |
---|
1045 | matrix monmat[1][ncols(mons)]=mons; |
---|
1046 | ideal dummy=monmat*motemp; |
---|
1047 | for(i=1;i<=size(sm);i++) |
---|
1048 | { |
---|
1049 | if(sm[i]-dummy[i]!=0) |
---|
1050 | { |
---|
1051 | ERROR("the second argument is not a submodule of the assumed structure"); |
---|
1052 | } |
---|
1053 | } |
---|
1054 | //---------------------------------------------------------------------- |
---|
1055 | // do the interreduction and convert back |
---|
1056 | //---------------------------------------------------------------------- |
---|
1057 | execute(tempstr); |
---|
1058 | def motemp=imap(r,motemp); |
---|
1059 | motemp=interred(motemp); |
---|
1060 | setring r; |
---|
1061 | kill motemp; |
---|
1062 | def motemp=imap(rtemp,motemp); |
---|
1063 | list ret=monmat,motemp,monmat*motemp; |
---|
1064 | for(i=1;i<=ncols(ret[2]);i++) |
---|
1065 | { |
---|
1066 | ret[2][i]=cleardenom(ret[2][i]); |
---|
1067 | } |
---|
1068 | for(i=1;i<=ncols(ret[3]);i++) |
---|
1069 | { |
---|
1070 | ret[3][i]=cleardenom(ret[3][i]); |
---|
1071 | } |
---|
1072 | return(ret); |
---|
1073 | } |
---|
1074 | example |
---|
1075 | { "EXAMPLE:"; echo=2; |
---|
1076 | ring r=0,(x,y,z),dp; |
---|
1077 | ideal i=x^2+x*y^2,x*y+x^2*y,z; |
---|
1078 | ideal j=x^2+x*y^2,x*y,z; |
---|
1079 | ideal mon=x^2,z,x*y; |
---|
1080 | intvec iv=1,3; |
---|
1081 | subrInterred(mon,i,iv); |
---|
1082 | subrInterred(mon,j,iv); |
---|
1083 | } |
---|
1084 | //////////////////////////////////////////////////////////////////////// |
---|
1085 | |
---|
1086 | |
---|
1087 | |
---|