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