Changeset b5cc09 in git
- Timestamp:
- Oct 5, 2015, 7:26:34 PM (9 years ago)
- Branches:
- (u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b4f17ed1d25f93d46dbe29e4b499baecc2fd51bb')
- Children:
- 7544b6110f6c5a9a1708916223da2c76ca5f8b4d
- Parents:
- 4ce1b50655a7e6cda593952e66da5382178e6f14
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/GND.lib
r4ce1b5 rb5cc09 2 2 version="version desingularization.lib 4.0.0.0 Mai_2015"; 3 3 category="Commutative Algebra"; 4 4 5 info=" 5 6 LIBRARY : GND.lib 6 7 AUTHOR : Adrian Popescu, popescu@mathematik.uni-kl.de 7 OVERVIEW : 8 9 OVERVIEW : 8 10 A method to compute the General Neron Desingularization in the frame of 9 11 one dimensional local domains 10 12 11 SEE ALSO : A. Popescu and D. Popescu, "A method to compute the General Neron Desingularization in the frame of one dimensional local domains", http://arxiv.org/abs/1508.05511 13 REFERENCES: 14 @* [1] A. Popescu, D. Popescu, \"A method to compute the General Neron Desingularization in the frame of one dimensional local domains\", arxiv.org/abs/1508.05511 15 16 KEYWORDS: desingularization; general neron. 12 17 13 18 PROCEDURES: … … 109 114 /////////////////////////////////////////////////////////////////////////////////////////// 110 115 111 static proc finde(ideal i,poly d)112 "USAGE: Finds smallest e>0 s.t. (0):(d^e) will become stationary"113 {114 ring BASE = basering;115 qring r = i;116 int e = 2;117 poly d = imap(BASE,d);118 ideal i1 = std(quotient(ideal(0),ideal(d^(e-1))));119 ideal i2 = std(quotient(ideal(0),ideal(d^(e))));120 while(reduce(i2,i1)!= 0)121 {122 e = e+1;123 i1 = i2;124 i2 = std(quotient(ideal(0),ideal(d^(e))));125 }126 setring(BASE);127 return(e-1);128 }129 130 ///////////////////////////////////////////////////////////////////////////////////////////131 116 132 117 static proc extractP(list wP, ideal I) … … 294 279 295 280 static proc inY(poly p,int nrx,int nry) 296 "USAGE: tests if poly is just in the Y variables"281 "USAGE: tests if poly p contains just the variables var(nrx+1),...,var(nry)" 297 282 { 298 283 int i; … … 313 298 314 299 static proc myjet(ideal p, int bound,string param,string variab) 315 "USAGE: computes the jet in the right ring (instead of switching on and on)"300 "USAGE: computes the jet in the right ring (instead of switching it on and on)" 316 301 { 317 302 def r = basering; 318 if(param != "")319 {320 execute("ring R2 = (0,"+param+"),("+variab+"),ds;");321 }322 else323 {324 execute("ring R2 = 0,("+variab+"),ds;");325 }303 list l; 304 l[1] = list(0,param,list(list("lp",1)),ideal(0)); 305 l[2] = variab; 306 l[3] = ringlist(r)[3]; 307 l[3][1][1] = "ds"; 308 l[4] = ideal(0); 309 def R2 = ring(l); 310 setring(R2); 326 311 ideal p = imap(r,p); 327 312 for(int i=1;i<=ncols(p);i++) … … 335 320 /////////////////////////////////////////////////////////////////////////////////////////// 336 321 337 proc invp(poly p, int bound, string param,stringvariab)322 proc invp(poly p, int bound,list param,list variab) 338 323 "USAGE: computes the inverse in Q(param)[variab]" 339 324 { 340 325 def r = basering; 341 if(param != "")342 {343 execute("ring R2 = (0,"+param+"),("+variab+"),ds;");344 }345 else346 {347 execute("ring R2 = 0,("+variab+"),ds;");348 }326 list l; 327 l[1] = list(0,param,list(list("lp",1)),ideal(0)); 328 l[2] = variab; 329 l[3] = ringlist(r)[3]; 330 l[3][1][1] = "ds"; 331 l[4] = ideal(0); 332 def R2 = ring(l); 333 setring(R2); 349 334 poly p = imap(r,p); 350 335 if(bound<=0) … … 356 341 ERROR("My inverse : p is 0"); 357 342 } 358 poly original;343 poly original,res; 359 344 original = p; 360 if( leadcoef(p) == 0)345 if(ord(p) != 0) 361 346 { 362 347 ERROR("My inverse : the power series is not a unit."); 363 348 } 364 349 poly q = 1/p[1]; 365 polyres = q;350 res = q; 366 351 p = q * (p[1] - jet(p,bound)); 367 352 poly s = p; … … 371 356 p = jet(p*s,bound); 372 357 } 358 //res = division(poly(1),p,bound); 359 373 360 //TEST 374 361 if(jet(original*res,bound) != poly(1)) … … 414 401 if(total2[i][n] <= n) 415 402 { 416 execute("resf[size(resf)+1] = intvec("+string(total2[i])+");"); 403 def dummy = total2[i]; 404 resf[size(resf)+1] = intvec(dummy[1..size(total2[i])]); 405 kill dummy; 417 406 } 418 407 if(total2[i][n] <= y) 419 408 { 420 execute("resy[size(resy)+1] = intvec("+string(total2[i])+");"); 409 def dummy = total2[i]; 410 resy[size(resy)+1] = intvec(dummy[1..size(total2[i])]); 411 kill dummy; 421 412 } 422 413 i = i+1; … … 439 430 /////////////////////////////////////////////////////////////////////////////////////////// 440 431 441 static proc mymin(a,b) 442 "USAGE: finds minimium between a and b" 443 { 444 if(a<b) 445 { 446 return(a); 447 } 448 else 449 { 450 return(b); 451 } 452 } 453 454 /////////////////////////////////////////////////////////////////////////////////////////// 455 456 static proc findnewvar(string Z) 432 433 proc findnewvar(string Z) 457 434 "USAGE: Finds new unused variable Z,Z1,Z2,..." 458 435 { 459 int i,j; 460 int exists = 0; 461 for(i=1;i<=nvars(basering) && exists==0;i++) 462 { 463 if(string(var(i)) == Z) 464 { 465 exists=1; 466 } 467 } 468 if(exists==0) 469 { 470 return(Z); 471 } 472 string newvar; 473 while(1) 474 { 475 j = j+1; 476 newvar = Z+"("+string(j)+")"; 477 exists=0; 478 for(i=1;i<=nvars(basering) && exists==0;i++) 479 { 480 if(string(var(i)) == newvar) 481 { 482 exists=1; 483 } 484 } 485 if(exists==0) 486 { 487 return(newvar); 488 } 489 } 490 } 491 492 /////////////////////////////////////////////////////////////////////////////////////////// 493 494 static proc mylcm(ideal id,string as) 436 string S = "," + charstr(basering) + "," + varstr(basering) + ","; 437 Z = "," + Z + ","; 438 if(find(S,Z) == 0) 439 { 440 return(Z[2..size(Z)-1]); 441 } 442 int i=1; 443 while(find(S,Z+"("+string(i)+")") == 1) 444 { 445 i++; 446 } 447 Z = string(Z[2..size(Z)-1])+"("+string(i)+")"; 448 return(Z); 449 } 450 451 /////////////////////////////////////////////////////////////////////////////////////////// 452 453 static proc mylcm(ideal id,list asl) 495 454 "USAGE: computes the lcm of these numbers" 496 455 { 497 456 int i; 498 457 def r1 = basering; 499 execute("ring r2 = 0,("+as+"),dp;"); 458 list l; 459 l[1] = 0; 460 l[2] = asl; 461 l[3] = list(list("dp",1),list("C",0)); 462 l[4] = ideal(0); 463 def r2 = ring(l); 464 setring(r2); 500 465 ideal id = imap(r1,id); 501 466 poly lcmi=1; … … 559 524 /////////////////////////////////////////////////////////////////////////////////////////// 560 525 526 static proc crl(list param, list variab, string order) 527 "USAGE: construct the desired ringlist" 528 { 529 list l; 530 if(size(param)==0) 531 { 532 l[1] = 0; 533 } 534 else 535 { 536 l[1] = list(0,param,list(list("lp",1)),ideal(0)); 537 } 538 l[2] = variab; 539 l[3] = list(list(order,1),list("C",0)); 540 l[4] = ideal(0); 541 return(l); 542 } 543 544 /////////////////////////////////////////////////////////////////////////////////////////// 545 561 546 proc desingularization(all, int nra, int nrx, int nry, ideal xid, ideal yid, ideal aid, ideal f, 562 547 list #) … … 586 571 } 587 572 string as,xs,ys; 573 list asl,xsl,ysl; 588 574 ideal a,x,y; 589 575 int i,j,k; … … 591 577 { 592 578 as = string(var(1)); 579 asl[1] = string(var(1)); 593 580 a[1] = var(1); 594 581 for(i=2;i<=nra;i++) … … 596 583 as = as+","+string(var(i)); 597 584 a[ncols(a)+1] = var(i); 585 asl[size(asl)+1] = string(var(i)); 598 586 } 599 587 } … … 603 591 xs = string(var(nra+1)); 604 592 x[1] = var(nra+1); 593 xsl[1] = string(var(nra+1)); 605 594 for(i=nra+2;i<=nra+nrx;i++) 606 595 { 607 596 xs = xs+","+string(var(i)); 608 597 x[ncols(x)+1] = var(i); 598 xsl[size(xsl)+1] = string(var(i)); 609 599 } 610 600 } … … 614 604 ys = string(var(nra+nrx+1)); 615 605 y[1] = var(nra+nrx+1); 606 ysl[1] = string(var(nra+nrx+1)); 616 607 for(i=nra+nrx+2;i<=nra+nrx+nry;i++) 617 608 { 618 609 ys = ys+","+string(var(i)); 619 610 y[ncols(y)+1] = var(i); 611 ysl[size(ysl)+1] = string(var(i)); 620 612 } 621 613 } 622 614 //ys; 623 execute("ring Abase = 0,("+xs+"),dp;"); 615 def Abase = ring(crl(list(),xsl,"dp")); 616 setring Abase; 624 617 qring A = std(imap(all,xid)); 625 execute("ring Bbase = (0,"+xs+"),("+ys+"),dp;"); 618 def Bbase = ring(crl(xsl,ysl,"dp")); 619 setring Bbase; 626 620 ideal yid = imap(all, yid); 627 621 i = ncols(yid); … … 637 631 } 638 632 qring B = yid; 639 execute("ring Bfalsebase = 0,("+xs+","+ys+"),dp;"); 633 def Bfalsebase = ring(crl(list(),xsl+ysl,"dp")); 634 setring(Bfalsebase); 640 635 qring Bfalse = std(imap(all,yid)); 641 636 if(nra!=0) 642 637 { 643 execute("ring Apreal = (0,"+as+"),("+xs+"),dp;"); 644 execute("ring Apbase = 0,("+as+","+xs+"),dp;"); 638 def Apreal = ring(crl(asl,xsl,"dp")); 639 def Apbase = ring(crl(list(),asl+xsl,"dp")); 640 setring Apbase; 645 641 qring Ap = std(imap(all,xid)+imap(all,aid)); 646 642 } 647 643 else 648 644 { 649 execute("ring Apbase = 0,("+xs+"),dp;"); 645 def Apbase = ring(crl(list(),xsl,"dp")); 646 setring Apbase; 650 647 qring Ap = std(imap(all,xid)); 651 648 } … … 734 731 // It should be enough to consider the minor I have found. 735 732 int found = 0; 736 //ideal mino = minor(jacob(yid),m ymin(nrows(jacob(yid)),ncols(jacob(yid))));733 //ideal mino = minor(jacob(yid),min(nrows(jacob(yid)),ncols(jacob(yid)))); 737 734 //mino;~; 738 735 //for(k = 1;size(mino);i++) … … 749 746 dz[2] = v(M) /dz[1]; 750 747 //z is actually 1/z so we have to construct it's invers 751 dz[2] = 1/leadcoef(dz[2])*invp(dz[2]/leadcoef(dz[2]),3*i,as ,xs);748 dz[2] = 1/leadcoef(dz[2])*invp(dz[2]/leadcoef(dz[2]),3*i,asl,xsl); 752 749 found = 1; 753 750 } … … 781 778 denz[size(denz)+1] = denominator(leadcoef(z[i])); 782 779 } 783 def den = mylcm(denz,as );780 def den = mylcm(denz,asl); 784 781 z = z * den; 785 782 int newAptest = 0; … … 788 785 newAptest = 1; 789 786 string aa = "a"; 790 def newpar = findnewvar(aa); 791 as = as + "," + newpar; 787 list newpar; 788 newpar[1] = string(findnewvar(aa)); 789 as = as + "," + findnewvar(aa); 792 790 nra = nra + 1; 793 execute("ring newApbase = 0,("+as+","+xs+"),dp;"); 791 def newApbase = ring(crl(list(),asl+xsl,"dp")); 792 setring newApbase; 794 793 poly p = var(nra)*imap(Apreal,den)-1; 795 794 qring newAp = std(imap(all,xid)+imap(all,aid)+p); … … 803 802 Apbase = newApbase; 804 803 Ap = newAp; 805 execute("ring newApreal = (0,"+as+"),("+xs+"),dp;"); 804 def newApreal = ring(crl(asl,xsl,"dp")); 805 setring newApreal; 806 806 def z = imap(Ap,z); 807 807 def d = imap(Apreal,d); 808 808 def den = imap(Apreal,den); 809 809 Apreal = newApreal; 810 execute("ring newaring = 0,"+newpar+",dp;"); 810 def newaring = ring(crl(list(),newpar,"dp")); 811 setring newaring; 811 812 } 812 813 else … … 821 822 822 823 string Z = "Z"; 823 def newvar = findnewvar(Z); 824 execute("ring newvarring = 0,"+newvar+",dp;"); 824 list newvar; 825 newvar[1] = string(findnewvar(Z)); 826 def newvarring = ring(crl(list(),newvar,"dp")); 825 827 def B1base = Bbase+newvarring; 826 828 def B1 = B+newvarring; … … 836 838 setring B1false; 837 839 //If P is constant (i.e. no Y), than I can work with d = d' = P = P' 838 execute("poly Pconst = reduce(imap(Bfalse,Pprim),ideal("+ys+"));"); 840 poly Pconst = reduce(imap(Bfalse,Pprim),imap(all,y)); 841 839 842 if(Pconst == imap(Bfalse,Pprim)) 840 843 { … … 907 910 setring Ap; 908 911 def lam = imap(Apreal,lam); 909 stringnewas;912 list newas; 910 913 int howmanyfinala; 911 914 if(size(lam) == 0) 912 915 { 913 newas = "";914 execute("ring Cbase = 0,("+xs+"),dp;");916 def Cbase = ring(crl(list(),xsl,"dp")); 917 setring Cbase; 915 918 } 916 919 else … … 918 921 ideal varlam = variables(lam); 919 922 howmanyfinala = ncols(varlam); 920 newas = string(varlam[1]);921 for(i=2;i<=ncols(varlam);i++)922 {923 newas = newas + "," + string(varlam[i]);924 }925 execute("ring Cbase = 0,("+newas+","+xs+"),dp;");923 for(i=1;i<=ncols(varlam);i++) 924 { 925 newas[size(newas)+1] = string(varlam[i]); 926 } 927 def Cbase = ring(crl(list(),newas+xsl,"dp")); 928 setring Cbase; 926 929 } 927 930 setring Ap; … … 983 986 matrix adj = adjoint(PM); 984 987 matrix G = imap(Bfalse,Plist)[2]*var(nvars(B1))^2*adj; 985 execute("ring seeGring = 0,"+string(ys[1])+string(nvars(B1))+",dp;");988 def seeGring = ring(crl(list(),list(string(ys[1])+string(nvars(B1))),"dp")); 986 989 def seeG = B1 + seeGring; 987 990 setring seeG; … … 1032 1035 // Start constructing E 1033 1036 1034 stringts,yz;1035 yz = ys + "," + string(ys[1]) + string(nry+1);1037 list ts,yz; 1038 yz = ysl + list(string(ys[1]) + string(nry+1)); 1036 1039 string t = "T"; 1037 1040 def newvart = findnewvar(t); 1038 1041 for(i=1;i<=nry+1;i++) 1039 1042 { 1040 ts = ts + newvart + string(i) + ",";1041 } 1042 ts = ts[1..size(ts)-1];1043 execute("ring Efalsebase = 0,("+varstr(D)+","+yz+","+ts+"),dp;");1043 ts[size(ts)+1] = newvart + string(i); 1044 } 1045 def Efalsebase = ring(crl(list(),newas+xsl+yz+ts,"dp")); 1046 setring Efalsebase; 1044 1047 qring Efalse = std(imap(D,dquot)); 1045 1048 def f = imap(Bfalse,Plist)[4]; … … 1067 1070 "This is cc";cc; 1068 1071 } 1069 execute("ring E = (0,"+varstr(D)+"),("+yz+","+ts+"),dp;"); 1072 def E = ring(crl(newas+xsl,yz+ts,"dp")); 1073 setring E; 1070 1074 def d = imap(D,d); 1071 1075 number ss = leadcoef(imap(D,ss)); … … 1075 1079 map mapvidjet = E,vidjet[nrx+1..ncols(vidjet)]; 1076 1080 G = mapvidjet(G); 1077 execute("matrix T[ "+string(nry+1)+"][1] = "+ts+";");1078 execute("ideal tid = "+ ts+";");1081 execute("matrix T[nry+1][1] = "+string(ts)+";"); 1082 execute("ideal tid = "+string(ts)+";"); 1079 1083 matrix ymy[nry+1][1]; 1080 1084 for(i=1;i<=nry+1;i++) … … 1303 1307 "And this is the map:"; 1304 1308 } 1305 execute("ideal themap = vidjet[nrx+1..size(vidjet)],"+t+",invssp;");1309 //execute("ideal themap = vidjet[nrx+1..size(vidjet)],"+t+",invssp;"); 1306 1310 invssp = 1/invssp; 1307 1311 if(debug) … … 1311 1315 string zz = "Z"; 1312 1316 def newvar2 = findnewvar(zz); 1313 execute("ring Efinalbase = (0,"+parstr(E)+"),("+varstr(E)+","+newvar2+"),dp;");1317 //execute("ring Efinalbase = (0,"+parstr(E)+"),("+varstr(E)+","+newvar2+"),dp;"); 1314 1318 def q = imap(E,e2); 1315 1319 poly teta = var(nvars(Efinalbase)) * imap(E,ss)*imap(E,ssp) - 1; … … 1360 1364 y4 = 1+a4*x2+(a2^2+a3)*x2^10+(a1+a2+a3+a4)*x1*x2^11; 1361 1365 1362 stringas,xs;1366 list as,xs; 1363 1367 if(nra != 0) 1364 1368 { 1365 as = string(var(1)); 1366 for( int i=2;i<=nra;i++) 1367 { 1368 as = as+","+string(var(i)); 1369 for(int i=1;i<=nra;i++) 1370 { 1371 as[size(as)+1] = string(var(i)); 1369 1372 } 1370 1373 } 1371 1374 if(nrx!=0) 1372 1375 { 1373 xs = string(var(nra+1)); 1374 for(int i=nra+2;i<=nra+nrx;i++) 1375 { 1376 xs = xs+","+string(var(i)); 1376 for(int i=nra+1;i<=nra+nrx;i++) 1377 { 1378 xs[size(xs)+1] = string(var(i)); 1377 1379 } 1378 1380 }
Note: See TracChangeset
for help on using the changeset viewer.