50 #define MAXPOINTS 10000 51 #define MAXINITELEMS 256 52 #define LIFT_COOR 50000 // siRand() % LIFT_COOR gives random lift value 53 #define SCALEDOWN 100.0 // lift value scale down for linear program 55 #define RVMULT 0.0001 // multiplicator for random shift vector 56 #define MAXRVVAL 50000 83 number
getDetAt(
const number* evpoint );
85 poly
getUDet(
const number* evpoint );
148 typedef struct onePoint * onePointP;
158 typedef struct _entry * entry;
178 inline onePointP operator[] (
const int index );
184 bool addPoint(
const onePointP vert );
190 bool addPoint(
const int * vert );
196 bool addPoint(
const Coord_t * vert );
199 bool removePoint(
const int indx );
204 bool mergeWithExp(
const onePointP vert );
209 bool mergeWithExp(
const int * vert );
212 void mergeWithPoly(
const poly
p );
215 void getRowMP(
const int indx,
int * vert );
218 int getExpPos(
const poly
p );
237 inline bool smaller(
int,
int );
240 inline bool larger(
int,
int );
245 inline bool checkMem();
262 ideal newtonPolytopesI(
const ideal gls );
268 bool inHull(poly
p, poly pointPoly,
int m,
int site);
297 void runMayanPyramid(
int dim );
316 bool storeMinkowskiSumPoint();
332 #if defined(mprDEBUG_PROT) || defined(mprDEBUG_ALL) 333 void print_mat(
mprfloat **a,
int maxrow,
int maxcol)
337 for (i = 1; i <= maxrow; i++)
340 for (j = 1; j <= maxcol; j++)
Print(
"% 7.2f, ", a[i][j]);
348 printf(
"Output matrix from LinProg");
349 for (i = 1; i <=
nrows; i++)
352 if (i == 1) printf(
" ");
353 else if (iposv[i-1] <=
N) printf(
"X%d", iposv[i-1]);
354 else printf(
"Y%d", iposv[i-1]-
N+1);
355 for (j = 1; j <=
ncols; j++) printf(
" %7.2f ",(
double)a[i][j]);
361 void print_exp(
const onePointP vert,
int n )
364 for ( i= 1; i <=
n; i++ )
366 Print(
" %d",vert->point[i] );
368 if ( i < n )
PrintS(
", ");
372 void print_matrix(
matrix omat )
377 for ( i= 1; i <=
MATROWS( omat ); i++ )
379 for ( j= 1; j <=
MATCOLS( omat ); j++ )
418 points = (onePointP *)
omAlloc( (count+1) *
sizeof(onePointP) );
419 for ( i= 0; i <=
max; i++ )
431 for ( i= 0; i <=
max; i++ )
452 (
max+1) *
sizeof(onePointP),
453 (2*
max + 1) *
sizeof(onePointP) );
454 for ( i=
max+1; i <=
max*2; i++ )
473 for ( i= 1; i <=
dim; i++ )
points[
num]->point[i]= vert->point[i];
495 for ( i= 0; i <
dim; i++ )
points[
num]->point[i+1]= vert[i];
518 for ( i= 1; i <=
num; i++ )
520 for ( j= 1; j <=
dim; j++ )
521 if (
points[i]->point[j] != vert->point[j] )
break;
522 if ( j >
dim )
break;
537 for ( i= 1; i <=
num; i++ )
539 for ( j= 1; j <=
dim; j++ )
541 if ( j >
dim )
break;
557 vert= (
int *)
omAlloc( (
dim+1) *
sizeof(int) );
563 for ( i= 1; i <=
num; i++ )
565 for ( j= 1; j <=
dim; j++ )
567 if ( j >
dim )
break;
586 vert= (
int *)
omAlloc( (
dim+1) *
sizeof(int) );
589 for ( i= 1; i <=
num; i++ )
591 for ( j= 1; j <=
dim; j++ )
593 if ( j >
dim )
break;
597 if ( i >
num )
return 0;
607 for ( i= 1; i <=
dim; i++ )
608 vert[i]= (
int)(
points[indx]->point[
i] -
points[indx]->rcPnt->point[
i]);
615 for ( i= 1; i <=
dim; i++ )
634 for ( i= 1; i <=
dim; i++ )
658 for ( i= 1; i <
num; i++ )
685 for(i = 1; i <
dim; i++)
690 for ( j=1; j <=
num; j++ )
693 for ( i=1; i <
dim; i++ )
695 sum += (int)
points[j]->point[i] * l[i];
702 for ( j=1; j <
dim; j++ )
Print(
" %d ",l[j] );
705 PrintS(
" lifted points: \n");
706 for ( j=1; j <=
num; j++ )
716 if ( !outerL )
omFreeSize( (
void *) l, (dim+1) *
sizeof(
int) );
739 pLP->LiPM[1][1] = +0.0;
740 pLP->LiPM[1][2] = +1.0;
741 pLP->LiPM[2][1] = +1.0;
742 pLP->LiPM[2][2] = -1.0;
744 for ( j=3; j <= pLP->n; j++)
746 pLP->LiPM[1][
j] = +0.0;
747 pLP->LiPM[2][
j] = -1.0;
750 for( i= 1; i <= n; i++) {
753 for( j= 1; j <=
m; j++ )
764 PrintS(
"Matrix of Linear Programming\n");
765 print_mat( pLP->LiPM, pLP->m+1,pLP->n);
772 return (pLP->icase == 0);
786 vert= (
int *)
omAlloc( (idelem+1) *
sizeof(int) );
789 for ( i= 0; i < idelem; i++ )
792 for( i= 0; i < idelem; i++ )
798 for( j= 1; j <=
m; j++) {
799 if( !inHull( (gls->m)[i], p, m, j ) )
802 Q[
i]->addPoint( vert );
815 omFreeSize( (
void *) vert, (idelem+1) *
sizeof(
int) );
819 for( i= 0; i < idelem; i++ )
821 Print(
" \\Conv(Qi[%d]): #%d\n", i,
Q[i]->
num );
822 for ( j=1; j <=
Q[
i]->num; j++ )
846 vert= (
int *)
omAlloc( (idelem+1) *
sizeof(int) );
849 for( i= 0; i < idelem; i++ )
854 for( j= 1; j <=
m; j++) {
855 if( !inHull( (gls->m)[i], p, m, j ) )
857 if ( (id->m)[
i] ==
NULL )
859 (
id->m)[i]=
pHead(p);
879 omFreeSize( (
void *) vert, (idelem+1) *
sizeof(
int) );
883 for( i= 0; i < idelem; i++ )
902 for ( i= 0; i <
MAXVARS+2; i++ ) acoords[i]= 0;
913 int i, ii,
j,
k, col, r;
919 numverts += Qi[
i]->num;
926 pLP->LiPM[1][1] = 0.0;
927 pLP->LiPM[1][2] = 1.0;
928 for( j=3; j<=cols; j++) pLP->LiPM[1][j] = 0.0;
930 for( i=0; i <= n; i++ )
932 pLP->LiPM[i+2][1] = 1.0;
933 pLP->LiPM[i+2][2] = 0.0;
935 for( i=1; i<=
dim; i++)
937 pLP->LiPM[n+2+
i][1] = (
mprfloat)(acoords_a[i-1]);
938 pLP->LiPM[n+2+
i][2] = -shift[
i];
943 for ( i= 0; i <= n; i++ )
946 for( k= 1; k <= Qi[ii]->num; k++ )
949 for ( r= 0; r <= n; r++ )
951 if ( r == i ) pLP->LiPM[r+2][col] = -1.0;
952 else pLP->LiPM[r+2][col] = 0.0;
954 for( r= 1; r <=
dim; r++ )
955 pLP->LiPM[r+n+2][col] = -(
mprfloat)((*Qi[ii])[k]->point[r]);
960 Werror(
"mayanPyramidAlg::vDistance:" 961 "setting up matrix for udist: col %d != cols %d",col,cols);
968 Print(
"vDistance LP, known koords dim=%d, constr %d, cols %d, acoords= ",
970 for( i= 0; i <
dim; i++ )
971 Print(
" %d",acoords_a[i]);
973 print_mat( pLP->LiPM, pLP->m+1, cols);
979 PrintS(
"LP returns matrix\n");
980 print_bmat( pLP->LiPM, pLP->m+1, cols+1-pLP->m, cols, pLP->iposv);
983 if( pLP->icase != 0 )
985 WerrorS(
"mayanPyramidAlg::vDistance:");
986 if( pLP->icase == 1 )
987 WerrorS(
" Unbounded v-distance: probably 1st v-coor=0");
988 else if( pLP->icase == -1 )
989 WerrorS(
" Infeasible v-distance");
995 return pLP->LiPM[1][1];
1000 int i,
j,
k, cols, cons;
1009 pLP->LiPM[1][1] = 0.0;
1010 for( i=2; i<=n+2; i++)
1012 pLP->LiPM[
i][1] = 1.0;
1013 pLP->LiPM[
i][2] = 0.0;
1018 for( i=0; i<=n; i++)
1021 for( j=1; j<= Qi[
i]->num; j++)
1024 pLP->LiPM[1][cols] = 0.0;
1025 for( k=2; k<=n+2; k++)
1027 if( k != la_cons_row) pLP->LiPM[
k][cols] = 0.0;
1028 else pLP->LiPM[
k][cols] = -1.0;
1030 for( k=1; k<=n; k++)
1031 pLP->LiPM[k+n+2][cols] = -(
mprfloat)((*Qi[
i])[j]->point[k]);
1035 for( i= 0; i <
dim; i++ )
1037 pLP->LiPM[i+n+3][1] = acoords[
i];
1038 pLP->LiPM[i+n+3][2] = 0.0;
1040 pLP->LiPM[dim+n+3][1] = 0.0;
1043 pLP->LiPM[1][2] = -1.0;
1044 pLP->LiPM[dim+n+3][2] = 1.0;
1047 Print(
"\nThats the matrix for minR, dim= %d, acoords= ",dim);
1048 for( i= 0; i <
dim; i++ )
1049 Print(
" %d",acoords[i]);
1051 print_mat( pLP->LiPM, cons+1, cols);
1061 if ( pLP->icase != 0 )
1064 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: minR: infeasible");
1065 else if( pLP->icase > 0)
1066 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: minR: unbounded");
1075 pLP->LiPM[1][1] = 0.0;
1076 for( i=2; i<=n+2; i++)
1078 pLP->LiPM[
i][1] = 1.0;
1079 pLP->LiPM[
i][2] = 0.0;
1083 for( i=0; i<=n; i++)
1086 for( j=1; j<=Qi[
i]->num; j++)
1089 pLP->LiPM[1][cols] = 0.0;
1090 for( k=2; k<=n+2; k++)
1092 if( k != la_cons_row) pLP->LiPM[
k][cols] = 0.0;
1093 else pLP->LiPM[
k][cols] = -1.0;
1095 for( k=1; k<=n; k++)
1096 pLP->LiPM[k+n+2][cols] = -(
mprfloat)((*Qi[
i])[j]->point[k]);
1100 for( i= 0; i <
dim; i++ )
1102 pLP->LiPM[i+n+3][1] = acoords[
i];
1103 pLP->LiPM[i+n+3][2] = 0.0;
1105 pLP->LiPM[dim+n+3][1] = 0.0;
1107 pLP->LiPM[1][2] = 1.0;
1108 pLP->LiPM[dim+n+3][2] = 1.0;
1111 Print(
"\nThats the matrix for maxR, dim= %d\n",dim);
1112 print_mat( pLP->LiPM, cons+1, cols);
1122 if ( pLP->icase != 0 )
1125 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: maxR: infeasible");
1126 else if( pLP->icase > 0)
1127 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: maxR: unbounded");
1133 Print(
" Range for dim=%d: [%d,%d]\n", dim, *minR, *maxR);
1145 dist= vDistance( &(acoords[0]), n );
1154 E->addPoint( &(acoords[0]) );
1170 mn_mx_MinkowskiSum( dim, &minR, &maxR );
1174 for( i=0; i <=
dim; i++)
Print(
"acoords[%d]=%d ",i,(
int)acoords[i]);
1175 Print(
":: [%d,%d]\n", minR, maxR);
1183 acoords[
dim] = minR;
1184 while( acoords[dim] <= maxR )
1186 if( !storeMinkowskiSumPoint() )
1195 acoords[
dim] = minR;
1196 while ( acoords[dim] <= maxR )
1198 if ( (acoords[dim] > minR) && (acoords[dim] <= maxR) )
1201 runMayanPyramid( dim + 1 );
1206 dist= vDistance( &(acoords[0]), dim + 1 );
1211 runMayanPyramid( dim + 1 );
1224 for ( i= 0; i <= nn; i++ )
1226 if ( (loffset < indx) && (indx <= pQ[
i]->
num + loffset) )
1232 else loffset+= pQ[
i]->
num;
1253 for ( i= 0; i <= n; i++ )
1256 for ( k= 1; k <=
size; k++ )
1264 for ( j = 0; j <= n; j++ )
1267 LP->LiPM[j+2][LP->n] = -1.0;
1269 LP->LiPM[j+2][LP->n] = 0.0;
1273 for ( j = 1; j <= n; j++ )
1275 LP->LiPM[j+n+2][LP->n] = - ( (
mprfloat) (*pQ[i])[
k]->point[
j] );
1280 for ( j = 0; j <= n; j++ ) LP->LiPM[j+2][1] = 1.0;
1281 for ( j= 1; j <= n; j++ )
1283 LP->LiPM[j+n+2][1]= (
mprfloat)(*E)[vert]->point[
j] - shift[
j];
1287 LP->LiPM[1][1] = 0.0;
1291 Print(
" n= %d, LP->m=M= %d, LP->n=N= %d\n",n,LP->m,LP->n);
1292 print_mat(LP->LiPM, LP->m+1, LP->n+1);
1299 if ( LP->icase < 0 )
1306 (*E)[vert]->point[E->
dim]= (int)(-LP->LiPM[1][1] *
SCALEDOWN);
1309 Print(
" simplx returned %d, Objective value = %f\n", LP->icase, LP->LiPM[1][1]);
1319 for ( i= 1; i < LP->m; i++ )
1321 if ( LP->iposv[i] > LP->iposv[i+1] )
1325 LP->iposv[
i]=LP->iposv[i+1];
1328 cd=LP->LiPM[i+1][1];
1329 LP->LiPM[i+1][1]=LP->LiPM[i+2][1];
1330 LP->LiPM[i+2][1]=
cd;
1339 print_bmat(LP->LiPM, LP->m + 1, LP->n+1-LP->m, LP->n+1, LP->iposv);
1340 PrintS(
" now split into sets\n");
1345 for ( i= 0; i <= E->
dim; i++ ) bucket[i]= 0;
1349 for ( i= 0; i < LP->m; i++ )
1352 if ( LP->LiPM[i+2][1] > 1e-12 )
1354 if ( !remapXiToPoint( LP->iposv[i+1], pQ, &(optSum[c].
set), &(optSum[c].
pnt) ) )
1356 Werror(
" resMatrixSparse::RC: Found bad solution in LP: %d!",LP->iposv[i+1]);
1357 WerrorS(
" resMatrixSparse::RC: remapXiToPoint failed!");
1360 bucket[optSum[c].
set]++;
1368 for ( i= 1; i < E->
dim; i++ )
1370 if ( bucket[c] >= bucket[i] )
1376 for ( i= onum - 1; i >= 0; i-- )
1378 if ( optSum[i].
set == c )
1382 (*E)[vert]->rc.set= c;
1383 (*E)[vert]->rc.pnt= optSum[
i].
pnt;
1384 (*E)[vert]->rcPnt= (*pQ[c])[optSum[i].
pnt];
1386 if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
1388 #ifdef mprDEBUG_PROT 1389 Print(
"\n Point E[%d] was <",vert);print_exp((*E)[vert],E->
dim-1);
Print(
">, bucket={");
1390 for ( j= 0; j < E->
dim; j++ )
1392 Print(
" %d",bucket[j]);
1394 PrintS(
" }\n optimal Sum: Qi ");
1395 for ( j= 0; j < LP->m; j++ )
1397 Print(
" [ %d, %d ]",optSum[j].
set,optSum[j].
pnt);
1399 Print(
" -> i= %d, j = %d\n",(*E)[vert]->rc.set,optSum[i].
pnt);
1407 return (
int)(-LP->LiPM[1][1] *
SCALEDOWN);
1422 int *epp_mon, *eexp;
1424 epp_mon= (
int *)
omAlloc( (n+2) *
sizeof(int) );
1442 for ( i= 1; i <= E->
num; i++ )
1448 rowp=
ppMult_qq( epp, (gls->m)[(*E)[i]->rc.set] );
1453 while ( iterp!=
NULL )
1460 Werror(
"resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
1461 i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
1467 if ( (*E)[i]->rc.set == linPolyS )
1474 if ( (*E)[i]->rc.set == linPolyS )
1479 (rmat->m)[i-1]= rowp;
1483 omFreeSize( (
void *) epp_mon, (n+2) *
sizeof(
int) );
1492 for ( i= 1; i <= numSet0; i++ )
1494 Print(
" row %d contains coeffs of f_%d\n",
IMATELEM(*uRPos,i,1),linPolyS);
1496 PrintS(
" Sparse Matrix done\n");
1512 for ( j= 1; j < i-1; j++ )
1533 for ( j= 1; j <= Q1->
num; j++ )
1535 for ( k= 1; k <= Q2->
num; k++ )
1537 for ( l= 1; l <=
dim; l++ )
1539 vert.
point[
l]= (*Q1)[
j]->point[
l] + (*Q2)[
k]->point[
l];
1558 for ( j= 1; j <= pQ[0]->
num; j++ ) vs->
addPoint( (*pQ[0])[j] );
1560 for ( j= 1; j < numq; j++ )
1563 vs= minkSumTwo( vs_old, pQ[j], dim );
1585 WerrorS(
"resMatrixSparse::resMatrixSparse: Too many variables!");
1604 LP =
new simplex( idelem+totverts*2+5, totverts+5 );
1608 shift[0]=0.005; shift[1]=0.003; shift[2]=0.008; shift[3]=0.005; shift[4]=0.002;
1609 shift[5]=0.1; shift[6]=0.3; shift[7]=0.2; shift[8]=0.4; shift[9]=0.2;
1613 #ifdef mprDEBUG_PROT 1614 PrintS(
" shift vector: ");
1615 for ( i= 1; i <=
idelem; i++ )
Print(
" %.12f ",(
double)shift[i]);
1631 #ifdef mprDEBUG_PROT 1635 PrintS(
"\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
1636 for ( pnt= 1; pnt <= E->
num; pnt++ )
1645 lift[0][1]=3; lift[0][2]=4; lift[0][3]=8; lift[0][4]=2;
1646 lift[1][1]=6; lift[1][2]=1; lift[1][3]=7; lift[1][4]=4;
1647 lift[2][1]=2; lift[2][2]=5; lift[2][3]=9; lift[2][4]=6;
1648 lift[3][1]=2; lift[3][2]=1; lift[3][3]=9; lift[3][4]=5;
1649 lift[4][1]=3; lift[4][2]=7; lift[4][3]=1; lift[4][4]=5;
1651 for ( i= 0; i <=
n; i++ ) Qi[i]->
lift( lift[i] );
1654 for ( i= 0; i <=
n; i++ ) Qi[i]->
lift();
1659 for ( pnt= 1; pnt <= E->
num; pnt++ )
1661 RC( Qi, E, pnt, shift );
1666 for ( pnt= k; pnt > 0; pnt-- )
1668 if ( (*E)[pnt]->rcPnt ==
NULL )
1676 #ifdef mprDEBUG_PROT 1677 PrintS(
" points which lie in a cell:\n");
1678 for ( pnt= 1; pnt <= E->
num; pnt++ )
1686 for ( i= 0; i <=
n; i++ ) Qi[i]->unlift();
1690 #ifdef mprDEBUG_PROT 1691 Print(
" points with a[ij] (%d):\n",E->
num);
1692 for ( pnt= 1; pnt <= E->
num; pnt++ )
1694 Print(
"Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->
dim );
1695 Print(
">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
1697 print_exp( (*E)[pnt]->rcPnt, E->
dim );
PrintS(
">\n");
1704 WerrorS(
"could not handle a degenerate situation: no inner points found");
1711 WerrorS(
"resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
1717 for ( i= 0; i <
idelem; i++ )
1745 for ( i= 1; i <=
numSet0; i++ )
1804 for ( i= 1; i <=
numSet0; i++ )
1812 for ( cp= 2; cp <=
idelem; cp++ )
1814 if ( !
nIsZero(evpoint[cp-1]) )
1864 for ( i= 1; i <=
numSet0; i++ )
1870 for ( cp= 2; cp <=
idelem; cp++ )
1872 if ( !
nIsZero(evpoint[cp-1]) )
1955 number
getDetAt(
const number* evpoint );
1970 void generateBaseData();
1975 void generateMonomData(
int deg,
intvec* polyDegs ,
intvec* iVO );
1980 void generateMonoms( poly
m,
int var,
int deg );
2019 poly getElem(
const int i );
2022 number getElemNum(
const int i );
2050 assume( 0 < i || i > numColVectorSize );
2059 assume( i >= 0 && i < numColVectorSize );
2060 return( numColVector[i] );
2103 numVectors *
sizeof( number ) );
2128 for ( i= 1; i <=
MATROWS(
m ); i++ )
2129 for ( j= 1; j <=
MATCOLS(
m ); j++ )
2141 for ( i= 0; i < (
currRing->N); i++ )
2167 for ( i=0; i < (
currRing->N); i++ )
2192 poly mon =
pCopy( mm );
2213 if ( var == (
currRing->N)+1 )
return;
2214 poly newm =
pCopy( mm );
2251 ideal pDegDiv=
idInit( polyDegs->
rows(), 1 );
2252 for ( k= 0; k < polyDegs->
rows(); k++ )
2255 pSetExp( p, k + 1, (*polyDegs)[k] );
2267 for ( k= 0; k <
IDELEMS(pDegDiv); k++ )
2278 for ( k= 0; k < iVO->
rows(); k++)
2291 for ( i= 0; i <
k; i++ )
2316 for ( i= 0; i < polyDegs->
rows(); i++ )
2319 for ( k= 0; k < polyDegs->
rows(); k++ )
2320 if ( i != k ) sub*= (*polyDegs)[
k];
2333 Print(
"// %s, S(%d), db ",
2363 for ( k= (
currRing->N) - 1; k >= 0; k-- )
2375 for ( k= 0; k < (
currRing->N); k++ )
2381 for ( k= 0; k < polyDegs.
rows(); k++ )
2382 sumDeg+= polyDegs[k];
2383 sumDeg-= polyDegs.
rows() - 1;
2405 while ( pi !=
NULL )
2437 while ( pi !=
NULL )
2495 for (j=1; j <= (
currRing->N); j++ )
2497 if (
MATELEM(resmat,numVectors-i,
2561 for ( i= 0; i < (
currRing->N); i++ )
2566 nCopy(evpoint[i]) );
2602 for ( i= 1; i <=
MATROWS( mat ); i++ )
2604 for ( j= 1; j <=
MATCOLS( mat ); j++ )
2654 #define MAXEVPOINT 1000000 // 0x7fffffff 2660 unsigned long over(
const unsigned long n ,
const unsigned long d )
2665 mpz_init(m);mpz_set_ui(m,1);
2666 mpz_init(md);mpz_set_ui(md,1);
2667 mpz_init(mn);mpz_set_ui(mn,1);
2674 mpz_tdiv_q(res,m,res);
2676 mpz_clear(m);mpz_clear(md);mpz_clear(mn);
2678 unsigned long result = mpz_get_ui(res);
2707 WerrorS(
"uResultant::uResultant: Unknown chosen resultant matrix type!");
2718 ideal newGls=
idCopy( igls );
2721 (
IDELEMS(igls) + 1) *
sizeof(poly) );
2730 for ( i=
IDELEMS(newGls)-1; i > 0; i-- )
2732 newGls->m[
i]= newGls->m[i-1];
2734 newGls->m[0]= linPoly;
2738 WerrorS(
"uResultant::extendIdeal: Unknown chosen resultant matrix type!");
2749 poly actlp, rootlp= newlp;
2751 for ( i= 1; i <= (
currRing->N); i++ )
2781 long mdg=
over(
n-1, tdg );
2784 long l=(long)
pow( (
double)(tdg+1),
n );
2786 #ifdef mprDEBUG_PROT 2787 Print(
"// total deg of D: tdg %ld\n",tdg);
2788 Print(
"// maximum number of terms in D: mdg: %ld\n",mdg);
2789 Print(
"// maximum number of terms in polynom of deg tdg: l %ld\n",l);
2794 presults= (number *)
omAlloc( mdg *
sizeof( number ) );
2795 for (i=0; i < mdg; i++) presults[i]=
nInit(0);
2797 number *pevpoint= (number *)
omAlloc(
n *
sizeof( number ) );
2798 number *pev= (number *)
omAlloc(
n *
sizeof( number ) );
2799 for (i=0; i <
n; i++) pev[i]=
nInit(0);
2801 mprPROTnl(
"// initial evaluation point: ");
2804 for (i=0; i <
n; i++)
2815 for ( i=0; i < mdg; i++ )
2817 for (j=0; j <
n; j++)
2820 nPower(pevpoint[j],i,&pev[j]);
2840 if ( subDetVal !=
NULL )
2843 for ( i= 0; i <= mdg; i++ )
2845 detdiv=
nDiv( ncpoly[i], subDetVal );
2854 for ( i=0; i < mdg; i++ )
2863 for ( i=0; i < mdg; i++ )
2865 if (
nEqual(ncpoly[i],nn) )
2875 for ( i= 0; i <
n; i++ ) exp[i]=0;
2882 for ( i=0; i <
l; i++ )
2891 for ( j= 0; j <
n; j++ )
pSetExp( p, j+1, exp[j] );
2895 for ( j= 1; j <
n; j++ )
pSetExp( p, j, exp[j] );
2899 if (result!=
NULL) result=
pAdd( result, p );
2906 for ( j= 0; j < n - 1; j++ )
2927 int loops= (matchUp?
n-2:
n-1);
2929 mprPROTnl(
"uResultant::interpolateDenseSP");
2937 presults= (number *)
omAlloc( (tdg + 1) *
sizeof( number ) );
2938 for ( i=0; i <= tdg; i++ ) presults[i]=
nInit(0);
2944 number *pevpoint= (number *)
omAlloc(
n *
sizeof( number ) );
2945 for (i=0; i <
n; i++) pevpoint[i]=
nInit(0);
2947 number *pev= (number *)
omAlloc( n *
sizeof( number ) );
2948 for (i=0; i <
n; i++) pev[i]=
nInit(0);
2954 for ( uvar= 0; uvar < loops; uvar++ )
2959 for (i=0; i <
n; i++)
2968 else if ( i <= uvar + 2 )
2980 for (i=0; i <
n; i++)
2991 if ( i == (uvar + 1) ) pevpoint[
i]=
nInit(-1);
2992 else pevpoint[
i]=
nInit(0);
2999 for (i=0; i <
n; i++)
3002 pev[
i]=
nCopy( pevpoint[i] );
3005 for ( i=0; i <= tdg; i++ )
3008 nPower(pevpoint[0],i,&pev[0]);
3024 if ( subDetVal !=
NULL )
3027 for ( i= 0; i <= tdg; i++ )
3029 detdiv=
nDiv( ncpoly[i], subDetVal );
3038 for ( i=0; i <= tdg; i++ )
3052 for ( i=0; i <
n; i++ )
nDelete( pev + i );
3053 omFreeSize( (
void *)pev, n *
sizeof( number ) );
3055 for ( i=0; i <= tdg; i++ )
nDelete( presults+i );
3056 omFreeSize( (
void *)presults, (tdg + 1) *
sizeof( number ) );
3066 int loops=(matchUp?
n-2:
n-1);
3068 if (loops==0) { loops=1;nn++;}
3078 number *pevpoint= (number *)
omAlloc( nn *
sizeof( number ) );
3079 for (i=0; i < nn; i++) pevpoint[i]=
nInit(0);
3084 for ( uvar= 0; uvar < loops; uvar++ )
3089 for (i=0; i <
n; i++)
3093 if ( i <= uvar + 2 )
3098 else pevpoint[
i]=
nInit(0);
3104 for (i=0; i <
n; i++)
3108 if ( i == (uvar + 1) ) pevpoint[
i]=
nInit(-1);
3109 else pevpoint[
i]=
nInit(0);
3116 number *ncpoly= (number *)
omAlloc( (tdg+1) *
sizeof(number) );
3123 for ( i= tdg; i >= 0; i-- )
3145 if ( subDetVal !=
NULL )
3148 for ( i= 0; i <= tdg; i++ )
3150 detdiv=
nDiv( ncpoly[i], subDetVal );
3168 for ( i=0; i <
n; i++ )
nDelete( pevpoint + i );
3169 omFreeSize( (
void *)pevpoint, n *
sizeof( number ) );
3196 int totverts,idelem;
3203 for( i=0; i < idelem; i++) totverts +=
pLength( (id->m)[i] );
3205 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)
#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'.
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
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)