attachment:f5.py of days10/CodingSprint

   1 """
   2 Faugere's F5
   3 
   4 These implementations are heavily inspired by John Perry's pseudocode
   5 and Singular implementation of these algorithms.
   6 
   7 See http://www.math.usm.edu/perry/Research/ for details.
   8 
   9 This implementation runs faster than the Singular script
  10 implementation but uses much more RAM for some reason.
  11 
  12 AUTHOR:
  13     -- 20081013 Martin Albrecht (initial version based on John Perry's pseudocode)
  14     -- 20081013 John Perry (loop from 0 to m-1 instead of m-1 to 0)
  15 
  16 EXAMPLE:
  17     sage: R.<x,y,z> = PolynomialRing(GF(29))
  18     sage: I =  R* [3*x^4*y + 18*x*y^4 + 4*x^3*y*z + 20*x*y^3*z + 3*x^2*z^3, \
  19                    3*x^3*y^2 + 7*x^2*y^3 + 24*y^2*z^3, 
  20                    12*x*y^4 + 17*x^4*z + 27*y^4*z + 11*x^3*z^2]
  21     sage: J = I.homogenize()
  22 
  23     sage: f5 = F5() # original F5
  24     sage: gb = f5(J)
  25     sage: f5.zero_reductions, len(gb)
  26     d 7
  27     d 9
  28     d 11
  29     d 6
  30     d 7
  31     d 8
  32     d 9
  33     d 10
  34     d 11
  35     verbose 0 (179: 283.py, top_reduction) reduction to zero.
  36     verbose 0 (179: 283.py, top_reduction) reduction to zero.
  37     verbose 0 (179: 283.py, top_reduction) reduction to zero.
  38     d 12
  39     d 13
  40     d 14
  41     d 16
  42     (3, 27)
  43 
  44     sage: f5 = F5R() # F5 with interreduced B
  45     sage: gb = f5(J)
  46     sage: f5.zero_reductions, len(gb)
  47     d 7
  48     d 9
  49     d 11
  50     d 6
  51     d 7
  52     d 8
  53     d 9
  54     d 10
  55     d 11
  56     verbose 0 (179: 283.py, top_reduction) reduction to zero.
  57     verbose 0 (179: 283.py, top_reduction) reduction to zero.
  58     verbose 0 (179: 283.py, top_reduction) reduction to zero.
  59     d 12
  60     d 13
  61     d 14
  62     d 16
  63     (3, 18)
  64 
  65     sage: f5 = F5C() # F5 with interreduced B and Gprev
  66     sage: gb = f5(J)
  67     sage: f5.zero_reductions, len(gb)
  68     d 7
  69     d 9
  70     d 11
  71     d 6
  72     d 7
  73     d 8
  74     d 9
  75     d 10
  76     d 11
  77     verbose 0 (179: 283.py, top_reduction) reduction to zero.
  78     verbose 0 (179: 283.py, top_reduction) reduction to zero.
  79     verbose 0 (179: 283.py, top_reduction) reduction to zero.
  80     d 12
  81     d 13
  82     d 14
  83     d 16
  84     (3, 18)
  85 """
  86 
  87 divides = lambda x,y: x.parent().monomial_divides(x,y)
  88 LCM = lambda f,g: f.parent().monomial_lcm(f,g)
  89 LM = lambda f: f.lm()
  90 LT = lambda f: f.lt()
  91 
  92 def compare_by_degree(f,g):
  93     if f.total_degree() > g.total_degree():
  94         return 1
  95     elif f.total_degree() < g.total_degree():
  96         return -1
  97     else:
  98         return cmp(f.lm(),g.lm())
  99 
 100 class F5:
 101     def __init__(self, F=None):
 102         if F is not None:
 103             self.Rules = [[] for _ in range(len(F))]
 104             self.L = [0]
 105             self.zero_reductions = 0
 106 
 107     def poly(self, i):
 108         return self.L[i][1]
 109 
 110     def sig(self, i):
 111         return self.L[i][0]
 112 
 113     def __call__(self, F):
 114         if isinstance(F, sage.rings.polynomial.multi_polynomial_ideal.MPolynomialIdeal):
 115             F = F.reduced_basis()
 116         else:
 117             F = Ideal(list(F)).reduced_basis()
 118         if not all(f.is_homogeneous() for f in F):
 119             F = Ideal(F).homogenize()
 120             F = F.gens()
 121         return self.basis(F)
 122 
 123     def basis(self, F):
 124         poly = self.poly
 125         incremental_basis = self.incremental_basis
 126 
 127         self.__init__(F)
 128 
 129         Rules = self.Rules
 130         L = self.L
 131 
 132         m = len(F)
 133         F = sorted(F, cmp=compare_by_degree)
 134         
 135         f0 = F[0]
 136         L[0] = (Signature(1, 0), f0*f0.lc()**(-1))
 137         Rules[0] = list()
 138 
 139         Gprev = set([0])
 140         B = [f0]
 141 
 142         for i in xrange(1,m):
 143             Gcurr = incremental_basis(F[i], i, B, Gprev)
 144             if any(poly(lambd) == 1 for lambd in Gcurr):
 145                 return set(1)
 146             Gprev = Gcurr
 147             B = [poly(l) for l in Gprev]
 148         return B
 149 
 150     def incremental_basis(self, f, i, B, Gprev):
 151         L = self.L
 152         critical_pair = self.critical_pair
 153         compute_spols = self.compute_spols
 154         reduction = self.reduction
 155         Rules = self.Rules
 156 
 157         L.append( (Signature(1,i), f*f.lc()**(-1)) )
 158         curr_idx = len(L) - 1
 159         Gcurr = Gprev.union([curr_idx])
 160         Rules[i] = list()
 161 
 162         P = reduce(lambda x,y: x.union(y), [critical_pair(curr_idx, j, i, Gprev) for j in Gprev], set())
 163         while len(P) != 0:
 164             d = min(t.degree() for (t,k,u,l,v) in P)
 165             print "d", d
 166             Pd = [(t,k,u,l,v) for (t,k,u,l,v) in P if t.degree() == d]
 167             P = P.difference(Pd)
 168             S = compute_spols(Pd)
 169             R = reduction(S, B, Gprev, Gcurr)
 170             for k in R:
 171                 P = reduce(lambda x,y: x.union(y), [critical_pair(j, k, i, Gprev) for j in Gcurr], P)
 172                 Gcurr.add(k)
 173         return Gcurr
 174 
 175     def critical_pair(self, k, l, i, Gprev):
 176         poly = self.poly
 177         sig = self.sig
 178         is_top_reducible = self.is_top_reducible
 179         is_rewritable = self.is_rewritable
 180 
 181         #print "crit_pair(%s,%s,%s,%s)"%(k, l, i, Gprev)
 182         #print self.L
 183         tk = poly(k).lt()
 184         tl = poly(l).lt()
 185         t = LCM(tk, tl)
 186         u0 = t//tk
 187         u1 = t//tl
 188         m0, e0 = sig(k)
 189         m1, e1 = sig(l)
 190         if e0 == e1 and u0*m0 == u1*m1:
 191             return set()
 192         #print "test1", e0, i, u0, m0
 193         #print "test2", e1, i, u1, m1
 194         if e0 == i and is_top_reducible(u0*m0, Gprev):
 195             #print "test1 done"
 196             return set()
 197         if e1 == i and is_top_reducible(u1*m1, Gprev):
 198             #print "test2 done"
 199             return set()
 200         if is_rewritable(u0, k) or is_rewritable(u1, l):
 201             #print "test3", is_rewritable(u0, k), is_rewritable(u1, l)
 202             return set()
 203         if u0 * sig(k) < u1 * sig(l):
 204             u0, u1 = u1, u0
 205             k, l = l, k
 206         return set([(t,k,u0,l,u1)])
 207         
 208     def compute_spols(self, P):
 209         poly = self.poly
 210         sig = self.sig
 211         spol = self.spol
 212         is_rewritable = self.is_rewritable
 213         add_rule = self.add_rule
 214 
 215         L = self.L
 216 
 217         S = list()
 218         P = sorted(P, key=lambda x: x[0])
 219         for (t,k,u,l,v) in P:
 220             if not is_rewritable(u,k) and not is_rewritable(v,l):
 221                 s = spol(poly(k), poly(l))
 222                 if s != 0:
 223                     L.append( (u * sig(k), s) )
 224                     add_rule(u * sig(k), len(L)-1)
 225                     S.append(len(L)-1)
 226         S = sorted(S, key=lambda x: sig(x))
 227         return S
 228 
 229     def spol(self, f, g):
 230         return LCM(LM(f),LM(g)) // LT(f) * f - LCM(LM(f),LM(g)) // LT(g) * g
 231 
 232     def reduction(self, S, B, Gprev, Gcurr):
 233         L = self.L
 234         sig = self.sig
 235         poly = self.poly
 236         top_reduction = self.top_reduction
 237 
 238         to_do = S
 239         completed = set()
 240         while len(to_do):
 241             k, to_do = to_do[0], to_do[1:]
 242             h = poly(k).reduce(B)
 243             L[k] = (sig(k), h)
 244             newly_completed, redo = top_reduction(k, Gprev, Gcurr.union(completed))
 245             completed = completed.union( newly_completed )
 246             for j in redo:
 247                 # insert j in to_do, sorted by increasing signature
 248                 to_do.append(j)
 249                 to_do.sort(key=lambda x: sig(x))
 250         return completed
 251 
 252     def top_reduction(self, k, Gprev, Gcurr):
 253         find_reductor = self.find_reductor
 254         add_rule = self.add_rule
 255         poly = self.poly
 256         sig = self.sig
 257         L = self.L
 258 
 259         if poly(k) == 0:
 260             verbose("reduction to zero.",level=0)
 261             self.zero_reductions += 1
 262             return set(),set()
 263         p = poly(k)
 264         J = find_reductor(k, Gprev, Gcurr)
 265         if J == set():
 266             L[k] = ( sig(k), p * p.lc()**(-1) )
 267             return set([k]),set()
 268         j = J.pop()
 269         q = poly(j)
 270         u = p.lt()//q.lt()
 271         p = p - u*q
 272         if p != 0:
 273             p = p * p.lc()**(-1)
 274         if u * sig(j) < sig(k):
 275             L[k] = (sig(k), p)
 276             return set(), set([k])
 277         else:
 278             L.append((u * sig(j), p))
 279             add_rule(u * sig(j), len(L)-1)
 280             return set(), set([k, len(L)-1])
 281 
 282     def find_reductor(self, k, Gprev, Gcurr):
 283         is_rewritable = self.is_rewritable
 284         is_top_reducible = self.is_top_reducible
 285         
 286         poly = self.poly
 287         sig = self.sig
 288         t = poly(k).lt()
 289         for j in Gcurr:
 290             tprime = poly(j).lt()
 291             if divides(tprime,t):
 292                 u = t // tprime
 293                 mj, ej = sig(j)
 294                 if u * sig(j) != sig(k) and not is_rewritable(u, j) \
 295                         and not is_top_reducible(u*mj, Gprev):
 296                     return set([j])
 297         return set()
 298                 
 299     def add_rule(self, s, k):
 300         self.Rules[s[1]].append( (s[0],k) )
 301 
 302     def is_rewritable(self, u, k):
 303         j = self.find_rewriting(u, k)
 304         return j != k
 305 
 306     def find_rewriting(self, u, k):
 307         Rules = self.Rules
 308         mk, v = self.sig(k)
 309         for ctr in reversed(xrange(len(Rules[v]))):
 310             mj, j = Rules[v][ctr]
 311             if divides(mj, u * mk):
 312                 return j
 313         return k
 314 
 315     def is_top_reducible(self, t, l):
 316         poly = self.poly
 317         for g in l:
 318             if divides(poly(g).lt(), t):
 319                 return True
 320         return False
 321 
 322 class F5R(F5):
 323     def basis(self, F):
 324         poly = self.poly
 325         incremental_basis = self.incremental_basis
 326 
 327         self.__init__(F)
 328 
 329         Rules = self.Rules
 330         L = self.L
 331 
 332         m = len(F)
 333         F = sorted(F, cmp=compare_by_degree)
 334         
 335         f0 = F[0]
 336         L[0] = (Signature(1, 0), f0*f0.lc()**(-1))
 337         Rules[0] = list()
 338 
 339         Gprev = set([0])
 340         B = [f0]
 341 
 342         for i in xrange(1,m):
 343             Gcurr = incremental_basis(F[i], i, B, Gprev)
 344             if any(poly(lambd) == 1 for lambd in Gcurr):
 345                 return set(1)
 346             Gprev = Gcurr
 347             B = Ideal([poly(l) for l in Gprev]).reduced_basis()
 348             
 349         return B
 350 
 351 class F5C(F5):
 352     def basis(self, F):
 353         poly = self.poly
 354 
 355         self.__init__(F)
 356 
 357         Rules = self.Rules
 358         L = self.L
 359 
 360         m = len(F)
 361         F = sorted(F, cmp=compare_by_degree)
 362         
 363         f0 = F[0]
 364         L[0] = (Signature(1, 0), f0*f0.lc()**(-1))
 365         Rules[0] = list()
 366 
 367         Gprev = set([0])
 368         B = set([f0])
 369 
 370         for i in xrange(1,m):
 371             Gcurr = self.incremental_basis(F[i], B, Gprev)
 372             if any(poly(lambd) == 1 for lambd in Gcurr):
 373                 return set(1)
 374             B = Ideal([poly(l) for l in Gcurr]).reduced_basis()
 375             if i != m-1:
 376                 Gprev = self.setup_reduced_basis(B)
 377         return B
 378 
 379     def incremental_basis(self, f, B, Gprev):
 380         L = self.L
 381         critical_pair = self.critical_pair
 382         compute_spols = self.compute_spols
 383         reduction = self.reduction
 384         Rules = self.Rules
 385 
 386         i = len(L)
 387         L.append( (Signature(1,i), f*f.lc()**(-1)) )
 388         Rules.append( list() )
 389         Gcurr = Gprev.union([i])
 390         P = reduce(lambda x,y: x.union(y), [critical_pair(i, j, i, Gprev) for j in Gprev], set())
 391         while len(P) != 0:
 392             d = min(t.degree() for (t,k,u,l,v) in P)
 393             print "d", d
 394             Pd = [(t,k,u,l,v) for (t,k,u,l,v) in P if t.degree() == d]
 395             P = P.difference(Pd)
 396             S = compute_spols(Pd)
 397             R = reduction(S, B, Gprev, Gcurr)
 398             for k in R:
 399                 P = reduce(lambda x,y: x.union(y), [critical_pair(j, k, i, Gprev) for j in Gcurr], P)
 400                 Gcurr.add(k)
 401         return Gcurr
 402 
 403     def setup_reduced_basis(self, B):
 404         add_rule = self.add_rule
 405         self.L = range(len(B))
 406         self.Rules = [[] for _ in range(len(B))]
 407 
 408         L = self.L
 409         Rules = self.Rules
 410         Gcurr = set()
 411 
 412         for i in range(len(B)):
 413             L[i] = (Signature(1,i), B[i])
 414             Gcurr.add( i )
 415             Rules[i] = []
 416             t = B[i].lt()
 417             for j in range(i+1, len(B)):
 418                 u = LCM(t, B[j].lt())//B[j].lt()
 419                 add_rule( Signature(u, j), -1 )
 420         return Gcurr
 421 
 422 from UserList import UserList
 423 
 424 class Signature(UserList):
 425      def __init__(self, monomial, index):
 426          UserList.__init__(self,[monomial, index])
 427          
 428      def __lt__(self, other):
 429          if self[1] < other[1]:
 430              return True
 431          elif self[1] > other[1]:
 432              return False
 433          else:
 434              return (self[0] < other[0])
 435 
 436      def __eq__(self, other):
 437          return self[0] == other[0] and self[1] == other[1]
 438    
 439      def __neq__(self, other):
 440          return self[0] != other[0] or self[1] != other[1]
 441    
 442      def __rmul__(self, other):
 443          if isinstance(self, Signature):
 444              return Signature(other * self[0], self[1])
 445          else:
 446              raise TypeError
 447 
 448 f5 = F5C()