Coverage for /usr/lib/python3/dist-packages/sympy/polys/fglmtools.py: 10%

1"""Implementation of matrix FGLM Groebner basis conversion algorithm. """ 

2 

3 

4from sympy.polys.monomials import monomial_mul, monomial_div 

5 

6def matrix_fglm(F, ring, O_to): 

7 """ 

8 Converts the reduced Groebner basis ``F`` of a zero-dimensional 

9 ideal w.r.t. ``O_from`` to a reduced Groebner basis 

10 w.r.t. ``O_to``. 

11 

12 References 

13 ========== 

14 

15 .. [1] J.C. Faugere, P. Gianni, D. Lazard, T. Mora (1994). Efficient 

16 Computation of Zero-dimensional Groebner Bases by Change of 

17 Ordering 

18 """ 

19 domain = ring.domain 

20 ngens = ring.ngens 

21 

22 ring_to = ring.clone(order=O_to) 

23 

24 old_basis = _basis(F, ring) 

25 M = _representing_matrices(old_basis, F, ring) 

26 

27 # V contains the normalforms (wrt O_from) of S 

28 S = [ring.zero_monom] 

29 V = [[domain.one] + [domain.zero] * (len(old_basis) - 1)] 

30 G = [] 

31 

32 L = [(i, 0) for i in range(ngens)] # (i, j) corresponds to x_i * S[j] 

33 L.sort(key=lambda k_l: O_to(_incr_k(S[k_l[1]], k_l[0])), reverse=True) 

34 t = L.pop() 

35 

36 P = _identity_matrix(len(old_basis), domain) 

37 

38 while True: 

39 s = len(S) 

40 v = _matrix_mul(M[t[0]], V[t[1]]) 

41 _lambda = _matrix_mul(P, v) 

42 

43 if all(_lambda[i] == domain.zero for i in range(s, len(old_basis))): 

44 # there is a linear combination of v by V 

45 lt = ring.term_new(_incr_k(S[t[1]], t[0]), domain.one) 

46 rest = ring.from_dict({S[i]: _lambda[i] for i in range(s)}) 

47 

48 g = (lt - rest).set_ring(ring_to) 

49 if g: 

50 G.append(g) 

51 else: 

52 # v is linearly independent from V 

53 P = _update(s, _lambda, P) 

54 S.append(_incr_k(S[t[1]], t[0])) 

55 V.append(v) 

56 

57 L.extend([(i, s) for i in range(ngens)]) 

58 L = list(set(L)) 

59 L.sort(key=lambda k_l: O_to(_incr_k(S[k_l[1]], k_l[0])), reverse=True) 

60 

61 L = [(k, l) for (k, l) in L if all(monomial_div(_incr_k(S[l], k), g.LM) is None for g in G)] 

62 

63 if not L: 

64 G = [ g.monic() for g in G ] 

65 return sorted(G, key=lambda g: O_to(g.LM), reverse=True) 

66 

67 t = L.pop() 

68 

69 

70def _incr_k(m, k): 

71 return tuple(list(m[:k]) + [m[k] + 1] + list(m[k + 1:])) 

72 

73 

74def _identity_matrix(n, domain): 

75 M = [[domain.zero]*n for _ in range(n)] 

76 

77 for i in range(n): 

78 M[i][i] = domain.one 

79 

80 return M 

81 

82 

83def _matrix_mul(M, v): 

84 return [sum([row[i] * v[i] for i in range(len(v))]) for row in M] 

85 

86 

87def _update(s, _lambda, P): 

88 """ 

89 Update ``P`` such that for the updated `P'` `P' v = e_{s}`. 

90 """ 

91 k = min([j for j in range(s, len(_lambda)) if _lambda[j] != 0]) 

92 

93 for r in range(len(_lambda)): 

94 if r != k: 

95 P[r] = [P[r][j] - (P[k][j] * _lambda[r]) / _lambda[k] for j in range(len(P[r]))] 

96 

97 P[k] = [P[k][j] / _lambda[k] for j in range(len(P[k]))] 

98 P[k], P[s] = P[s], P[k] 

99 

100 return P 

101 

102 

103def _representing_matrices(basis, G, ring): 

104 r""" 

105 Compute the matrices corresponding to the linear maps `m \mapsto 

106 x_i m` for all variables `x_i`. 

107 """ 

108 domain = ring.domain 

109 u = ring.ngens-1 

110 

111 def var(i): 

112 return tuple([0] * i + [1] + [0] * (u - i)) 

113 

114 def representing_matrix(m): 

115 M = [[domain.zero] * len(basis) for _ in range(len(basis))] 

116 

117 for i, v in enumerate(basis): 

118 r = ring.term_new(monomial_mul(m, v), domain.one).rem(G) 

119 

120 for monom, coeff in r.terms(): 

121 j = basis.index(monom) 

122 M[j][i] = coeff 

123 

124 return M 

125 

126 return [representing_matrix(var(i)) for i in range(u + 1)] 

127 

128 

129def _basis(G, ring): 

130 r""" 

131 Computes a list of monomials which are not divisible by the leading 

132 monomials wrt to ``O`` of ``G``. These monomials are a basis of 

133 `K[X_1, \ldots, X_n]/(G)`. 

134 """ 

135 order = ring.order 

136 

137 leading_monomials = [g.LM for g in G] 

138 candidates = [ring.zero_monom] 

139 basis = [] 

140 

141 while candidates: 

142 t = candidates.pop() 

143 basis.append(t) 

144 

145 new_candidates = [_incr_k(t, k) for k in range(ring.ngens) 

146 if all(monomial_div(_incr_k(t, k), lmg) is None 

147 for lmg in leading_monomials)] 

148 candidates.extend(new_candidates) 

149 candidates.sort(key=order, reverse=True) 

150 

151 basis = list(set(basis)) 

152 

153 return sorted(basis, key=order)