Coverage for /usr/lib/python3/dist-packages/sympy/polys/heuristicgcd.py: 6%

1"""Heuristic polynomial GCD algorithm (HEUGCD). """ 

2 

3from .polyerrors import HeuristicGCDFailed 

4 

5HEU_GCD_MAX = 6 

6 

7def heugcd(f, g): 

8 """ 

9 Heuristic polynomial GCD in ``Z[X]``. 

10 

11 Given univariate polynomials ``f`` and ``g`` in ``Z[X]``, returns 

12 their GCD and cofactors, i.e. polynomials ``h``, ``cff`` and ``cfg`` 

13 such that:: 

14 

15 h = gcd(f, g), cff = quo(f, h) and cfg = quo(g, h) 

16 

17 The algorithm is purely heuristic which means it may fail to compute 

18 the GCD. This will be signaled by raising an exception. In this case 

19 you will need to switch to another GCD method. 

20 

21 The algorithm computes the polynomial GCD by evaluating polynomials 

22 ``f`` and ``g`` at certain points and computing (fast) integer GCD 

23 of those evaluations. The polynomial GCD is recovered from the integer 

24 image by interpolation. The evaluation process reduces f and g variable 

25 by variable into a large integer. The final step is to verify if the 

26 interpolated polynomial is the correct GCD. This gives cofactors of 

27 the input polynomials as a side effect. 

28 

29 Examples 

30 ======== 

31 

32 >>> from sympy.polys.heuristicgcd import heugcd 

33 >>> from sympy.polys import ring, ZZ 

34 

35 >>> R, x,y, = ring("x,y", ZZ) 

36 

37 >>> f = x**2 + 2*x*y + y**2 

38 >>> g = x**2 + x*y 

39 

40 >>> h, cff, cfg = heugcd(f, g) 

41 >>> h, cff, cfg 

42 (x + y, x + y, x) 

43 

44 >>> cff*h == f 

45 True 

46 >>> cfg*h == g 

47 True 

48 

49 References 

50 ========== 

51 

52 .. [1] [Liao95]_ 

53 

54 """ 

55 assert f.ring == g.ring and f.ring.domain.is_ZZ 

56 

57 ring = f.ring 

58 x0 = ring.gens[0] 

59 domain = ring.domain 

60 

61 gcd, f, g = f.extract_ground(g) 

62 

63 f_norm = f.max_norm() 

64 g_norm = g.max_norm() 

65 

66 B = domain(2*min(f_norm, g_norm) + 29) 

67 

68 x = max(min(B, 99*domain.sqrt(B)), 

69 2*min(f_norm // abs(f.LC), 

70 g_norm // abs(g.LC)) + 4) 

71 

72 for i in range(0, HEU_GCD_MAX): 

73 ff = f.evaluate(x0, x) 

74 gg = g.evaluate(x0, x) 

75 

76 if ff and gg: 

77 if ring.ngens == 1: 

78 h, cff, cfg = domain.cofactors(ff, gg) 

79 else: 

80 h, cff, cfg = heugcd(ff, gg) 

81 

82 h = _gcd_interpolate(h, x, ring) 

83 h = h.primitive()[1] 

84 

85 cff_, r = f.div(h) 

86 

87 if not r: 

88 cfg_, r = g.div(h) 

89 

90 if not r: 

91 h = h.mul_ground(gcd) 

92 return h, cff_, cfg_ 

93 

94 cff = _gcd_interpolate(cff, x, ring) 

95 

96 h, r = f.div(cff) 

97 

98 if not r: 

99 cfg_, r = g.div(h) 

100 

101 if not r: 

102 h = h.mul_ground(gcd) 

103 return h, cff, cfg_ 

104 

105 cfg = _gcd_interpolate(cfg, x, ring) 

106 

107 h, r = g.div(cfg) 

108 

109 if not r: 

110 cff_, r = f.div(h) 

111 

112 if not r: 

113 h = h.mul_ground(gcd) 

114 return h, cff_, cfg 

115 

116 x = 73794*x * domain.sqrt(domain.sqrt(x)) // 27011 

117 

118 raise HeuristicGCDFailed('no luck') 

119 

120def _gcd_interpolate(h, x, ring): 

121 """Interpolate polynomial GCD from integer GCD. """ 

122 f, i = ring.zero, 0 

123 

124 # TODO: don't expose poly repr implementation details 

125 if ring.ngens == 1: 

126 while h: 

127 g = h % x 

128 if g > x // 2: g -= x 

129 h = (h - g) // x 

130 

131 # f += X**i*g 

132 if g: 

133 f[(i,)] = g 

134 i += 1 

135 else: 

136 while h: 

137 g = h.trunc_ground(x) 

138 h = (h - g).quo_ground(x) 

139 

140 # f += X**i*g 

141 if g: 

142 for monom, coeff in g.iterterms(): 

143 f[(i,) + monom] = coeff 

144 i += 1 

145 

146 if f.LC < 0: 

147 return -f 

148 else: 

149 return f