"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/poly/groebner.py" (19 Nov 2012, 4802 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "groebner.py" see the Fossies "Dox" file reference documentation.

    1 """
    2 Groebner basis
    3 """
    4 
    5 from __future__ import division
    6 import nzmath.compatibility
    7 import nzmath.ring as ring
    8 
    9 
   10 def s_polynomial(f, g, order):
   11     """
   12     Return S-polynomial of f and g with respect to the order.
   13 
   14     S(f, g) = (lc(g)*T/lb(f))*f - (lc(f)*T/lb(g))*g,
   15     where T = lcm(lb(f), lb(g)).
   16     """
   17     f_lb, f_lc = order.leading_term(f) # term = (base, coeff)
   18     g_lb, g_lc = order.leading_term(g)
   19     t = f_lb.lcm(g_lb)
   20     return f.term_mul((t - f_lb, g_lc)) - g.term_mul((t - g_lb, f_lc))
   21 
   22 def step_reduce(reducee, reducer, order):
   23     """
   24     Return the reduced polynomial of reducee by reducer.  That is, if
   25     one of reducee's bases is divisible by the leading base of reducer
   26     with respect to the order, the returned polynomial is the result
   27     of canceling out the term.
   28     """
   29     lb, lc = order.leading_term(reducer) # term = (base, coeff)
   30     for b, c in reducee.iterterms():
   31         if b.lcm(lb) == b:
   32             return reducee + reducer.term_mul((b - lb, -c / lc)), True
   33     return reducee, False
   34 
   35 def reduce_closure(f, reducers, order):
   36     """
   37     Return normalized form of f with respect to reducer, a set of
   38     polynomials, and order.
   39     """
   40     while True:
   41         for reducer in reducers:
   42             f, reduced = step_reduce(f, reducer, order)
   43             if reduced:
   44                 break
   45         else:
   46             return f
   47 
   48 def buchberger(generating, order):
   49     """
   50     Return a Groebner basis of the ideal generated by given generating
   51     set of polynomials with respect to the order.
   52 
   53     Be careful, this implementation is very naive.
   54     """
   55     pairs = []
   56     for f in generating:
   57         for g in generating:
   58             if f is g:
   59                 continue
   60             pairs.append((f, g))
   61     groebner = list(generating)
   62 
   63     # main loop
   64     while pairs:
   65         f, g = pairs.pop()
   66         h = reduce_closure(s_polynomial(f, g, order), groebner, order)
   67         if h:
   68             for f in groebner:
   69                 pairs.append((f, h))
   70             groebner.append(h)
   71 
   72     # finish
   73     return groebner
   74 
   75 def reduce_groebner(gbasis, order):
   76     """
   77     Return the reduced Groebner basis constructed from a Groebner
   78     basis.
   79 
   80     1) lb(f) divides lb(g) => g is not in reduced Groebner basis
   81     2) monic
   82     """
   83     reduced_basis = []
   84     lb_rel = dict([(order.leading_term(g)[0], g) for g in gbasis])
   85     lbs = sorted([order.leading_term(g)[0] for g in gbasis])[::-1]
   86     for i, lbi in enumerate(lbs):
   87         for lbj in lbs[(len(lbs) - 1):i:-1]:
   88             if lbi == lbj.lcm(lbi):
   89                 # divisor found
   90                 break
   91         else:
   92             g = lb_rel[lbi]
   93             if g[lbi] != ring.getRing(g[lbi]).one:
   94                 # make it monic
   95                 g = g.scalar_mul(ring.inverse(g[lbi]))
   96             reduced_basis.append(g)
   97     return reduced_basis
   98 
   99 def normal_strategy(generating, order):
  100     """
  101     Return a Groebner basis of the ideal generated by given generating
  102     set of polynomials with respect to the order.
  103 
  104     This function uses the 'normal strategy'.
  105     """
  106     groebner = list(generating) # omit reduction
  107     pairs, lcms = [], []
  108     treat = set()
  109     for i, f in enumerate(groebner):
  110         lb_f = order.leading_term(f)[0]
  111         for j, g in enumerate(groebner):
  112             if j <= i:
  113                 continue
  114             lb_g = order.leading_term(g)[0]
  115             lcm_f_g = lb_f.lcm(lb_g)
  116             if lcm_f_g == lb_f + lb_g: # disjoint
  117                 treat.add((i, j))
  118                 treat.add((j, i))
  119             else:
  120                 # keep lcms sorted, and so pairs in paralell
  121                 k = order.bisect(lcms, lcm_f_g)
  122                 pairs.insert(k, (i, j))
  123                 lcms.insert(k, lcm_f_g)
  124 
  125     # main loop
  126     while pairs:
  127         i, j = pairs.pop(0)
  128         lcm_f_g = lcms.pop(0)
  129         treat.add((i, j))
  130         for p in range(len(groebner)):
  131             if (i, p) in treat and (p, j) in treat:
  132                 pivot = order.leading_term(groebner[p])[0]
  133                 if pivot.lcm(lcm_f_g) == lcm_f_g:
  134                     break
  135         else:
  136             f, g = groebner[i], groebner[j]
  137             h = reduce_closure(s_polynomial(f, g, order), groebner, order)
  138             if h:
  139                 lb_h = order.leading_term(h)[0]
  140                 hindex = len(groebner)
  141                 for i, f in enumerate(groebner):
  142                     lb_f = order.leading_term(f)[0]
  143                     lcm_f_h = lb_f.lcm(lb_h)
  144                     if lcm_f_h == lb_f + lb_h: # disjoint
  145                         treat.add((i, hindex))
  146                     else:
  147                         k = order.bisect(lcms, lcm_f_h)
  148                         pairs.insert(k, (i, hindex))
  149                         lcms.insert(k, lcm_f_h)
  150                 groebner.append(h)
  151 
  152     # finish
  153     return groebner