Skip to content

Commit

Permalink
Tidying up the optimization approach
Browse files Browse the repository at this point in the history
Added some comments and remove prints.
  • Loading branch information
slap committed Sep 21, 2023
1 parent ba53b22 commit 80f2055
Showing 1 changed file with 42 additions and 20 deletions.
62 changes: 42 additions & 20 deletions Singular/LIB/integralbasis.lib
Original file line number Diff line number Diff line change
Expand Up @@ -1187,6 +1187,12 @@ proc in_array(list l, def elem)
return(found);
}

// We compute the Puiseux segments using Hensel lifting.
// The segments are order according to the initial exponent from
// smaller to larger.
// The classes of Puiseux expansions are grouped matching the
// Puiseux segments.

proc getSegments(poly f, list classes, list slopes, int globOrder)
{
int i, j;
Expand All @@ -1209,6 +1215,7 @@ proc getSegments(poly f, list classes, list slopes, int globOrder)
}
int locOrder = maxOrder;

// If loc == 1, we remove the component outside the origin.
if(loc == 1)
{
list comps;
Expand All @@ -1230,13 +1237,13 @@ proc getSegments(poly f, list classes, list slopes, int globOrder)
fSegment;
for(i = 1; i <= size(fSegment); i++)
{
// We compute the newton poly of each polynomial to find out the
// corresponding slope
list l = newtonpoly(fSegment[i]);
"l(", i , ") = ", l;

slopeSegment = number(l[2][1]) / number(l[1][2]);

// We look for all classes with same slope as each segment
classesNew = list();

// We look for all classes with same slope as this segment
for(j = 1; j <= size(classes); j++)
{
if(slopes[j][1] * l[1][2] == slopes[j][2] * l[2][1]) // same slope
Expand All @@ -1245,24 +1252,23 @@ proc getSegments(poly f, list classes, list slopes, int globOrder)
classesNew[size(classesNew) + 1] = classes[j];
}
}

// We put everything together
out[i] = list(fSegment[i], slopeSegment, classesNew);
}
}

// We reorder according to the slope
out = sortSegments(out);


"output of getSegments: "; out;
//~;
return(out);
}

// This is used so that the puiseux expansions are always given in
// the same order. This procedure is also used in Puiseux library.
proc sortSegments(list segment);
{
//"START - integralbasis.lib - sortFactors - 1";
//"START - integralbasis.lib - sortSegments - 1";
int i, j;
int n = size(segment);
list segTemp;
Expand All @@ -1282,8 +1288,8 @@ proc sortSegments(list segment);
}


// We compute the local integral at f by exhaustive search
// See "Direct approach in the combinatorial paper
// We compute the local integral basis at f by exhaustive search
// See "Direct approach in the combinatorial paper.
proc ibSegment(poly f, list classes, int degExpand)
{
"ibSegment begins";
Expand All @@ -1292,7 +1298,8 @@ proc ibSegment(poly f, list classes, int degExpand)
list polySet;
intvec degPolySet;
intvec vy = (0,1);
// We get all the possible polynomials, includign the full polys

// We get all the possible polynomials, including the full polys
list bF = buildFactors(classes);

int degF = deg(f, vy);
Expand All @@ -1301,13 +1308,15 @@ proc ibSegment(poly f, list classes, int degExpand)
poly pGround;
for(i = 1; i <= size(classes); i++)
{
// We put everything in a list, with no repeated elements
for(j = 1; j <= size(bF[1][i]); j++)
{
if(in_array(polySet, bF[1][i][j]) == 0)
{
polySet[size(polySet)+1] = bF[1][i][j];
}
}
// We compute the full polynomial up to the integrality exponent
pGround = buildPolyGroundXRootClass(classes[i], degExpand);
if(in_array(polySet, pGround) == 0)
{
Expand All @@ -1331,17 +1340,20 @@ proc ibSegment(poly f, list classes, int degExpand)
ordExpAtPoly[i][j] = ordAtPol(polySet[i], classes[j][1]);
}
}
"polySet = "; polySet;
"degPolySet = "; degPolySet;
"ordExpAtPoly = "; ordExpAtPoly;
~;
//"polySet = "; polySet;
//"degPolySet = "; degPolySet;
//"ordExpAtPoly = "; ordExpAtPoly;
//~;

// We compute the integral basis by exhaustive search
list ib = ibExhaustive(polySet, degPolySet, ordExpAtPoly, degF);

return(ib);
}

// For each degree up to the degree of f, with compute the best
// possible product of factors by checking exhaustively all possible
// combinations such that the product has the correct degree.
proc ibExhaustive(list polySet, intvec degPolySet, list ordExpAtPoly, int degF)
{
int i;
Expand All @@ -1356,7 +1368,12 @@ proc ibExhaustive(list polySet, intvec degPolySet, list ordExpAtPoly, int degF)
{
// The best element of degree d;
oMax = 0;

// We generate the list of all possible combinations
c = summandsFac(degPolySet, d);

// We compute the valuation of each combination,
// and we keep the largest one.
for(i = 1; i <= size(c); i++)
{
o = getValuation(c[i], ordExpAtPoly);
Expand All @@ -1376,8 +1393,8 @@ proc ibExhaustive(list polySet, intvec degPolySet, list ordExpAtPoly, int degF)
ib[d];
"---";
}
"ib computation finished!";
ib;
//"ib computation finished!";
//ib;
return(ib);
}

Expand Down Expand Up @@ -1407,7 +1424,10 @@ proc getValuation(c, ordExpAtPoly)
return(oBest);
}


// List of intvecs suchs that the sum of each coordinate times the degree
// given in degPolySet equals degTotal
// This is used to compute a product of factors such that the total
// degree is equal to degTotal
proc summandsFac(intvec degPolySet, int degTotal)
{
intvec degsT;
Expand Down Expand Up @@ -1436,7 +1456,7 @@ proc summandsFac(intvec degPolySet, int degTotal)
return(l);
}


//We implement the optimization approach from scratch
proc optiBase(list ibSeg, list segmentSlope)
{
if(size(ibSeg) == 1) {return(ibSeg[1]);}
Expand Down Expand Up @@ -1554,6 +1574,8 @@ static proc ibAt0(poly f, int locBasis, int optimize)
int degExpand = int(maxExp) + 1;
dbprint(dbg, "---- Degree of expansions: ", degExpand);


// Optimization algorithm from scratch
if(optimize == 1)
{
// 1) Compute Puiseux sets (we use segments for easier implementation)
Expand All @@ -1577,7 +1599,7 @@ static proc ibAt0(poly f, int locBasis, int optimize)
int intExp = int(ib[size(ib)][1]);
IOut[1] = var(1)^intExp;

// We convery to ib shape
// We convert to ib shape
for(i = 1; i <= size(ib); i++)
{
IOut[i + 1] = ib[i][3] * var(1)^(intExp - int(ib[i][1]));
Expand Down

0 comments on commit 80f2055

Please sign in to comment.