1 /*++ 2 Copyright (c) 2011 Microsoft Corporation 3 4 Module Name: 5 6 upolynomial.h 7 8 Abstract: 9 10 Goodies for creating and handling univariate polynomials. 11 12 A dense representation is much better for Root isolation algorithms, 13 encoding algebraic numbers, factorization, etc. 14 15 We also use integers as the coefficients of univariate polynomials. 16 17 Author: 18 19 Leonardo (leonardo) 2011-11-29 20 21 Notes: 22 23 --*/ 24 #pragma once 25 26 #include "util/mpzzp.h" 27 #include "util/rational.h" 28 #include "math/polynomial/polynomial.h" 29 #include "util/z3_exception.h" 30 #include "util/mpbq.h" 31 #include "util/rlimit.h" 32 #define FACTOR_VERBOSE_LVL 1000 33 34 namespace upolynomial { 35 36 typedef polynomial::factor_params factor_params; 37 38 // It is used only for signing cancellation. 39 class upolynomial_exception : public default_exception { 40 public: upolynomial_exception(char const * msg)41 upolynomial_exception(char const * msg):default_exception(msg) {} 42 }; 43 44 typedef mpz numeral; 45 typedef mpzzp_manager numeral_manager; 46 typedef mpzzp_manager zp_numeral_manager; 47 typedef unsynch_mpz_manager z_numeral_manager; 48 typedef svector<numeral> numeral_vector; 49 50 class core_manager { 51 public: 52 typedef _scoped_numeral_vector<numeral_manager> scoped_numeral_vector; 53 typedef _scoped_numeral<numeral_manager> scoped_numeral; 54 /** 55 \brief Convenient vector of polynomials that manages its own memory and keeps the degree, of each polynomial. 56 Polynomial is c*f_1^k_1*...*f_n^k_n. 57 */ 58 class factors { 59 private: 60 vector<numeral_vector> m_factors; 61 svector<unsigned> m_degrees; 62 core_manager & m_upm; 63 numeral m_constant; 64 unsigned m_total_factors; 65 unsigned m_total_degree; 66 public: 67 factors(core_manager & upm); 68 ~factors(); 69 upm()70 core_manager & upm() const { return m_upm; } pm()71 core_manager & pm() const { return m_upm; } 72 numeral_manager & nm() const; 73 distinct_factors()74 unsigned distinct_factors() const { return m_factors.size(); } total_factors()75 unsigned total_factors() const { return m_total_factors; } 76 void clear(); reset()77 void reset() { clear(); } 78 79 numeral_vector const & operator[](unsigned i) const { return m_factors[i]; } 80 get_constant()81 numeral const & get_constant() const { return m_constant; } 82 void set_constant(numeral const & constant); 83 get_degree()84 unsigned get_degree() const { return m_total_degree; } get_degree(unsigned i)85 unsigned get_degree(unsigned i) const { return m_degrees[i]; } 86 void set_degree(unsigned i, unsigned degree); 87 void push_back(numeral_vector const & p, unsigned degree); 88 // push p to vectors and kill it 89 void push_back_swap(numeral_vector & p, unsigned degree); 90 91 void swap_factor(unsigned i, numeral_vector & p); 92 void swap(factors & other); 93 void multiply(numeral_vector & out) const; 94 95 void display(std::ostream & out) const; 96 97 friend std::ostream & operator<<(std::ostream & out, factors const & fs) { 98 fs.display(out); 99 return out; 100 } 101 }; 102 103 protected: 104 reslimit& m_limit; 105 numeral_manager m_manager; 106 numeral_vector m_basic_tmp; 107 numeral_vector m_div_tmp1; 108 numeral_vector m_div_tmp2; 109 numeral_vector m_exact_div_tmp; 110 numeral_vector m_gcd_tmp1; 111 numeral_vector m_gcd_tmp2; 112 numeral_vector m_CRA_tmp; 113 #define UPOLYNOMIAL_MGCD_TMPS 6 114 numeral_vector m_mgcd_tmp[UPOLYNOMIAL_MGCD_TMPS]; 115 numeral_vector m_sqf_tmp1; 116 numeral_vector m_sqf_tmp2; 117 numeral_vector m_pw_tmp; 118 is_alias(numeral const * p,numeral_vector & buffer)119 static bool is_alias(numeral const * p, numeral_vector & buffer) { return buffer.c_ptr() != nullptr && buffer.c_ptr() == p; } 120 void neg_core(unsigned sz1, numeral const * p1, numeral_vector & buffer); 121 void add_core(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & buffer); 122 void sub_core(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & buffer); 123 void mul_core(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & buffer); 124 125 void flip_sign_if_lm_neg(numeral_vector & buffer); 126 127 void mod_gcd(unsigned sz_u, numeral const * u, unsigned sz_v, numeral const * v, numeral_vector & result); 128 void CRA_combine_images(numeral_vector const & q, numeral const & p, numeral_vector & C, numeral & bound); 129 130 public: 131 core_manager(reslimit& lim, z_numeral_manager & m); 132 ~core_manager(); 133 zm()134 z_numeral_manager & zm() const { return m_manager.m(); } m()135 numeral_manager & m() const { return const_cast<core_manager*>(this)->m_manager; } lim()136 reslimit& lim() const { return m_limit; } 137 /** 138 \brief Return true if Z_p[X] 139 */ modular()140 bool modular() const { return m().modular(); } field()141 bool field() const { return m().field(); } 142 /** 143 \brief Return p in Z_p[X] 144 \pre modular 145 */ p()146 numeral const & p() const { return m().p(); } 147 /** 148 \brief Set manager as Z[X] 149 */ set_z()150 void set_z() { m().set_z(); } 151 /** 152 \brief Set manager as Z_p[X] 153 */ set_zp(numeral const & p)154 void set_zp(numeral const & p) { m().set_zp(p); } set_zp(uint64_t p)155 void set_zp(uint64_t p) { m().set_zp(p); } 156 157 void checkpoint(); 158 159 160 /** 161 \brief set p size to 0. That is, p is the zero polynomial after this operation. 162 */ 163 void reset(numeral_vector & p); 164 165 /** 166 \brief Remove zero leading coefficients. 167 After applying this method, we have that p is empty() or p[p.size()-1] is not zero. 168 */ 169 void trim(numeral_vector & p); 170 171 void set_size(unsigned sz, numeral_vector & buffer); 172 173 /** 174 \brief Return the actual degree of p. 175 */ degree(numeral_vector const & p)176 unsigned degree(numeral_vector const & p) { 177 unsigned sz = p.size(); 178 return sz == 0 ? 0 : sz - 1; 179 } 180 181 /** 182 \brief Return true if p is the zero polynomial. 183 */ is_zero(numeral_vector & p)184 bool is_zero(numeral_vector & p) { trim(p); return p.empty(); } 185 186 /** 187 \brief Return true if p is a constant polynomial 188 */ is_const(numeral_vector & p)189 bool is_const(numeral_vector & p) { trim(p); return p.size() <= 1; } 190 191 /** 192 \brief Copy p to buffer. 193 */ 194 void set(unsigned sz, numeral const * p, numeral_vector & buffer); set(numeral_vector & target,numeral_vector const & source)195 void set(numeral_vector & target, numeral_vector const & source) { set(source.size(), source.c_ptr(), target); } 196 197 /** 198 \brief Copy p to buffer. 199 200 Coefficients of p must be integer. 201 */ 202 void set(unsigned sz, rational const * p, numeral_vector & buffer); 203 204 /** 205 \brief Compute the primitive part and the content of f (pp can alias f). 206 */ 207 void get_primitive_and_content(unsigned f_sz, numeral const * f, numeral_vector & pp, numeral & cont); get_primitive_and_content(numeral_vector const & f,numeral_vector & pp,numeral & cont)208 void get_primitive_and_content(numeral_vector const & f, numeral_vector & pp, numeral & cont) { 209 get_primitive_and_content(f.size(), f.c_ptr(), pp, cont); 210 } get_primitive(numeral_vector const & f,numeral_vector & pp)211 void get_primitive(numeral_vector const & f, numeral_vector & pp) { 212 scoped_numeral cont(m()); 213 get_primitive_and_content(f.size(), f.c_ptr(), pp, cont); 214 } 215 216 /** 217 \brief p := -p 218 */ 219 void neg(unsigned sz, numeral * p); neg(numeral_vector & p)220 void neg(numeral_vector & p) { neg(p.size(), p.c_ptr()); } 221 222 /** 223 \brief buffer := -p 224 */ 225 void neg(unsigned sz, numeral const * p, numeral_vector & buffer); neg(numeral_vector const & p,numeral_vector & p_neg)226 void neg(numeral_vector const & p, numeral_vector & p_neg) { neg(p.size(), p.c_ptr(), p_neg); } 227 228 /** 229 \brief buffer := p1 + p2 230 */ 231 void add(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & buffer); add(numeral_vector const & a,numeral_vector const & b,numeral_vector & c)232 void add(numeral_vector const & a, numeral_vector const & b, numeral_vector & c) { add(a.size(), a.c_ptr(), b.size(), b.c_ptr(), c); } 233 234 /** 235 \brief buffer := p1 - p2 236 */ 237 void sub(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & buffer); sub(numeral_vector const & a,numeral_vector const & b,numeral_vector & c)238 void sub(numeral_vector const & a, numeral_vector const & b, numeral_vector & c) { sub(a.size(), a.c_ptr(), b.size(), b.c_ptr(), c); } 239 240 /** 241 \brief buffer := p1 * p2 242 */ 243 void mul(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & buffer); mul(numeral_vector const & a,numeral_vector const & b,numeral_vector & c)244 void mul(numeral_vector const & a, numeral_vector const & b, numeral_vector & c) { mul(a.size(), a.c_ptr(), b.size(), b.c_ptr(), c); } 245 246 /** 247 \brief r := p^k 248 */ 249 void pw(unsigned sz, numeral const * p, unsigned k, numeral_vector & r); 250 251 /** 252 \brief buffer := dp/dx 253 */ 254 void derivative(unsigned sz1, numeral const * p, numeral_vector & buffer); derivative(numeral_vector const & p,numeral_vector & d_p)255 void derivative(numeral_vector const & p, numeral_vector & d_p) { derivative(p.size(), p.c_ptr(), d_p); } 256 257 /** 258 \brief Divide coefficients of p by their GCD 259 */ 260 void normalize(unsigned sz, numeral * p); 261 262 /** 263 \brief Divide coefficients of p by their GCD 264 */ 265 void normalize(numeral_vector & p); 266 267 /** 268 \brief Divide the coefficients of p by b. 269 This method assumes that every coefficient of p is a multiple of b, and b != 0. 270 */ 271 void div(unsigned sz, numeral * p, numeral const & b); div(numeral_vector & p,numeral const & b)272 void div(numeral_vector & p, numeral const & b) { div(p.size(), p.c_ptr(), b); } 273 274 /** 275 \brief Multiply the coefficients of p by b. 276 277 This method assume b != 0. 278 */ 279 void mul(unsigned sz, numeral * p, numeral const & b); 280 281 /** 282 \brief Multiply the coefficients of p by b. 283 If b == 0, it is equivalent to a reset() 284 */ 285 void mul(numeral_vector & p, numeral const & b); 286 287 /** 288 \brief Similar to div_rem but p1 and p2 must not be q and r. 289 */ 290 void div_rem_core(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, unsigned & d, numeral_vector & q, numeral_vector & r); 291 292 /** 293 \brief If numeral is a field, then 294 return q and r s.t. p1 = q*p2 + r 295 And degree(r) < degree(p2). 296 297 If numeral is not a field, then 298 return q and r s.t. (b_m)^d * p1 = q * p2 + r 299 where b_m is the leading coefficient of p2 and d <= sz1 - sz2 + 1 300 if sz1 >= sz2. 301 302 The output value d is irrelevant if numeral is a field. 303 */ 304 void div_rem(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, unsigned & d, numeral_vector & q, numeral_vector & r); 305 306 /** 307 \see div_rem 308 */ div_rem(unsigned sz1,numeral const * p1,unsigned sz2,numeral const * p2,numeral_vector & q,numeral_vector & r)309 void div_rem(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & q, numeral_vector & r) { 310 unsigned d = 0; 311 div_rem(sz1, p1, sz2, p2, d, q, r); 312 } 313 div_rem(numeral_vector const & p1,numeral_vector const & p2,numeral_vector & q,numeral_vector & r)314 void div_rem(numeral_vector const & p1, numeral_vector const & p2, numeral_vector & q, numeral_vector & r) { 315 div_rem(p1.size(), p1.c_ptr(), p2.size(), p2.c_ptr(), q, r); 316 } 317 318 /** 319 \see div_rem 320 */ 321 void div(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & q); 322 323 /** 324 \see div_rem 325 */ 326 void rem(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, unsigned & d, numeral_vector & r); 327 328 /** 329 \see div_rem 330 */ rem(unsigned sz1,numeral const * p1,unsigned sz2,numeral const * p2,numeral_vector & r)331 void rem(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & r) { 332 unsigned d = 0; 333 rem(sz1, p1, sz2, p2, d, r); 334 } 335 336 /** 337 \brief Signed pseudo-remainder. 338 Alias for rem(sz1, p1, sz2, p2, r); neg(r); 339 */ 340 void srem(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & r); 341 342 /** 343 \brief Return true if p2 divides p1. 344 */ 345 bool divides(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2); divides(numeral_vector const & p1,numeral_vector const & p2)346 bool divides(numeral_vector const & p1, numeral_vector const & p2) { return divides(p1.size(), p1.c_ptr(), p2.size(), p2.c_ptr()); } 347 348 /** 349 \brief Return true if p2 divides p1, and store the quotient in q. 350 */ 351 bool exact_div(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & q); exact_div(numeral_vector const & p1,numeral_vector const & p2,numeral_vector & q)352 bool exact_div(numeral_vector const & p1, numeral_vector const & p2, numeral_vector & q) { 353 return exact_div(p1.size(), p1.c_ptr(), p2.size(), p2.c_ptr(), q); 354 } 355 356 /** 357 \brief Assuming that we can, make the polynomial monic by dividing with the leading coefficient. It 358 puts the leading coefficient into lc, and it's inverse into lc_inv. 359 */ 360 void mk_monic(unsigned sz, numeral * p, numeral & lc, numeral & lc_inv); mk_monic(unsigned sz,numeral * p,numeral & lc)361 void mk_monic(unsigned sz, numeral * p, numeral & lc) { numeral lc_inv; mk_monic(sz, p, lc, lc_inv); m().del(lc_inv); } mk_monic(unsigned sz,numeral * p)362 void mk_monic(unsigned sz, numeral * p) { numeral lc, lc_inv; mk_monic(sz, p, lc, lc_inv); m().del(lc); m().del(lc_inv); } mk_monic(numeral_vector & p)363 void mk_monic(numeral_vector & p) { mk_monic(p.size(), p.c_ptr()); } 364 365 /** 366 \brief g := gcd(p1, p2) 367 If in a field the coefficients don't matter, so we also make sure that D is monic. 368 */ 369 void gcd(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & g); 370 void euclid_gcd(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & g); 371 void subresultant_gcd(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, numeral_vector & g); gcd(numeral_vector const & p1,numeral_vector const & p2,numeral_vector & g)372 void gcd(numeral_vector const & p1, numeral_vector const & p2, numeral_vector & g) { 373 gcd(p1.size(), p1.c_ptr(), p2.size(), p2.c_ptr(), g); 374 } subresultant_gcd(numeral_vector const & p1,numeral_vector const & p2,numeral_vector & g)375 void subresultant_gcd(numeral_vector const & p1, numeral_vector const & p2, numeral_vector & g) { 376 subresultant_gcd(p1.size(), p1.c_ptr(), p2.size(), p2.c_ptr(), g); 377 } 378 379 /** 380 \brief g := square free part of p 381 */ 382 void square_free(unsigned sz, numeral const * p, numeral_vector & g); 383 384 /** 385 \brief Return true if p is a square-free polynomial. 386 */ 387 bool is_square_free(unsigned sz, numeral const * p); 388 389 /** 390 \brief Return true if p is a square-free polynomial. 391 */ is_square_free(numeral_vector const & p)392 bool is_square_free(numeral_vector const & p) { 393 return is_square_free(p.size(), p.c_ptr()); 394 } 395 396 /** 397 \brief Convert a multi-variate polynomial (that is) actually representing an univariate polynomial 398 into a vector of coefficients. 399 */ 400 template<typename polynomial_ref> to_numeral_vector(polynomial_ref const & p,numeral_vector & r)401 void to_numeral_vector(polynomial_ref const & p, numeral_vector & r) { 402 typename polynomial_ref::manager & pm = p.m(); 403 SASSERT(pm.is_univariate(p)); 404 polynomial_ref np(pm); 405 np = pm.normalize(p); 406 unsigned sz = pm.size(p); 407 unsigned deg = pm.total_degree(p); 408 r.reserve(deg+1); 409 for (unsigned i = 0; i <= deg; i++) { 410 m().reset(r[i]); 411 } 412 for (unsigned i = 0; i < sz; i++) { 413 unsigned k = pm.total_degree(pm.get_monomial(p, i)); 414 SASSERT(k <= deg); 415 m().set(r[k], pm.coeff(p, i)); 416 } 417 set_size(deg+1, r); 418 } 419 420 /** 421 \brief Convert a multi-variate polynomial in [x, y1, ..., yn] to a univariate polynomial in just x 422 by removing everything multivariate. 423 */ 424 template<typename polynomial_ref> to_numeral_vector(polynomial_ref const & p,polynomial::var x,numeral_vector & r)425 void to_numeral_vector(polynomial_ref const & p, polynomial::var x, numeral_vector & r) { 426 typename polynomial_ref::manager & pm = p.m(); 427 polynomial_ref np(pm); 428 np = pm.normalize(p); 429 unsigned sz = pm.size(p); 430 unsigned deg = pm.degree(p, x); 431 r.reserve(deg+1); 432 for (unsigned i = 0; i <= deg; i++) { 433 m().reset(r[i]); 434 } 435 for (unsigned i = 0; i < sz; i++) { 436 typename polynomial::monomial * mon = pm.get_monomial(p, i); 437 if (pm.size(mon) == 0) { 438 m().set(r[0], pm.coeff(p, i)); 439 } else if (pm.size(mon) == 1 && pm.get_var(mon, 0) == x) { 440 unsigned m_deg_x = pm.degree(mon, 0); 441 m().set(r[m_deg_x], pm.coeff(p, i)); 442 } 443 } 444 set_size(deg+1, r); 445 } 446 447 /** 448 \brief Extended GCD 449 This method assumes that numeral is a field. 450 It determines U, V, D such that 451 A*U + B*V = D and D is the GCD of A and B. 452 Since in a field the coefficients don't matter, we also make sure that D is monic. 453 */ 454 void ext_gcd(unsigned szA, numeral const * A, unsigned szB, numeral const * B, numeral_vector & U, numeral_vector & V, numeral_vector & D); ext_gcd(numeral_vector const & A,numeral_vector const & B,numeral_vector & U,numeral_vector & V,numeral_vector & D)455 void ext_gcd(numeral_vector const & A, numeral_vector const & B, numeral_vector & U, numeral_vector & V, numeral_vector & D) { 456 ext_gcd(A.size(), A.c_ptr(), B.size(), B.c_ptr(), U, V, D); 457 } 458 459 bool eq(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2); eq(numeral_vector const & p1,numeral_vector const & p2)460 bool eq(numeral_vector const & p1, numeral_vector const & p2) { return eq(p1.size(), p1.c_ptr(), p2.size(), p2.c_ptr()); } 461 462 std::ostream& display(std::ostream & out, unsigned sz, numeral const * p, char const * var_name = "x", bool use_star = false) const; 463 std::ostream& display(std::ostream & out, numeral_vector const & p, char const * var_name = "x") const { return display(out, p.size(), p.c_ptr(), var_name); } display_star(std::ostream & out,unsigned sz,numeral const * p)464 std::ostream& display_star(std::ostream & out, unsigned sz, numeral const * p) { return display(out, sz, p, "x", true); } display_star(std::ostream & out,numeral_vector const & p)465 std::ostream& display_star(std::ostream & out, numeral_vector const & p) { return display_star(out, p.size(), p.c_ptr()); } 466 467 std::ostream& display_smt2(std::ostream & out, unsigned sz, numeral const * p, char const * var_name = "x") const; 468 std::ostream& display_smt2(std::ostream & out, numeral_vector const & p, char const * var_name = "x") const { 469 return display_smt2(out, p.size(), p.c_ptr(), var_name); 470 } 471 }; 472 473 class scoped_set_z { 474 core_manager & m; 475 bool m_modular; 476 core_manager::scoped_numeral m_p; 477 public: scoped_set_z(core_manager & _m)478 scoped_set_z(core_manager & _m):m(_m), m_modular(m.modular()), m_p(m.m()) { m_p = m.p(); m.set_z(); } ~scoped_set_z()479 ~scoped_set_z() { if (m_modular) m.set_zp(m_p); } 480 }; 481 482 class scoped_set_zp { 483 core_manager & m; 484 bool m_modular; 485 core_manager::scoped_numeral m_p; 486 public: scoped_set_zp(core_manager & _m,numeral const & p)487 scoped_set_zp(core_manager & _m, numeral const & p):m(_m), m_modular(m.modular()), m_p(m.m()) { m_p = m.p(); m.set_zp(p); } scoped_set_zp(core_manager & _m,uint64_t p)488 scoped_set_zp(core_manager & _m, uint64_t p):m(_m), m_modular(m.modular()), m_p(m.m()) { m_p = m.p(); m.set_zp(p); } ~scoped_set_zp()489 ~scoped_set_zp() { if (m_modular) m.set_zp(m_p); else m.set_z(); } 490 }; 491 492 class manager; 493 494 typedef core_manager z_manager; 495 typedef core_manager zp_manager; 496 497 typedef z_manager::factors factors; 498 typedef zp_manager::factors zp_factors; 499 500 typedef svector<numeral> numeral_vector; 501 502 class scoped_numeral_vector : public _scoped_numeral_vector<numeral_manager> { 503 public: scoped_numeral_vector(numeral_manager & m)504 scoped_numeral_vector(numeral_manager & m):_scoped_numeral_vector<numeral_manager>(m) {} 505 scoped_numeral_vector(manager & m); 506 }; 507 508 class upolynomial_sequence { 509 numeral_vector m_seq_coeffs; // coefficients of all polynomials in the sequence 510 unsigned_vector m_begins; // start position (in m_seq_coeffs) of each polynomial in the sequence 511 unsigned_vector m_szs; // size of each polynomial in the sequence 512 friend class manager; 513 public: 514 /** 515 \brief Add a new polynomial to the sequence. 516 The contents of p is erased. 517 */ 518 void push(unsigned sz, numeral * p); 519 520 /** 521 \brief Add a new polynomial to the sequence. 522 The contents of p is preserved. 523 */ 524 void push(numeral_manager & m, unsigned sz, numeral const * p); 525 526 /** 527 \brief Return the number of polynomials in the sequence. 528 */ size()529 unsigned size() const { return m_szs.size(); } 530 531 /** 532 \brief Return the vector of coefficients for the i-th polynomial in the sequence. 533 */ coeffs(unsigned i)534 numeral const * coeffs(unsigned i) const { return m_seq_coeffs.c_ptr() + m_begins[i]; } 535 536 /** 537 \brief Return the size of the i-th polynomial in the sequence. 538 */ size(unsigned i)539 unsigned size(unsigned i) const { return m_szs[i]; } 540 }; 541 542 class scoped_upolynomial_sequence : public upolynomial_sequence { 543 manager & m_manager; 544 public: scoped_upolynomial_sequence(manager & m)545 scoped_upolynomial_sequence(manager & m):m_manager(m) {} 546 ~scoped_upolynomial_sequence(); 547 }; 548 549 class manager : public core_manager { 550 numeral_vector m_db_tmp; 551 numeral_vector m_dbab_tmp1; 552 numeral_vector m_dbab_tmp2; 553 numeral_vector m_tr_tmp; 554 numeral_vector m_push_tmp; 555 556 sign sign_of(numeral const & c); 557 struct drs_frame; 558 void pop_top_frame(numeral_vector & p_stack, svector<drs_frame> & frame_stack); 559 void push_child_frames(unsigned sz, numeral const * p, numeral_vector & p_stack, svector<drs_frame> & frame_stack); 560 void add_isolating_interval(svector<drs_frame> const & frame_stack, mpbq_manager & bqm, mpbq_vector & lowers, mpbq_vector & uppers); 561 void add_root(svector<drs_frame> const & frame_stack, mpbq_manager & bqm, mpbq_vector & roots); 562 void drs_isolate_0_1_roots(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers); 563 void drs_isolate_roots(unsigned sz, numeral * p, numeral & U, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers); 564 void drs_isolate_roots(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers); 565 void sqf_nz_isolate_roots(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers); 566 void sturm_seq_core(upolynomial_sequence & seq); 567 enum location { PLUS_INF, MINUS_INF, ZERO, MPBQ }; 568 template<location loc> 569 unsigned sign_variations_at_core(upolynomial_sequence const & seq, mpbq const & b); 570 571 void flip_sign(factors & r); 572 void flip_factor_sign_if_lm_neg(numeral_vector & p, factors & r, unsigned k); 573 void factor_2_sqf_pp(numeral_vector & p, factors & r, unsigned k); 574 bool factor_sqf_pp(numeral_vector & p, factors & r, unsigned k, factor_params const & params); 575 bool factor_core(unsigned sz, numeral const * p, factors & r, factor_params const & params); 576 577 public: manager(reslimit & lim,z_numeral_manager & m)578 manager(reslimit& lim, z_numeral_manager & m):core_manager(lim, m) {} 579 ~manager(); 580 reset(numeral_vector & p)581 void reset(numeral_vector & p) { core_manager::reset(p); } 582 583 void reset(upolynomial_sequence & seq); 584 585 /** 586 \brief Return true if 0 is a root of p. 587 */ has_zero_roots(unsigned sz,numeral const * p)588 bool has_zero_roots(unsigned sz, numeral const * p) { SASSERT(sz > 0); return m().is_zero(p[0]); } 589 590 /** 591 \brief Store in buffer a polynomial that has the same roots of p but the zero roots. 592 We have that: 593 forall u, p(u) = 0 and u != 0 implies buffer(u) = 0 594 forall u, buffer(u) = 0 implies p(u) = 0 595 596 This method assumes p is not the zero polynomial 597 */ 598 void remove_zero_roots(unsigned sz, numeral const * p, numeral_vector & buffer); 599 600 /** 601 \brief Return true if 1/2 is a root of p. 602 */ 603 bool has_one_half_root(unsigned sz, numeral const * p); 604 605 /** 606 \brief Store in buffer a polynomial that has the same roots of p, but a 1/2 root is removed. 607 608 This method assumes that 1/2 is a root of p. 609 */ 610 void remove_one_half_root(unsigned sz, numeral const * p, numeral_vector & buffer); 611 612 /** 613 \brief Return the number of sign changes in the coefficients of p. 614 Zero coefficients are ignored. 615 */ 616 unsigned sign_changes(unsigned sz, numeral const * p); 617 618 /** 619 \brief Return the descartes bound for the number of roots of p in the interval (0, +oo) 620 621 Result: 622 0 - p has no roots in (0,1) 623 1 - p has one root in (0,1) 624 >1 - p has more than one root in (0,1) 625 */ 626 unsigned descartes_bound(unsigned sz, numeral const * p); 627 628 /** 629 \brief Return the descartes bound for the number of roots of p in the interval (0, 1) 630 631 \see descartes_bound 632 */ 633 unsigned descartes_bound_0_1(unsigned sz, numeral const * p); 634 635 /** 636 \brief Return the descartes bound for the number of roots of p in the interval (a, b) 637 638 \see descartes_bound 639 */ 640 unsigned descartes_bound_a_b(unsigned sz, numeral const * p, mpbq_manager & m, mpbq const & a, mpbq const & b); 641 642 /** 643 \brief p(x) := p(x+1) 644 */ 645 void translate(unsigned sz, numeral * p); translate(unsigned sz,numeral const * p,numeral_vector & buffer)646 void translate(unsigned sz, numeral const * p, numeral_vector & buffer) { set(sz, p, buffer); translate(sz, buffer.c_ptr()); } 647 648 /** 649 \brief p(x) := p(x+2^k) 650 */ 651 void translate_k(unsigned sz, numeral * p, unsigned k); translate_k(unsigned sz,numeral const * p,unsigned k,numeral_vector & buffer)652 void translate_k(unsigned sz, numeral const * p, unsigned k, numeral_vector & buffer) { set(sz, p, buffer); translate_k(sz, buffer.c_ptr(), k); } 653 654 /** 655 \brief p(x) := p(x+c) 656 */ 657 void translate_z(unsigned sz, numeral * p, numeral const & c); translate_z(unsigned sz,numeral const * p,numeral const & c,numeral_vector & buffer)658 void translate_z(unsigned sz, numeral const * p, numeral const & c, numeral_vector & buffer) { set(sz, p, buffer); translate_z(sz, buffer.c_ptr(), c); } 659 660 /** 661 \brief p(x) := p(x+b) where b = c/2^k 662 buffer := (2^k)^n * p(x + c/(2^k)) 663 */ 664 void translate_bq(unsigned sz, numeral * p, mpbq const & b); translate_bq(unsigned sz,numeral const * p,mpbq const & b,numeral_vector & buffer)665 void translate_bq(unsigned sz, numeral const * p, mpbq const & b, numeral_vector & buffer) { set(sz, p, buffer); translate_bq(sz, buffer.c_ptr(), b); } 666 667 /** 668 \brief p(x) := p(x+b) where b = c/d 669 buffer := d^n * p(x + c/d) 670 */ 671 void translate_q(unsigned sz, numeral * p, mpq const & b); translate_q(unsigned sz,numeral const * p,mpq const & b,numeral_vector & buffer)672 void translate_q(unsigned sz, numeral const * p, mpq const & b, numeral_vector & buffer) { set(sz, p, buffer); translate_q(sz, buffer.c_ptr(), b); } 673 674 /** 675 \brief p(x) := 2^n*p(x/2) where n = sz-1 676 */ 677 void compose_2n_p_x_div_2(unsigned sz, numeral * p); 678 679 /** 680 \brief p(x) := (2^k)^n * p(x/(2^k)) 681 */ 682 void compose_2kn_p_x_div_2k(unsigned sz, numeral * p, unsigned k); 683 684 /** 685 \brief p(x) := p(2^k * x) 686 687 If u is a root of old(p), then u/2^k is a root of p 688 */ 689 void compose_p_2k_x(unsigned sz, numeral * p, unsigned k); 690 691 /** 692 \brief p(x) := p(b * x) 693 694 If u is a root of old(p), then u/b is a root of p 695 */ 696 void compose_p_b_x(unsigned sz, numeral * p, numeral const & b); 697 698 /** 699 \brief p(x) := p(b * x) 700 701 If u is a root of old(p), then u/b is a root of p 702 703 Let b be of the form c/(2^k), then this operation is equivalent to: 704 (2^k)^n*p(c*x/(2^k)) 705 706 Let old(p) be of the form: 707 a_n * x^n + a_{n-1}*x^{n-1} + ... + a_1 * x + a_0 708 709 Then p is of the form: 710 a_n * c^n * x^n + a_{n-1} * c^{n-1} * 2^k * x^{n-1} + ... + a_1 * c * (2^k)^(n-1) * x + a_0 711 */ 712 void compose_p_b_x(unsigned sz, numeral * p, mpbq const & b); 713 714 /** 715 \brief p(x) := p(q*x) 716 */ 717 void compose_p_q_x(unsigned sz, numeral * p, mpq const & q); 718 719 /** 720 \brief p(x) := a^n * p(x/a) 721 */ 722 void compose_an_p_x_div_a(unsigned sz, numeral * p, numeral const & a); 723 724 /** 725 \brief p(x) := p(-x) 726 */ 727 void p_minus_x(unsigned sz, numeral * p); 728 729 /** 730 \brief p(x) := x^n * p(1/x) 731 */ 732 void p_1_div_x(unsigned sz, numeral * p); 733 734 /** 735 \brief Evaluate the sign of p(b) 736 */ 737 sign eval_sign_at(unsigned sz, numeral const * p, mpbq const & b); 738 739 /** 740 \brief Evaluate the sign of p(b) 741 */ 742 sign eval_sign_at(unsigned sz, numeral const * p, mpq const & b); 743 744 /** 745 \brief Evaluate the sign of p(b) 746 */ 747 sign eval_sign_at(unsigned sz, numeral const * p, mpz const & b); 748 749 /** 750 \brief Evaluate the sign of p(0) 751 */ 752 sign eval_sign_at_zero(unsigned sz, numeral const * p); 753 754 /** 755 \brief Evaluate the sign of p(+oo) 756 */ 757 sign eval_sign_at_plus_inf(unsigned sz, numeral const * p); 758 759 /** 760 \brief Evaluate the sign of p(-oo) 761 */ 762 sign eval_sign_at_minus_inf(unsigned sz, numeral const * p); 763 764 /** 765 \brief Evaluate the sign variations in the polynomial sequence at -oo 766 */ 767 unsigned sign_variations_at_minus_inf(upolynomial_sequence const & seq); 768 769 /** 770 \brief Evaluate the sign variations in the polynomial sequence at +oo 771 */ 772 unsigned sign_variations_at_plus_inf(upolynomial_sequence const & seq); 773 774 /** 775 \brief Evaluate the sign variations in the polynomial sequence at 0 776 */ 777 unsigned sign_variations_at_zero(upolynomial_sequence const & seq); 778 779 /** 780 \brief Evaluate the sign variations in the polynomial sequence at b 781 */ 782 unsigned sign_variations_at(upolynomial_sequence const & seq, mpbq const & b); 783 784 /** 785 \brief Return an upper bound U for all roots of p. 786 U is a positive value. 787 We have that if u is a root of p, then |u| < U 788 */ 789 void root_upper_bound(unsigned sz, numeral const * p, numeral & U); 790 791 unsigned knuth_positive_root_upper_bound(unsigned sz, numeral const * p); 792 unsigned knuth_negative_root_upper_bound(unsigned sz, numeral const * p); 793 794 /** 795 \brief Return k s.t. for any nonzero root alpha of p(x): 796 |alpha| > 1/2^k 797 */ 798 unsigned nonzero_root_lower_bound(unsigned sz, numeral const * p); 799 800 /** 801 \brief Isolate roots of a square free polynomial p. 802 The result is stored in three vectors: roots, lowers and uppers. 803 The vector roots contains actual roots of p. 804 The vectors lowers and uppers have the same size, and 805 For all i in [0, lowers.size()), we have that there is only and only one root of p in the interval (lowers[i], uppers[i]). 806 Every root of p in roots or in an interval (lowers[i], uppers[i]) 807 808 The total number of roots of p is roots.size() + lowers.size() 809 810 \pre p is not the zero polynomial, that is, sz > 0 811 */ 812 void sqf_isolate_roots(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers); 813 814 /** 815 \brief Isolate roots of an arbitrary polynomial p. 816 817 \see sqf_isolate_roots. 818 819 \pre p is not the zero polynomial, that is, sz > 0 820 */ 821 void isolate_roots(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers); 822 823 void drs_isolate_roots(unsigned sz, numeral * p, unsigned neg_k, unsigned pos_k, 824 mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers); 825 826 void sturm_isolate_roots_core(unsigned sz, numeral * p, unsigned neg_k, unsigned pos_k, 827 mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers); 828 829 void sturm_isolate_roots(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers); 830 831 /** 832 \brief Compute the sturm sequence for p1 and p2. 833 */ 834 void sturm_seq(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, upolynomial_sequence & seq); 835 836 /** 837 \brief Compute the sturm sequence for p and p'. 838 */ 839 void sturm_seq(unsigned sz, numeral const * p, upolynomial_sequence & seq); 840 841 /** 842 \brief Compute the sturm tarski sequence for p1 and p1'*p2. 843 */ 844 void sturm_tarski_seq(unsigned sz1, numeral const * p1, unsigned sz2, numeral const * p2, upolynomial_sequence & seq); 845 846 /** 847 \brief Compute the Fourier sequence for p. 848 */ 849 void fourier_seq(unsigned sz, numeral const * p, upolynomial_sequence & seq); 850 851 /** 852 \brief Convert an isolating interval into a refinable one. 853 See comments in upolynomial.cpp. 854 */ 855 bool isolating2refinable(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq & a, mpbq & b); 856 857 // 858 // Interval refinement procedures 859 // They all assume p is square free and (a, b) is a refinable isolating interval. 860 // 861 // Return TRUE, if interval was squeezed, and new interval is stored in (a,b). 862 // Return FALSE, if the actual root was found, it is stored in a. 863 // 864 // See upolynomial.cpp for additional comments 865 bool refine_core(unsigned sz, numeral const * p, sign sign_a, mpbq_manager & bqm, mpbq & a, mpbq & b); 866 867 bool refine(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq & a, mpbq & b); 868 869 bool refine_core(unsigned sz, numeral const * p, sign sign_a, mpbq_manager & bqm, mpbq & a, mpbq & b, unsigned prec_k); 870 871 bool refine(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq & a, mpbq & b, unsigned prec_k); 872 ///////////////////// 873 874 /** 875 \brief Convert a isolating (refinable) rational interval into a 876 isolating refinable binary rational interval. 877 878 Return TRUE, if interval was found and the result is stored in (c, d). 879 Return FALSE, if the actual root was found, it is stored in c. 880 */ 881 bool convert_q2bq_interval(unsigned sz, numeral const * p, mpq const & a, mpq const & b, mpbq_manager & bqm, mpbq & c, mpbq & d); 882 883 /** 884 \brief Given a polynomial p, and a lower bound l. Return 885 the root id i. That is, the first root u > l is the i-th root of p. 886 */ 887 unsigned get_root_id(unsigned sz, numeral const * p, mpbq const & l); 888 889 /** 890 \brief Make sure that isolating interval (a, b) for p does not contain zero. 891 892 Return TRUE, if updated (a, b) does not contain zero. 893 Return FALSE, if zero is a root of p 894 */ 895 bool normalize_interval_core(unsigned sz, numeral const * p, int sign_a, mpbq_manager & m, mpbq & a, mpbq & b); 896 897 /** 898 \brief Similar to normalize_interval_core, but sign_a does not need to be provided. 899 */ 900 bool normalize_interval(unsigned sz, numeral const * p, mpbq_manager & m, mpbq & a, mpbq & b); 901 902 /** 903 \brief Return true if all irreducible factors were found. 904 That is, if the result if false, there is no guarantee that the factors in r are irreducible. 905 This can happen when limits (e.g., on the search space size) are set in params. 906 */ 907 bool factor(unsigned sz, numeral const * p, factors & r, factor_params const & params = factor_params()); 908 bool factor(numeral_vector const & p, factors & r, factor_params const & params = factor_params()) { return factor(p.size(), p.c_ptr(), r, params); } 909 910 std::ostream& display(std::ostream & out, unsigned sz, numeral const * p, char const * var_name = "x", bool use_star = false) const { 911 return core_manager::display(out, sz, p, var_name); 912 } 913 std::ostream& display(std::ostream & out, numeral_vector const & p, char const * var_name = "x") const { 914 return core_manager::display(out, p, var_name); 915 } 916 std::ostream& display(std::ostream & out, upolynomial_sequence const & seq, char const * var_name = "x") const; 917 }; 918 919 }; 920 921