Changeset b7e650a in git


Ignore:
Timestamp:
Jul 7, 2017, 11:02:55 PM (7 years ago)
Author:
bendooru <bendooru@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
680ab1e30905d8a758193c8630cc2bb24c0f4da0
Parents:
fe9a47c708125057131952a3c852748ff4f27ba0
git-author:
bendooru <bendooru@users.noreply.github.com>2017-07-07 23:02:55+02:00
git-committer:
bendooru <bendooru@users.noreply.github.com>2017-07-07 23:07:00+02:00
Message:
FGLM for Groebner bases

Uses `fglm` for Groebner basis computations in `rootIsolation`, which
drops the requirement of calling `rootIsolation` in a ring with
elimination ordering. Rings in examples were changed accordingly.

We also only perform fraction 'rounding' if the ring does not use
floating point numbers.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/rootisolation.lib

    rfe9a47 rb7e650a  
    9797{
    9898    "EXAMPLE:"; echo = 2;
    99     ring R = 0,x,lp;
     99    ring R = 0,x,dp;
    100100    interval I = bounds(1);
    101101    I;
    … …  
    199199{
    200200    "EXAMPLE:"; echo = 2;
    201     ring R = 0,(x,y),lp;
     201    ring R = 0,(x,y),dp;
    202202
    203203    box B = list(bounds(0,1),
    … …  
    265265{
    266266    "EXAMPLE:"; echo = 2;
    267     ring R = 0,x(1..5),lp;
     267    ring R = 0,x(1..5),dp;
    268268
    269269    ivmat A = ivmatInit(4, 5); A;
    … …  
    280280{
    281281    "EXAMPLE:"; echo = 2;
    282     ring R = 0,x,lp;
     282    ring R = 0,x,dp;
    283283
    284284    ivmat A = ivmatInit(2,3);
    … …  
    319319{
    320320    "EXAMPLE:"; echo = 2;
    321     ring R = 0,x,lp;
     321    ring R = 0,x,dp;
    322322    ivmat A = ivmatInit(2,2);             A;
    323323    A = ivmatSet(A, 1, 2, bounds(1, 2));  A;
    … …  
    339339{
    340340    "EXAMPLE:"; echo = 2;
    341     ring R = 0,x,lp;
     341    ring R = 0,x,dp;
    342342    ivmat A = diagMatrix(2, bounds(1, 2)); A;
    343343}
    … …  
    352352example
    353353{
    354     ring R = 0,(x,y),lp;
     354    "EXAMPLE:"; echo = 2;
     355    ring R = 0,(x,y),dp;
    355356    unitMatrix(2);
    356357    unitMatrix(3);
    … …  
    463464{
    464465    "EXAMPLE:"; echo = 2;
    465     ring R = 0,(x,y),lp;
     466    ring R = 0,(x,y),dp;
    466467
    467468    ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;
    … …  
    521522
    522523proc evalJacobianAtBox(ideal I, box B)
    523 "USAGE: @code{evalJacobianAtBox(I, B)}; @code{I ideal B box}
     524"USAGE: @code{evalJacobianAtBox(I, B)}; @code{I ideal, B box}
    524525RETURN: Jacobian matrix of @code{I} where polynomials are evaluated at @code{B}
    525526PURPOSE: evaluates each polynomial of the Jacobian matrix of @code{I} using
    … …  
    546547{
    547548    "EXAMPLE:"; echo = 2;
    548     ring R = 0,(x,y),lp;
     549    ring R = 0,(x,y),dp;
    549550    ideal I = 2x2-xy+2y2-2, 2x2-3xy+3y2-2;
    550551
    … …  
    560561        0 if test is inconclusive;
    561562        box is intersection of Newton step and supplied box if applicable
    562 NOTE:   rounding is performed on fractions obtained by intersecting to prevent
    563         the size of denominators and numerators from increasing dramatically
     563NOTE:   rounding is performed on fractions obtained by matrix inversion to
     564        prevent the size of denominators and numerators from increasing
     565        dramatically
    564566EXAMPLE: example testPolyBox; tests the above for intersection of ellipses."
    565567{
    566568    int N = nvars(basering);
     569    int isreal = find(charstr(basering), "Float(");
    567570    int i;
    568571
    … …  
    614617        }
    615618
    616         // in this case, fB may have horrible fractions, so try to simplify
    617         //B = Bint;
    618 
    619         // attempt simplification of fractions
    620         // add check if we work over reals
    621         list bb;
    622         for (i = 1; i <= N; i++)
    623         {
    624             lo = B[i][1];
    625             up = B[i][2];
    626 
    627             // modify numerators of B to tighten box
    628             if (lo < Bint[i][1])
     619        if (isreal)
     620        {
     621            // fraction simplification over real basering is pointless
     622            B = Bint;
     623        }
     624        else
     625        {
     626            // attempt simplification of fractions
     627            // add check if we work over reals
     628            list bb;
     629            for (i = 1; i <= N; i++)
    629630            {
    630                 n = denominator(lo);
    631                 // floor
    632                 lo = round(Bint[i][1]*n - 1/2)/n;
     631                lo = B[i][1];
     632                up = B[i][2];
     633
     634                // modify numerators of B to tighten box
     635                if (lo < Bint[i][1])
     636                {
     637                    n = denominator(lo);
     638                    // floor
     639                    lo = round(Bint[i][1]*n - 1/2)/n;
     640                }
     641                if (up > Bint[i][2])
     642                {
     643                    n = denominator(up);
     644                    // ceil
     645                    up = round(Bint[i][2]*n + 1/2)/n;
     646                }
     647
     648                // make sure box does not grow
     649                if (lo >= B[i][1] && up <= B[i][2])
     650                {
     651                    bb[i] = bounds(lo, up);
     652                }
     653                else
     654                {
     655                    bb[i] = Bint[i];
     656                }
    633657            }
    634             if (up > Bint[i][2])
    635             {
    636                 n = denominator(up);
    637                 // ceil
    638                 up = round(Bint[i][2]*n + 1/2)/n;
    639             }
    640 
    641             // make sure box does not grow
    642             if (lo >= B[i][1] && up <= B[i][2])
    643             {
    644                 bb[i] = bounds(lo, up);
    645             }
    646             else
    647             {
    648                 bb[i] = Bint[i];
    649             }
    650         }
    651 
    652         B = bb;
     658
     659            B = bb;
     660        }
    653661
    654662        if (laststep) { return(1, B); }
    … …  
    661669{
    662670    "EXAMPLE:"; echo = 2;
    663     ring R = 0,(x,y),lp;
     671    ring R = 0,(x,y),dp;
    664672    ideal I = 2x2-xy+2y2-2, 2x2-3xy+3y2-2;
    665673
    … …  
    690698{
    691699    "EXAMPLE:"; echo = 2;
    692     ring R = 0,(x,y),lp;
     700    ring R = 0,(x,y),dp;
    693701    interval I1 = bounds(0, 1); I1;
    694702    interval I2 = bounds(0, 1); I2;
    … …  
    704712
    705713proc rootIsolationNoPreprocessing(ideal I, def start, number eps)
    706 "USAGE:  @code{rootIsolationNoPreprocessing(I, B, eps)}; @code{I ideal, B box/list of boxes,
    707         eps number};
     714"USAGE:  @code{rootIsolationNoPreprocessing(I, B, eps)}; @code{I ideal,
     715        B box/list of boxes, eps number};
    708716ASSUME: @code{I} is a zero-dimensional radical ideal
    709 RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which may contain an
    710         element of V(@code{I}), i.e. a root and @code{L2} contains boxes which contain
    711         exactly one element of V(@code{I})
     717RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which
     718        may contain an element of V(@code{I}), i.e. a root and @code{L2}
     719        contains boxes which contain exactly one element of V(@code{I})
    712720PURPOSE: Given input box(es) @code{start} we try to find all roots of @code{I}
    713721        lying in @code{start} by computing boxes that contain exactly one root.
    … …  
    798806    "EXAMPLE:"; echo = 2;
    799807
    800     ring R = 0,(x,y),lp;
     808    ring R = 0,(x,y),dp;
    801809    ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;  // V(I) has four elements
    802810
    … …  
    848856{
    849857    "EXAMPLE:"; echo = 2;
    850     ring R = 0,(x,y),lp;
     858    ring R = 0,(x,y),dp;
    851859
    852860    ideal I = x+1,y-2;
    … …  
    861869//im Moment geht das nur mit eingegebener eliminationsordnung
    862870proc rootIsolation(ideal I, box start, number eps)
    863 "USAGE:  @code{rootIsolation(I, start, eps)}; @code{I ideal, start box, eps number}
     871"USAGE:  @code{rootIsolation(I, start, eps)}; @code{I ideal, start box,
     872        eps number}
    864873ASSUME: @code{I} is a zero-dimensional radical ideal
    865         and @code{basering} is defined with an elimination ordering.
    866 RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which may contain an
    867         element of V(@code{I}), i.e. a root and @code{L2} contains boxes which contain
    868         exactly one element of V(@code{I})
    869 PURPOSE: same as @code{rootIsolationNoPreprocessing}, but speeds up computation by
    870         preprocessing starting box
     874RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which
     875        may contain an element of V(@code{I}), i.e. a root and @code{L2}
     876        contains boxes which contain exactly one element of V(@code{I})
     877PURPOSE: same as @code{rootIsolationNoPreprocessing}, but speeds up computation
     878        by preprocessing starting box
    871879THEORY: As every root of @code{I} is a root of the polynomials @code{I[i]}, we
    872880        use Groebner elimination to find univariate polynomials for every
    873         variable which have these roots as well. Applying root isolation to
    874         these univariate polynomials then provides smaller starting boxes which
    875         speed up computations in the multivariate case.
     881        variable which have these roots as well. Using that @code{I} is
     882        zero-dimensional these Groebner bases may be quickly computed using
     883        FGLM. Applying root isolation to these univariate polynomials then
     884        provides smaller starting boxes which speed up computations in the
     885        multivariate case.
     886NOTE:   This algorithm and some procedures used therein perform Groebner basis
     887        computations in @code{basering}. It is thus advised to define @code{I}
     888        w.r.t. a fast monomial ordering.
    876889EXAMPLE: example rootIsolation; for intersection of two ellipses"
    877890{
    … …  
    907920    }
    908921
    909     // construct single variable ring from basering
     922    // compute reduced GB in (hopefully) fast ordering
     923    option(redSB);
     924    ideal fastGB = groebner(I);
     925
    910926    ring rSource = basering;
    911927    list rList = ringlist(rSource);
    912     // delete all but first variable
     928    // first construct univariate ring
     929    int numords = size(rList[3]);
     930    // remove all but first variables
    913931    rList[2] = list(rList[2][1]);
    914     // adjust weight vector
    915     rList[3][1][2] = intvec(1);
    916 
    917     // compute reduced groebner bases
    918     option(redSB);
    919 
    920     ideal gbUnivarPolys, gbTemp, Ipermutation;
    921     intvec perm;
    922 
    923     // permute every variable to the end once and compute Groebner basis to
    924     // find univariate polynomial in that variable
    925     for (i = 1; i <= N; i++)
    926     {
    927         perm = 1..N;
    928         perm[i] = N;
    929         perm[N] = i;
    930 
    931         Ipermutation = fetch(rSource, I, perm);
    932         gbTemp = groebner(Ipermutation);
    933 
    934         // first poly in GB now only depends on last variable
    935         gbUnivarPolys[i] = gbTemp[1];
    936     }
    937 
    938     // do computations in univariate ring
     932    // change ordering accordingly (keep last block)
     933    rList[3] = list( list(rList[3][1][1], intvec(1)), rList[3][numords] );
     934
     935    // construct and switch to univariate ring
    939936    ring rUnivar = ring(rList);
    940     setring rUnivar;
    941 
    942     // gbUnivarPoly only contains var(N),
    943     // make sure variables are properly mapped
    944     ideal gbUnivarPolys = fetch(rSource, gbUnivarPolys, intvec(0:(N-1),1));
     937
     938    // some necessary variables
     939    ideal gbUnivar;
     940    // maps var(N) in rElimOrd(i) to var(1) in rUnivar
     941    intvec varmap = (0:(N-1)),1;
    945942    number eps = fetch(rSource, eps);
    946943    list univarResult, startBoxesPerDim;
    947944
     945    // now prepare lp-ring
     946    setring rSource;
     947    rList = ringlist(rSource);
     948    // set ordering to lp
     949    rList[3] = list( list( "lp", 1:N ), list( "C", 0 ) );
     950    // save variable order for later
     951    list varstrs = rList[2];
     952
    948953    for (i = 1; i <= N; i++)
    949954    {
    950         univarResult = rootIsolationNoPreprocessing(ideal(gbUnivarPolys[i]),
     955        // permute variables in ringlist,
     956        // create and switch to ring with elimination ordering
     957        rList[2][i] = varstrs[N];
     958        rList[2][N] = varstrs[i];
     959        ring rElimOrd = ring(rList);
     960        // get GB in elimination ordering, note that GB[1] only depends on
     961        // var(i) of rSource, which is var(N) in rElimOrd
     962        ideal GB = fglm(rSource, fastGB);
     963        GB = GB[1];
     964
     965        setring rUnivar;
     966        gbUnivar = fetch(rElimOrd, GB, varmap);
     967        // clean up ring and its elements
     968        kill rElimOrd;
     969
     970        // get boxes containing roots of univariate polynomial
     971        univarResult = rootIsolationNoPreprocessing(gbUnivar,
    951972            box(list(start[i])), eps);
    952973        // maybe result[1] is not empty, so take both
    953974        startBoxesPerDim[i] = univarResult[1] + univarResult[2];
    954975        // debug:
    955         print(string("Sieved variable ", varstr(rSource, i), " to ",
     976        print(string("Sieved variable ", varstrs[i], " to ",
    956977            size(startBoxesPerDim[i]), " intervals."));
     978
     979        // reset ring variable order
     980        setring rSource;
     981        rList[2] = varstrs;
    957982    }
    958983
    … …  
    973998    list startBoxes, sbTemp;
    974999
     1000    // prepare a list of lists
    9751001    for (i = 1; i <= numBoxes; i++)
    9761002    {
    … …  
    9971023    }
    9981024
    999     // switch back for final computation and so that boxes have proper size
    1000     setring rSource;
     1025    // since we're back in rSource box(...) will return box of proper size
    10011026    for (i = 1; i <= size(sbTemp); i++)
    10021027    {
    … …  
    10101035    "EXAMPLE:"; echo = 2;
    10111036
    1012     ring R = 0,(x,y),lp;
     1037    ring R = 0,(x,y),dp;
    10131038    ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;  // V(I) has four elements
    10141039
    … …  
    10161041    box B = list(i, i);
    10171042
    1018     list result = rootIsolation(I, B, 1/512);
    1019     size(result[1]);
    1020     size(result[2]);
     1043    list result = rootIsolation(I, B, 0);
    10211044
    10221045    result;
Note: See TracChangeset for help on using the changeset viewer.