Changeset e212bd in git
 Timestamp:
 Nov 20, 2013, 4:07:38 PM (9 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '48f1dd268d0ff74ef2f7dccbf02545425002ddcc')
 Children:
 fd1101c51f85808be5f4f270758360b46a425a4f
 Parents:
 867d59cf46048ef5ad1c520db05b0b6d5f71d305
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

kernel/hilb.cc
r867d59 re212bd 17 17 #include <kernel/hutil.h> 18 18 #include <kernel/stairc.h> 19 //ADICHANGES: 20 #include <kernel/ideals.h> 21 #include <kernel/kstd1.h> 22 #include<gmp.h> 23 #include<vector> 24 19 25 20 26 static int **Qpol; 21 27 static int *Q0, *Ql; 22 28 static int hLength; 29 //ADICHANGES 30 static int NNN; 31 static mpz_t ec; 32 static mpz_ptr hilbertcoef = (mpz_ptr)omAlloc(NNN*sizeof(mpz_t)); 33 std::vector<int> hilbertpowers; 34 std::vector<int> hilbertpowerssorted; 35 23 36 24 37 static int hMinModulweight(intvec *modulweight) … … 213 226 } 214 227 } 215 228 //  ADICHANGES  229 //!!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!! 230 void PutInVector(int degvalue, int oldposition, int where) 231 { 232 hilbertpowerssorted.resize(hilbertpowerssorted.size()+2); 233 int i; 234 for(i=hilbertpowerssorted.size()1;i>=where+2;i) 235 { 236 hilbertpowerssorted[i] = hilbertpowerssorted[i2]; 237 } 238 hilbertpowerssorted[where] = degvalue; 239 hilbertpowerssorted[where+1] = oldposition; 240 } 241 242 void SortPowerVec() 243 { 244 int i,j; 245 int test; 246 bool flag; 247 hilbertpowerssorted.resize(2); 248 hilbertpowerssorted[0] = hilbertpowers[0]; 249 hilbertpowerssorted[1] = 0; 250 for(i=1;i<hilbertpowers.size();i++) 251 { 252 flag=TRUE; 253 if(hilbertpowers[i]<=hilbertpowerssorted[0]) 254 { 255 flag=FALSE; 256 PutInVector(hilbertpowers[i], i, 0); 257 } 258 if(hilbertpowers[i]>=hilbertpowerssorted[hilbertpowerssorted.size()2]) 259 { 260 flag=FALSE; 261 PutInVector(hilbertpowers[i], i, hilbertpowerssorted.size()2); 262 } 263 if(flag==TRUE) 264 { 265 for(j=2;(j<=hilbertpowerssorted.size()4)&&(flag);j=j+2) 266 { 267 if(hilbertpowers[i]<=hilbertpowerssorted[j]) 268 { 269 flag=FALSE; 270 PutInVector(hilbertpowers[i], i, j); 271 } 272 } 273 } 274 //printf("\nTEST\n"); 275 //for(test=0;test<hilbertpowerssorted.size();test++) 276 // { 277 // printf(" %d ", hilbertpowerssorted[test]); 278 // } 279 280 } 281 } 282 283 int DegMon(poly p, ring tailRing) 284 { 285 int i,deg; 286 deg = 0; 287 for(i=0;i<=tailRing>N;i++) 288 { 289 deg = deg + p_GetExp(p, i, tailRing); 290 } 291 return(deg); 292 } 293 294 poly SearchP(ideal I, ring tailRing) 295 { 296 int i,j,exp; 297 poly res; 298 for(i=IDELEMS(I)1;i>=0;i) 299 { 300 if(DegMon(I>m[i], tailRing)>=2) 301 { 302 res = p_ISet(1, tailRing); 303 res = p_Copy(I>m[i], tailRing); 304 for(j=0;j<=tailRing>N;j++) 305 { 306 exp = p_GetExp(I>m[i], j, tailRing); 307 if(exp > 0) 308 { 309 p_SetExp(res, j, exp  1, tailRing); 310 return(res); 311 } 312 } 313 } 314 } 315 return(NULL); 316 } 317 318 poly ChoosePVar (ideal I, ring tailRing) 319 { 320 bool flag=TRUE; 321 int i,j; 322 poly res; 323 for(i=1;i<=tailRing>N;i++) 324 { 325 flag=TRUE; 326 for(j=IDELEMS(I)1;(j>=0)&&(flag);j) 327 { 328 if(p_GetExp(I>m[j], i, tailRing)>0) 329 { 330 flag=FALSE; 331 } 332 } 333 334 if(flag == TRUE) 335 { 336 res = p_ISet(1, tailRing); 337 p_SetExp(res, i, 1, tailRing); 338 //idPrint(I); 339 //pWrite(res); 340 return(res); 341 } 342 } 343 return(NULL); //i.e. it is the maximal ideal 344 } 345 346 ideal idSimplify(ideal I, ring tailRing) 347 { 348 int i,j; 349 bool flag; 350 /*std::vector<int> var; 351 for(i=0;i<=tailRing>N;i++) 352 { 353 if(p_GetExp(I>[IDELEMS(I)1], tailRing)>0) 354 { 355 var.resize(var.size()+1); 356 var[var.size()1] = i; 357 } 358 } 359 */ 360 for(i=IDELEMS(I)2;i>=0;i) 361 { 362 flag=TRUE; 363 for(j=0;(j<=tailRing>N)&&(flag);j++) 364 { 365 if(p_GetExp(I>m[i], j, tailRing) < p_GetExp(I>m[IDELEMS(I)1], j, tailRing)) 366 { 367 flag=FALSE; 368 } 369 } 370 if(flag) 371 { 372 I>m[i]=NULL; 373 } 374 } 375 idSkipZeroes(I); 376 return(I); 377 } 378 379 380 void eulerchar (ideal I, ring tailRing) 381 { 382 //gmp_printf("\nEuler char: %Zd\n", ec); 383 mpz_t dummy; 384 if( idElem(I) == 0) 385 { 386 mpz_init(dummy); 387 //change: era 1 in loc de 0 388 mpz_set_si(dummy, 1); 389 mpz_sub(ec, ec, dummy); 390 return; 391 } 392 if( idElem(I) == 1) 393 { 394 if(!p_IsOne(I>m[0], tailRing)) 395 { 396 //aici nu era nimic 397 //mpz_set_si(dummy, 1); 398 //mpz_add(ec, ec, dummy); 399 return; 400 } 401 else 402 { 403 mpz_init(dummy); 404 mpz_set_si(dummy, 1); 405 mpz_sub(ec, ec, dummy); 406 return; 407 } 408 } 409 ideal p = idInit(1,1); 410 p>m[0] = SearchP(I, tailRing); 411 //idPrint(I); 412 //printf("\nSearchP founded: \n");pWrite(p>m[0]); 413 if(p>m[0] == NULL) 414 { 415 mpz_init(dummy); 416 mpz_set_si(dummy, IDELEMS(I)1); 417 mpz_add(ec, ec, dummy); 418 return; 419 } 420 ideal Ip = idQuot(I,p,TRUE,TRUE); 421 ideal Iplusp = id_Add(I,p,tailRing); 422 Iplusp=idSimplify(Iplusp, tailRing); 423 //mpz_t i,j; 424 //i = eulerchar(Ip, ec, tailRing); 425 //j = eulerchar(Iplusp, ec, tailRing); 426 //mpz_add(ec, i,j); 427 eulerchar(Ip, tailRing); 428 eulerchar(Iplusp, tailRing); 429 return; 430 } 431 432 poly SqFree (ideal I, ring tailRing) 433 { 434 int i,j; 435 bool flag=TRUE; 436 poly notsqrfree = NULL; 437 for(i=IDELEMS(I)1;(i>=0)&&(flag);i) 438 { 439 for(j=0;(j<=tailRing>N)&&(flag);j++) 440 { 441 if(p_GetExp(I>m[i],j,tailRing)>1) 442 { 443 flag=FALSE; 444 notsqrfree = p_ISet(1,tailRing); 445 p_SetExp(notsqrfree,j,1,tailRing); 446 } 447 } 448 } 449 return(notsqrfree); 450 } 451 452 void rouneslice(ideal I, ideal S, poly q, ring tailRing, poly x) 453 { 454 int i,j; 455 int dummy; 456 mpz_t dummympz; 457 poly m; 458 ideal p, koszsimp; 459 I = idMinBase(I); 460 461 //Work on it 462 // PRUNING OF S  463 S = idMinBase(S); 464 for(i=IDELEMS(S)1;i>=0;i) 465 { 466 if(kNF(I,NULL,S>m[i],NULL,NULL)==NULL) 467 { 468 S>m[i]=NULL; 469 } 470 } 471 idSkipZeroes(S); 472 // 473 for(i=IDELEMS(I)1;i>=0;i) 474 { 475 m = p_Copy(I>m[i],tailRing); 476 for(j=0;j<=tailRing>N;j++) 477 { 478 dummy = p_GetExp(m,j,tailRing); 479 if(dummy > 0) 480 { 481 p_SetExp(m,j,dummy1,tailRing); 482 } 483 } 484 if(kNF(S,NULL,m,NULL,NULL)==NULL) 485 { 486 I>m[i]=NULL; 487 //printf("\n Check if m is in S: \n");pWrite(m);idPrint(S); 488 //printf("\n Deleted, since pi(m) is in S\n");pWrite(m); 489 } 490 } 491 idSkipZeroes(I); 492 493 // MORE PRUNING OF S  494 //strictly divide??? 495 // 496 497 //printf("Ideal I: \n");idPrint(I); 498 //printf("Ideal S: \n");idPrint(S); 499 //printf("Poly q: \n");pWrite(q); 500 if(idElem(I)==0) 501 { 502 //I = 0 503 //printf("\nx does not divide lcm(I)"); 504 //pWrite(x);pWrite(m);idPrint(I); 505 printf("\nEmpty Set: ");pWrite(q); 506 return; 507 } 508 m = p_ISet(1,tailRing); 509 for(i=0;i<=tailRing>N;i++) 510 { 511 dummy=0; 512 for(j=IDELEMS(I)1;j>=0;j) 513 { 514 if(p_GetExp(I>m[j],i,tailRing) > dummy) 515 { 516 dummy = p_GetExp(I>m[j],i,tailRing); 517 } 518 } 519 p_SetExp(m,i,dummy,tailRing); 520 } 521 p_Setm(m,tailRing); 522 if(p_DivisibleBy(x,m,tailRing)==FALSE) 523 { 524 //printf("\nx does not divide lcm(I)"); 525 //pWrite(x);pWrite(m);idPrint(I); 526 //printf("\nEmpty set");pWrite(q); 527 return; 528 } 529 m = SqFree(I, tailRing); 530 koszsimp = idInit(IDELEMS(I),1); 531 if(m==NULL) 532 { 533 printf("\n Corner: "); 534 pWrite(q); 535 printf("\n With the facets of the simplex:\n"); 536 for(i=IDELEMS(I)1;i>=0;i) 537 { 538 m = p_ISet(1,tailRing); 539 //m = p_Divide(x,I>m[i],tailRing); 540 for(j=tailRing>N; j>=0; j) 541 { 542 p_SetExp(m,j, p_GetExp(x,j,tailRing) p_GetExp(I>m[i],j,tailRing),tailRing); 543 } 544 printf(" "); 545 p_Setm(m, tailRing); 546 p_Write(m, tailRing); 547 koszsimp>m[i] = p_Copy(m, tailRing); 548 } 549 printf("\n Euler Characteristic = "); 550 mpz_init(dummympz); 551 eulerchar(koszsimp, tailRing); 552 //gmp_printf("dummy %Zd \n", dummympz); 553 gmp_printf("ec %Zd \n", ec); 554 mpz_set(dummympz, ec); 555 mpz_sub(ec, ec, ec); 556 hilbertcoef = (mpz_ptr)omRealloc(hilbertcoef, (NNN+1)*sizeof(mpz_t)); 557 mpz_init(&hilbertcoef[NNN]); 558 mpz_set(&hilbertcoef[NNN], dummympz); 559 hilbertpowers.resize(NNN+1); 560 hilbertpowers[NNN] = DegMon(q, tailRing); 561 NNN++; 562 return; 563 } 564 if(IDELEMS(S)!=1) 565 { 566 m = SearchP(S, tailRing); 567 if(m == NULL) 568 { 569 m = ChoosePVar(S, tailRing); 570 } 571 } 572 p = idInit(1,1); 573 p>m[0] = m; 574 printf("Poly p: \n");pWrite(m); 575 ideal Ip = idQuot(I,p,TRUE,TRUE); 576 ideal Sp = idQuot(S,p,TRUE,TRUE); 577 ideal Splusp = id_Add(S,p,tailRing); 578 Splusp = idSimplify(Splusp, tailRing); 579 poly pq = pp_Mult_mm(q,m,tailRing); 580 rouneslice(Ip, Sp, pq, tailRing, x); 581 rouneslice(I, Splusp, q, tailRing, x); 582 return; 583 } 584 585 void slicehilb(ideal I, ring tailRing) 586 { 587 printf("Adi changes are here: \n"); 588 int i; 589 ideal S = idInit(0,1); 590 poly q = p_ISet(1,tailRing); 591 ideal X = idInit(1,1); 592 X>m[0]=p_One(tailRing); 593 for(i=1;i<=tailRing>N;i++) 594 { 595 p_SetExp(X>m[0],i,1,tailRing); 596 } 597 p_Setm(X>m[0],tailRing); 598 I = id_Mult(I,X,tailRing); 599 printf("\nRouneSlice\n"); 600 rouneslice(I,S,q,tailRing,X>m[0]); 601 for(i=0;i<NNN;i++) 602 { 603 gmp_printf(" %Zd ", &hilbertcoef[i]); 604 } 605 printf("\n"); 606 for(i=0;i<NNN;i++) 607 { 608 printf(" %d ", hilbertpowers[i]); 609 } 610 SortPowerVec(); 611 printf("\n"); 612 for(i=0;i<hilbertpowerssorted.size();i++) 613 { 614 printf(" %d ", hilbertpowerssorted[i]); 615 } 616 mpz_t coefhilb; 617 mpz_t dummy; 618 mpz_init(coefhilb); 619 mpz_init(dummy); 620 printf("\n// 1 t^0"); 621 for(i=0;i<hilbertpowerssorted.size();i=i+2) 622 { 623 mpz_set(dummy, &hilbertcoef[hilbertpowerssorted[i+1]]); 624 mpz_add(coefhilb, coefhilb, dummy); 625 while((hilbertpowerssorted[i]==hilbertpowerssorted[i+2]) && (i<=hilbertpowerssorted.size()2)) 626 { 627 i=i+2; 628 mpz_add(coefhilb, coefhilb, &hilbertcoef[hilbertpowerssorted[i+1]]); 629 } 630 if(mpz_sgn(coefhilb)!=0) 631 { 632 gmp_printf("\n// %Zd t^%d",coefhilb,hilbertpowerssorted[i]); 633 } 634 mpz_sub(coefhilb, coefhilb, coefhilb); 635 } 636 omFreeSize(hilbertcoef, NNN*sizeof(mpz_t)); 637 printf("\n\n"); 638 } 639 640 //  END OF CHANGES  216 641 static intvec * hSeries(ideal S, intvec *modulweight, 217 642 int /*notstc*/, intvec *wdegree, ideal Q, ring tailRing) 218 643 { 644 slicehilb(S,tailRing); // ADICHANGES 219 645 intvec *work, *hseries1=NULL; 220 646 int mc; … … 472 898 } 473 899 900 901
Note: See TracChangeset
for help on using the changeset viewer.