Coverage for /usr/lib/python3/dist-packages/sympy/polys/specialpolys.py: 23%

1"""Functions for generating interesting polynomials, e.g. for benchmarking. """ 

2 

3 

4from sympy.core import Add, Mul, Symbol, sympify, Dummy, symbols 

5from sympy.core.containers import Tuple 

6from sympy.core.singleton import S 

7from sympy.ntheory import nextprime 

8from sympy.polys.densearith import ( 

9 dmp_add_term, dmp_neg, dmp_mul, dmp_sqr 

10) 

11from sympy.polys.densebasic import ( 

12 dmp_zero, dmp_one, dmp_ground, 

13 dup_from_raw_dict, dmp_raise, dup_random 

14) 

15from sympy.polys.domains import ZZ 

16from sympy.polys.factortools import dup_zz_cyclotomic_poly 

17from sympy.polys.polyclasses import DMP 

18from sympy.polys.polytools import Poly, PurePoly 

19from sympy.polys.polyutils import _analyze_gens 

20from sympy.utilities import subsets, public, filldedent 

21 

22 

23@public 

24def swinnerton_dyer_poly(n, x=None, polys=False): 

25 """Generates n-th Swinnerton-Dyer polynomial in `x`. 

26 

27 Parameters 

28 ---------- 

29 n : int 

30 `n` decides the order of polynomial 

31 x : optional 

32 polys : bool, optional 

33 ``polys=True`` returns an expression, otherwise 

34 (default) returns an expression. 

35 """ 

36 if n <= 0: 

37 raise ValueError( 

38 "Cannot generate Swinnerton-Dyer polynomial of order %s" % n) 

39 

40 if x is not None: 

41 sympify(x) 

42 else: 

43 x = Dummy('x') 

44 

45 if n > 3: 

46 from sympy.functions.elementary.miscellaneous import sqrt 

47 from .numberfields import minimal_polynomial 

48 p = 2 

49 a = [sqrt(2)] 

50 for i in range(2, n + 1): 

51 p = nextprime(p) 

52 a.append(sqrt(p)) 

53 return minimal_polynomial(Add(*a), x, polys=polys) 

54 

55 if n == 1: 

56 ex = x**2 - 2 

57 elif n == 2: 

58 ex = x**4 - 10*x**2 + 1 

59 elif n == 3: 

60 ex = x**8 - 40*x**6 + 352*x**4 - 960*x**2 + 576 

61 

62 return PurePoly(ex, x) if polys else ex 

63 

64 

65@public 

66def cyclotomic_poly(n, x=None, polys=False): 

67 """Generates cyclotomic polynomial of order `n` in `x`. 

68 

69 Parameters 

70 ---------- 

71 n : int 

72 `n` decides the order of polynomial 

73 x : optional 

74 polys : bool, optional 

75 ``polys=True`` returns an expression, otherwise 

76 (default) returns an expression. 

77 """ 

78 if n <= 0: 

79 raise ValueError( 

80 "Cannot generate cyclotomic polynomial of order %s" % n) 

81 

82 poly = DMP(dup_zz_cyclotomic_poly(int(n), ZZ), ZZ) 

83 

84 if x is not None: 

85 poly = Poly.new(poly, x) 

86 else: 

87 poly = PurePoly.new(poly, Dummy('x')) 

88 

89 return poly if polys else poly.as_expr() 

90 

91 

92@public 

93def symmetric_poly(n, *gens, polys=False): 

94 """ 

95 Generates symmetric polynomial of order `n`. 

96 

97 Parameters 

98 ========== 

99 

100 polys: bool, optional (default: False) 

101 Returns a Poly object when ``polys=True``, otherwise 

102 (default) returns an expression. 

103 """ 

104 gens = _analyze_gens(gens) 

105 

106 if n < 0 or n > len(gens) or not gens: 

107 raise ValueError("Cannot generate symmetric polynomial of order %s for %s" % (n, gens)) 

108 elif not n: 

109 poly = S.One 

110 else: 

111 poly = Add(*[Mul(*s) for s in subsets(gens, int(n))]) 

112 

113 return Poly(poly, *gens) if polys else poly 

114 

115 

116@public 

117def random_poly(x, n, inf, sup, domain=ZZ, polys=False): 

118 """Generates a polynomial of degree ``n`` with coefficients in 

119 ``[inf, sup]``. 

120 

121 Parameters 

122 ---------- 

123 x 

124 `x` is the independent term of polynomial 

125 n : int 

126 `n` decides the order of polynomial 

127 inf 

128 Lower limit of range in which coefficients lie 

129 sup 

130 Upper limit of range in which coefficients lie 

131 domain : optional 

132 Decides what ring the coefficients are supposed 

133 to belong. Default is set to Integers. 

134 polys : bool, optional 

135 ``polys=True`` returns an expression, otherwise 

136 (default) returns an expression. 

137 """ 

138 poly = Poly(dup_random(n, inf, sup, domain), x, domain=domain) 

139 

140 return poly if polys else poly.as_expr() 

141 

142 

143@public 

144def interpolating_poly(n, x, X='x', Y='y'): 

145 """Construct Lagrange interpolating polynomial for ``n`` 

146 data points. If a sequence of values are given for ``X`` and ``Y`` 

147 then the first ``n`` values will be used. 

148 """ 

149 ok = getattr(x, 'free_symbols', None) 

150 

151 if isinstance(X, str): 

152 X = symbols("%s:%s" % (X, n)) 

153 elif ok and ok & Tuple(*X).free_symbols: 

154 ok = False 

155 

156 if isinstance(Y, str): 

157 Y = symbols("%s:%s" % (Y, n)) 

158 elif ok and ok & Tuple(*Y).free_symbols: 

159 ok = False 

160 

161 if not ok: 

162 raise ValueError(filldedent(''' 

163 Expecting symbol for x that does not appear in X or Y. 

164 Use `interpolate(list(zip(X, Y)), x)` instead.''')) 

165 

166 coeffs = [] 

167 numert = Mul(*[x - X[i] for i in range(n)]) 

168 

169 for i in range(n): 

170 numer = numert/(x - X[i]) 

171 denom = Mul(*[(X[i] - X[j]) for j in range(n) if i != j]) 

172 coeffs.append(numer/denom) 

173 

174 return Add(*[coeff*y for coeff, y in zip(coeffs, Y)]) 

175 

176 

177def fateman_poly_F_1(n): 

178 """Fateman's GCD benchmark: trivial GCD """ 

179 Y = [Symbol('y_' + str(i)) for i in range(n + 1)] 

180 

181 y_0, y_1 = Y[0], Y[1] 

182 

183 u = y_0 + Add(*Y[1:]) 

184 v = y_0**2 + Add(*[y**2 for y in Y[1:]]) 

185 

186 F = ((u + 1)*(u + 2)).as_poly(*Y) 

187 G = ((v + 1)*(-3*y_1*y_0**2 + y_1**2 - 1)).as_poly(*Y) 

188 

189 H = Poly(1, *Y) 

190 

191 return F, G, H 

192 

193 

194def dmp_fateman_poly_F_1(n, K): 

195 """Fateman's GCD benchmark: trivial GCD """ 

196 u = [K(1), K(0)] 

197 

198 for i in range(n): 

199 u = [dmp_one(i, K), u] 

200 

201 v = [K(1), K(0), K(0)] 

202 

203 for i in range(0, n): 

204 v = [dmp_one(i, K), dmp_zero(i), v] 

205 

206 m = n - 1 

207 

208 U = dmp_add_term(u, dmp_ground(K(1), m), 0, n, K) 

209 V = dmp_add_term(u, dmp_ground(K(2), m), 0, n, K) 

210 

211 f = [[-K(3), K(0)], [], [K(1), K(0), -K(1)]] 

212 

213 W = dmp_add_term(v, dmp_ground(K(1), m), 0, n, K) 

214 Y = dmp_raise(f, m, 1, K) 

215 

216 F = dmp_mul(U, V, n, K) 

217 G = dmp_mul(W, Y, n, K) 

218 

219 H = dmp_one(n, K) 

220 

221 return F, G, H 

222 

223 

224def fateman_poly_F_2(n): 

225 """Fateman's GCD benchmark: linearly dense quartic inputs """ 

226 Y = [Symbol('y_' + str(i)) for i in range(n + 1)] 

227 

228 y_0 = Y[0] 

229 

230 u = Add(*Y[1:]) 

231 

232 H = Poly((y_0 + u + 1)**2, *Y) 

233 

234 F = Poly((y_0 - u - 2)**2, *Y) 

235 G = Poly((y_0 + u + 2)**2, *Y) 

236 

237 return H*F, H*G, H 

238 

239 

240def dmp_fateman_poly_F_2(n, K): 

241 """Fateman's GCD benchmark: linearly dense quartic inputs """ 

242 u = [K(1), K(0)] 

243 

244 for i in range(n - 1): 

245 u = [dmp_one(i, K), u] 

246 

247 m = n - 1 

248 

249 v = dmp_add_term(u, dmp_ground(K(2), m - 1), 0, n, K) 

250 

251 f = dmp_sqr([dmp_one(m, K), dmp_neg(v, m, K)], n, K) 

252 g = dmp_sqr([dmp_one(m, K), v], n, K) 

253 

254 v = dmp_add_term(u, dmp_one(m - 1, K), 0, n, K) 

255 

256 h = dmp_sqr([dmp_one(m, K), v], n, K) 

257 

258 return dmp_mul(f, h, n, K), dmp_mul(g, h, n, K), h 

259 

260 

261def fateman_poly_F_3(n): 

262 """Fateman's GCD benchmark: sparse inputs (deg f ~ vars f) """ 

263 Y = [Symbol('y_' + str(i)) for i in range(n + 1)] 

264 

265 y_0 = Y[0] 

266 

267 u = Add(*[y**(n + 1) for y in Y[1:]]) 

268 

269 H = Poly((y_0**(n + 1) + u + 1)**2, *Y) 

270 

271 F = Poly((y_0**(n + 1) - u - 2)**2, *Y) 

272 G = Poly((y_0**(n + 1) + u + 2)**2, *Y) 

273 

274 return H*F, H*G, H 

275 

276 

277def dmp_fateman_poly_F_3(n, K): 

278 """Fateman's GCD benchmark: sparse inputs (deg f ~ vars f) """ 

279 u = dup_from_raw_dict({n + 1: K.one}, K) 

280 

281 for i in range(0, n - 1): 

282 u = dmp_add_term([u], dmp_one(i, K), n + 1, i + 1, K) 

283 

284 v = dmp_add_term(u, dmp_ground(K(2), n - 2), 0, n, K) 

285 

286 f = dmp_sqr( 

287 dmp_add_term([dmp_neg(v, n - 1, K)], dmp_one(n - 1, K), n + 1, n, K), n, K) 

288 g = dmp_sqr(dmp_add_term([v], dmp_one(n - 1, K), n + 1, n, K), n, K) 

289 

290 v = dmp_add_term(u, dmp_one(n - 2, K), 0, n - 1, K) 

291 

292 h = dmp_sqr(dmp_add_term([v], dmp_one(n - 1, K), n + 1, n, K), n, K) 

293 

294 return dmp_mul(f, h, n, K), dmp_mul(g, h, n, K), h 

295 

296# A few useful polynomials from Wang's paper ('78). 

297 

298from sympy.polys.rings import ring 

299 

300def _f_0(): 

301 R, x, y, z = ring("x,y,z", ZZ) 

302 return x**2*y*z**2 + 2*x**2*y*z + 3*x**2*y + 2*x**2 + 3*x + 4*y**2*z**2 + 5*y**2*z + 6*y**2 + y*z**2 + 2*y*z + y + 1 

303 

304def _f_1(): 

305 R, x, y, z = ring("x,y,z", ZZ) 

306 return x**3*y*z + x**2*y**2*z**2 + x**2*y**2 + 20*x**2*y*z + 30*x**2*y + x**2*z**2 + 10*x**2*z + x*y**3*z + 30*x*y**2*z + 20*x*y**2 + x*y*z**3 + 10*x*y*z**2 + x*y*z + 610*x*y + 20*x*z**2 + 230*x*z + 300*x + y**2*z**2 + 10*y**2*z + 30*y*z**2 + 320*y*z + 200*y + 600*z + 6000 

307 

308def _f_2(): 

309 R, x, y, z = ring("x,y,z", ZZ) 

310 return x**5*y**3 + x**5*y**2*z + x**5*y*z**2 + x**5*z**3 + x**3*y**2 + x**3*y*z + 90*x**3*y + 90*x**3*z + x**2*y**2*z - 11*x**2*y**2 + x**2*z**3 - 11*x**2*z**2 + y*z - 11*y + 90*z - 990 

311 

312def _f_3(): 

313 R, x, y, z = ring("x,y,z", ZZ) 

314 return x**5*y**2 + x**4*z**4 + x**4 + x**3*y**3*z + x**3*z + x**2*y**4 + x**2*y**3*z**3 + x**2*y*z**5 + x**2*y*z + x*y**2*z**4 + x*y**2 + x*y*z**7 + x*y*z**3 + x*y*z**2 + y**2*z + y*z**4 

315 

316def _f_4(): 

317 R, x, y, z = ring("x,y,z", ZZ) 

318 return -x**9*y**8*z - x**8*y**5*z**3 - x**7*y**12*z**2 - 5*x**7*y**8 - x**6*y**9*z**4 + x**6*y**7*z**3 + 3*x**6*y**7*z - 5*x**6*y**5*z**2 - x**6*y**4*z**3 + x**5*y**4*z**5 + 3*x**5*y**4*z**3 - x**5*y*z**5 + x**4*y**11*z**4 + 3*x**4*y**11*z**2 - x**4*y**8*z**4 + 5*x**4*y**7*z**2 + 15*x**4*y**7 - 5*x**4*y**4*z**2 + x**3*y**8*z**6 + 3*x**3*y**8*z**4 - x**3*y**5*z**6 + 5*x**3*y**4*z**4 + 15*x**3*y**4*z**2 + x**3*y**3*z**5 + 3*x**3*y**3*z**3 - 5*x**3*y*z**4 + x**2*z**7 + 3*x**2*z**5 + x*y**7*z**6 + 3*x*y**7*z**4 + 5*x*y**3*z**4 + 15*x*y**3*z**2 + y**4*z**8 + 3*y**4*z**6 + 5*z**6 + 15*z**4 

319 

320def _f_5(): 

321 R, x, y, z = ring("x,y,z", ZZ) 

322 return -x**3 - 3*x**2*y + 3*x**2*z - 3*x*y**2 + 6*x*y*z - 3*x*z**2 - y**3 + 3*y**2*z - 3*y*z**2 + z**3 

323 

324def _f_6(): 

325 R, x, y, z, t = ring("x,y,z,t", ZZ) 

326 return 2115*x**4*y + 45*x**3*z**3*t**2 - 45*x**3*t**2 - 423*x*y**4 - 47*x*y**3 + 141*x*y*z**3 + 94*x*y*z*t - 9*y**3*z**3*t**2 + 9*y**3*t**2 - y**2*z**3*t**2 + y**2*t**2 + 3*z**6*t**2 + 2*z**4*t**3 - 3*z**3*t**2 - 2*z*t**3 

327 

328def _w_1(): 

329 R, x, y, z = ring("x,y,z", ZZ) 

330 return 4*x**6*y**4*z**2 + 4*x**6*y**3*z**3 - 4*x**6*y**2*z**4 - 4*x**6*y*z**5 + x**5*y**4*z**3 + 12*x**5*y**3*z - x**5*y**2*z**5 + 12*x**5*y**2*z**2 - 12*x**5*y*z**3 - 12*x**5*z**4 + 8*x**4*y**4 + 6*x**4*y**3*z**2 + 8*x**4*y**3*z - 4*x**4*y**2*z**4 + 4*x**4*y**2*z**3 - 8*x**4*y**2*z**2 - 4*x**4*y*z**5 - 2*x**4*y*z**4 - 8*x**4*y*z**3 + 2*x**3*y**4*z + x**3*y**3*z**3 - x**3*y**2*z**5 - 2*x**3*y**2*z**3 + 9*x**3*y**2*z - 12*x**3*y*z**3 + 12*x**3*y*z**2 - 12*x**3*z**4 + 3*x**3*z**3 + 6*x**2*y**3 - 6*x**2*y**2*z**2 + 8*x**2*y**2*z - 2*x**2*y*z**4 - 8*x**2*y*z**3 + 2*x**2*y*z**2 + 2*x*y**3*z - 2*x*y**2*z**3 - 3*x*y*z + 3*x*z**3 - 2*y**2 + 2*y*z**2 

331 

332def _w_2(): 

333 R, x, y = ring("x,y", ZZ) 

334 return 24*x**8*y**3 + 48*x**8*y**2 + 24*x**7*y**5 - 72*x**7*y**2 + 25*x**6*y**4 + 2*x**6*y**3 + 4*x**6*y + 8*x**6 + x**5*y**6 + x**5*y**3 - 12*x**5 + x**4*y**5 - x**4*y**4 - 2*x**4*y**3 + 292*x**4*y**2 - x**3*y**6 + 3*x**3*y**3 - x**2*y**5 + 12*x**2*y**3 + 48*x**2 - 12*y**3 

335 

336def f_polys(): 

337 return _f_0(), _f_1(), _f_2(), _f_3(), _f_4(), _f_5(), _f_6() 

338 

339def w_polys(): 

340 return _w_1(), _w_2()