48 #define MAXPOINTS 10000 49 #define MAXINITELEMS 256 50 #define LIFT_COOR 50000 // siRand() % LIFT_COOR gives random lift value 51 #define SCALEDOWN 100.0 // lift value scale down for linear program 53 #define RVMULT 0.0001 // multiplicator for random shift vector 54 #define MAXRVVAL 50000 81 number
getDetAt(
const number* evpoint );
83 poly
getUDet(
const number* evpoint );
146 typedef struct onePoint * onePointP;
156 typedef struct _entry * entry;
176 inline onePointP operator[] (
const int index );
182 bool addPoint(
const onePointP vert );
188 bool addPoint(
const int * vert );
194 bool addPoint(
const Coord_t * vert );
197 bool removePoint(
const int indx );
202 bool mergeWithExp(
const onePointP vert );
207 bool mergeWithExp(
const int * vert );
210 void mergeWithPoly(
const poly
p );
213 void getRowMP(
const int indx,
int * vert );
216 int getExpPos(
const poly
p );
235 inline bool smaller(
int,
int );
238 inline bool larger(
int,
int );
243 inline bool checkMem();
260 ideal newtonPolytopesI(
const ideal gls );
266 bool inHull(poly
p, poly pointPoly,
int m,
int site);
295 void runMayanPyramid(
int dim );
314 bool storeMinkowskiSumPoint();
330 #if defined(mprDEBUG_PROT) || defined(mprDEBUG_ALL) 331 void print_mat(
mprfloat **a,
int maxrow,
int maxcol)
335 for (i = 1; i <= maxrow; i++)
338 for (j = 1; j <= maxcol; j++)
Print(
"% 7.2f, ", a[i][j]);
346 printf(
"Output matrix from LinProg");
347 for (i = 1; i <=
nrows; i++)
350 if (i == 1) printf(
" ");
351 else if (iposv[i-1] <=
N) printf(
"X%d", iposv[i-1]);
352 else printf(
"Y%d", iposv[i-1]-
N+1);
353 for (j = 1; j <=
ncols; j++) printf(
" %7.2f ",(
double)a[i][j]);
359 void print_exp(
const onePointP vert,
int n )
362 for ( i= 1; i <=
n; i++ )
364 Print(
" %d",vert->point[i] );
366 if ( i < n )
PrintS(
", ");
370 void print_matrix(
matrix omat )
375 for ( i= 1; i <=
MATROWS( omat ); i++ )
377 for ( j= 1; j <=
MATCOLS( omat ); j++ )
416 points = (onePointP *)
omAlloc( (count+1) *
sizeof(onePointP) );
417 for ( i= 0; i <=
max; i++ )
429 for ( i= 0; i <=
max; i++ )
450 (
max+1) *
sizeof(onePointP),
451 (2*
max + 1) *
sizeof(onePointP) );
452 for ( i=
max+1; i <=
max*2; i++ )
471 for ( i= 1; i <=
dim; i++ )
points[
num]->point[i]= vert->point[i];
493 for ( i= 0; i <
dim; i++ )
points[
num]->point[i+1]= vert[i];
516 for ( i= 1; i <=
num; i++ )
518 for ( j= 1; j <=
dim; j++ )
519 if (
points[i]->point[j] != vert->point[j] )
break;
520 if ( j >
dim )
break;
535 for ( i= 1; i <=
num; i++ )
537 for ( j= 1; j <=
dim; j++ )
539 if ( j >
dim )
break;
555 vert= (
int *)
omAlloc( (
dim+1) *
sizeof(int) );
561 for ( i= 1; i <=
num; i++ )
563 for ( j= 1; j <=
dim; j++ )
565 if ( j >
dim )
break;
584 vert= (
int *)
omAlloc( (
dim+1) *
sizeof(int) );
587 for ( i= 1; i <=
num; i++ )
589 for ( j= 1; j <=
dim; j++ )
591 if ( j >
dim )
break;
595 if ( i >
num )
return 0;
605 for ( i= 1; i <=
dim; i++ )
606 vert[i]= (
int)(
points[indx]->point[
i] -
points[indx]->rcPnt->point[
i]);
613 for ( i= 1; i <=
dim; i++ )
632 for ( i= 1; i <=
dim; i++ )
656 for ( i= 1; i <
num; i++ )
683 for(i = 1; i <
dim; i++)
688 for ( j=1; j <=
num; j++ )
691 for ( i=1; i <
dim; i++ )
693 sum += (int)
points[j]->point[i] * l[i];
700 for ( j=1; j <
dim; j++ )
Print(
" %d ",l[j] );
703 PrintS(
" lifted points: \n");
704 for ( j=1; j <=
num; j++ )
714 if ( !outerL )
omFreeSize( (
void *) l, (dim+1) *
sizeof(
int) );
737 pLP->LiPM[1][1] = +0.0;
738 pLP->LiPM[1][2] = +1.0;
739 pLP->LiPM[2][1] = +1.0;
740 pLP->LiPM[2][2] = -1.0;
742 for ( j=3; j <= pLP->n; j++)
744 pLP->LiPM[1][
j] = +0.0;
745 pLP->LiPM[2][
j] = -1.0;
748 for( i= 1; i <= n; i++) {
751 for( j= 1; j <=
m; j++ )
762 PrintS(
"Matrix of Linear Programming\n");
763 print_mat( pLP->LiPM, pLP->m+1,pLP->n);
770 return (pLP->icase == 0);
784 vert= (
int *)
omAlloc( (idelem+1) *
sizeof(int) );
787 for ( i= 0; i < idelem; i++ )
790 for( i= 0; i < idelem; i++ )
796 for( j= 1; j <=
m; j++) {
797 if( !inHull( (gls->m)[i], p, m, j ) )
800 Q[
i]->addPoint( vert );
813 omFreeSize( (
void *) vert, (idelem+1) *
sizeof(
int) );
817 for( i= 0; i < idelem; i++ )
819 Print(
" \\Conv(Qi[%d]): #%d\n", i,
Q[i]->
num );
820 for ( j=1; j <=
Q[
i]->num; j++ )
844 vert= (
int *)
omAlloc( (idelem+1) *
sizeof(int) );
847 for( i= 0; i < idelem; i++ )
852 for( j= 1; j <=
m; j++) {
853 if( !inHull( (gls->m)[i], p, m, j ) )
855 if ( (id->m)[
i] ==
NULL )
857 (
id->m)[i]=
pHead(p);
877 omFreeSize( (
void *) vert, (idelem+1) *
sizeof(
int) );
881 for( i= 0; i < idelem; i++ )
900 for ( i= 0; i <
MAXVARS+2; i++ ) acoords[i]= 0;
911 int i, ii,
j,
k, col, r;
917 numverts += Qi[
i]->num;
924 pLP->LiPM[1][1] = 0.0;
925 pLP->LiPM[1][2] = 1.0;
926 for( j=3; j<=cols; j++) pLP->LiPM[1][j] = 0.0;
928 for( i=0; i <= n; i++ )
930 pLP->LiPM[i+2][1] = 1.0;
931 pLP->LiPM[i+2][2] = 0.0;
933 for( i=1; i<=
dim; i++)
935 pLP->LiPM[n+2+
i][1] = (
mprfloat)(acoords_a[i-1]);
936 pLP->LiPM[n+2+
i][2] = -shift[
i];
941 for ( i= 0; i <= n; i++ )
944 for( k= 1; k <= Qi[ii]->num; k++ )
947 for ( r= 0; r <= n; r++ )
949 if ( r == i ) pLP->LiPM[r+2][col] = -1.0;
950 else pLP->LiPM[r+2][col] = 0.0;
952 for( r= 1; r <=
dim; r++ )
953 pLP->LiPM[r+n+2][col] = -(
mprfloat)((*Qi[ii])[k]->point[r]);
958 Werror(
"mayanPyramidAlg::vDistance:" 959 "setting up matrix for udist: col %d != cols %d",col,cols);
966 Print(
"vDistance LP, known koords dim=%d, constr %d, cols %d, acoords= ",
968 for( i= 0; i <
dim; i++ )
969 Print(
" %d",acoords_a[i]);
971 print_mat( pLP->LiPM, pLP->m+1, cols);
977 PrintS(
"LP returns matrix\n");
978 print_bmat( pLP->LiPM, pLP->m+1, cols+1-pLP->m, cols, pLP->iposv);
981 if( pLP->icase != 0 )
983 WerrorS(
"mayanPyramidAlg::vDistance:");
984 if( pLP->icase == 1 )
985 WerrorS(
" Unbounded v-distance: probably 1st v-coor=0");
986 else if( pLP->icase == -1 )
987 WerrorS(
" Infeasible v-distance");
993 return pLP->LiPM[1][1];
998 int i,
j,
k, cols, cons;
1007 pLP->LiPM[1][1] = 0.0;
1008 for( i=2; i<=n+2; i++)
1010 pLP->LiPM[
i][1] = 1.0;
1011 pLP->LiPM[
i][2] = 0.0;
1016 for( i=0; i<=n; i++)
1019 for( j=1; j<= Qi[
i]->num; j++)
1022 pLP->LiPM[1][cols] = 0.0;
1023 for( k=2; k<=n+2; k++)
1025 if( k != la_cons_row) pLP->LiPM[
k][cols] = 0.0;
1026 else pLP->LiPM[
k][cols] = -1.0;
1028 for( k=1; k<=n; k++)
1029 pLP->LiPM[k+n+2][cols] = -(
mprfloat)((*Qi[
i])[j]->point[k]);
1033 for( i= 0; i <
dim; i++ )
1035 pLP->LiPM[i+n+3][1] = acoords[
i];
1036 pLP->LiPM[i+n+3][2] = 0.0;
1038 pLP->LiPM[dim+n+3][1] = 0.0;
1041 pLP->LiPM[1][2] = -1.0;
1042 pLP->LiPM[dim+n+3][2] = 1.0;
1045 Print(
"\nThats the matrix for minR, dim= %d, acoords= ",dim);
1046 for( i= 0; i <
dim; i++ )
1047 Print(
" %d",acoords[i]);
1049 print_mat( pLP->LiPM, cons+1, cols);
1059 if ( pLP->icase != 0 )
1062 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: minR: infeasible");
1063 else if( pLP->icase > 0)
1064 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: minR: unbounded");
1073 pLP->LiPM[1][1] = 0.0;
1074 for( i=2; i<=n+2; i++)
1076 pLP->LiPM[
i][1] = 1.0;
1077 pLP->LiPM[
i][2] = 0.0;
1081 for( i=0; i<=n; i++)
1084 for( j=1; j<=Qi[
i]->num; j++)
1087 pLP->LiPM[1][cols] = 0.0;
1088 for( k=2; k<=n+2; k++)
1090 if( k != la_cons_row) pLP->LiPM[
k][cols] = 0.0;
1091 else pLP->LiPM[
k][cols] = -1.0;
1093 for( k=1; k<=n; k++)
1094 pLP->LiPM[k+n+2][cols] = -(
mprfloat)((*Qi[
i])[j]->point[k]);
1098 for( i= 0; i <
dim; i++ )
1100 pLP->LiPM[i+n+3][1] = acoords[
i];
1101 pLP->LiPM[i+n+3][2] = 0.0;
1103 pLP->LiPM[dim+n+3][1] = 0.0;
1105 pLP->LiPM[1][2] = 1.0;
1106 pLP->LiPM[dim+n+3][2] = 1.0;
1109 Print(
"\nThats the matrix for maxR, dim= %d\n",dim);
1110 print_mat( pLP->LiPM, cons+1, cols);
1120 if ( pLP->icase != 0 )
1123 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: maxR: infeasible");
1124 else if( pLP->icase > 0)
1125 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: maxR: unbounded");
1131 Print(
" Range for dim=%d: [%d,%d]\n", dim, *minR, *maxR);
1143 dist= vDistance( &(acoords[0]), n );
1152 E->addPoint( &(acoords[0]) );
1168 mn_mx_MinkowskiSum( dim, &minR, &maxR );
1172 for( i=0; i <=
dim; i++)
Print(
"acoords[%d]=%d ",i,(
int)acoords[i]);
1173 Print(
":: [%d,%d]\n", minR, maxR);
1181 acoords[
dim] = minR;
1182 while( acoords[dim] <= maxR )
1184 if( !storeMinkowskiSumPoint() )
1193 acoords[
dim] = minR;
1194 while ( acoords[dim] <= maxR )
1196 if ( (acoords[dim] > minR) && (acoords[dim] <= maxR) )
1199 runMayanPyramid( dim + 1 );
1204 dist= vDistance( &(acoords[0]), dim + 1 );
1209 runMayanPyramid( dim + 1 );
1222 for ( i= 0; i <= nn; i++ )
1224 if ( (loffset < indx) && (indx <= pQ[
i]->
num + loffset) )
1230 else loffset+= pQ[
i]->
num;
1251 for ( i= 0; i <= n; i++ )
1254 for ( k= 1; k <=
size; k++ )
1262 for ( j = 0; j <= n; j++ )
1265 LP->LiPM[j+2][LP->n] = -1.0;
1267 LP->LiPM[j+2][LP->n] = 0.0;
1271 for ( j = 1; j <= n; j++ )
1273 LP->LiPM[j+n+2][LP->n] = - ( (
mprfloat) (*pQ[i])[
k]->point[
j] );
1278 for ( j = 0; j <= n; j++ ) LP->LiPM[j+2][1] = 1.0;
1279 for ( j= 1; j <= n; j++ )
1281 LP->LiPM[j+n+2][1]= (
mprfloat)(*E)[vert]->point[
j] - shift[
j];
1285 LP->LiPM[1][1] = 0.0;
1289 Print(
" n= %d, LP->m=M= %d, LP->n=N= %d\n",n,LP->m,LP->n);
1290 print_mat(LP->LiPM, LP->m+1, LP->n+1);
1297 if ( LP->icase < 0 )
1304 (*E)[vert]->point[E->
dim]= (int)(-LP->LiPM[1][1] *
SCALEDOWN);
1307 Print(
" simplx returned %d, Objective value = %f\n", LP->icase, LP->LiPM[1][1]);
1317 for ( i= 1; i < LP->m; i++ )
1319 if ( LP->iposv[i] > LP->iposv[i+1] )
1323 LP->iposv[
i]=LP->iposv[i+1];
1326 cd=LP->LiPM[i+1][1];
1327 LP->LiPM[i+1][1]=LP->LiPM[i+2][1];
1328 LP->LiPM[i+2][1]=
cd;
1337 print_bmat(LP->LiPM, LP->m + 1, LP->n+1-LP->m, LP->n+1, LP->iposv);
1338 PrintS(
" now split into sets\n");
1343 for ( i= 0; i <= E->
dim; i++ ) bucket[i]= 0;
1347 for ( i= 0; i < LP->m; i++ )
1350 if ( LP->LiPM[i+2][1] > 1e-12 )
1352 if ( !remapXiToPoint( LP->iposv[i+1], pQ, &(optSum[c].
set), &(optSum[c].
pnt) ) )
1354 Werror(
" resMatrixSparse::RC: Found bad solution in LP: %d!",LP->iposv[i+1]);
1355 WerrorS(
" resMatrixSparse::RC: remapXiToPoint failed!");
1358 bucket[optSum[c].
set]++;
1366 for ( i= 1; i < E->
dim; i++ )
1368 if ( bucket[c] >= bucket[i] )
1374 for ( i= onum - 1; i >= 0; i-- )
1376 if ( optSum[i].
set == c )
1380 (*E)[vert]->rc.set= c;
1381 (*E)[vert]->rc.pnt= optSum[
i].
pnt;
1382 (*E)[vert]->rcPnt= (*pQ[c])[optSum[i].
pnt];
1384 if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
1386 #ifdef mprDEBUG_PROT 1387 Print(
"\n Point E[%d] was <",vert);print_exp((*E)[vert],E->
dim-1);
Print(
">, bucket={");
1388 for ( j= 0; j < E->
dim; j++ )
1390 Print(
" %d",bucket[j]);
1392 PrintS(
" }\n optimal Sum: Qi ");
1393 for ( j= 0; j < LP->m; j++ )
1395 Print(
" [ %d, %d ]",optSum[j].
set,optSum[j].
pnt);
1397 Print(
" -> i= %d, j = %d\n",(*E)[vert]->rc.set,optSum[i].
pnt);
1405 return (
int)(-LP->LiPM[1][1] *
SCALEDOWN);
1420 int *epp_mon, *eexp;
1422 epp_mon= (
int *)
omAlloc( (n+2) *
sizeof(int) );
1440 for ( i= 1; i <= E->
num; i++ )
1446 rowp=
ppMult_qq( epp, (gls->m)[(*E)[i]->rc.set] );
1451 while ( iterp!=
NULL )
1458 Werror(
"resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
1459 i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
1465 if ( (*E)[i]->rc.set == linPolyS )
1472 if ( (*E)[i]->rc.set == linPolyS )
1477 (rmat->m)[i-1]= rowp;
1481 omFreeSize( (
void *) epp_mon, (n+2) *
sizeof(
int) );
1490 for ( i= 1; i <= numSet0; i++ )
1492 Print(
" row %d contains coeffs of f_%d\n",
IMATELEM(*uRPos,i,1),linPolyS);
1494 PrintS(
" Sparse Matrix done\n");
1510 for ( j= 1; j < i-1; j++ )
1531 for ( j= 1; j <= Q1->
num; j++ )
1533 for ( k= 1; k <= Q2->
num; k++ )
1535 for ( l= 1; l <=
dim; l++ )
1537 vert.
point[
l]= (*Q1)[
j]->point[
l] + (*Q2)[
k]->point[
l];
1556 for ( j= 1; j <= pQ[0]->
num; j++ ) vs->
addPoint( (*pQ[0])[j] );
1558 for ( j= 1; j < numq; j++ )
1561 vs= minkSumTwo( vs_old, pQ[j], dim );
1583 WerrorS(
"resMatrixSparse::resMatrixSparse: Too many variables!");
1602 LP =
new simplex( idelem+totverts*2+5, totverts+5 );
1606 shift[0]=0.005; shift[1]=0.003; shift[2]=0.008; shift[3]=0.005; shift[4]=0.002;
1607 shift[5]=0.1; shift[6]=0.3; shift[7]=0.2; shift[8]=0.4; shift[9]=0.2;
1611 #ifdef mprDEBUG_PROT 1612 PrintS(
" shift vector: ");
1613 for ( i= 1; i <=
idelem; i++ )
Print(
" %.12f ",(
double)shift[i]);
1629 #ifdef mprDEBUG_PROT 1633 PrintS(
"\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
1634 for ( pnt= 1; pnt <= E->
num; pnt++ )
1643 lift[0][1]=3; lift[0][2]=4; lift[0][3]=8; lift[0][4]=2;
1644 lift[1][1]=6; lift[1][2]=1; lift[1][3]=7; lift[1][4]=4;
1645 lift[2][1]=2; lift[2][2]=5; lift[2][3]=9; lift[2][4]=6;
1646 lift[3][1]=2; lift[3][2]=1; lift[3][3]=9; lift[3][4]=5;
1647 lift[4][1]=3; lift[4][2]=7; lift[4][3]=1; lift[4][4]=5;
1649 for ( i= 0; i <=
n; i++ ) Qi[i]->
lift( lift[i] );
1652 for ( i= 0; i <=
n; i++ ) Qi[i]->
lift();
1657 for ( pnt= 1; pnt <= E->
num; pnt++ )
1659 RC( Qi, E, pnt, shift );
1664 for ( pnt= k; pnt > 0; pnt-- )
1666 if ( (*E)[pnt]->rcPnt ==
NULL )
1674 #ifdef mprDEBUG_PROT 1675 PrintS(
" points which lie in a cell:\n");
1676 for ( pnt= 1; pnt <= E->
num; pnt++ )
1684 for ( i= 0; i <=
n; i++ ) Qi[i]->unlift();
1688 #ifdef mprDEBUG_PROT 1689 Print(
" points with a[ij] (%d):\n",E->
num);
1690 for ( pnt= 1; pnt <= E->
num; pnt++ )
1692 Print(
"Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->
dim );
1693 Print(
">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
1695 print_exp( (*E)[pnt]->rcPnt, E->
dim );
PrintS(
">\n");
1702 WerrorS(
"could not handle a degenerate situation: no inner points found");
1709 WerrorS(
"resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
1715 for ( i= 0; i <
idelem; i++ )
1743 for ( i= 1; i <=
numSet0; i++ )
1802 for ( i= 1; i <=
numSet0; i++ )
1810 for ( cp= 2; cp <=
idelem; cp++ )
1812 if ( !
nIsZero(evpoint[cp-1]) )
1862 for ( i= 1; i <=
numSet0; i++ )
1868 for ( cp= 2; cp <=
idelem; cp++ )
1870 if ( !
nIsZero(evpoint[cp-1]) )
1953 number
getDetAt(
const number* evpoint );
1968 void generateBaseData();
1973 void generateMonomData(
int deg,
intvec* polyDegs ,
intvec* iVO );
1978 void generateMonoms( poly
m,
int var,
int deg );
2017 poly getElem(
const int i );
2020 number getElemNum(
const int i );
2048 assume( 0 < i || i > numColVectorSize );
2057 assume( i >= 0 && i < numColVectorSize );
2058 return( numColVector[i] );
2101 numVectors *
sizeof( number ) );
2126 for ( i= 1; i <=
MATROWS(
m ); i++ )
2127 for ( j= 1; j <=
MATCOLS(
m ); j++ )
2139 for ( i= 0; i < (
currRing->N); i++ )
2165 for ( i=0; i < (
currRing->N); i++ )
2190 poly mon =
pCopy( mm );
2211 if ( var == (
currRing->N)+1 )
return;
2212 poly newm =
pCopy( mm );
2249 ideal pDegDiv=
idInit( polyDegs->
rows(), 1 );
2250 for ( k= 0; k < polyDegs->
rows(); k++ )
2253 pSetExp( p, k + 1, (*polyDegs)[k] );
2265 for ( k= 0; k <
IDELEMS(pDegDiv); k++ )
2276 for ( k= 0; k < iVO->
rows(); k++)
2289 for ( i= 0; i <
k; i++ )
2314 for ( i= 0; i < polyDegs->
rows(); i++ )
2317 for ( k= 0; k < polyDegs->
rows(); k++ )
2318 if ( i != k ) sub*= (*polyDegs)[
k];
2331 Print(
"// %s, S(%d), db ",
2361 for ( k= (
currRing->N) - 1; k >= 0; k-- )
2373 for ( k= 0; k < (
currRing->N); k++ )
2379 for ( k= 0; k < polyDegs.
rows(); k++ )
2380 sumDeg+= polyDegs[k];
2381 sumDeg-= polyDegs.
rows() - 1;
2403 while ( pi !=
NULL )
2435 while ( pi !=
NULL )
2493 for (j=1; j <= (
currRing->N); j++ )
2495 if (
MATELEM(resmat,numVectors-i,
2559 for ( i= 0; i < (
currRing->N); i++ )
2564 nCopy(evpoint[i]) );
2600 for ( i= 1; i <=
MATROWS( mat ); i++ )
2602 for ( j= 1; j <=
MATCOLS( mat ); j++ )
2652 #define MAXEVPOINT 1000000 // 0x7fffffff 2658 unsigned long over(
const unsigned long n ,
const unsigned long d )
2663 mpz_init(m);mpz_set_ui(m,1);
2664 mpz_init(md);mpz_set_ui(md,1);
2665 mpz_init(mn);mpz_set_ui(mn,1);
2672 mpz_tdiv_q(res,m,res);
2674 mpz_clear(m);mpz_clear(md);mpz_clear(mn);
2676 unsigned long result = mpz_get_ui(res);
2705 WerrorS(
"uResultant::uResultant: Unknown chosen resultant matrix type!");
2716 ideal newGls=
idCopy( igls );
2719 (
IDELEMS(igls) + 1) *
sizeof(poly) );
2728 for ( i=
IDELEMS(newGls)-1; i > 0; i-- )
2730 newGls->m[
i]= newGls->m[i-1];
2732 newGls->m[0]= linPoly;
2736 WerrorS(
"uResultant::extendIdeal: Unknown chosen resultant matrix type!");
2747 poly actlp, rootlp= newlp;
2749 for ( i= 1; i <= (
currRing->N); i++ )
2779 long mdg=
over(
n-1, tdg );
2782 long l=(long)
pow( (
double)(tdg+1),
n );
2784 #ifdef mprDEBUG_PROT 2785 Print(
"// total deg of D: tdg %ld\n",tdg);
2786 Print(
"// maximum number of terms in D: mdg: %ld\n",mdg);
2787 Print(
"// maximum number of terms in polynom of deg tdg: l %ld\n",l);
2792 presults= (number *)
omAlloc( mdg *
sizeof( number ) );
2793 for (i=0; i < mdg; i++) presults[i]=
nInit(0);
2795 number *pevpoint= (number *)
omAlloc(
n *
sizeof( number ) );
2796 number *pev= (number *)
omAlloc(
n *
sizeof( number ) );
2797 for (i=0; i <
n; i++) pev[i]=
nInit(0);
2799 mprPROTnl(
"// initial evaluation point: ");
2802 for (i=0; i <
n; i++)
2813 for ( i=0; i < mdg; i++ )
2815 for (j=0; j <
n; j++)
2818 nPower(pevpoint[j],i,&pev[j]);
2838 if ( subDetVal !=
NULL )
2841 for ( i= 0; i <= mdg; i++ )
2843 detdiv=
nDiv( ncpoly[i], subDetVal );
2852 for ( i=0; i < mdg; i++ )
2861 for ( i=0; i < mdg; i++ )
2863 if (
nEqual(ncpoly[i],nn) )
2873 for ( i= 0; i <
n; i++ ) exp[i]=0;
2880 for ( i=0; i <
l; i++ )
2889 for ( j= 0; j <
n; j++ )
pSetExp( p, j+1, exp[j] );
2893 for ( j= 1; j <
n; j++ )
pSetExp( p, j, exp[j] );
2897 if (result!=
NULL) result=
pAdd( result, p );
2904 for ( j= 0; j < n - 1; j++ )
2925 int loops= (matchUp?
n-2:
n-1);
2927 mprPROTnl(
"uResultant::interpolateDenseSP");
2935 presults= (number *)
omAlloc( (tdg + 1) *
sizeof( number ) );
2936 for ( i=0; i <= tdg; i++ ) presults[i]=
nInit(0);
2942 number *pevpoint= (number *)
omAlloc(
n *
sizeof( number ) );
2943 for (i=0; i <
n; i++) pevpoint[i]=
nInit(0);
2945 number *pev= (number *)
omAlloc( n *
sizeof( number ) );
2946 for (i=0; i <
n; i++) pev[i]=
nInit(0);
2952 for ( uvar= 0; uvar < loops; uvar++ )
2957 for (i=0; i <
n; i++)
2966 else if ( i <= uvar + 2 )
2978 for (i=0; i <
n; i++)
2989 if ( i == (uvar + 1) ) pevpoint[
i]=
nInit(-1);
2990 else pevpoint[
i]=
nInit(0);
2997 for (i=0; i <
n; i++)
3000 pev[
i]=
nCopy( pevpoint[i] );
3003 for ( i=0; i <= tdg; i++ )
3006 nPower(pevpoint[0],i,&pev[0]);
3022 if ( subDetVal !=
NULL )
3025 for ( i= 0; i <= tdg; i++ )
3027 detdiv=
nDiv( ncpoly[i], subDetVal );
3036 for ( i=0; i <= tdg; i++ )
3050 for ( i=0; i <
n; i++ )
nDelete( pev + i );
3051 omFreeSize( (
void *)pev, n *
sizeof( number ) );
3053 for ( i=0; i <= tdg; i++ )
nDelete( presults+i );
3054 omFreeSize( (
void *)presults, (tdg + 1) *
sizeof( number ) );
3064 int loops=(matchUp?
n-2:
n-1);
3066 if (loops==0) { loops=1;nn++;}
3076 number *pevpoint= (number *)
omAlloc( nn *
sizeof( number ) );
3077 for (i=0; i < nn; i++) pevpoint[i]=
nInit(0);
3082 for ( uvar= 0; uvar < loops; uvar++ )
3087 for (i=0; i <
n; i++)
3091 if ( i <= uvar + 2 )
3096 else pevpoint[
i]=
nInit(0);
3102 for (i=0; i <
n; i++)
3106 if ( i == (uvar + 1) ) pevpoint[
i]=
nInit(-1);
3107 else pevpoint[
i]=
nInit(0);
3114 number *ncpoly= (number *)
omAlloc( (tdg+1) *
sizeof(number) );
3121 for ( i= tdg; i >= 0; i-- )
3143 if ( subDetVal !=
NULL )
3146 for ( i= 0; i <= tdg; i++ )
3148 detdiv=
nDiv( ncpoly[i], subDetVal );
3166 for ( i=0; i <
n; i++ )
nDelete( pevpoint + i );
3167 omFreeSize( (
void *)pevpoint, n *
sizeof( number ) );
3194 int totverts,idelem;
3201 for( i=0; i < idelem; i++) totverts +=
pLength( (id->m)[i] );
3203 LP =
new simplex( idelem+totverts*2+5, totverts+5 );
int status int void size_t count
#define pSetmComp(p)
TODO:
#define mprSTICKYPROT(msg)
complex root finder for univariate polynomials based on laguers algorithm
pointSet(const int _dim, const int _index=0, const int count=MAXINITELEMS)
int RC(pointSet **pQ, pointSet *E, int vert, mprfloat shift[])
Row Content Function Finds the largest i such that F[i] is a point, F[i]= a[ij] in A[i] for some j...
long isReduced(const nmod_mat_t M)
void randomVector(const int dim, mprfloat shift[])
number * interpolateDense(const number *q)
Solves the Vandermode linear system {i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
vandermonde system solver for interpolating polynomials from their values
CanonicalForm cd(bCommonDen(FF))
#define idDelete(H)
delete an ideal
gmp_float exp(const gmp_float &a)
Linear Programming / Linear Optimization using Simplex - Algorithm.
Compatiblity layer for legacy polynomial operations (over currRing)
#define nPower(a, b, res)
#define mprPROTL(msg, intval)
poly sm_CallDet(ideal I, const ring R)
bool larger(int, int)
points[a] > points[b] ?
#define mprSTICKYPROT2(msg, arg)
static void p_GetExpV(poly p, int *ev, const ring r)
#define omFreeSize(addr, size)
bool checkMem()
Checks, if more mem is needed ( i.e.
#define omfreeSize(addr, size)
poly linearPoly(const resMatType rmt)
virtual poly getUDet(const number *)
void getRowMP(const int indx, int *vert)
void mn_mx_MinkowskiSum(int dim, Coord_t *minR, Coord_t *maxR)
LP for finding min/max coord in MinkowskiSum, given previous coors.
poly singclap_det(const matrix m, const ring s)
mprfloat vDistance(Coord_t *acoords, int dim)
Compute v-distance via Linear Programing Linear Program finds the v-distance of the point in accords[...
void WerrorS(const char *s)
ideal loNewtonPolytope(const ideal id)
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Base class for sparse and dense u-Resultant computation.
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
resMatrixDense(const ideal _gls, const int special=SNONE)
_gls: system of multivariate polynoms special: -1 -> resMatrixDense is a symbolic matrix 0...
resVector * resVectorList
virtual number getDetAt(const number *)
#define nPrint(a)
only for debug, over any initalized currRing
static FORCE_INLINE number n_Param(const int iParameter, const coeffs r)
return the (iParameter^th) parameter as a NEW number NOTE: parameter numbering: 1..n_NumberOfParameters(...)
poly getElem(const int i)
index von 0 ...
ideal lift(const ideal J, const ring r, const ideal inI, const ring s)
int getExpPos(const poly p)
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
#define omReallocSize(addr, o_size, size)
#define pGetExp(p, i)
Exponent.
number * numColVector
holds the column vector if (elementOfS != linPolyS)
poly monomAt(poly p, int i)
number getElemNum(const int i)
index von 0 ...
ideal newtonPolytopesI(const ideal gls)
static FORCE_INLINE long n_Int(number &n, const coeffs r)
conversion of n to an int; 0 if not possible in Z/pZ: the representing int lying in (-p/2 ...
const CanonicalForm CFMap CFMap & N
static int max(int a, int b)
pointSet * getInnerPoints(pointSet **_q_i, mprfloat _shift[])
Drive Mayan Pyramid Algorithm.
#define pLmDivisibleByNoComp(a, b)
like pLmDivisibleBy, does not check components
static long pTotaldegree(poly p)
for(int i=0;i<=n;i++) degsf[i]
#define pLmInit(p)
like pInit, except that expvector is initialized to that of p, p must be != NULL
void lift(int *l=NULL)
Lifts the point set using sufficiently generic linear lifting homogeneous forms l[1]..l[dim] in Z.
void mergeWithPoly(const poly p)
convexHull(simplex *_pLP)
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
int nextPrime(const int p)
poly getUDet(const number *evpoint)
void PrintS(const char *s)
mayanPyramidAlg(simplex *_pLP)
int elementOfS
number of the set S mon is element of
virtual ideal getSubMatrix()
static unsigned pLength(poly a)
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL ...
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
pointSet * minkSumTwo(pointSet *Q1, pointSet *Q2, int dim)
static int index(p_Length length, p_Ord ord)
int numColVectorSize
size of numColVector
matrix mpNew(int r, int c)
create a r x c zero-matrix
ideal idInit(int idsize, int rank)
initialise an ideal / module
onePointP operator[](const int index)
REvaluation E(1, terms.length(), IntRandom(25))
bool addPoint(const onePointP vert)
Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
bool removePoint(const int indx)
resMatrixSparse(const ideal _gls, const int special=SNONE)
void generateBaseData()
Generate the "matrix" M.
pointSet * minkSumAll(pointSet **pQ, int numq, int dim)
ideal getSubMatrix()
Returns the submatrix M' of M in an usable presentation.
bool smaller(int, int)
points[a] < points[b] ?
void createMatrix()
Creates quadratic matrix M of size numVectors for later use.
bool remapXiToPoint(const int indx, pointSet **pQ, int *set, int *vtx)
#define pInit()
allocates a new monomial and initializes everything to 0
poly interpolateDense(const number subDetVal=NULL)
number getSubDet()
Evaluates the determinant of the submatrix M'.
bool mergeWithExp(const onePointP vert)
Adds point to pointSet, iff pointSet point = .
#define mprPROTN(msg, nval)
number getDetAt(const number *evpoint)
Fills in resMat[][] with evpoint[] and gets determinant uRPos[i][1]: row of matrix uRPos[i][idelem+1]...
ideal getMatrix()
Returns the matrix M in an usable presentation.
void sort(CFArray &A, int l=0)
quick sort A
pointSet ** newtonPolytopesP(const ideal gls)
Computes the point sets of the convex hulls of the supports given by the polynoms in gls...
void generateMonomData(int deg, intvec *polyDegs, intvec *iVO)
Generates needed set of monoms, split them into sets S0, ...
int * numColParNr
holds the index of u0, u1, ..., un, if (elementOfS == linPolyS) the size is given by (currRing->N) ...
int createMatrix(pointSet *E)
create coeff matrix uRPos[i][1]: row of matrix uRPos[i][idelem+1]: col of u(0) uRPos[i][2..idelem]: col of u(1) .
uResultant(const ideal _gls, const resMatType _rmt=sparseResMat, BOOLEAN extIdeal=true)
Rational pow(const Rational &a, int e)
#define IMATELEM(M, I, J)
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
unsigned long over(const unsigned long n, const unsigned long d)
ideal id_Matrix2Module(matrix mat, const ring R)
converts mat to module, destroys mat
resVector * getMVector(const int i)
column vector of matrix, index von 0 ...
void Werror(const char *fmt,...)
#define mprPROTNnl(msg, nval)
ideal extendIdeal(const ideal gls, poly linPoly, const resMatType rmt)
virtual number getSubDet()
void generateMonoms(poly m, int var, int deg)
Recursively generate all homogeneous monoms of (currRing->N) variables of degree deg.
bool inHull(poly p, poly pointPoly, int m, int site)
Returns true iff the support of poly pointPoly is inside the convex hull of all points given by the s...
number getDetAt(const number *evpoint)
Evaluate the determinant of the matrix M at the point evpoint where the ui's are replaced by the comp...
void runMayanPyramid(int dim)
Recursive Mayan Pyramid algorithm for directly computing MinkowskiSum lattice points for (n+1)-fold M...
bool storeMinkowskiSumPoint()
Stores point in E->points[pt], iff v-distance != 0 Returns true iff point was stored, else flase.
#define pCopy(p)
return a copy of the poly
#define MATELEM(mat, i, j)
1-based access to matrix