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