Changeset 204092 in git
- Timestamp:
- Jul 3, 2012, 7:51:13 PM (12 years ago)
- Branches:
- (u'spielwiese', '4a9821a93ffdc22a6696668bd4f6b8c9de3e6c5f')
- Children:
- 7088f18ba4c47a3e87c604a59adbacd8ef35e979
- Parents:
- cd5fefcfcb3762114c087ec176782c753f2a14d8
- git-author:
- Oleksandr Motsak <motsak@mathematik.uni-kl.de>2012-07-03 19:51:13+02:00
- git-committer:
- Oleksandr Motsak <motsak@mathematik.uni-kl.de>2014-05-07 04:41:45+02:00
- Location:
- dyn_modules/syzextra
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
dyn_modules/syzextra/mod_main.cc
rcd5fefc r204092 34 34 35 35 #include "singularxx_defs.h" 36 36 37 #include "DebugPrint.h" 37 38 #include "myNF.h" … … 51 52 #include <string.h> 52 53 53 54 #ifdef _GNU_SOURCE55 #define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp, r)56 #else57 #define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp)58 #endif59 54 60 55 … … 64 59 extern ring rAssure_InducedSchreyerOrdering(const ring r, BOOLEAN complete, int sign); 65 60 extern int rGetISPos(const int p, const ring r); 61 62 // USING_NAMESPACE_SINGULARXX; 66 63 67 64 USING_NAMESPACE( SINGULARXXNAME :: DEBUG ) … … 509 506 510 507 511 #ifdef _GNU_SOURCE 512 static int cmp_c_ds(const void *p1, const void *p2, void *R) 513 { 514 #else 515 static int cmp_c_ds(const void *p1, const void *p2) 516 { 517 void *R = currRing; 518 #endif 519 520 const int YES = 1; 521 const int NO = -1; 522 523 const ring r = (const ring) R; // TODO/NOTE: the structure is known: C, lp!!! 524 525 assume( r == currRing ); 526 527 const poly a = *(const poly*)p1; 528 const poly b = *(const poly*)p2; 529 530 assume( a != NULL ); 531 assume( b != NULL ); 532 533 assume( p_LmTest(a, r) ); 534 assume( p_LmTest(b, r) ); 535 536 537 const signed long iCompDiff = p_GetComp(a, r) - p_GetComp(b, r); 538 539 // TODO: test this!!!!!!!!!!!!!!!! 540 541 //return -( compare (c, qsorts) ) 542 543 const int __DEBUG__ = 0; 544 545 #ifndef NDEBUG 546 if( __DEBUG__ ) 547 { 548 PrintS("cmp_c_ds: a, b: \np1: "); dPrint(a, r, r, 2); 549 PrintS("b: "); dPrint(b, r, r, 2); 550 PrintLn(); 551 } 552 #endif 553 554 555 if( iCompDiff > 0 ) 556 return YES; 557 558 if( iCompDiff < 0 ) 559 return NO; 560 561 assume( iCompDiff == 0 ); 562 563 const signed long iDegDiff = p_Totaldegree(a, r) - p_Totaldegree(b, r); 564 565 if( iDegDiff > 0 ) 566 return YES; 567 568 if( iDegDiff < 0 ) 569 return NO; 570 571 assume( iDegDiff == 0 ); 572 573 #ifndef NDEBUG 574 if( __DEBUG__ ) 575 { 576 PrintS("cmp_c_ds: a & b have the same comp & deg! "); PrintLn(); 577 } 578 #endif 579 580 for (int v = rVar(r); v > 0; v--) 581 { 582 assume( v > 0 ); 583 assume( v <= rVar(r) ); 584 585 const signed int d = p_GetExp(a, v, r) - p_GetExp(b, v, r); 586 587 if( d > 0 ) 588 return YES; 589 590 if( d < 0 ) 591 return NO; 592 593 assume( d == 0 ); 594 } 595 596 return 0; 597 } 598 599 600 601 static ideal ComputeLeadingSyzygyTerms(const ideal& id, const ring r) 602 { 603 // 1. set of components S? 604 // 2. for each component c from S: set of indices of leading terms 605 // with this component? 606 // 3. short exp. vectors for each leading term? 607 608 const int size = IDELEMS(id); 609 610 if( size < 2 ) 611 { 612 const ideal newid = idInit(1, 0); newid->m[0] = NULL; // zero ideal... 613 return newid; 614 } 615 616 617 // TODO/NOTE: input is supposed to be (reverse-) sorted wrt "(c,ds)"!?? 618 619 // components should come in groups: count elements in each group 620 // && estimate the real size!!! 621 622 623 // use just a vector instead??? 624 const ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case! 625 626 int k = 0; 627 628 for (int j = 0; j < size; j++) 629 { 630 const poly p = id->m[j]; 631 assume( p != NULL ); 632 const int c = p_GetComp(p, r); 633 634 for (int i = j - 1; i >= 0; i--) 635 { 636 const poly pp = id->m[i]; 637 assume( pp != NULL ); 638 const int cc = p_GetComp(pp, r); 639 640 if( c != cc ) 641 continue; 642 643 const poly m = p_Init(r); // p_New??? 644 645 // m = LCM(p, pp) / p! // TODO: optimize: knowing the ring structure: (C/lp)! 646 for (int v = rVar(r); v > 0; v--) 647 { 648 assume( v > 0 ); 649 assume( v <= rVar(r) ); 650 651 const short e1 = p_GetExp(p , v, r); 652 const short e2 = p_GetExp(pp, v, r); 653 654 if( e1 >= e2 ) 655 p_SetExp(m, v, 0, r); 656 else 657 p_SetExp(m, v, e2 - e1, r); 658 659 } 660 661 assume( (j > i) && (i >= 0) ); 662 663 p_SetComp(m, j + 1, r); 664 pNext(m) = NULL; 665 p_SetCoeff0(m, n_Init(1, r->cf), r); // for later... 666 667 p_Setm(m, r); // should not do anything!!! 668 669 newid->m[k++] = m; 670 } 671 } 672 673 // if( __DEBUG__ && FALSE ) 674 // { 675 // PrintS("ComputeLeadingSyzygyTerms::Temp0: \n"); 676 // dPrint(newid, r, r, 1); 677 // } 678 679 // the rest of newid is assumed to be zeroes... 680 681 // simplify(newid, 2 + 32)?? 682 // sort(newid, "C,ds")[1]??? 683 id_DelDiv(newid, r); // #define SIMPL_LMDIV 32 684 685 // if( __DEBUG__ && FALSE ) 686 // { 687 // PrintS("ComputeLeadingSyzygyTerms::Temp1: \n"); 688 // dPrint(newid, r, r, 1); 689 // } 690 691 idSkipZeroes(newid); // #define SIMPL_NULL 2 692 693 // if( __DEBUG__ ) 694 // { 695 // PrintS("ComputeLeadingSyzygyTerms::Output: \n"); 696 // dPrint(newid, r, r, 1); 697 // } 698 699 const int sizeNew = IDELEMS(newid); 700 701 if( sizeNew >= 2 ) 702 qsort_my(newid->m, sizeNew, sizeof(poly), r, cmp_c_ds); 703 704 newid->rank = id_RankFreeModule(newid, r); 705 706 return newid; 707 } 708 709 710 static BOOLEAN ComputeLeadingSyzygyTerms(leftv res, leftv h) 508 static BOOLEAN _ComputeLeadingSyzygyTerms(leftv res, leftv h) 711 509 { 712 510 const ring r = currRing; … … 767 565 /// change the input inplace!!! 768 566 // TODO: use a ring with >_{c, ds}!??? 769 static BOOLEAN Sort_c_ds(leftv res, leftv h)567 static BOOLEAN _Sort_c_ds(leftv res, leftv h) 770 568 { 771 569 NoReturn(res); … … 817 615 id_Test(id, r); 818 616 819 const int size = IDELEMS(id); 820 821 const ideal newid = id; // id_Copy(id, r); // copy??? 822 823 if( size >= 2 ) 824 qsort_my(newid->m, size, sizeof(poly), r, cmp_c_ds); 825 826 // res->data = newid; 617 Sort_c_ds(id, r); // NOT A COPY! inplace sorting!!! 618 619 // res->data = id; 827 620 // res->rtyp = h->Typ(); 828 621 … … 830 623 { 831 624 PrintS("Sort_c_ds::Output: \n"); 832 dPrint( newid, r, r, 1);625 dPrint(id, r, r, 1); 833 626 } 834 627 628 // NOTE: nothing is to be returned!!! 835 629 return FALSE; 836 630 } … … 841 635 842 636 843 844 static ideal Compute2LeadingSyzygyTerms(const ideal& id, const ring r) 845 { 637 static BOOLEAN _Compute2LeadingSyzygyTerms(leftv res, leftv h) 638 { 639 const ring r = currRing; 640 NoReturn(res); 641 642 if( h == NULL ) 643 { 644 WarnS("Compute2LeadingSyzygyTerms needs an argument..."); 645 return TRUE; 646 } 647 648 assume( h != NULL ); 649 846 650 #ifndef NDEBUG 847 651 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE))); … … 850 654 #endif 851 655 852 const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));853 const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)0)));854 855 856 857 // 1. set of components S?858 // 2. for each component c from S: set of indices of leading terms859 // with this component?860 // 3. short exp. vectors for each leading term?861 862 const int size = IDELEMS(id);863 864 if( size < 2 )865 {866 const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module...867 return newid;868 }869 870 871 // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!??872 873 // components should come in groups: count elements in each group874 // && estimate the real size!!!875 876 877 // use just a vector instead???878 ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!879 880 int k = 0;881 882 for (int j = 0; j < size; j++)883 {884 const poly p = id->m[j];885 assume( p != NULL );886 const int c = p_GetComp(p, r);887 888 for (int i = j - 1; i >= 0; i--)889 {890 const poly pp = id->m[i];891 assume( pp != NULL );892 const int cc = p_GetComp(pp, r);893 894 if( c != cc )895 continue;896 897 // allocate memory & zero it out!898 const poly m = p_Init(r); const poly mm = p_Init(r);899 900 901 // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp!902 // TODO: optimize: knowing the ring structure: (C/lp)!903 904 for (int v = rVar(r); v > 0; v--)905 {906 assume( v > 0 );907 assume( v <= rVar(r) );908 909 const short e1 = p_GetExp(p , v, r);910 const short e2 = p_GetExp(pp, v, r);911 912 if( e1 >= e2 )913 p_SetExp(mm, v, e1 - e2, r); // p_SetExp(m, v, 0, r);914 else915 p_SetExp(m, v, e2 - e1, r); // p_SetExp(mm, v, 0, r);916 917 }918 919 assume( (j > i) && (i >= 0) );920 921 p_SetComp(m, j + 1, r);922 p_SetComp(mm, i + 1, r);923 924 const number& lc1 = p_GetCoeff(p , r);925 const number& lc2 = p_GetCoeff(pp, r);926 927 number g = n_Lcm( lc1, lc2, r );928 929 p_SetCoeff0(m , n_Div(g, lc1, r), r);930 p_SetCoeff0(mm, n_Neg(n_Div(g, lc2, r), r), r);931 932 n_Delete(&g, r);933 934 p_Setm(m, r); // should not do anything!!!935 p_Setm(mm, r); // should not do anything!!!936 937 pNext(m) = mm; // pNext(mm) = NULL;938 939 newid->m[k++] = m;940 }941 }942 943 // if( __DEBUG__ && FALSE )944 // {945 // PrintS("Compute2LeadingSyzygyTerms::Temp0: \n");946 // dPrint(newid, r, r, 1);947 // }948 949 if( !__TAILREDSYZ__ )950 {951 // simplify(newid, 2 + 32)??952 // sort(newid, "C,ds")[1]???953 id_DelDiv(newid, r); // #define SIMPL_LMDIV 32954 955 // if( __DEBUG__ && FALSE )956 // {957 // PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n");958 // dPrint(newid, r, r, 1);959 // }960 }961 else962 {963 // option(redSB); option(redTail);964 // TEST_OPT_REDSB965 // TEST_OPT_REDTAIL966 assume( r == currRing );967 BITSET _save_test = test;968 test |= (Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB));969 970 intvec* w=new intvec(IDELEMS(newid));971 ideal tmp = kStd(newid, currQuotient, isHomog, &w);972 delete w;973 974 test = _save_test;975 976 id_Delete(&newid, r);977 newid = tmp;978 979 // if( __DEBUG__ && FALSE )980 // {981 // PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n");982 // dPrint(newid, r, r, 1);983 // }984 985 }986 987 idSkipZeroes(newid);988 989 const int sizeNew = IDELEMS(newid);990 991 if( sizeNew >= 2 )992 qsort_my(newid->m, sizeNew, sizeof(poly), r, cmp_c_ds);993 994 newid->rank = id_RankFreeModule(newid, r);995 996 return newid;997 }998 999 1000 1001 static BOOLEAN Compute2LeadingSyzygyTerms(leftv res, leftv h)1002 {1003 const ring r = currRing;1004 NoReturn(res);1005 1006 if( h == NULL )1007 {1008 WarnS("Compute2LeadingSyzygyTerms needs an argument...");1009 return TRUE;1010 }1011 1012 assume( h != NULL );1013 1014 #ifndef NDEBUG1015 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));1016 #else1017 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));1018 #endif1019 1020 656 const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)0))); 1021 657 … … 1047 683 1048 684 1049 /// return a new term: leading coeff * leading monomial of p1050 /// with 0 leading component!1051 static inline poly leadmonom(const poly p, const ring r)1052 {1053 poly m = NULL;1054 1055 if( p != NULL )1056 {1057 assume( p != NULL );1058 assume( p_LmTest(p, r) );1059 1060 m = p_LmInit(p, r);1061 p_SetCoeff0(m, n_Copy(p_GetCoeff(p, r), r), r);1062 1063 p_SetComp(m, 0, r);1064 p_Setm(m, r);1065 1066 assume( p_GetComp(m, r) == 0 );1067 assume( m != NULL );1068 assume( pNext(m) == NULL );1069 assume( p_LmTest(m, r) );1070 }1071 1072 return m;1073 }1074 1075 1076 1077 static poly FindReducer(poly product, poly syzterm,1078 ideal L, ideal LS,1079 const ring r)1080 {1081 assume( product != NULL );1082 assume( L != NULL );1083 1084 int c = 0;1085 1086 if (syzterm != NULL)1087 c = p_GetComp(syzterm, r) - 1;1088 1089 assume( c >= 0 && c < IDELEMS(L) );1090 1091 #ifndef NDEBUG1092 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));1093 #else1094 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));1095 #endif1096 1097 long debug = __DEBUG__;1098 const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)debug)));1099 1100 if (__SYZCHECK__ && syzterm != NULL)1101 {1102 const poly m = L->m[c];1103 1104 assume( m != NULL ); assume( pNext(m) == NULL );1105 1106 poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);1107 assume( p_EqualPolys(lm, product, r) );1108 1109 // def @@c = leadcomp(syzterm); int @@r = int(@@c);1110 // def @@product = leadmonomial(syzterm) * L[@@r];1111 1112 p_Delete(&lm, r);1113 }1114 1115 // looking for an appropriate diviser q = L[k]...1116 for( int k = IDELEMS(L)-1; k>= 0; k-- )1117 {1118 const poly p = L->m[k];1119 1120 // ... which divides the product, looking for the _1st_ appropriate one!1121 if( !p_LmDivisibleBy(p, product, r) )1122 continue;1123 1124 1125 const poly q = p_New(r);1126 pNext(q) = NULL;1127 p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))1128 p_SetComp(q, k + 1, r);1129 p_Setm(q, r);1130 1131 // cannot allow something like: a*gen(i) - a*gen(i)1132 if (syzterm != NULL && (k == c))1133 if (p_ExpVectorEqual(syzterm, q, r))1134 {1135 if( __DEBUG__ )1136 {1137 Print("FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");1138 dPrint(syzterm, r, r, 1);1139 }1140 1141 p_LmFree(q, r);1142 continue;1143 }1144 1145 // while the complement (the fraction) is not reducible by leading syzygies1146 if( LS != NULL )1147 {1148 BOOLEAN ok = TRUE;1149 1150 for(int kk = IDELEMS(LS)-1; kk>= 0; kk-- )1151 {1152 const poly pp = LS->m[kk];1153 1154 if( p_LmDivisibleBy(pp, q, r) )1155 {1156 1157 if( __DEBUG__ )1158 {1159 Print("FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", kk+1);1160 dPrint(pp, r, r, 1);1161 }1162 1163 ok = FALSE; // q in <LS> :((1164 break;1165 }1166 }1167 1168 if(!ok)1169 {1170 p_LmFree(q, r);1171 continue;1172 }1173 }1174 1175 p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);1176 return q;1177 1178 }1179 1180 1181 return NULL;1182 }1183 1184 685 1185 686 /// TODO: save shortcut (syz: |-.->) LM(LM(m) * "t") -> syz? 1186 687 /// proc SSFindReducer(def product, def syzterm, def L, def T, list #) 1187 static BOOLEAN FindReducer(leftv res, leftv h)688 static BOOLEAN _FindReducer(leftv res, leftv h) 1188 689 { 1189 690 const char* usage = "`FindReducer(<poly/vector>, <vector/0>, <ideal/module>[,<module>])` expected"; … … 1278 779 } 1279 780 1280 static poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,1281 ideal L, ideal T, ideal LS, const ring r);1282 1283 static poly TraverseTail(poly multiplier, poly tail,1284 ideal L, ideal T, ideal LS,1285 const ring r)1286 {1287 assume( multiplier != NULL );1288 1289 assume( L != NULL );1290 assume( T != NULL );1291 1292 poly s = NULL;1293 1294 // iterate over the tail1295 for(poly p = tail; p != NULL; p = pNext(p))1296 s = p_Add_q(s, ReduceTerm(multiplier, p, NULL, L, T, LS, r), r);1297 1298 return s;1299 }1300 1301 1302 static poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,1303 ideal L, ideal T, ideal LS, const ring r)1304 {1305 assume( multiplier != NULL );1306 assume( term4reduction != NULL );1307 1308 1309 assume( L != NULL );1310 assume( T != NULL );1311 1312 // assume(r == currRing); // ?1313 1314 // simple implementation with FindReducer:1315 poly s = NULL;1316 1317 if(1)1318 {1319 // NOTE: only LT(term4reduction) should be used in the following:1320 poly product = pp_Mult_mm(multiplier, term4reduction, r);1321 s = FindReducer(product, syztermCheck, L, LS, r);1322 p_Delete(&product, r);1323 }1324 1325 if( s == NULL ) // No Reducer?1326 return s;1327 1328 poly b = leadmonom(s, r);1329 1330 const int c = p_GetComp(s, r) - 1;1331 assume( c >= 0 && c < IDELEMS(T) );1332 1333 const poly tail = T->m[c];1334 1335 if( tail != NULL )1336 s = p_Add_q(s, TraverseTail(b, tail, L, T, LS, r), r);1337 1338 return s;1339 }1340 1341 1342 static poly SchreyerSyzygyNF(poly syz_lead, poly syz_2, ideal L, ideal T, ideal LS, const ring r)1343 {1344 assume( syz_lead != NULL );1345 assume( syz_2 != NULL );1346 1347 assume( L != NULL );1348 assume( T != NULL );1349 1350 assume( IDELEMS(L) == IDELEMS(T) );1351 1352 int c = p_GetComp(syz_lead, r) - 1;1353 1354 assume( c >= 0 && c < IDELEMS(T) );1355 1356 poly p = leadmonom(syz_lead, r); // :(1357 poly spoly = pp_Mult_qq(p, T->m[c], r);1358 p_Delete(&p, r);1359 1360 1361 c = p_GetComp(syz_2, r) - 1;1362 assume( c >= 0 && c < IDELEMS(T) );1363 1364 p = leadmonom(syz_2, r); // :(1365 spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);1366 p_Delete(&p, r);1367 1368 poly tail = p_Copy(syz_2, r); // TODO: use bucket!?1369 1370 while (spoly != NULL)1371 {1372 poly t = FindReducer(spoly, NULL, L, LS, r);1373 1374 p_LmDelete(&spoly, r);1375 1376 if( t != NULL )1377 {1378 p = leadmonom(t, r); // :(1379 c = p_GetComp(t, r) - 1;1380 1381 assume( c >= 0 && c < IDELEMS(T) );1382 1383 spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);1384 1385 p_Delete(&p, r);1386 1387 tail = p_Add_q(tail, t, r);1388 }1389 }1390 1391 return tail;1392 }1393 1394 781 // proc SchreyerSyzygyNF(vector syz_lead, vector syz_2, def L, def T, list #) 1395 static BOOLEAN SchreyerSyzygyNF(leftv res, leftv h)782 static BOOLEAN _SchreyerSyzygyNF(leftv res, leftv h) 1396 783 { 1397 784 const char* usage = "`SchreyerSyzygyNF(<vector>, <vector>, <ideal/module>, <ideal/module>[,<module>])` expected"; … … 1503 890 /// TODO: save shortcut (syz: |-.->) LM(m) * "t" -> ? 1504 891 /// proc SSReduceTerm(poly m, def t, def syzterm, def L, def T, list #) 1505 static BOOLEAN ReduceTerm(leftv res, leftv h)892 static BOOLEAN _ReduceTerm(leftv res, leftv h) 1506 893 { 1507 894 const char* usage = "`ReduceTerm(<poly>, <poly/vector>, <vector/0>, <ideal/module>, <ideal/module>[,<module>])` expected"; … … 1648 1035 // TODO: store m * @tail -.-^-.-^-.--> ? 1649 1036 // proc SSTraverseTail(poly m, def @tail, def L, def T, list #) 1650 static BOOLEAN TraverseTail(leftv res, leftv h)1037 static BOOLEAN _TraverseTail(leftv res, leftv h) 1651 1038 { 1652 1039 const char* usage = "`TraverseTail(<poly>, <poly/vector>, <ideal/module>, <ideal/module>[,<module>])` expected"; … … 1750 1137 1751 1138 1752 1753 static void ComputeSyzygy(const ideal L, const ideal T, ideal& LL, ideal& TT, const ring R) 1754 { 1755 assume( R == currRing ); // For attributes :-/ 1756 1757 assume( IDELEMS(L) == IDELEMS(T) ); 1758 1139 // module (N, LL, TT) = SSComputeSyzygy(L, T); 1140 // Compute Syz(L ++ T) = N = LL ++ TT 1141 // proc SSComputeSyzygy(def L, def T) 1142 static BOOLEAN _ComputeSyzygy(leftv res, leftv h) 1143 { 1144 const char* usage = "`ComputeSyzygy(<ideal/module>, <ideal/module>])` expected"; 1145 const ring r = currRing; 1146 1147 NoReturn(res); 1148 1759 1149 #ifndef NDEBUG 1760 1150 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE))); … … 1763 1153 #endif 1764 1154 1765 const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)1)));1766 const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)1)));1767 const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)__DEBUG__)));1768 1769 const BOOLEAN __HYBRIDNF__ = (BOOLEAN)((long)(atGet(currRingHdl,"HYBRIDNF",INT_CMD, (void*)0)));1770 1771 1772 if( __LEAD2SYZ__ )1773 LL = Compute2LeadingSyzygyTerms(L, R); // 2 terms!1774 else1775 LL = ComputeLeadingSyzygyTerms(L, R); // 1 term!1776 1777 const int size = IDELEMS(LL);1778 1779 TT = idInit(size, 0);1780 1781 if( size == 1 && LL->m[0] == NULL )1782 return;1783 1784 1785 ideal LS = NULL;1786 1787 if( __TAILREDSYZ__ )1788 LS = LL;1789 1790 for( int k = size - 1; k >= 0; k-- )1791 {1792 const poly a = LL->m[k]; assume( a != NULL );1793 1794 const int r = p_GetComp(a, R) - 1;1795 1796 assume( r >= 0 && r < IDELEMS(T) );1797 assume( r >= 0 && r < IDELEMS(L) );1798 1799 poly aa = leadmonom(a, R); assume( aa != NULL); // :(1800 1801 poly a2 = pNext(a);1802 1803 if( a2 != NULL )1804 {1805 TT->m[k] = a2; pNext(a) = NULL;1806 }1807 1808 if( ! __HYBRIDNF__ )1809 {1810 poly t = TraverseTail(aa, T->m[r], L, T, LS, R);1811 1812 if( a2 != NULL )1813 {1814 assume( __LEAD2SYZ__ );1815 1816 const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(1817 1818 assume( r2 >= 0 && r2 < IDELEMS(T) );1819 1820 TT->m[k] = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, T->m[r2], L, T, LS, R), R), R);1821 1822 p_Delete(&aa2, R);1823 } else1824 TT->m[k] = p_Add_q(t, ReduceTerm(aa, L->m[r], a, L, T, LS, R), R);1825 1826 } else1827 {1828 if( a2 == NULL )1829 {1830 aa = p_Mult_mm(aa, L->m[r], R);1831 a2 = FindReducer(aa, a, L, LS, R);1832 }1833 assume( a2 != NULL );1834 1835 TT->m[k] = SchreyerSyzygyNF(a, a2, L, T, LS, R); // will copy a2 :(1836 1837 p_Delete(&a2, R);1838 }1839 p_Delete(&aa, R);1840 }1841 1842 TT->rank = id_RankFreeModule(TT, R);1843 }1844 1845 1846 1847 // module (N, LL, TT) = SSComputeSyzygy(L, T);1848 // Compute Syz(L ++ T) = N = LL ++ TT1849 // proc SSComputeSyzygy(def L, def T)1850 static BOOLEAN ComputeSyzygy(leftv res, leftv h)1851 {1852 const char* usage = "`ComputeSyzygy(<ideal/module>, <ideal/module>])` expected";1853 const ring r = currRing;1854 1855 NoReturn(res);1856 1857 #ifndef NDEBUG1858 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));1859 #else1860 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));1861 #endif1862 1863 1155 if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL)) 1864 1156 { … … 1915 1207 } 1916 1208 1917 1918 1919 1920 1921 1922 1923 1209 /// Get leading term without a module component 1924 static BOOLEAN leadmonom(leftv res, leftv h)1210 static BOOLEAN _leadmonom(leftv res, leftv h) 1925 1211 { 1926 1212 NoReturn(res); … … 1931 1217 const poly p = (poly)(h->Data()); 1932 1218 1933 const poly m = leadmonom(p, r); 1934 1219 res->data = reinterpret_cast<void *>( leadmonom(p, r) ); 1935 1220 res->rtyp = POLY_CMD; 1936 res->data = reinterpret_cast<void *>(m);1937 1221 1938 1222 return FALSE; … … 2530 1814 2531 1815 ADD(psModulFunctions, currPack->libname, "DetailedPrint", FALSE, DetailedPrint); 2532 ADD(psModulFunctions, currPack->libname, "leadmonomial", FALSE, leadmonom);1816 ADD(psModulFunctions, currPack->libname, "leadmonomial", FALSE, _leadmonom); 2533 1817 ADD(psModulFunctions, currPack->libname, "leadcomp", FALSE, leadcomp); 2534 1818 ADD(psModulFunctions, currPack->libname, "leadrawexp", FALSE, leadrawexp); … … 2551 1835 ADD(psModulFunctions, currPack->libname, "Tail", FALSE, Tail); 2552 1836 2553 ADD(psModulFunctions, currPack->libname, "ComputeLeadingSyzygyTerms", FALSE, ComputeLeadingSyzygyTerms);2554 ADD(psModulFunctions, currPack->libname, "Compute2LeadingSyzygyTerms", FALSE, Compute2LeadingSyzygyTerms);2555 2556 ADD(psModulFunctions, currPack->libname, "Sort_c_ds", FALSE, Sort_c_ds);2557 ADD(psModulFunctions, currPack->libname, "FindReducer", FALSE, FindReducer);2558 2559 2560 ADD(psModulFunctions, currPack->libname, "ReduceTerm", FALSE, ReduceTerm);2561 ADD(psModulFunctions, currPack->libname, "TraverseTail", FALSE, TraverseTail);1837 ADD(psModulFunctions, currPack->libname, "ComputeLeadingSyzygyTerms", FALSE, _ComputeLeadingSyzygyTerms); 1838 ADD(psModulFunctions, currPack->libname, "Compute2LeadingSyzygyTerms", FALSE, _Compute2LeadingSyzygyTerms); 1839 1840 ADD(psModulFunctions, currPack->libname, "Sort_c_ds", FALSE, _Sort_c_ds); 1841 ADD(psModulFunctions, currPack->libname, "FindReducer", FALSE, _FindReducer); 1842 1843 1844 ADD(psModulFunctions, currPack->libname, "ReduceTerm", FALSE, _ReduceTerm); 1845 ADD(psModulFunctions, currPack->libname, "TraverseTail", FALSE, _TraverseTail); 2562 1846 2563 1847 2564 ADD(psModulFunctions, currPack->libname, "SchreyerSyzygyNF", FALSE, SchreyerSyzygyNF);2565 ADD(psModulFunctions, currPack->libname, "ComputeSyzygy", FALSE, ComputeSyzygy);1848 ADD(psModulFunctions, currPack->libname, "SchreyerSyzygyNF", FALSE, _SchreyerSyzygyNF); 1849 ADD(psModulFunctions, currPack->libname, "ComputeSyzygy", FALSE, _ComputeSyzygy); 2566 1850 2567 1851 // ADD(psModulFunctions, currPack->libname, "GetAMData", FALSE, GetAMData); -
dyn_modules/syzextra/syzextra.cc
rcd5fefc r204092 19 19 #include "syzextra.h" 20 20 21 #include "DebugPrint.h" 22 21 23 #include <omalloc/omalloc.h> 24 25 #include <misc/intvec.h> 26 #include <misc/options.h> 27 28 #include <coeffs/coeffs.h> 29 22 30 #include <polys/monomials/p_polys.h> 23 31 #include <polys/monomials/ring.h> 32 33 #include <kernel/kstd1.h> 34 #include <kernel/polys.h> 35 #include <kernel/syz.h> 24 36 #include <kernel/ideals.h> 25 37 26 38 39 40 #include <Singular/tok.h> 41 #include <Singular/ipid.h> 42 #include <Singular/lists.h> 43 #include <Singular/attrib.h> 44 45 #include <Singular/ipid.h> 46 #include <Singular/ipshell.h> // For iiAddCproc 47 48 #include <stdio.h> 49 #include <stdlib.h> 50 #include <string.h> 51 52 // USING_NAMESPACE_SINGULARXX; 53 USING_NAMESPACE( SINGULARXXNAME :: DEBUG ) 54 55 27 56 BEGIN_NAMESPACE_SINGULARXX BEGIN_NAMESPACE(SYZEXTRA) 28 57 29 58 BEGIN_NAMESPACE(SORT_c_ds) 59 60 61 #ifdef _GNU_SOURCE 62 static int cmp_c_ds(const void *p1, const void *p2, void *R) 63 { 64 #else 65 static int cmp_c_ds(const void *p1, const void *p2) 66 { 67 void *R = currRing; 68 #endif 69 70 const int YES = 1; 71 const int NO = -1; 72 73 const ring r = (const ring) R; // TODO/NOTE: the structure is known: C, lp!!! 74 75 assume( r == currRing ); 76 77 const poly a = *(const poly*)p1; 78 const poly b = *(const poly*)p2; 79 80 assume( a != NULL ); 81 assume( b != NULL ); 82 83 assume( p_LmTest(a, r) ); 84 assume( p_LmTest(b, r) ); 85 86 87 const signed long iCompDiff = p_GetComp(a, r) - p_GetComp(b, r); 88 89 // TODO: test this!!!!!!!!!!!!!!!! 90 91 //return -( compare (c, qsorts) ) 92 93 const int __DEBUG__ = 0; 94 95 #ifndef NDEBUG 96 if( __DEBUG__ ) 97 { 98 PrintS("cmp_c_ds: a, b: \np1: "); dPrint(a, r, r, 2); 99 PrintS("b: "); dPrint(b, r, r, 2); 100 PrintLn(); 101 } 102 #endif 103 104 105 if( iCompDiff > 0 ) 106 return YES; 107 108 if( iCompDiff < 0 ) 109 return NO; 110 111 assume( iCompDiff == 0 ); 112 113 const signed long iDegDiff = p_Totaldegree(a, r) - p_Totaldegree(b, r); 114 115 if( iDegDiff > 0 ) 116 return YES; 117 118 if( iDegDiff < 0 ) 119 return NO; 120 121 assume( iDegDiff == 0 ); 122 123 #ifndef NDEBUG 124 if( __DEBUG__ ) 125 { 126 PrintS("cmp_c_ds: a & b have the same comp & deg! "); PrintLn(); 127 } 128 #endif 129 130 for (int v = rVar(r); v > 0; v--) 131 { 132 assume( v > 0 ); 133 assume( v <= rVar(r) ); 134 135 const signed int d = p_GetExp(a, v, r) - p_GetExp(b, v, r); 136 137 if( d > 0 ) 138 return YES; 139 140 if( d < 0 ) 141 return NO; 142 143 assume( d == 0 ); 144 } 145 146 return 0; 147 } 148 149 END_NAMESPACE 150 /* namespace SORT_c_ds */ 151 152 /// return a new term: leading coeff * leading monomial of p 153 /// with 0 leading component! 154 poly leadmonom(const poly p, const ring r) 155 { 156 poly m = NULL; 157 158 if( p != NULL ) 159 { 160 assume( p != NULL ); 161 assume( p_LmTest(p, r) ); 162 163 m = p_LmInit(p, r); 164 p_SetCoeff0(m, n_Copy(p_GetCoeff(p, r), r), r); 165 166 p_SetComp(m, 0, r); 167 p_Setm(m, r); 168 169 assume( p_GetComp(m, r) == 0 ); 170 assume( m != NULL ); 171 assume( pNext(m) == NULL ); 172 assume( p_LmTest(m, r) ); 173 } 174 175 return m; 176 } 177 178 179 30 180 poly p_Tail(const poly p, const ring r) 31 181 { … … 54 204 55 205 206 void Sort_c_ds(const ideal id, const ring r) 207 { 208 const int sizeNew = IDELEMS(id); 209 210 #ifdef _GNU_SOURCE 211 #define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp, r) 212 #else 213 #define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp) 214 #endif 215 216 if( sizeNew >= 2 ) 217 qsort_my(id->m, sizeNew, sizeof(poly), r, SORT_c_ds::cmp_c_ds); 218 219 #undef qsort_my 220 221 id->rank = id_RankFreeModule(id, r); 222 } 223 224 225 226 ideal ComputeLeadingSyzygyTerms(const ideal& id, const ring r) 227 { 228 // 1. set of components S? 229 // 2. for each component c from S: set of indices of leading terms 230 // with this component? 231 // 3. short exp. vectors for each leading term? 232 233 const int size = IDELEMS(id); 234 235 if( size < 2 ) 236 { 237 const ideal newid = idInit(1, 0); newid->m[0] = NULL; // zero ideal... 238 return newid; 239 } 240 241 242 // TODO/NOTE: input is supposed to be (reverse-) sorted wrt "(c,ds)"!?? 243 244 // components should come in groups: count elements in each group 245 // && estimate the real size!!! 246 247 248 // use just a vector instead??? 249 const ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case! 250 251 int k = 0; 252 253 for (int j = 0; j < size; j++) 254 { 255 const poly p = id->m[j]; 256 assume( p != NULL ); 257 const int c = p_GetComp(p, r); 258 259 for (int i = j - 1; i >= 0; i--) 260 { 261 const poly pp = id->m[i]; 262 assume( pp != NULL ); 263 const int cc = p_GetComp(pp, r); 264 265 if( c != cc ) 266 continue; 267 268 const poly m = p_Init(r); // p_New??? 269 270 // m = LCM(p, pp) / p! // TODO: optimize: knowing the ring structure: (C/lp)! 271 for (int v = rVar(r); v > 0; v--) 272 { 273 assume( v > 0 ); 274 assume( v <= rVar(r) ); 275 276 const short e1 = p_GetExp(p , v, r); 277 const short e2 = p_GetExp(pp, v, r); 278 279 if( e1 >= e2 ) 280 p_SetExp(m, v, 0, r); 281 else 282 p_SetExp(m, v, e2 - e1, r); 283 284 } 285 286 assume( (j > i) && (i >= 0) ); 287 288 p_SetComp(m, j + 1, r); 289 pNext(m) = NULL; 290 p_SetCoeff0(m, n_Init(1, r->cf), r); // for later... 291 292 p_Setm(m, r); // should not do anything!!! 293 294 newid->m[k++] = m; 295 } 296 } 297 298 // if( __DEBUG__ && FALSE ) 299 // { 300 // PrintS("ComputeLeadingSyzygyTerms::Temp0: \n"); 301 // dPrint(newid, r, r, 1); 302 // } 303 304 // the rest of newid is assumed to be zeroes... 305 306 // simplify(newid, 2 + 32)?? 307 // sort(newid, "C,ds")[1]??? 308 id_DelDiv(newid, r); // #define SIMPL_LMDIV 32 309 310 // if( __DEBUG__ && FALSE ) 311 // { 312 // PrintS("ComputeLeadingSyzygyTerms::Temp1: \n"); 313 // dPrint(newid, r, r, 1); 314 // } 315 316 idSkipZeroes(newid); // #define SIMPL_NULL 2 317 318 // if( __DEBUG__ ) 319 // { 320 // PrintS("ComputeLeadingSyzygyTerms::Output: \n"); 321 // dPrint(newid, r, r, 1); 322 // } 323 324 Sort_c_ds(newid, r); 325 326 return newid; 327 } 328 329 330 331 332 ideal Compute2LeadingSyzygyTerms(const ideal& id, const ring r) 333 { 334 #ifndef NDEBUG 335 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE))); 336 #else 337 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE))); 338 #endif 339 340 const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0))); 341 const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)0))); 342 343 344 345 // 1. set of components S? 346 // 2. for each component c from S: set of indices of leading terms 347 // with this component? 348 // 3. short exp. vectors for each leading term? 349 350 const int size = IDELEMS(id); 351 352 if( size < 2 ) 353 { 354 const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module... 355 return newid; 356 } 357 358 359 // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!?? 360 361 // components should come in groups: count elements in each group 362 // && estimate the real size!!! 363 364 365 // use just a vector instead??? 366 ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case! 367 368 int k = 0; 369 370 for (int j = 0; j < size; j++) 371 { 372 const poly p = id->m[j]; 373 assume( p != NULL ); 374 const int c = p_GetComp(p, r); 375 376 for (int i = j - 1; i >= 0; i--) 377 { 378 const poly pp = id->m[i]; 379 assume( pp != NULL ); 380 const int cc = p_GetComp(pp, r); 381 382 if( c != cc ) 383 continue; 384 385 // allocate memory & zero it out! 386 const poly m = p_Init(r); const poly mm = p_Init(r); 387 388 389 // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp! 390 // TODO: optimize: knowing the ring structure: (C/lp)! 391 392 for (int v = rVar(r); v > 0; v--) 393 { 394 assume( v > 0 ); 395 assume( v <= rVar(r) ); 396 397 const short e1 = p_GetExp(p , v, r); 398 const short e2 = p_GetExp(pp, v, r); 399 400 if( e1 >= e2 ) 401 p_SetExp(mm, v, e1 - e2, r); // p_SetExp(m, v, 0, r); 402 else 403 p_SetExp(m, v, e2 - e1, r); // p_SetExp(mm, v, 0, r); 404 405 } 406 407 assume( (j > i) && (i >= 0) ); 408 409 p_SetComp(m, j + 1, r); 410 p_SetComp(mm, i + 1, r); 411 412 const number& lc1 = p_GetCoeff(p , r); 413 const number& lc2 = p_GetCoeff(pp, r); 414 415 number g = n_Lcm( lc1, lc2, r ); 416 417 p_SetCoeff0(m , n_Div(g, lc1, r), r); 418 p_SetCoeff0(mm, n_Neg(n_Div(g, lc2, r), r), r); 419 420 n_Delete(&g, r); 421 422 p_Setm(m, r); // should not do anything!!! 423 p_Setm(mm, r); // should not do anything!!! 424 425 pNext(m) = mm; // pNext(mm) = NULL; 426 427 newid->m[k++] = m; 428 } 429 } 430 431 // if( __DEBUG__ && FALSE ) 432 // { 433 // PrintS("Compute2LeadingSyzygyTerms::Temp0: \n"); 434 // dPrint(newid, r, r, 1); 435 // } 436 437 if( !__TAILREDSYZ__ ) 438 { 439 // simplify(newid, 2 + 32)?? 440 // sort(newid, "C,ds")[1]??? 441 id_DelDiv(newid, r); // #define SIMPL_LMDIV 32 442 443 // if( __DEBUG__ && FALSE ) 444 // { 445 // PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n"); 446 // dPrint(newid, r, r, 1); 447 // } 448 } 449 else 450 { 451 // option(redSB); option(redTail); 452 // TEST_OPT_REDSB 453 // TEST_OPT_REDTAIL 454 assume( r == currRing ); 455 BITSET _save_test = test; 456 test |= (Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB)); 457 458 intvec* w=new intvec(IDELEMS(newid)); 459 ideal tmp = kStd(newid, currQuotient, isHomog, &w); 460 delete w; 461 462 test = _save_test; 463 464 id_Delete(&newid, r); 465 newid = tmp; 466 467 // if( __DEBUG__ && FALSE ) 468 // { 469 // PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n"); 470 // dPrint(newid, r, r, 1); 471 // } 472 473 } 474 475 idSkipZeroes(newid); 476 477 Sort_c_ds(newid, r); 478 479 return newid; 480 } 481 482 poly FindReducer(poly product, poly syzterm, 483 ideal L, ideal LS, 484 const ring r) 485 { 486 assume( product != NULL ); 487 assume( L != NULL ); 488 489 int c = 0; 490 491 if (syzterm != NULL) 492 c = p_GetComp(syzterm, r) - 1; 493 494 assume( c >= 0 && c < IDELEMS(L) ); 495 496 #ifndef NDEBUG 497 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE))); 498 #else 499 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE))); 500 #endif 501 502 long debug = __DEBUG__; 503 const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)debug))); 504 505 if (__SYZCHECK__ && syzterm != NULL) 506 { 507 const poly m = L->m[c]; 508 509 assume( m != NULL ); assume( pNext(m) == NULL ); 510 511 poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r); 512 assume( p_EqualPolys(lm, product, r) ); 513 514 // def @@c = leadcomp(syzterm); int @@r = int(@@c); 515 // def @@product = leadmonomial(syzterm) * L[@@r]; 516 517 p_Delete(&lm, r); 518 } 519 520 // looking for an appropriate diviser q = L[k]... 521 for( int k = IDELEMS(L)-1; k>= 0; k-- ) 522 { 523 const poly p = L->m[k]; 524 525 // ... which divides the product, looking for the _1st_ appropriate one! 526 if( !p_LmDivisibleBy(p, product, r) ) 527 continue; 528 529 530 const poly q = p_New(r); 531 pNext(q) = NULL; 532 p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k])) 533 p_SetComp(q, k + 1, r); 534 p_Setm(q, r); 535 536 // cannot allow something like: a*gen(i) - a*gen(i) 537 if (syzterm != NULL && (k == c)) 538 if (p_ExpVectorEqual(syzterm, q, r)) 539 { 540 if( __DEBUG__ ) 541 { 542 Print("FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: "); 543 dPrint(syzterm, r, r, 1); 544 } 545 546 p_LmFree(q, r); 547 continue; 548 } 549 550 // while the complement (the fraction) is not reducible by leading syzygies 551 if( LS != NULL ) 552 { 553 BOOLEAN ok = TRUE; 554 555 for(int kk = IDELEMS(LS)-1; kk>= 0; kk-- ) 556 { 557 const poly pp = LS->m[kk]; 558 559 if( p_LmDivisibleBy(pp, q, r) ) 560 { 561 562 if( __DEBUG__ ) 563 { 564 Print("FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", kk+1); 565 dPrint(pp, r, r, 1); 566 } 567 568 ok = FALSE; // q in <LS> :(( 569 break; 570 } 571 } 572 573 if(!ok) 574 { 575 p_LmFree(q, r); 576 continue; 577 } 578 } 579 580 p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r); 581 return q; 582 583 } 584 585 586 return NULL; 587 } 588 589 poly TraverseTail(poly multiplier, poly tail, 590 ideal L, ideal T, ideal LS, 591 const ring r) 592 { 593 assume( multiplier != NULL ); 594 595 assume( L != NULL ); 596 assume( T != NULL ); 597 598 poly s = NULL; 599 600 // iterate over the tail 601 for(poly p = tail; p != NULL; p = pNext(p)) 602 s = p_Add_q(s, ReduceTerm(multiplier, p, NULL, L, T, LS, r), r); 603 604 return s; 605 } 606 607 608 poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck, 609 ideal L, ideal T, ideal LS, const ring r) 610 { 611 assume( multiplier != NULL ); 612 assume( term4reduction != NULL ); 613 614 615 assume( L != NULL ); 616 assume( T != NULL ); 617 618 // assume(r == currRing); // ? 619 620 // simple implementation with FindReducer: 621 poly s = NULL; 622 623 if(1) 624 { 625 // NOTE: only LT(term4reduction) should be used in the following: 626 poly product = pp_Mult_mm(multiplier, term4reduction, r); 627 s = FindReducer(product, syztermCheck, L, LS, r); 628 p_Delete(&product, r); 629 } 630 631 if( s == NULL ) // No Reducer? 632 return s; 633 634 poly b = leadmonom(s, r); 635 636 const int c = p_GetComp(s, r) - 1; 637 assume( c >= 0 && c < IDELEMS(T) ); 638 639 const poly tail = T->m[c]; 640 641 if( tail != NULL ) 642 s = p_Add_q(s, TraverseTail(b, tail, L, T, LS, r), r); 643 644 return s; 645 } 646 647 648 poly SchreyerSyzygyNF(poly syz_lead, poly syz_2, ideal L, ideal T, ideal LS, const ring r) 649 { 650 assume( syz_lead != NULL ); 651 assume( syz_2 != NULL ); 652 653 assume( L != NULL ); 654 assume( T != NULL ); 655 656 assume( IDELEMS(L) == IDELEMS(T) ); 657 658 int c = p_GetComp(syz_lead, r) - 1; 659 660 assume( c >= 0 && c < IDELEMS(T) ); 661 662 poly p = leadmonom(syz_lead, r); // :( 663 poly spoly = pp_Mult_qq(p, T->m[c], r); 664 p_Delete(&p, r); 665 666 667 c = p_GetComp(syz_2, r) - 1; 668 assume( c >= 0 && c < IDELEMS(T) ); 669 670 p = leadmonom(syz_2, r); // :( 671 spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r); 672 p_Delete(&p, r); 673 674 poly tail = p_Copy(syz_2, r); // TODO: use bucket!? 675 676 while (spoly != NULL) 677 { 678 poly t = FindReducer(spoly, NULL, L, LS, r); 679 680 p_LmDelete(&spoly, r); 681 682 if( t != NULL ) 683 { 684 p = leadmonom(t, r); // :( 685 c = p_GetComp(t, r) - 1; 686 687 assume( c >= 0 && c < IDELEMS(T) ); 688 689 spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r); 690 691 p_Delete(&p, r); 692 693 tail = p_Add_q(tail, t, r); 694 } 695 } 696 697 return tail; 698 } 699 700 701 void ComputeSyzygy(const ideal L, const ideal T, ideal& LL, ideal& TT, const ring R) 702 { 703 assume( R == currRing ); // For attributes :-/ 704 705 assume( IDELEMS(L) == IDELEMS(T) ); 706 707 #ifndef NDEBUG 708 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE))); 709 #else 710 const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE))); 711 #endif 712 713 const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)1))); 714 const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)1))); 715 const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)__DEBUG__))); 716 717 const BOOLEAN __HYBRIDNF__ = (BOOLEAN)((long)(atGet(currRingHdl,"HYBRIDNF",INT_CMD, (void*)0))); 718 719 720 if( __LEAD2SYZ__ ) 721 LL = Compute2LeadingSyzygyTerms(L, R); // 2 terms! 722 else 723 LL = ComputeLeadingSyzygyTerms(L, R); // 1 term! 724 725 const int size = IDELEMS(LL); 726 727 TT = idInit(size, 0); 728 729 if( size == 1 && LL->m[0] == NULL ) 730 return; 731 732 733 ideal LS = NULL; 734 735 if( __TAILREDSYZ__ ) 736 LS = LL; 737 738 for( int k = size - 1; k >= 0; k-- ) 739 { 740 const poly a = LL->m[k]; assume( a != NULL ); 741 742 const int r = p_GetComp(a, R) - 1; 743 744 assume( r >= 0 && r < IDELEMS(T) ); 745 assume( r >= 0 && r < IDELEMS(L) ); 746 747 poly aa = leadmonom(a, R); assume( aa != NULL); // :( 748 749 poly a2 = pNext(a); 750 751 if( a2 != NULL ) 752 { 753 TT->m[k] = a2; pNext(a) = NULL; 754 } 755 756 if( ! __HYBRIDNF__ ) 757 { 758 poly t = TraverseTail(aa, T->m[r], L, T, LS, R); 759 760 if( a2 != NULL ) 761 { 762 assume( __LEAD2SYZ__ ); 763 764 const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :( 765 766 assume( r2 >= 0 && r2 < IDELEMS(T) ); 767 768 TT->m[k] = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, T->m[r2], L, T, LS, R), R), R); 769 770 p_Delete(&aa2, R); 771 } else 772 TT->m[k] = p_Add_q(t, ReduceTerm(aa, L->m[r], a, L, T, LS, R), R); 773 774 } else 775 { 776 if( a2 == NULL ) 777 { 778 aa = p_Mult_mm(aa, L->m[r], R); 779 a2 = FindReducer(aa, a, L, LS, R); 780 } 781 assume( a2 != NULL ); 782 783 TT->m[k] = SchreyerSyzygyNF(a, a2, L, T, LS, R); // will copy a2 :( 784 785 p_Delete(&a2, R); 786 } 787 p_Delete(&aa, R); 788 } 789 790 TT->rank = id_RankFreeModule(TT, R); 791 } 792 793 794 56 795 57 796 -
dyn_modules/syzextra/syzextra.h
rcd5fefc r204092 33 33 BEGIN_NAMESPACE_SINGULARXX BEGIN_NAMESPACE(SYZEXTRA) 34 34 35 poly leadmonom(const poly p, const ring r); 35 36 36 37 /// return the tail of a given polynomial or vector … … 47 48 48 49 50 /// inplace sorting of the module (ideal) id wrt >_(c,ds) 51 void Sort_c_ds(const ideal id, const ring r); 52 53 ideal ComputeLeadingSyzygyTerms(const ideal& id, const ring r); 54 55 ideal Compute2LeadingSyzygyTerms(const ideal& id, const ring r); 56 57 58 poly FindReducer(poly product, poly syzterm, 59 ideal L, ideal LS, 60 const ring r); 61 62 63 poly TraverseTail(poly multiplier, poly tail, 64 ideal L, ideal T, ideal LS, 65 const ring r); 66 67 poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck, 68 ideal L, ideal T, ideal LS, const ring r); 69 70 71 72 poly SchreyerSyzygyNF(poly syz_lead, poly syz_2, ideal L, ideal T, ideal LS, const ring r); 73 74 75 void ComputeSyzygy(const ideal L, const ideal T, ideal& LL, ideal& TT, const ring R); 76 49 77 50 78 END_NAMESPACE END_NAMESPACE_SINGULARXX
Note: See TracChangeset
for help on using the changeset viewer.