51 #define MAXPOINTS 10000 52 #define MAXINITELEMS 256 53 #define LIFT_COOR 50000 // siRand() % LIFT_COOR gives random lift value 54 #define SCALEDOWN 100.0 // lift value scale down for linear program 56 #define RVMULT 0.0001 // multiplicator for random shift vector 57 #define MAXRVVAL 50000 84 number
getDetAt(
const number* evpoint );
149 typedef struct onePoint * onePointP;
159 typedef struct _entry * entry;
179 inline onePointP operator[] (
const int index );
185 bool addPoint(
const onePointP vert );
191 bool addPoint(
const int * vert );
197 bool addPoint(
const Coord_t * vert );
200 bool removePoint(
const int indx );
205 bool mergeWithExp(
const onePointP vert );
210 bool mergeWithExp(
const int * vert );
213 void mergeWithPoly(
const poly p );
216 void getRowMP(
const int indx,
int * vert );
219 int getExpPos(
const poly p );
238 inline bool smaller(
int,
int );
241 inline bool larger(
int,
int );
246 inline bool checkMem();
263 ideal newtonPolytopesI(
const ideal gls );
269 bool inHull(
poly p,
poly pointPoly,
int m,
int site);
298 void runMayanPyramid(
int dim );
317 bool storeMinkowskiSumPoint();
333 #if defined(mprDEBUG_PROT) || defined(mprDEBUG_ALL) 334 void print_mat(
mprfloat **
a,
int maxrow,
int maxcol)
338 for (i = 1; i <= maxrow; i++)
341 for (j = 1; j <= maxcol; j++)
Print(
"% 7.2f, ", a[i][j]);
349 printf(
"Output matrix from LinProg");
350 for (i = 1; i <=
nrows; i++)
353 if (i == 1) printf(
" ");
354 else if (iposv[i-1] <=
N) printf(
"X%d", iposv[i-1]);
355 else printf(
"Y%d", iposv[i-1]-
N+1);
356 for (j = 1; j <=
ncols; j++) printf(
" %7.2f ",(
double)
a[i][j]);
362 void print_exp(
const onePointP vert,
int n )
365 for ( i= 1; i <=
n; i++ )
367 Print(
" %d",vert->point[i] );
369 if ( i < n )
PrintS(
", ");
373 void print_matrix(
matrix omat )
378 for ( i= 1; i <=
MATROWS( omat ); i++ )
380 for ( j= 1; j <=
MATCOLS( omat ); j++ )
419 points = (onePointP *)
omAlloc( (count+1) *
sizeof(onePointP) );
420 for ( i= 0; i <=
max; i++ )
432 for ( i= 0; i <=
max; i++ )
453 (
max+1) *
sizeof(onePointP),
454 (2*
max + 1) *
sizeof(onePointP) );
455 for ( i=
max+1; i <=
max*2; i++ )
474 for ( i= 1; i <=
dim; i++ )
points[
num]->point[i]= vert->point[i];
496 for ( i= 0; i <
dim; i++ )
points[
num]->point[i+1]= vert[i];
519 for ( i= 1; i <=
num; i++ )
521 for ( j= 1; j <=
dim; j++ )
522 if (
points[i]->point[j] != vert->point[j] )
break;
523 if ( j >
dim )
break;
538 for ( i= 1; i <=
num; i++ )
540 for ( j= 1; j <=
dim; j++ )
542 if ( j >
dim )
break;
558 vert= (
int *)
omAlloc( (
dim+1) *
sizeof(int) );
564 for ( i= 1; i <=
num; i++ )
566 for ( j= 1; j <=
dim; j++ )
568 if ( j >
dim )
break;
587 vert= (
int *)
omAlloc( (
dim+1) *
sizeof(int) );
590 for ( i= 1; i <=
num; i++ )
592 for ( j= 1; j <=
dim; j++ )
594 if ( j >
dim )
break;
598 if ( i >
num )
return 0;
608 for ( i= 1; i <=
dim; i++ )
609 vert[i]= (
int)(
points[indx]->point[
i] -
points[indx]->rcPnt->point[
i]);
616 for ( i= 1; i <=
dim; i++ )
635 for ( i= 1; i <=
dim; i++ )
659 for ( i= 1; i <
num; i++ )
686 for(i = 1; i <
dim; i++)
691 for ( j=1; j <=
num; j++ )
694 for ( i=1; i <
dim; i++ )
696 sum += (int)
points[j]->point[i] * l[i];
703 for ( j=1; j <
dim; j++ )
Print(
" %d ",l[j] );
706 PrintS(
" lifted points: \n");
707 for ( j=1; j <=
num; j++ )
717 if ( !outerL )
omFreeSize( (
void *) l, (dim+1) *
sizeof(
int) );
740 pLP->LiPM[1][1] = +0.0;
741 pLP->LiPM[1][2] = +1.0;
742 pLP->LiPM[2][1] = +1.0;
743 pLP->LiPM[2][2] = -1.0;
745 for ( j=3; j <= pLP->n; j++)
747 pLP->LiPM[1][
j] = +0.0;
748 pLP->LiPM[2][
j] = -1.0;
751 for( i= 1; i <= n; i++) {
754 for( j= 1; j <=
m; j++ )
765 PrintS(
"Matrix of Linear Programming\n");
766 print_mat( pLP->LiPM, pLP->m+1,pLP->n);
773 return (pLP->icase == 0);
787 vert= (
int *)
omAlloc( (idelem+1) *
sizeof(int) );
790 for ( i= 0; i < idelem; i++ )
793 for( i= 0; i < idelem; i++ )
799 for( j= 1; j <=
m; j++) {
800 if( !inHull( (gls->m)[i], p, m, j ) )
803 Q[
i]->addPoint( vert );
816 omFreeSize( (
void *) vert, (idelem+1) *
sizeof(
int) );
820 for( i= 0; i < idelem; i++ )
822 Print(
" \\Conv(Qi[%d]): #%d\n", i,
Q[i]->
num );
823 for ( j=1; j <=
Q[
i]->num; j++ )
847 vert= (
int *)
omAlloc( (idelem+1) *
sizeof(int) );
850 for( i= 0; i < idelem; i++ )
855 for( j= 1; j <=
m; j++) {
856 if( !inHull( (gls->m)[i], p, m, j ) )
858 if ( (id->m)[
i] ==
NULL )
860 (
id->m)[i]=
pHead(p);
880 omFreeSize( (
void *) vert, (idelem+1) *
sizeof(
int) );
884 for( i= 0; i < idelem; i++ )
903 for ( i= 0; i <
MAXVARS+2; i++ ) acoords[i]= 0;
914 int i, ii,
j,
k, col,
r;
920 numverts += Qi[
i]->num;
927 pLP->LiPM[1][1] = 0.0;
928 pLP->LiPM[1][2] = 1.0;
929 for( j=3; j<=cols; j++) pLP->LiPM[1][j] = 0.0;
931 for( i=0; i <= n; i++ )
933 pLP->LiPM[i+2][1] = 1.0;
934 pLP->LiPM[i+2][2] = 0.0;
936 for( i=1; i<=
dim; i++)
938 pLP->LiPM[n+2+
i][1] = (
mprfloat)(acoords_a[i-1]);
939 pLP->LiPM[n+2+
i][2] = -shift[
i];
944 for ( i= 0; i <= n; i++ )
947 for( k= 1; k <= Qi[ii]->num; k++ )
950 for ( r= 0; r <= n; r++ )
952 if ( r == i ) pLP->LiPM[r+2][col] = -1.0;
953 else pLP->LiPM[r+2][col] = 0.0;
955 for( r= 1; r <=
dim; r++ )
956 pLP->LiPM[r+n+2][col] = -(
mprfloat)((*Qi[ii])[k]->point[r]);
961 Werror(
"mayanPyramidAlg::vDistance:" 962 "setting up matrix for udist: col %d != cols %d",col,cols);
969 Print(
"vDistance LP, known koords dim=%d, constr %d, cols %d, acoords= ",
971 for( i= 0; i <
dim; i++ )
972 Print(
" %d",acoords_a[i]);
974 print_mat( pLP->LiPM, pLP->m+1, cols);
980 PrintS(
"LP returns matrix\n");
981 print_bmat( pLP->LiPM, pLP->m+1, cols+1-pLP->m, cols, pLP->iposv);
984 if( pLP->icase != 0 )
986 WerrorS(
"mayanPyramidAlg::vDistance:");
987 if( pLP->icase == 1 )
988 WerrorS(
" Unbounded v-distance: probably 1st v-coor=0");
989 else if( pLP->icase == -1 )
990 WerrorS(
" Infeasible v-distance");
996 return pLP->LiPM[1][1];
1001 int i,
j,
k, cols, cons;
1010 pLP->LiPM[1][1] = 0.0;
1011 for( i=2; i<=n+2; i++)
1013 pLP->LiPM[
i][1] = 1.0;
1014 pLP->LiPM[
i][2] = 0.0;
1019 for( i=0; i<=n; i++)
1022 for( j=1; j<= Qi[
i]->num; j++)
1025 pLP->LiPM[1][cols] = 0.0;
1026 for( k=2; k<=n+2; k++)
1028 if( k != la_cons_row) pLP->LiPM[
k][cols] = 0.0;
1029 else pLP->LiPM[
k][cols] = -1.0;
1031 for( k=1; k<=n; k++)
1032 pLP->LiPM[k+n+2][cols] = -(
mprfloat)((*Qi[
i])[j]->point[k]);
1036 for( i= 0; i <
dim; i++ )
1038 pLP->LiPM[i+n+3][1] = acoords[
i];
1039 pLP->LiPM[i+n+3][2] = 0.0;
1041 pLP->LiPM[dim+n+3][1] = 0.0;
1044 pLP->LiPM[1][2] = -1.0;
1045 pLP->LiPM[dim+n+3][2] = 1.0;
1048 Print(
"\nThats the matrix for minR, dim= %d, acoords= ",dim);
1049 for( i= 0; i <
dim; i++ )
1050 Print(
" %d",acoords[i]);
1052 print_mat( pLP->LiPM, cons+1, cols);
1062 if ( pLP->icase != 0 )
1065 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: minR: infeasible");
1066 else if( pLP->icase > 0)
1067 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: minR: unbounded");
1076 pLP->LiPM[1][1] = 0.0;
1077 for( i=2; i<=n+2; i++)
1079 pLP->LiPM[
i][1] = 1.0;
1080 pLP->LiPM[
i][2] = 0.0;
1084 for( i=0; i<=n; i++)
1087 for( j=1; j<=Qi[
i]->num; j++)
1090 pLP->LiPM[1][cols] = 0.0;
1091 for( k=2; k<=n+2; k++)
1093 if( k != la_cons_row) pLP->LiPM[
k][cols] = 0.0;
1094 else pLP->LiPM[
k][cols] = -1.0;
1096 for( k=1; k<=n; k++)
1097 pLP->LiPM[k+n+2][cols] = -(
mprfloat)((*Qi[
i])[j]->point[k]);
1101 for( i= 0; i <
dim; i++ )
1103 pLP->LiPM[i+n+3][1] = acoords[
i];
1104 pLP->LiPM[i+n+3][2] = 0.0;
1106 pLP->LiPM[dim+n+3][1] = 0.0;
1108 pLP->LiPM[1][2] = 1.0;
1109 pLP->LiPM[dim+n+3][2] = 1.0;
1112 Print(
"\nThats the matrix for maxR, dim= %d\n",dim);
1113 print_mat( pLP->LiPM, cons+1, cols);
1123 if ( pLP->icase != 0 )
1126 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: maxR: infeasible");
1127 else if( pLP->icase > 0)
1128 WerrorS(
" mn_mx_MinkowskiSum: LinearProgram: maxR: unbounded");
1134 Print(
" Range for dim=%d: [%d,%d]\n", dim, *minR, *maxR);
1146 dist= vDistance( &(acoords[0]), n );
1155 E->addPoint( &(acoords[0]) );
1171 mn_mx_MinkowskiSum( dim, &minR, &maxR );
1175 for( i=0; i <=
dim; i++)
Print(
"acoords[%d]=%d ",i,(
int)acoords[i]);
1176 Print(
":: [%d,%d]\n", minR, maxR);
1184 acoords[
dim] = minR;
1185 while( acoords[dim] <= maxR )
1187 if( !storeMinkowskiSumPoint() )
1196 acoords[
dim] = minR;
1197 while ( acoords[dim] <= maxR )
1199 if ( (acoords[dim] > minR) && (acoords[dim] <= maxR) )
1202 runMayanPyramid( dim + 1 );
1207 dist= vDistance( &(acoords[0]), dim + 1 );
1212 runMayanPyramid( dim + 1 );
1225 for ( i= 0; i <= nn; i++ )
1227 if ( (loffset < indx) && (indx <= pQ[
i]->
num + loffset) )
1233 else loffset+= pQ[
i]->
num;
1254 for ( i= 0; i <= n; i++ )
1257 for ( k= 1; k <=
size; k++ )
1265 for ( j = 0; j <= n; j++ )
1268 LP->LiPM[j+2][LP->n] = -1.0;
1270 LP->LiPM[j+2][LP->n] = 0.0;
1274 for ( j = 1; j <= n; j++ )
1276 LP->LiPM[j+n+2][LP->n] = - ( (
mprfloat) (*pQ[i])[
k]->point[
j] );
1281 for ( j = 0; j <= n; j++ ) LP->LiPM[j+2][1] = 1.0;
1282 for ( j= 1; j <= n; j++ )
1284 LP->LiPM[j+n+2][1]= (
mprfloat)(*E)[vert]->point[
j] - shift[
j];
1288 LP->LiPM[1][1] = 0.0;
1292 Print(
" n= %d, LP->m=M= %d, LP->n=N= %d\n",n,LP->m,LP->n);
1293 print_mat(LP->LiPM, LP->m+1, LP->n+1);
1300 if ( LP->icase < 0 )
1307 (*E)[vert]->point[E->
dim]= (int)(-LP->LiPM[1][1] *
SCALEDOWN);
1310 Print(
" simplx returned %d, Objective value = %f\n", LP->icase, LP->LiPM[1][1]);
1320 for ( i= 1; i < LP->m; i++ )
1322 if ( LP->iposv[i] > LP->iposv[i+1] )
1326 LP->iposv[
i]=LP->iposv[i+1];
1329 cd=LP->LiPM[i+1][1];
1330 LP->LiPM[i+1][1]=LP->LiPM[i+2][1];
1331 LP->LiPM[i+2][1]=
cd;
1340 print_bmat(LP->LiPM, LP->m + 1, LP->n+1-LP->m, LP->n+1, LP->iposv);
1341 PrintS(
" now split into sets\n");
1346 for ( i= 0; i <= E->
dim; i++ ) bucket[i]= 0;
1350 for ( i= 0; i < LP->m; i++ )
1353 if ( LP->LiPM[i+2][1] > 1e-12 )
1355 if ( !remapXiToPoint( LP->iposv[i+1], pQ, &(optSum[c].
set), &(optSum[c].
pnt) ) )
1357 Werror(
" resMatrixSparse::RC: Found bad solution in LP: %d!",LP->iposv[i+1]);
1358 WerrorS(
" resMatrixSparse::RC: remapXiToPoint faild!");
1361 bucket[optSum[c].
set]++;
1369 for ( i= 1; i < E->
dim; i++ )
1371 if ( bucket[c] >= bucket[i] )
1377 for ( i= onum - 1; i >= 0; i-- )
1379 if ( optSum[i].
set == c )
1383 (*E)[vert]->rc.set= c;
1384 (*E)[vert]->rc.pnt= optSum[
i].
pnt;
1385 (*E)[vert]->rcPnt= (*pQ[c])[optSum[i].
pnt];
1387 if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
1389 #ifdef mprDEBUG_PROT 1390 Print(
"\n Point E[%d] was <",vert);print_exp((*E)[vert],E->
dim-1);
Print(
">, bucket={");
1391 for ( j= 0; j < E->
dim; j++ )
1393 Print(
" %d",bucket[j]);
1395 PrintS(
" }\n optimal Sum: Qi ");
1396 for ( j= 0; j < LP->m; j++ )
1398 Print(
" [ %d, %d ]",optSum[j].
set,optSum[j].
pnt);
1400 Print(
" -> i= %d, j = %d\n",(*E)[vert]->rc.set,optSum[i].
pnt);
1408 return (
int)(-LP->LiPM[1][1] *
SCALEDOWN);
1423 int *epp_mon, *eexp;
1425 epp_mon= (
int *)
omAlloc( (n+2) *
sizeof(int) );
1443 for ( i= 1; i <= E->
num; i++ )
1449 rowp=
ppMult_qq( epp, (gls->m)[(*E)[i]->rc.set] );
1454 while ( iterp!=
NULL )
1461 Werror(
"resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
1462 i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
1468 if ( (*E)[i]->rc.set == linPolyS )
1475 if ( (*E)[i]->rc.set == linPolyS )
1480 (rmat->m)[i-1]= rowp;
1484 omFreeSize( (
void *) epp_mon, (n+2) *
sizeof(
int) );
1493 for ( i= 1; i <= numSet0; i++ )
1495 Print(
" row %d contains coeffs of f_%d\n",
IMATELEM(*uRPos,i,1),linPolyS);
1497 PrintS(
" Sparse Matrix done\n");
1513 for ( j= 1; j < i-1; j++ )
1534 for ( j= 1; j <= Q1->
num; j++ )
1536 for ( k= 1; k <= Q2->
num; k++ )
1538 for ( l= 1; l <=
dim; l++ )
1540 vert.
point[
l]= (*Q1)[
j]->point[
l] + (*Q2)[
k]->point[
l];
1559 for ( j= 1; j <= pQ[0]->
num; j++ ) vs->
addPoint( (*pQ[0])[j] );
1561 for ( j= 1; j < numq; j++ )
1564 vs= minkSumTwo( vs_old, pQ[j], dim );
1586 WerrorS(
"resMatrixSparse::resMatrixSparse: Too many variables!");
1605 LP =
new simplex( idelem+totverts*2+5, totverts+5 );
1609 shift[0]=0.005; shift[1]=0.003; shift[2]=0.008; shift[3]=0.005; shift[4]=0.002;
1610 shift[5]=0.1; shift[6]=0.3; shift[7]=0.2; shift[8]=0.4; shift[9]=0.2;
1614 #ifdef mprDEBUG_PROT 1615 PrintS(
" shift vector: ");
1616 for ( i= 1; i <=
idelem; i++ )
Print(
" %.12f ",(
double)shift[i]);
1632 #ifdef mprDEBUG_PROT 1636 PrintS(
"\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
1637 for ( pnt= 1; pnt <= E->
num; pnt++ )
1646 lift[0][1]=3; lift[0][2]=4; lift[0][3]=8; lift[0][4]=2;
1647 lift[1][1]=6; lift[1][2]=1; lift[1][3]=7; lift[1][4]=4;
1648 lift[2][1]=2; lift[2][2]=5; lift[2][3]=9; lift[2][4]=6;
1649 lift[3][1]=2; lift[3][2]=1; lift[3][3]=9; lift[3][4]=5;
1650 lift[4][1]=3; lift[4][2]=7; lift[4][3]=1; lift[4][4]=5;
1652 for ( i= 0; i <=
n; i++ ) Qi[i]->
lift( lift[i] );
1655 for ( i= 0; i <=
n; i++ ) Qi[i]->
lift();
1660 for ( pnt= 1; pnt <= E->
num; pnt++ )
1662 RC( Qi, E, pnt, shift );
1667 for ( pnt= k; pnt > 0; pnt-- )
1669 if ( (*E)[pnt]->rcPnt ==
NULL )
1677 #ifdef mprDEBUG_PROT 1678 PrintS(
" points which lie in a cell:\n");
1679 for ( pnt= 1; pnt <= E->
num; pnt++ )
1687 for ( i= 0; i <=
n; i++ ) Qi[i]->unlift();
1691 #ifdef mprDEBUG_PROT 1692 Print(
" points with a[ij] (%d):\n",E->
num);
1693 for ( pnt= 1; pnt <= E->
num; pnt++ )
1695 Print(
"Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->
dim );
1696 Print(
">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
1698 print_exp( (*E)[pnt]->rcPnt, E->
dim );
PrintS(
">\n");
1705 WerrorS(
"could not handle a degenerate situation: no inner points found");
1712 WerrorS(
"resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
1718 for ( i= 0; i <
idelem; i++ )
1746 for ( i= 1; i <=
numSet0; i++ )
1805 for ( i= 1; i <=
numSet0; i++ )
1813 for ( cp= 2; cp <=
idelem; cp++ )
1815 if ( !
nIsZero(evpoint[cp-1]) )
1865 for ( i= 1; i <=
numSet0; i++ )
1871 for ( cp= 2; cp <=
idelem; cp++ )
1873 if ( !
nIsZero(evpoint[cp-1]) )
1956 number
getDetAt(
const number* evpoint );
1971 void generateBaseData();
1976 void generateMonomData(
int deg,
intvec* polyDegs ,
intvec* iVO );
1981 void generateMonoms(
poly m,
int var,
int deg );
2020 poly getElem(
const int i );
2023 number getElemNum(
const int i );
2051 assume( 0 < i || i > numColVectorSize );
2060 assume( i >= 0 && i < numColVectorSize );
2061 return( numColVector[i] );
2104 numVectors *
sizeof( number ) );
2129 for ( i= 1; i <=
MATROWS(
m ); i++ )
2130 for ( j= 1; j <=
MATCOLS(
m ); j++ )
2142 for ( i= 0; i < (
currRing->N); i++ )
2168 for ( i=0; i < (
currRing->N); i++ )
2214 if ( var == (
currRing->N)+1 )
return;
2252 ideal pDegDiv=
idInit( polyDegs->
rows(), 1 );
2253 for ( k= 0; k < polyDegs->
rows(); k++ )
2256 pSetExp( p, k + 1, (*polyDegs)[k] );
2268 for ( k= 0; k <
IDELEMS(pDegDiv); k++ )
2279 for ( k= 0; k < iVO->
rows(); k++)
2292 for ( i= 0; i <
k; i++ )
2317 for ( i= 0; i < polyDegs->
rows(); i++ )
2320 for ( k= 0; k < polyDegs->
rows(); k++ )
2321 if ( i != k ) sub*= (*polyDegs)[
k];
2334 Print(
"// %s, S(%d), db ",
2364 for ( k= (
currRing->N) - 1; k >= 0; k-- )
2376 for ( k= 0; k < (
currRing->N); k++ )
2382 for ( k= 0; k < polyDegs.
rows(); k++ )
2383 sumDeg+= polyDegs[k];
2384 sumDeg-= polyDegs.
rows() - 1;
2406 while ( pi !=
NULL )
2438 while ( pi !=
NULL )
2496 for (j=1; j <= (
currRing->N); j++ )
2498 if (
MATELEM(resmat,numVectors-i,
2562 for ( i= 0; i < (
currRing->N); i++ )
2565 nCopy(evpoint[i]) );
2601 for ( i= 1; i <=
MATROWS( mat ); i++ )
2603 for ( j= 1; j <=
MATCOLS( mat ); j++ )
2653 #define MAXEVPOINT 1000000 // 0x7fffffff 2659 unsigned long over(
const unsigned long n ,
const unsigned long d )
2664 mpz_init(m);mpz_set_ui(m,1);
2665 mpz_init(md);mpz_set_ui(md,1);
2666 mpz_init(mn);mpz_set_ui(mn,1);
2673 mpz_tdiv_q(res,m,res);
2675 mpz_clear(m);mpz_clear(md);mpz_clear(mn);
2677 unsigned long result = mpz_get_ui(res);
2706 WerrorS(
"uResultant::uResultant: Unknown chosen resultant matrix type!");
2717 ideal newGls=
idCopy( igls );
2729 for ( i=
IDELEMS(newGls)-1; i > 0; i-- )
2731 newGls->m[
i]= newGls->m[i-1];
2733 newGls->m[0]= linPoly;
2737 WerrorS(
"uResultant::extendIdeal: Unknown chosen resultant matrix type!");
2748 poly actlp, rootlp= newlp;
2750 for ( i= 1; i <= (
currRing->N); i++ )
2780 long mdg=
over(
n-1, tdg );
2783 long l=(long)
pow( (
double)(tdg+1),
n );
2785 #ifdef mprDEBUG_PROT 2786 Print(
"// total deg of D: tdg %ld\n",tdg);
2787 Print(
"// maximum number of terms in D: mdg: %ld\n",mdg);
2788 Print(
"// maximum number of terms in polynom of deg tdg: l %ld\n",l);
2793 presults= (number *)
omAlloc( mdg *
sizeof( number ) );
2794 for (i=0; i < mdg; i++) presults[i]=
nInit(0);
2796 number *pevpoint= (number *)
omAlloc(
n *
sizeof( number ) );
2797 number *pev= (number *)
omAlloc(
n *
sizeof( number ) );
2798 for (i=0; i <
n; i++) pev[i]=
nInit(0);
2800 mprPROTnl(
"// initial evaluation point: ");
2803 for (i=0; i <
n; i++)
2814 for ( i=0; i < mdg; i++ )
2816 for (j=0; j <
n; j++)
2819 nPower(pevpoint[j],i,&pev[j]);
2839 if ( subDetVal !=
NULL )
2842 for ( i= 0; i <= mdg; i++ )
2844 detdiv=
nDiv( ncpoly[i], subDetVal );
2853 for ( i=0; i < mdg; i++ )
2862 for ( i=0; i < mdg; i++ )
2864 if (
nEqual(ncpoly[i],nn) )
2874 for ( i= 0; i <
n; i++ ) exp[i]=0;
2881 for ( i=0; i <
l; i++ )
2890 for ( j= 0; j <
n; j++ )
pSetExp( p, j+1, exp[j] );
2894 for ( j= 1; j <
n; j++ )
pSetExp( p, j, exp[j] );
2898 if (result!=
NULL) result=
pAdd( result, p );
2905 for ( j= 0; j < n - 1; j++ )
2926 int loops= (matchUp?
n-2:
n-1);
2928 mprPROTnl(
"uResultant::interpolateDenseSP");
2936 presults= (number *)
omAlloc( (tdg + 1) *
sizeof( number ) );
2937 for ( i=0; i <= tdg; i++ ) presults[i]=
nInit(0);
2943 number *pevpoint= (number *)
omAlloc(
n *
sizeof( number ) );
2944 for (i=0; i <
n; i++) pevpoint[i]=
nInit(0);
2946 number *pev= (number *)
omAlloc( n *
sizeof( number ) );
2947 for (i=0; i <
n; i++) pev[i]=
nInit(0);
2953 for ( uvar= 0; uvar < loops; uvar++ )
2958 for (i=0; i <
n; i++)
2967 else if ( i <= uvar + 2 )
2979 for (i=0; i <
n; i++)
2990 if ( i == (uvar + 1) ) pevpoint[
i]=
nInit(-1);
2991 else pevpoint[
i]=
nInit(0);
2998 for (i=0; i <
n; i++)
3001 pev[
i]=
nCopy( pevpoint[i] );
3004 for ( i=0; i <= tdg; i++ )
3007 nPower(pevpoint[0],i,&pev[0]);
3023 if ( subDetVal !=
NULL )
3026 for ( i= 0; i <= tdg; i++ )
3028 detdiv=
nDiv( ncpoly[i], subDetVal );
3037 for ( i=0; i <= tdg; i++ )
3051 for ( i=0; i <
n; i++ )
nDelete( pev + i );
3052 omFreeSize( (
void *)pev, n *
sizeof( number ) );
3054 for ( i=0; i <= tdg; i++ )
nDelete( presults+i );
3055 omFreeSize( (
void *)presults, (tdg + 1) *
sizeof( number ) );
3065 int loops=(matchUp?
n-2:
n-1);
3067 if (loops==0) { loops=1;nn++;}
3077 number *pevpoint= (number *)
omAlloc( nn *
sizeof( number ) );
3078 for (i=0; i < nn; i++) pevpoint[i]=
nInit(0);
3083 for ( uvar= 0; uvar < loops; uvar++ )
3088 for (i=0; i <
n; i++)
3092 if ( i <= uvar + 2 )
3097 else pevpoint[
i]=
nInit(0);
3103 for (i=0; i <
n; i++)
3107 if ( i == (uvar + 1) ) pevpoint[
i]=
nInit(-1);
3108 else pevpoint[
i]=
nInit(0);
3115 number *ncpoly= (number *)
omAlloc( (tdg+1) *
sizeof(number) );
3122 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
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)
#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)
static int pLength(poly a)
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)
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
#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
for(int i=0;i< R->ExpL_Size;i++) Print("%09lx "
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)
#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()
#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)
#define pGetExpV(p, e)
Gets a copy of (resp. set) the exponent vector, where e is assumed to point to (r->N +1)*sizeof(long)...
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)
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)