mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-08-22 02:57:50 +00:00
Code Activity

Centralize and document TRACE tags using X-macros (#7657)

Browse source 
* Introduce X-macro-based trace tag definition
- Created trace_tags.def to centralize TRACE tag definitions
- Each tag includes a symbolic name and description
- Set up enum class TraceTag for type-safe usage in TRACE macros

* Add script to generate Markdown documentation from trace_tags.def
- Python script parses trace_tags.def and outputs trace_tags.md

* Refactor TRACE_NEW to prepend TraceTag and pass enum to is_trace_enabled

* trace: improve trace tag handling system with hierarchical tagging

- Introduce hierarchical tag-class structure: enabling a tag class activates all child tags
- Unify TRACE, STRACE, SCTRACE, and CTRACE under enum TraceTag
- Implement initial version of trace_tag.def using X(tag, tag_class, description)
  (class names and descriptions to be refined in a future update)

* trace: replace all string-based TRACE tags with enum TraceTag
- Migrated all TRACE, STRACE, SCTRACE, and CTRACE macros to use enum TraceTag values instead of raw string literals

* trace : add cstring header

* trace : Add Markdown documentation generation from trace_tags.def via mk_api_doc.py

* trace : rename macro parameter 'class' to 'tag_class' and remove Unicode comment in trace_tags.h.

* trace : Add TODO comment for future implementation of tag_class activation

* trace : Disable code related to tag_class until implementation is ready (#7663).
This commit is contained in:
LeeYoungJoon 2025-05-28 22:31:25 +09:00 • committed by GitHub
parent d766292dab
commit 0a93ff515d
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
583 changed files with 8698 additions and 7299 deletions
Download patch file Download diff file Expand all files Collapse all files

162
src/math/polynomial/algebraic_numbers.cpp
View file

@ -287,7 +287,7 @@ namespace algebraic_numbers {
return true;
if (a.to_algebraic()->m_not_rational)
return false; // we know for sure a is not a rational
TRACE("algebraic_bug", tout << "is_rational(a):\n"; display_root(tout, a); tout << "\n"; display_interval(tout, a); tout << "\n";);
TRACE(algebraic_bug, tout << "is_rational(a):\n"; display_root(tout, a); tout << "\n"; display_interval(tout, a); tout << "\n";);
algebraic_cell * c = a.to_algebraic();
save_intervals saved_a(*this, a);
mpz & a_n = c->m_p[c->m_p_sz - 1];
@ -299,7 +299,7 @@ namespace algebraic_numbers {
unsigned k = qm().log2(abs_a_n);
k++;
TRACE("algebraic_bug", tout << "abs(an): " << qm().to_string(abs_a_n) << ", k: " << k << "\n";);
TRACE(algebraic_bug, tout << "abs(an): " << qm().to_string(abs_a_n) << ", k: " << k << "\n";);
// make sure the isolating interval size is less than 1/2^k
if (!refine_until_prec(a, k)) {
@ -307,7 +307,7 @@ namespace algebraic_numbers {
return true;
}
TRACE("algebraic_bug", tout << "interval after refinement: "; display_interval(tout, a); tout << "\n";);
TRACE(algebraic_bug, tout << "interval after refinement: "; display_interval(tout, a); tout << "\n";);
// Find unique candidate rational in the isolating interval
scoped_mpbq a_n_lower(bqm());
@ -507,7 +507,7 @@ namespace algebraic_numbers {
}
SASSERT(acell_inv(*a.to_algebraic()));
}
TRACE("algebraic", tout << "a: "; display_root(tout, a); tout << "\n";);
TRACE(algebraic, tout << "a: "; display_root(tout, a); tout << "\n";);
}
void set(numeral & a, numeral const & b) {
@ -552,7 +552,7 @@ namespace algebraic_numbers {
scoped_upoly & up_sqf = m_isolate_tmp3;
up_sqf.reset();
upm().square_free(up.size(), up.data(), up_sqf);
TRACE("algebraic", upm().display(tout, up_sqf.size(), up_sqf.data()); tout << "\n";);
TRACE(algebraic, upm().display(tout, up_sqf.size(), up_sqf.data()); tout << "\n";);
r.push_back(up_sqf, 1);
return false;
}
@ -582,7 +582,7 @@ namespace algebraic_numbers {
(void)c_lt_a;
// (a <= b & b <= c) => a <= c
// b < a or c < b or !(c < a)
CTRACE("algebraic_bug",
CTRACE(algebraic_bug,
(!b_lt_a && !c_lt_b && c_lt_a),
display_root(tout << "a ", a) << "\n";
display_root(tout << "b ", b) << "\n";
@ -601,7 +601,7 @@ namespace algebraic_numbers {
}
void isolate_roots(scoped_upoly const & up, numeral_vector & roots) {
TRACE("algebraic", upm().display(tout, up); tout << "\n";);
TRACE(algebraic, upm().display(tout, up); tout << "\n";);
if (up.empty())
return; // ignore the zero polynomial
factors & fs = m_isolate_factors;
@ -622,12 +622,12 @@ namespace algebraic_numbers {
upolynomial::numeral_vector const & f = fs[i];
// polynomial f contains the non zero roots
unsigned d = upm().degree(f);
TRACE("algebraic", tout << "factor " << i << " degree: " << d << "\n";);
TRACE(algebraic, tout << "factor " << i << " degree: " << d << "\n";);
if (d == 0)
continue; // found all roots of f
scoped_mpq r(qm());
if (d == 1) {
TRACE("algebraic", tout << "linear polynomial...\n";);
TRACE(algebraic, tout << "linear polynomial...\n";);
// f is a linear polynomial ax + b
// set r <- -b/a
qm().set(r, f[0]);
@ -640,7 +640,7 @@ namespace algebraic_numbers {
upm().sqf_isolate_roots(f.size(), f.data(), bqm(), m_isolate_roots, m_isolate_lowers, m_isolate_uppers);
// collect rational/basic roots
unsigned sz = m_isolate_roots.size();
TRACE("algebraic", tout << "isolated roots: " << sz << "\n";);
TRACE(algebraic, tout << "isolated roots: " << sz << "\n";);
for (unsigned i = 0; i < sz; i++) {
to_mpq(qm(), m_isolate_roots[i], r);
roots.push_back(numeral(mk_basic_cell(r)));
@ -670,7 +670,7 @@ namespace algebraic_numbers {
void isolate_roots(polynomial_ref const & p, numeral_vector & roots) {
SASSERT(is_univariate(p));
TRACE("algebraic", tout << "isolating roots of: " << p << "\n";);
TRACE(algebraic, tout << "isolating roots of: " << p << "\n";);
if (::is_zero(p))
return; // ignore the zero polynomial
scoped_upoly & up = m_isolate_tmp1;
@ -688,7 +688,7 @@ namespace algebraic_numbers {
scoped_numeral_vector roots(m_wrapper);
isolate_roots(up, roots);
unsigned num_roots = roots.size();
TRACE("algebraic", tout << "num-roots: " << num_roots << "\n";
TRACE(algebraic, tout << "num-roots: " << num_roots << "\n";
for (unsigned i = 0; i < num_roots; i++) {
display_interval(tout, roots[i]);
tout << "\n";
@ -701,7 +701,7 @@ namespace algebraic_numbers {
void mk_root(polynomial_ref const & p, unsigned i, numeral & r) {
SASSERT(i != 0);
SASSERT(is_univariate(p));
TRACE("algebraic", tout << "isolating roots of: " << p << "\n";);
TRACE(algebraic, tout << "isolating roots of: " << p << "\n";);
scoped_upoly & up = m_isolate_tmp1;
upm().to_numeral_vector(p, up);
mk_root(up, i, r);
@ -711,7 +711,7 @@ namespace algebraic_numbers {
SASSERT(i != 0);
scoped_upoly & up = m_isolate_tmp1;
sexpr2upolynomial(upm(), p, up);
TRACE("algebraic", tout << "mk_root " << i << "\n"; upm().display(tout, up); tout << "\n";);
TRACE(algebraic, tout << "mk_root " << i << "\n"; upm().display(tout, up); tout << "\n";);
mk_root(up, i, r);
}
@ -985,13 +985,13 @@ namespace algebraic_numbers {
\pre lV and uV are the sign variations (in seq) for r_i.lower() and r_i.upper()
*/
void set_core(numeral & c, scoped_upoly & p, mpbqi & r_i, upolynomial::scoped_upolynomial_sequence & seq, int lV, int uV, bool minimal) {
TRACE("algebraic", tout << "set_core p: "; upm().display(tout, p); tout << "\n";);
TRACE(algebraic, tout << "set_core p: "; upm().display(tout, p); tout << "\n";);
if (bqim().contains_zero(r_i)) {
if (upm().has_zero_roots(p.size(), p.data())) {
// zero is a root of p, and r_i is an isolating interval containing zero,
// then c is zero
reset(c);
TRACE("algebraic", tout << "resetting\nresult: "; display_root(tout, c); tout << "\n";);
TRACE(algebraic, tout << "resetting\nresult: "; display_root(tout, c); tout << "\n";);
return;
}
int zV = upm().sign_variations_at_zero(seq);
@ -1024,7 +1024,7 @@ namespace algebraic_numbers {
set(c, r);
}
else {
TRACE("algebraic", tout << "set_core...\n";);
TRACE(algebraic, tout << "set_core...\n";);
set(c, nz_p.size(), nz_p.data(), r_i.lower(), r_i.upper(), minimal);
}
}
@ -1052,14 +1052,14 @@ namespace algebraic_numbers {
scoped_upoly p(upm());
scoped_upoly f(upm());
mk_poly(cell_a, cell_b, p);
TRACE("anum_mk_binary", tout << "a: "; display_root(tout, a); tout << "\nb: "; display_root(tout, b); tout << "\np: ";
TRACE(anum_mk_binary, tout << "a: "; display_root(tout, a); tout << "\nb: "; display_root(tout, b); tout << "\np: ";
upm().display(tout, p); tout << "\n";);
factors fs(upm());
bool full_fact = factor(p, fs);
unsigned num_fs = fs.distinct_factors();
scoped_ptr_vector<typename upolynomial::scoped_upolynomial_sequence> seqs;
for (unsigned i = 0; i < num_fs; i++) {
TRACE("anum_mk_binary", tout << "factor " << i << "\n"; upm().display(tout, fs[i]); tout << "\n";);
TRACE(anum_mk_binary, tout << "factor " << i << "\n"; upm().display(tout, fs[i]); tout << "\n";);
typename upolynomial::scoped_upolynomial_sequence * seq = alloc(typename upolynomial::scoped_upolynomial_sequence, upm());
upm().sturm_seq(fs[i].size(), fs[i].data(), *seq);
seqs.push_back(seq);
@ -1082,11 +1082,11 @@ namespace algebraic_numbers {
for (unsigned i = 0; i < num_fs; i++) {
if (seqs[i] == nullptr)
continue; // sequence was discarded because it does not contain the root.
TRACE("anum_mk_binary", tout << "sequence " << i << "\n"; upm().display(tout, *(seqs[i])); tout << "\n";);
TRACE(anum_mk_binary, tout << "sequence " << i << "\n"; upm().display(tout, *(seqs[i])); tout << "\n";);
int lV = upm().sign_variations_at(*(seqs[i]), r_i.lower());
int uV = upm().sign_variations_at(*(seqs[i]), r_i.upper());
int V = lV - uV;
TRACE("algebraic", tout << "r_i: "; bqim().display(tout, r_i); tout << "\n";
TRACE(algebraic, tout << "r_i: "; bqim().display(tout, r_i); tout << "\n";
tout << "lV: " << lV << ", uV: " << uV << "\n";
tout << "a.m_interval: "; bqim().display(tout, cell_a->m_interval); tout << "\n";
tout << "b.m_interval: "; bqim().display(tout, cell_b->m_interval); tout << "\n";
@ -1108,7 +1108,7 @@ namespace algebraic_numbers {
if (num_rem == 1 && target_i != UINT_MAX) {
// found isolating interval
TRACE("anum_mk_binary", tout << "target_i: " << target_i << "\n";);
TRACE(anum_mk_binary, tout << "target_i: " << target_i << "\n";);
saved_a.restore_if_too_small();
saved_b.restore_if_too_small();
upm().set(fs[target_i].size(), fs[target_i].data(), f);
@ -1163,7 +1163,7 @@ namespace algebraic_numbers {
int lV = upm().sign_variations_at(*(seqs[i]), r_i.lower());
int uV = upm().sign_variations_at(*(seqs[i]), r_i.upper());
int V = lV - uV;
TRACE("algebraic", tout << "r_i: "; bqim().display(tout, r_i); tout << "\n";
TRACE(algebraic, tout << "r_i: "; bqim().display(tout, r_i); tout << "\n";
tout << "lV: " << lV << ", uV: " << uV << "\n";
tout << "a.m_interval: "; bqim().display(tout, cell_a->m_interval); tout << "\n";
);
@ -1319,7 +1319,7 @@ namespace algebraic_numbers {
if (qm().root(a_val, k, r_a_val)) {
// the result is rational
TRACE("root_core", tout << "r_a_val: " << r_a_val << " a_val: "; qm().display(tout, a_val); tout << "\n";);
TRACE(root_core, tout << "r_a_val: " << r_a_val << " a_val: "; qm().display(tout, a_val); tout << "\n";);
set(b, r_a_val);
return;
}
@ -1360,7 +1360,7 @@ namespace algebraic_numbers {
bqm().add(upper, mpz(1), upper); // make sure (a_val)^{1/k} < upper
}
SASSERT(bqm().lt(lower, upper));
TRACE("algebraic", tout << "root_core:\n"; upm().display(tout, p.size(), p.data()); tout << "\n";);
TRACE(algebraic, tout << "root_core:\n"; upm().display(tout, p.size(), p.data()); tout << "\n";);
// p is not necessarily a minimal polynomial.
// So, we set the m_minimal flag to false. TODO: try to factor.
set(b, p.size(), p.data(), lower, upper, false);
@ -1434,7 +1434,7 @@ namespace algebraic_numbers {
template<bool IsAdd>
void add(algebraic_cell * a, basic_cell * b, numeral & c) {
TRACE("algebraic", tout << "adding algebraic and basic cells:\n";
TRACE(algebraic, tout << "adding algebraic and basic cells:\n";
tout << "a: "; upm().display(tout, a->m_p_sz, a->m_p); tout << " "; bqim().display(tout, a->m_interval); tout << "\n";
tout << "b: "; qm().display(tout, b->m_value); tout << "\n";);
scoped_mpq nbv(qm());
@ -1457,18 +1457,18 @@ namespace algebraic_numbers {
scoped_mpq iu(qm());
to_mpq(qm(), i.lower(), il);
to_mpq(qm(), i.upper(), iu);
TRACE("algebraic",
TRACE(algebraic,
tout << "nbv: " << nbv << "\n";
tout << "il: " << il << ", iu: " << iu << "\n";);
qm().add(il, nbv, il);
qm().add(iu, nbv, iu);
// (il, iu) is an isolating refinable (rational) interval for the new polynomial.
if (!upm().convert_q2bq_interval(m_add_tmp.size(), m_add_tmp.data(), il, iu, bqm(), l, u)) {
TRACE("algebraic", tout << "conversion failed\n");
TRACE(algebraic, tout << "conversion failed\n");
}
}
TRACE("algebraic",
TRACE(algebraic,
upm().display(tout, m_add_tmp.size(), m_add_tmp.data());
tout << ", l: " << l << ", u: " << u << "\n";
tout << "l_sign: " << upm().eval_sign_at(m_add_tmp.size(), m_add_tmp.data(), l) << "\n";
@ -1541,7 +1541,7 @@ namespace algebraic_numbers {
}
void mul(algebraic_cell * a, basic_cell * b, numeral & c) {
TRACE("algebraic", tout << "mult algebraic and basic cells:\n";
TRACE(algebraic, tout << "mult algebraic and basic cells:\n";
tout << "a: "; upm().display(tout, a->m_p_sz, a->m_p); tout << " "; bqim().display(tout, a->m_interval); tout << "\n";
tout << "b: "; qm().display(tout, b->m_value); tout << "\n";);
SASSERT(upm().eval_sign_at(a->m_p_sz, a->m_p, lower(a)) == -upm().eval_sign_at(a->m_p_sz, a->m_p, upper(a)));
@ -1568,7 +1568,7 @@ namespace algebraic_numbers {
scoped_mpq iu(qm());
to_mpq(qm(), i.lower(), il);
to_mpq(qm(), i.upper(), iu);
TRACE("algebraic",
TRACE(algebraic,
tout << "nbv: " << nbv << "\n";
tout << "il: " << il << ", iu: " << iu << "\n";);
qm().mul(il, nbv, il);
@ -1577,10 +1577,10 @@ namespace algebraic_numbers {
qm().swap(il, iu);
// (il, iu) is an isolating refinable (rational) interval for the new polynomial.
if (!upm().convert_q2bq_interval(mulp.size(), mulp.data(), il, iu, bqm(), l, u)) {
TRACE("algebraic", tout << "conversion failed\n");
TRACE(algebraic, tout << "conversion failed\n");
}
}
TRACE("algebraic",
TRACE(algebraic,
upm().display(tout, mulp.size(), mulp.data());
tout << ", l: " << l << ", u: " << u << "\n";
tout << "l_sign: " << upm().eval_sign_at(mulp.size(), mulp.data(), l) << "\n";
@ -1679,7 +1679,7 @@ namespace algebraic_numbers {
qm().inv(a.to_basic()->m_value);
}
else {
TRACE("algebraic_bug", tout << "before inv: "; display_root(tout, a); tout << "\n"; display_interval(tout, a); tout << "\n";);
TRACE(algebraic_bug, tout << "before inv: "; display_root(tout, a); tout << "\n"; display_interval(tout, a); tout << "\n";);
algebraic_cell * cell_a = a.to_algebraic();
upm().p_1_div_x(cell_a->m_p_sz, cell_a->m_p);
// convert binary rational bounds into rational bounds
@ -1690,13 +1690,13 @@ namespace algebraic_numbers {
qm().inv(inv_lower);
qm().inv(inv_upper);
qm().swap(inv_lower, inv_upper);
TRACE("algebraic_bug", tout << "inv new_bounds: " << qm().to_string(inv_lower) << ", " << qm().to_string(inv_upper) << "\n";);
TRACE(algebraic_bug, tout << "inv new_bounds: " << qm().to_string(inv_lower) << ", " << qm().to_string(inv_upper) << "\n";);
// convert isolating interval back as a binary rational bound
if (!upm().convert_q2bq_interval(cell_a->m_p_sz, cell_a->m_p, inv_lower, inv_upper, bqm(), lower(cell_a), upper(cell_a))) {
TRACE("algebraic_bug", tout << "root isolation failed\n");
TRACE(algebraic_bug, tout << "root isolation failed\n");
throw algebraic_exception("inversion of algebraic number failed");
}
TRACE("algebraic_bug", tout << "after inv: "; display_root(tout, a); tout << "\n"; display_interval(tout, a); tout << "\n";);
TRACE(algebraic_bug, tout << "after inv: "; display_root(tout, a); tout << "\n"; display_interval(tout, a); tout << "\n";);
update_sign_lower(cell_a);
SASSERT(acell_inv(*cell_a));
}
@ -1792,7 +1792,7 @@ namespace algebraic_numbers {
return sign_zero;
}
TRACE("algebraic", tout << "comparing\n";
TRACE(algebraic, tout << "comparing\n";
tout << "a: "; upm().display(tout, cell_a->m_p_sz, cell_a->m_p); tout << "\n"; bqim().display(tout, cell_a->m_interval);
tout << "\ncell_a->m_minimal: " << cell_a->m_minimal << "\n";
tout << "b: "; upm().display(tout, cell_b->m_p_sz, cell_b->m_p); tout << "\n"; bqim().display(tout, cell_b->m_interval);
@ -1893,7 +1893,7 @@ namespace algebraic_numbers {
unsigned V1 = upm().sign_variations_at(seq, a_lower);
unsigned V2 = upm().sign_variations_at(seq, a_upper);
int V = V1 - V2;
TRACE("algebraic", tout << "comparing using sturm\n";
TRACE(algebraic, tout << "comparing using sturm\n";
display_interval(tout, a) << "\n";
display_interval(tout, b) << "\n";
tout << "V: " << V << " V1 " << V1 << " V2 " << V2
@ -1943,7 +1943,7 @@ namespace algebraic_numbers {
}
::sign compare(numeral & a, numeral & b) {
TRACE("algebraic", tout << "comparing: "; display_interval(tout, a); tout << " "; display_interval(tout, b); tout << "\n";);
TRACE(algebraic, tout << "comparing: "; display_interval(tout, a); tout << " "; display_interval(tout, b); tout << "\n";);
if (a.is_basic()) {
if (b.is_basic())
return compare(basic_value(a), basic_value(b));
@ -2058,7 +2058,7 @@ namespace algebraic_numbers {
mpq const & operator()(polynomial::var x) const override {
anum const & v = m_x2v(x);
SASSERT(v.is_basic());
TRACE("var2basic", tout << "getting value of x" << x << " -> " << m().to_string(m_imp.basic_value(v)) << "\n";);
TRACE(var2basic, tout << "getting value of x" << x << " -> " << m().to_string(m_imp.basic_value(v)) << "\n";);
return m_imp.basic_value(v);
}
};
@ -2082,7 +2082,7 @@ namespace algebraic_numbers {
polynomial::var_vector m_eval_sign_vars;
sign eval_sign_at(polynomial_ref const & p, polynomial::var2anum const & x2v) {
polynomial::manager & ext_pm = p.m();
TRACE("anum_eval_sign", tout << "evaluating sign of: " << p << "\n";);
TRACE(anum_eval_sign, tout << "evaluating sign of: " << p << "\n";);
while (true) {
bool restart = false;
// Optimistic: maybe x2v contains only rational values
@ -2090,7 +2090,7 @@ namespace algebraic_numbers {
opt_var2basic x2v_basic(*this, x2v);
scoped_mpq r(qm());
ext_pm.eval(p, x2v_basic, r);
TRACE("anum_eval_sign", tout << "all variables are assigned to rationals, value of p: " << r << "\n";);
TRACE(anum_eval_sign, tout << "all variables are assigned to rationals, value of p: " << r << "\n";);
return ::to_sign(qm().sign(r));
}
catch (const opt_var2basic::failed &) {
@ -2101,7 +2101,7 @@ namespace algebraic_numbers {
polynomial_ref p_prime(ext_pm);
var2basic x2v_basic(*this, x2v);
p_prime = ext_pm.substitute(p, x2v_basic);
TRACE("anum_eval_sign", tout << "p after eliminating rationals: " << p_prime << "\n";);
TRACE(anum_eval_sign, tout << "p after eliminating rationals: " << p_prime << "\n";);
if (ext_pm.is_zero(p_prime)) {
// polynomial vanished after substituting rational values.
@ -2125,7 +2125,7 @@ namespace algebraic_numbers {
while (true) {
checkpoint();
ext_pm.eval(p_prime, x2v_interval, ri);
TRACE("anum_eval_sign", tout << "evaluating using intervals: " << ri << "\n";);
TRACE(anum_eval_sign, tout << "evaluating using intervals: " << ri << "\n";);
if (!bqim().contains_zero(ri)) {
return bqim().is_pos(ri) ? sign_pos : sign_neg;
}
@ -2144,7 +2144,7 @@ namespace algebraic_numbers {
restart = true;
break;
}
TRACE("anum_eval_sign", tout << "refined algebraic interval\n";);
TRACE(anum_eval_sign, tout << "refined algebraic interval\n";);
SASSERT(!v.is_basic());
refined = true;
}
@ -2158,7 +2158,7 @@ namespace algebraic_numbers {
if (restart) {
// Some non-basic value became basic.
// So, restarting the whole process
TRACE("anum_eval_sign", tout << "restarting some algebraic_cell became basic\n";);
TRACE(anum_eval_sign, tout << "restarting some algebraic_cell became basic\n";);
continue;
}
@ -2174,7 +2174,7 @@ namespace algebraic_numbers {
// Evaluating sign using algebraic arithmetic
scoped_anum ra(m_wrapper);
ext_pm.eval(p_prime, x2v, ra);
TRACE("anum_eval_sign", tout << "value of p as algebraic number " << ra << "\n";);
TRACE(anum_eval_sign, tout << "value of p as algebraic number " << ra << "\n";);
if (is_zero(ra))
return 0;
return is_pos(ra) ? 1 : -1;
@ -2234,7 +2234,7 @@ namespace algebraic_numbers {
scoped_upoly & _R = m_eval_sign_tmp;
upm().to_numeral_vector(R, _R);
unsigned k = upm().nonzero_root_lower_bound(_R.size(), _R.data());
TRACE("anum_eval_sign", tout << "R: " << R << "\nk: " << k << "\nri: "<< ri << "\n";);
TRACE(anum_eval_sign, tout << "R: " << R << "\nk: " << k << "\nri: "<< ri << "\n";);
scoped_mpbq mL(bqm()), L(bqm());
bqm().set(mL, -1);
bqm().set(L, 1);
@ -2261,7 +2261,7 @@ namespace algebraic_numbers {
while (!restart) {
checkpoint();
ext_pm.eval(p_prime, x2v_interval, ri);
TRACE("anum_eval_sign", tout << "evaluating using intervals: " << ri << "\n";
TRACE(anum_eval_sign, tout << "evaluating using intervals: " << ri << "\n";
tout << "zero lower bound is: " << L << "\n";);
if (!bqim().contains_zero(ri)) {
return bqim().is_pos(ri) ? sign_pos : sign_neg;
@ -2279,7 +2279,7 @@ namespace algebraic_numbers {
restart = true;
break;
}
TRACE("anum_eval_sign", tout << "refined algebraic interval\n";);
TRACE(anum_eval_sign, tout << "refined algebraic interval\n";);
SASSERT(!v.is_basic());
}
}
@ -2332,7 +2332,7 @@ namespace algebraic_numbers {
// Remove from roots any solution r such that p does not evaluate to 0 at x2v extended with x->r.
void filter_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, polynomial::var x, numeral_vector & roots) {
TRACE("isolate_roots", tout << "before filtering roots, x: x" << x << "\n";
TRACE(isolate_roots, tout << "before filtering roots, x: x" << x << "\n";
for (unsigned i = 0; i < roots.size(); i++) {
display_root(tout, roots[i]); tout << "\n";
});
@ -2342,9 +2342,9 @@ namespace algebraic_numbers {
for (unsigned i = 0; i < sz; i++) {
checkpoint();
ext_var2num ext_x2v(m_wrapper, x2v, x, roots[i]);
TRACE("isolate_roots", tout << "filter_roots i: " << i << ", ext_x2v: x" << x << " -> "; display_root(tout, roots[i]); tout << "\n";);
TRACE(isolate_roots, tout << "filter_roots i: " << i << ", ext_x2v: x" << x << " -> "; display_root(tout, roots[i]); tout << "\n";);
sign sign = eval_sign_at(p, ext_x2v);
TRACE("isolate_roots", tout << "filter_roots i: " << i << ", result sign: " << sign << "\n";);
TRACE(isolate_roots, tout << "filter_roots i: " << i << ", result sign: " << sign << "\n";);
if (sign != 0)
continue;
if (i != j)
@ -2355,7 +2355,7 @@ namespace algebraic_numbers {
del(roots[i]);
roots.shrink(j);
TRACE("isolate_roots", tout << "after filtering roots:\n";
TRACE(isolate_roots, tout << "after filtering roots:\n";
for (unsigned i = 0; i < roots.size(); i++) {
display_root(tout, roots[i]); tout << "\n";
});
@ -2381,16 +2381,16 @@ namespace algebraic_numbers {
polynomial::var_vector m_isolate_roots_vars;
void isolate_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, numeral_vector & roots, bool nested_call = false) {
TRACE("isolate_roots", tout << "isolating roots of: " << p << "\n";);
TRACE(isolate_roots, tout << "isolating roots of: " << p << "\n";);
SASSERT(roots.empty());
polynomial::manager & ext_pm = p.m();
if (ext_pm.is_zero(p) || ext_pm.is_const(p)) {
TRACE("isolate_roots", tout << "p is zero or the constant polynomial\n";);
TRACE(isolate_roots, tout << "p is zero or the constant polynomial\n";);
return;
}
if (ext_pm.is_univariate(p)) {
TRACE("isolate_roots", tout << "p is univariate, using univariate procedure\n";);
TRACE(isolate_roots, tout << "p is univariate, using univariate procedure\n";);
isolate_roots(p, roots);
return;
}
@ -2399,10 +2399,10 @@ namespace algebraic_numbers {
polynomial_ref p_prime(ext_pm);
var2basic x2v_basic(*this, x2v);
p_prime = ext_pm.substitute(p, x2v_basic);
TRACE("isolate_roots", tout << "p after applying (rational fragment of) x2v:\n" << p_prime << "\n";);
TRACE(isolate_roots, tout << "p after applying (rational fragment of) x2v:\n" << p_prime << "\n";);
if (ext_pm.is_zero(p_prime) || ext_pm.is_const(p_prime)) {
TRACE("isolate_roots", tout << "p is zero or the constant polynomial after applying (rational fragment of) x2v\n";);
TRACE(isolate_roots, tout << "p is zero or the constant polynomial after applying (rational fragment of) x2v\n";);
return;
}
@ -2413,7 +2413,7 @@ namespace algebraic_numbers {
// So, the polynomial does not have any roots
return;
}
TRACE("isolate_roots", tout << "p is univariate after applying (rational fragment of) x2v... using univariate procedure\n";);
TRACE(isolate_roots, tout << "p is univariate after applying (rational fragment of) x2v... using univariate procedure\n";);
isolate_roots(p_prime, roots);
return;
}
@ -2425,7 +2425,7 @@ namespace algebraic_numbers {
// sort variables by the degree of the values
std::stable_sort(xs.begin(), xs.end(), var_degree_lt(*this, x2v));
TRACE("isolate_roots", tout << "there are " << (xs.size() - 1) << " variables assigned to nonbasic numbers...\n";);
TRACE(isolate_roots, tout << "there are " << (xs.size() - 1) << " variables assigned to nonbasic numbers...\n";);
// last variables is the one not assigned by x2v, or the unassigned variable vanished
polynomial::var x = xs.back();
@ -2453,7 +2453,7 @@ namespace algebraic_numbers {
algebraic_cell * c = v.to_algebraic();
p_y = ext_pm.to_polynomial(c->m_p_sz, c->m_p, y);
ext_pm.resultant(q, p_y, y, q);
TRACE("isolate_roots", tout << "resultant loop i: " << i << ", y: x" << y << "\np_y: " << p_y << "\n";
TRACE(isolate_roots, tout << "resultant loop i: " << i << ", y: x" << y << "\np_y: " << p_y << "\n";
tout << "q: " << q << "\n";);
if (ext_pm.is_zero(q)) {
SASSERT(!nested_call);
@ -2462,7 +2462,7 @@ namespace algebraic_numbers {
}
if (ext_pm.is_zero(q)) {
TRACE("isolate_roots", tout << "q vanished\n";);
TRACE(isolate_roots, tout << "q vanished\n";);
// q may vanish at some of the other roots of the polynomial defining the values.
// To decide if p_prime vanishes at x2v or not, we start evaluating each coefficient of p_prime at x2v
// until we find one that is not zero at x2v.
@ -2473,7 +2473,7 @@ namespace algebraic_numbers {
SASSERT(n > 0);
if (n == 1) {
// p_prime is linear on p, so we just evaluate the coefficients...
TRACE("isolate_roots", tout << "p is linear after applying (rational fragment) of x2v\n";);
TRACE(isolate_roots, tout << "p is linear after applying (rational fragment) of x2v\n";);
polynomial_ref c0(ext_pm);
polynomial_ref c1(ext_pm);
c0 = ext_pm.coeff(p_prime, x, 0);
@ -2484,14 +2484,14 @@ namespace algebraic_numbers {
ext_pm.eval(c1, x2v, a1);
// the root must be - a0/a1 if a1 != 0
if (is_zero(a1)) {
TRACE("isolate_roots", tout << "coefficient of degree 1 vanished, so p does not have roots at x2v\n";);
TRACE(isolate_roots, tout << "coefficient of degree 1 vanished, so p does not have roots at x2v\n";);
// p_prime does not have any root
return;
}
roots.push_back(anum());
div(a0, a1, roots[0]);
neg(roots[0]);
TRACE("isolate_roots", tout << "after trivial solving p has only one root:\n"; display_root(tout, roots[0]); tout << "\n";);
TRACE(isolate_roots, tout << "after trivial solving p has only one root:\n"; display_root(tout, roots[0]); tout << "\n";);
}
else {
polynomial_ref c(ext_pm);
@ -2506,7 +2506,7 @@ namespace algebraic_numbers {
if (i == 0) {
// all coefficients of x vanished, so
// the polynomial has no roots
TRACE("isolate_roots", tout << "all coefficients vanished... polynomial does not have roots\n";);
TRACE(isolate_roots, tout << "all coefficients vanished... polynomial does not have roots\n";);
return;
}
SASSERT(!is_zero(a));
@ -2527,7 +2527,7 @@ namespace algebraic_numbers {
xi_p = pm().mk_polynomial(x, i);
z_p = pm().mk_polynomial(z);
q2 = z_p*xi_p + q2;
TRACE("isolate_roots", tout << "invoking isolate_roots with q2:\n" << q2 << "\n";
TRACE(isolate_roots, tout << "invoking isolate_roots with q2:\n" << q2 << "\n";
tout << "z: x" << z << " -> "; display_root(tout, a); tout << "\n";);
// extend x2p with z->a
ext_var2num ext_x2v(m_wrapper, x2v, z, a);
@ -2536,7 +2536,7 @@ namespace algebraic_numbers {
}
else if (ext_pm.is_const(q)) {
// q does not have any roots, so p_prime also does not have roots at x2v.
TRACE("isolate_roots", tout << "q is the constant polynomial, so p does not contain any roots at x2v\n";);
TRACE(isolate_roots, tout << "q is the constant polynomial, so p does not contain any roots at x2v\n";);
}
else {
SASSERT(is_univariate(q));
@ -2619,7 +2619,7 @@ namespace algebraic_numbers {
// Select a numeral between prev and curr.
// Pre: prev < curr
void select(numeral & prev, numeral & curr, numeral & result) {
TRACE("algebraic_select",
TRACE(algebraic_select,
tout << "prev: "; display_interval(tout, prev); tout << "\n";
tout << "curr: "; display_interval(tout, curr); tout << "\n";);
SASSERT(lt(prev, curr));
@ -2674,7 +2674,7 @@ namespace algebraic_numbers {
signs.push_back(s);
}
else {
TRACE("isolate_roots_bug", tout << "p: " << p << "\n";
TRACE(isolate_roots_bug, tout << "p: " << p << "\n";
polynomial::var_vector xs;
p.m().vars(p, xs);
for (unsigned i = 0; i < xs.size(); i++) {
@ -2694,7 +2694,7 @@ namespace algebraic_numbers {
scoped_anum w(m_wrapper);
int_lt(roots[0], w);
TRACE("isolate_roots_bug", tout << "w: "; display_root(tout, w); tout << "\n";);
TRACE(isolate_roots_bug, tout << "w: "; display_root(tout, w); tout << "\n";);
{
ext2_var2num ext_x2v(m_wrapper, x2v, w);
auto s = eval_sign_at(p, ext_x2v);
@ -2996,15 +2996,15 @@ namespace algebraic_numbers {
}
void manager::power(numeral const & a, unsigned k, numeral & b) {
TRACE("anum_detail", display_root(tout, a); tout << "^" << k << "\n";);
TRACE(anum_detail, display_root(tout, a); tout << "^" << k << "\n";);
m_imp->power(const_cast<numeral&>(a), k, b);
TRACE("anum_detail", tout << "^ result: "; display_root(tout, b); tout << "\n";);
TRACE(anum_detail, tout << "^ result: "; display_root(tout, b); tout << "\n";);
}
void manager::add(numeral const & a, numeral const & b, numeral & c) {
TRACE("anum_detail", display_root(tout, a); tout << " + "; display_root(tout, b); tout << "\n";);
TRACE(anum_detail, display_root(tout, a); tout << " + "; display_root(tout, b); tout << "\n";);
m_imp->add(const_cast<numeral&>(a), const_cast<numeral&>(b), c);
TRACE("anum_detail", tout << "+ result: "; display_root(tout, c); tout << "\n";);
TRACE(anum_detail, tout << "+ result: "; display_root(tout, c); tout << "\n";);
}
void manager::add(numeral const & a, mpz const & b, numeral & c) {
@ -3014,15 +3014,15 @@ namespace algebraic_numbers {
}
void manager::sub(numeral const & a, numeral const & b, numeral & c) {
TRACE("anum_detail", display_root(tout, a); tout << " - "; display_root(tout, b); tout << "\n";);
TRACE(anum_detail, display_root(tout, a); tout << " - "; display_root(tout, b); tout << "\n";);
m_imp->sub(const_cast<numeral&>(a), const_cast<numeral&>(b), c);
TRACE("anum_detail", tout << "- result: "; display_root(tout, c); tout << "\n";);
TRACE(anum_detail, tout << "- result: "; display_root(tout, c); tout << "\n";);
}
void manager::mul(numeral const & a, numeral const & b, numeral & c) {
TRACE("anum_detail", display_root(tout, a); tout << " * "; display_root(tout, b); tout << "\n";);
TRACE(anum_detail, display_root(tout, a); tout << " * "; display_root(tout, b); tout << "\n";);
m_imp->mul(const_cast<numeral&>(a), const_cast<numeral&>(b), c);
TRACE("anum_detail", tout << "* result: "; display_root(tout, c); tout << "\n";);
TRACE(anum_detail, tout << "* result: "; display_root(tout, c); tout << "\n";);
}
void manager::div(numeral const & a, numeral const & b, numeral & c) {

4
src/math/polynomial/linear_eq_solver.h
View file

@ -86,7 +86,7 @@ public:
// Return false if the matrix is singular
bool solve(numeral * xs) {
for (unsigned k = 0; k < n; k++) {
TRACE("linear_eq_solver", tout << "iteration " << k << "\n"; display(tout););
TRACE(linear_eq_solver, tout << "iteration " << k << "\n"; display(tout););
// find pivot
unsigned i = k;
for (; i < n; i++) {
@ -120,7 +120,7 @@ public:
unsigned k = n;
while (k > 0) {
--k;
TRACE("linear_eq_solver", tout << "iteration " << k << "\n"; display(tout););
TRACE(linear_eq_solver, tout << "iteration " << k << "\n"; display(tout););
SASSERT(m.is_one(A[k][k]));
// save result
m.set(xs[k], b[k]);

278
src/math/polynomial/polynomial.cpp
View file

File diff suppressed because it is too large Load diff

142
src/math/polynomial/upolynomial.cpp
View file

@ -386,7 +386,7 @@ namespace upolynomial {
#ifdef Z3DEBUG
scoped_numeral tmp(m());
m().mul(g, p[i], tmp);
CTRACE("div_bug", !m().eq(tmp, old_p_i), tout << "old(p[i]): " << m().to_string(old_p_i) << ", g: " << m().to_string(g) << ", p[i]: " <<
CTRACE(div_bug, !m().eq(tmp, old_p_i), tout << "old(p[i]): " << m().to_string(old_p_i) << ", g: " << m().to_string(g) << ", p[i]: " <<
m().to_string(p[i]) << ", tmp: " << m().to_string(tmp) << "\n";);
SASSERT(tmp == old_p_i);
#endif
@ -403,7 +403,7 @@ namespace upolynomial {
if (m().is_one(b))
return;
for (unsigned i = 0; i < sz; i++) {
CTRACE("upolynomial", !m().divides(b, p[i]), tout << "b: " << m().to_string(b) << ", p[i]: " << m().to_string(p[i]) << "\n";);
CTRACE(upolynomial, !m().divides(b, p[i]), tout << "b: " << m().to_string(b) << ", p[i]: " << m().to_string(p[i]) << "\n";);
SASSERT(m().divides(b, p[i]));
m().div(p[i], b, p[i]);
}
@ -527,10 +527,10 @@ namespace upolynomial {
SASSERT(!m().is_zero(b_n));
scoped_numeral a_m(m());
while (m_limit.inc()) {
TRACE("rem_bug", tout << "rem loop, p2:\n"; display(tout, sz2, p2); tout << "\nbuffer:\n"; display(tout, buffer); tout << "\n";);
TRACE(rem_bug, tout << "rem loop, p2:\n"; display(tout, sz2, p2); tout << "\nbuffer:\n"; display(tout, buffer); tout << "\n";);
sz1 = buffer.size();
if (sz1 < sz2) {
TRACE("rem_bug", tout << "finished\n";);
TRACE(rem_bug, tout << "finished\n";);
return;
}
unsigned m_n = sz1 - sz2;
@ -546,7 +546,7 @@ namespace upolynomial {
// p2: b_n * x^n + b_{n-1} * x^{n-1} + ... + b_0
d++;
m().set(a_m, buffer[sz1 - 1]);
TRACE("rem_bug", tout << "a_m: " << m().to_string(a_m) << ", b_n: " << m().to_string(b_n) << "\n";);
TRACE(rem_bug, tout << "a_m: " << m().to_string(a_m) << ", b_n: " << m().to_string(b_n) << "\n";);
// don't need to update position sz1 - 1, since it will become 0
for (unsigned i = 0; i < sz1 - 1; i++) {
m().mul(buffer[i], b_n, buffer[i]);
@ -619,11 +619,11 @@ namespace upolynomial {
_r.reserve(deg+1);
numeral_vector & _p1 = m_div_tmp1;
// std::cerr << "dividing with "; display(std::cerr, _p1); std::cerr << std::endl;
TRACE("factor_bug", tout << "sz1: " << sz1 << " p1: " << p1 << ", _p1.c_ptr(): " << _p1.data() << ", _p1.size(): " << _p1.size() << "\n";);
TRACE(factor_bug, tout << "sz1: " << sz1 << " p1: " << p1 << ", _p1.c_ptr(): " << _p1.data() << ", _p1.size(): " << _p1.size() << "\n";);
set(sz1, p1, _p1);
SASSERT(_p1.size() == sz1);
while (true) {
TRACE("upolynomial", tout << "exact_div loop...\n"; display(tout, _p1); tout << "\n"; display(tout, _r); tout << "\n";);
TRACE(upolynomial, tout << "exact_div loop...\n"; display(tout, _p1); tout << "\n"; display(tout, _r); tout << "\n";);
// std::cerr << "dividing with "; display(std::cerr, _p1); std::cerr << std::endl;
if (sz1 == 0) {
set_size(deg+1, _r);
@ -672,12 +672,12 @@ namespace upolynomial {
// inv2 is the inverse of b2 mod b1
m().m().mod(inv1, b2, inv1);
m().m().mod(inv2, b1, inv2);
TRACE("CRA", tout << "inv1: " << inv1 << ", inv2: " << inv2 << "\n";);
TRACE(CRA, tout << "inv1: " << inv1 << ", inv2: " << inv2 << "\n";);
scoped_numeral a1(m());
scoped_numeral a2(m());
m().mul(b2, inv2, a1); // a1 is the multiplicator for coefficients of C1
m().mul(b1, inv1, a2); // a2 is the multiplicator for coefficients of C2
TRACE("CRA", tout << "a1: " << a1 << ", a2: " << a2 << "\n";);
TRACE(CRA, tout << "a1: " << a1 << ", a2: " << a2 << "\n";);
// new bound
scoped_numeral new_bound(m());
m().mul(b1, b2, new_bound);
@ -687,7 +687,7 @@ namespace upolynomial {
m().div(new_bound, 2, upper);
m().set(lower, upper);
m().neg(lower);
TRACE("CRA", tout << "lower: " << lower << ", upper: " << upper << "\n";);
TRACE(CRA, tout << "lower: " << lower << ", upper: " << upper << "\n";);
#define ADD(A1, A2) { \
m().mul(A1, a1, tmp1); \
@ -721,7 +721,7 @@ namespace upolynomial {
void core_manager::mod_gcd(unsigned sz_u, numeral const * u,
unsigned sz_v, numeral const * v,
numeral_vector & result) {
TRACE("mgcd", tout << "u: "; display_star(tout, sz_u, u); tout << "\nv: "; display_star(tout, sz_v, v); tout << "\n";);
TRACE(mgcd, tout << "u: "; display_star(tout, sz_u, u); tout << "\nv: "; display_star(tout, sz_v, v); tout << "\n";);
SASSERT(sz_u > 0 && sz_v > 0);
SASSERT(!m().modular());
scoped_numeral c_u(m()), c_v(m());
@ -747,17 +747,17 @@ namespace upolynomial {
for (unsigned i = 0; i < NUM_BIG_PRIMES; i++) {
m().set(p, polynomial::g_big_primes[i]);
TRACE("mgcd", tout << "trying prime: " << p << "\n";);
TRACE(mgcd, tout << "trying prime: " << p << "\n";);
{
scoped_set_zp setZp(*this, p);
set(pp_u.size(), pp_u.data(), u_Zp);
set(pp_v.size(), pp_v.data(), v_Zp);
if (degree(u_Zp) < d_u) {
TRACE("mgcd", tout << "bad prime, leading coefficient vanished\n";);
TRACE(mgcd, tout << "bad prime, leading coefficient vanished\n";);
continue; // bad prime
}
if (degree(v_Zp) < d_v) {
TRACE("mgcd", tout << "bad prime, leading coefficient vanished\n";);
TRACE(mgcd, tout << "bad prime, leading coefficient vanished\n";);
continue; // bad prime
}
euclid_gcd(u_Zp.size(), u_Zp.data(), v_Zp.size(), v_Zp.data(), q);
@ -767,9 +767,9 @@ namespace upolynomial {
m().set(c, lc_g);
mul(q, c);
}
TRACE("mgcd", tout << "new q:\n"; display_star(tout, q); tout << "\n";);
TRACE(mgcd, tout << "new q:\n"; display_star(tout, q); tout << "\n";);
if (is_const(q)) {
TRACE("mgcd", tout << "done, modular gcd is one\n";);
TRACE(mgcd, tout << "done, modular gcd is one\n";);
reset(result);
result.push_back(numeral());
m().set(result.back(), c_g);
@ -781,27 +781,27 @@ namespace upolynomial {
}
else if (q.size() < C.size() || m().m().is_even(p) || m().m().is_even(bound)) {
// discard accumulated image, it was affected by unlucky primes
TRACE("mgcd", tout << "discarding image\n";);
TRACE(mgcd, tout << "discarding image\n";);
set(q.size(), q.data(), C);
m().set(bound, p);
}
else {
CRA_combine_images(q, p, C, bound);
TRACE("mgcd", tout << "new combined:\n"; display_star(tout, C); tout << "\n";);
TRACE(mgcd, tout << "new combined:\n"; display_star(tout, C); tout << "\n";);
}
numeral_vector & candidate = q;
get_primitive(C, candidate);
TRACE("mgcd", tout << "candidate:\n"; display_star(tout, candidate); tout << "\n";);
TRACE(mgcd, tout << "candidate:\n"; display_star(tout, candidate); tout << "\n";);
SASSERT(candidate.size() > 0);
numeral const & lc_candidate = candidate[candidate.size() - 1];
if (m().divides(lc_candidate, lc_g) &&
divides(pp_u, candidate) &&
divides(pp_v, candidate)) {
TRACE("mgcd", tout << "found GCD\n";);
TRACE(mgcd, tout << "found GCD\n";);
mul(candidate, c_g);
flip_sign_if_lm_neg(candidate);
candidate.swap(result);
TRACE("mgcd", tout << "r: "; display_star(tout, result); tout << "\n";);
TRACE(mgcd, tout << "r: "; display_star(tout, result); tout << "\n";);
return;
}
}
@ -828,10 +828,10 @@ namespace upolynomial {
numeral_vector & R = buffer;
set(sz1, p1, A);
set(sz2, p2, B);
TRACE("upolynomial", tout << "sz1: " << sz1 << ", p1: " << p1 << ", sz2: " << sz2 << ", p2: " << p2 << "\nB.size(): " << B.size() <<
TRACE(upolynomial, tout << "sz1: " << sz1 << ", p1: " << p1 << ", sz2: " << sz2 << ", p2: " << p2 << "\nB.size(): " << B.size() <<
", B.c_ptr(): " << B.data() << "\n";);
while (m_limit.inc()) {
TRACE("upolynomial", tout << "A: "; display(tout, A); tout <<"\nB: "; display(tout, B); tout << "\n";);
TRACE(upolynomial, tout << "A: "; display(tout, A); tout <<"\nB: "; display(tout, B); tout << "\n";);
if (B.empty()) {
normalize(A);
buffer.swap(A);
@ -842,7 +842,7 @@ namespace upolynomial {
else {
flip_sign_if_lm_neg(buffer);
}
TRACE("upolynomial", tout << "GCD\n"; display(tout, sz1, p1); tout << "\n"; display(tout, sz2, p2); tout << "\n--->\n";
TRACE(upolynomial, tout << "GCD\n"; display(tout, sz1, p1); tout << "\n"; display(tout, sz2, p2); tout << "\n--->\n";
display(tout, buffer); tout << "\n";);
return;
}
@ -899,7 +899,7 @@ namespace upolynomial {
A.swap(B);
while (true) {
SASSERT(A.size() >= B.size());
TRACE("upolynomial", tout << "A: "; display(tout, A); tout <<"\nB: "; display(tout, B); tout << "\n";
TRACE(upolynomial, tout << "A: "; display(tout, A); tout <<"\nB: "; display(tout, B); tout << "\n";
tout << "g: " << m().to_string(g) << ", h: " << m().to_string(h) << "\n";);
if (B.empty()) {
normalize(A);
@ -911,7 +911,7 @@ namespace upolynomial {
else {
flip_sign_if_lm_neg(buffer);
}
TRACE("upolynomial", tout << "subresultant GCD\n"; display(tout, sz1, p1); tout << "\n"; display(tout, sz2, p2); tout << "\n--->\n";
TRACE(upolynomial, tout << "subresultant GCD\n"; display(tout, sz1, p1); tout << "\n"; display(tout, sz2, p2); tout << "\n--->\n";
display(tout, buffer); tout << "\n";);
return;
}
@ -927,7 +927,7 @@ namespace upolynomial {
mul(R, aux);
}
d = pseudo_div_d;
TRACE("upolynomial", tout << "R: "; display(tout, R); tout << "\nd: " << d << "\n";);
TRACE(upolynomial, tout << "R: "; display(tout, R); tout << "\nd: " << d << "\n";);
// aux <- g*h^d
m().power(h, d, aux);
m().mul(g, aux, aux);
@ -1407,14 +1407,14 @@ namespace upolynomial {
// Basic idea: apply descartes_bound_0_1 to p2(x) where
// p1(x) = p(x+a)
// p2(x) = p1((b-a)*x)
TRACE("upolynomial", tout << "pos interval... " << bqm.to_string(a) << ", " << bqm.to_string(b) << "\n"; display(tout, sz, p); tout << "\n";);
TRACE(upolynomial, tout << "pos interval... " << bqm.to_string(a) << ", " << bqm.to_string(b) << "\n"; display(tout, sz, p); tout << "\n";);
numeral_vector & p_aux = m_dbab_tmp1;
translate_bq(sz, p, a, p_aux);
TRACE("upolynomial", tout << "after translation\n"; display(tout, p_aux); tout << "\n";);
TRACE(upolynomial, tout << "after translation\n"; display(tout, p_aux); tout << "\n";);
scoped_mpbq b_a(bqm);
bqm.sub(b, a, b_a);
compose_p_b_x(p_aux.size(), p_aux.data(), b_a);
TRACE("upolynomial", tout << "after composition: " << bqm.to_string(b_a) << "\n"; display(tout, p_aux); tout << "\n";);
TRACE(upolynomial, tout << "after composition: " << bqm.to_string(b_a) << "\n"; display(tout, p_aux); tout << "\n";);
unsigned result = descartes_bound_0_1(p_aux.size(), p_aux.data());
return result;
}
@ -1538,7 +1538,7 @@ namespace upolynomial {
return;
// Step 1
compose_2kn_p_x_div_2k(sz, p, b.k());
TRACE("upolynomial", tout << "after compose 2kn_p_x_div_2k\n"; display(tout, sz, p); tout << "\n";);
TRACE(upolynomial, tout << "after compose 2kn_p_x_div_2k\n"; display(tout, sz, p); tout << "\n";);
// Step 2
numeral const & c = b.numerator();
unsigned n = sz - 1;
@ -1551,7 +1551,7 @@ namespace upolynomial {
}
m().mul2k(p[n], b.k());
}
TRACE("upolynomial", tout << "after special translation\n"; display(tout, sz, p); tout << "\n";);
TRACE(upolynomial, tout << "after special translation\n"; display(tout, sz, p); tout << "\n";);
}
// Similar to translate_bq but for rationals
@ -1560,7 +1560,7 @@ namespace upolynomial {
return;
// Step 1
compose_an_p_x_div_a(sz, p, b.denominator());
TRACE("upolynomial", tout << "after compose_an_p_x_div_a\n"; display(tout, sz, p); tout << "\n";);
TRACE(upolynomial, tout << "after compose_an_p_x_div_a\n"; display(tout, sz, p); tout << "\n";);
// Step 2
numeral const & c = b.numerator();
unsigned n = sz - 1;
@ -1573,7 +1573,7 @@ namespace upolynomial {
}
m().mul(p[n], b.denominator(), p[n]);
}
TRACE("upolynomial", tout << "after special translation\n"; display(tout, sz, p); tout << "\n";);
TRACE(upolynomial, tout << "after special translation\n"; display(tout, sz, p); tout << "\n";);
}
// p(x) := 2^n*p(x/2) where n = sz-1
@ -2146,7 +2146,7 @@ namespace upolynomial {
unsigned idx = frame_stack.size() - 1;
while (idx != UINT_MAX) {
drs_frame const & fr = frame_stack[idx];
TRACE("upolynomial",
TRACE(upolynomial,
tout << "normalizing...\n";
tout << "idx: " << idx << ", left: " << fr.m_left << ", l: " << bqm.to_string(l) << ", u: " << bqm.to_string(u) << "\n";);
if (fr.m_left) {
@ -2161,7 +2161,7 @@ namespace upolynomial {
}
idx = fr.m_parent_idx;
}
TRACE("upolynomial", tout << "adding normalized interval (" << bqm.to_string(l) << ", " << bqm.to_string(u) << ")\n";);
TRACE(upolynomial, tout << "adding normalized interval (" << bqm.to_string(l) << ", " << bqm.to_string(u) << ")\n";);
lowers.push_back(mpbq());
uppers.push_back(mpbq());
swap(lowers.back(), l);
@ -2187,32 +2187,32 @@ namespace upolynomial {
}
idx = fr.m_parent_idx;
}
TRACE("upolynomial", tout << "adding normalized root: " << bqm.to_string(u) << "\n";);
TRACE(upolynomial, tout << "adding normalized root: " << bqm.to_string(u) << "\n";);
roots.push_back(mpbq());
swap(roots.back(), u);
}
// Isolate roots in the interval (0, 1)
void manager::drs_isolate_0_1_roots(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers) {
TRACE("upolynomial", tout << "isolating (0,1) roots of:\n"; display(tout, sz, p); tout << "\n";);
TRACE(upolynomial, tout << "isolating (0,1) roots of:\n"; display(tout, sz, p); tout << "\n";);
unsigned k = descartes_bound_0_1(sz, p);
// easy cases...
if (k == 0) {
TRACE("upolynomial", tout << "polynomial does not have any roots\n";);
TRACE(upolynomial, tout << "polynomial does not have any roots\n";);
return;
}
if (k == 1) {
TRACE("upolynomial", tout << "polynomial has one root in (0, 1)\n";);
TRACE(upolynomial, tout << "polynomial has one root in (0, 1)\n";);
lowers.push_back(mpbq(0));
uppers.push_back(mpbq(1));
return;
}
TRACE("upolynomial", tout << "polynomial has more than one root in (0, 1), starting search...\n";);
TRACE(upolynomial, tout << "polynomial has more than one root in (0, 1), starting search...\n";);
scoped_numeral_vector q(m());
scoped_numeral_vector p_stack(m());
svector<drs_frame> frame_stack;
if (has_one_half_root(sz, p)) {
TRACE("upolynomial", tout << "polynomial has a 1/2 root\n";);
TRACE(upolynomial, tout << "polynomial has a 1/2 root\n";);
roots.push_back(mpbq(1, 1));
remove_one_half_root(sz, p, q);
push_child_frames(q.size(), q.data(), p_stack, frame_stack);
@ -2227,7 +2227,7 @@ namespace upolynomial {
unsigned sz = fr.m_size;
SASSERT(sz <= p_stack.size());
numeral const * p = p_stack.data() + p_stack.size() - sz;
TRACE("upolynomial", tout << "processing frame #" << frame_stack.size() - 1 << "\n"
TRACE(upolynomial, tout << "processing frame #" << frame_stack.size() - 1 << "\n"
<< "first: " << fr.m_first << ", left: " << fr.m_left << ", sz: " << fr.m_size << ", parent_idx: ";
if (fr.m_parent_idx == UINT_MAX) tout << "<null>"; else tout << fr.m_parent_idx;
tout << "\np: "; display(tout, sz, p); tout << "\n";);
@ -2238,19 +2238,19 @@ namespace upolynomial {
fr.m_first = false;
unsigned k = descartes_bound_0_1(sz, p);
if (k == 0) {
TRACE("upolynomial", tout << "(0, 1) does not have roots\n";);
TRACE(upolynomial, tout << "(0, 1) does not have roots\n";);
pop_top_frame(p_stack, frame_stack);
continue;
}
if (k == 1) {
TRACE("upolynomial", tout << "(0, 1) is isolating interval\n";);
TRACE(upolynomial, tout << "(0, 1) is isolating interval\n";);
add_isolating_interval(frame_stack, bqm, lowers, uppers);
pop_top_frame(p_stack, frame_stack);
continue;
}
TRACE("upolynomial", tout << "polynomial has more than one root in (0, 1) creating child frames...\n";);
TRACE(upolynomial, tout << "polynomial has more than one root in (0, 1) creating child frames...\n";);
if (has_one_half_root(sz, p)) {
TRACE("upolynomial", tout << "1/2 is a root\n";);
TRACE(upolynomial, tout << "1/2 is a root\n";);
add_root(frame_stack, bqm, roots);
remove_one_half_root(sz, p, q);
push_child_frames(q.size(), q.data(), p_stack, frame_stack);
@ -2295,7 +2295,7 @@ namespace upolynomial {
// p(x) := p(2^{pos_k} * x)
// Since the desired positive roots of p(x) are in (0, 2^pos_k),
TRACE("upolynomial", tout << "searching at (0, 1)\n";);
TRACE(upolynomial, tout << "searching at (0, 1)\n";);
unsigned old_roots_sz = roots.size();
unsigned old_lowers_sz = lowers.size();
drs_isolate_0_1_roots(sz, aux_p.data(), bqm, roots, lowers, uppers);
@ -2308,7 +2308,7 @@ namespace upolynomial {
// p(x) := p(-x)
p_minus_x(sz, p);
compose_p_2k_x(sz, p, neg_k);
TRACE("upolynomial", tout << "searching at (-1, 0) using:\n"; display(tout, sz, p); tout << "\n";);
TRACE(upolynomial, tout << "searching at (-1, 0) using:\n"; display(tout, sz, p); tout << "\n";);
old_roots_sz = roots.size();
old_lowers_sz = lowers.size();
drs_isolate_0_1_roots(sz, p, bqm, roots, lowers, uppers);
@ -2328,7 +2328,7 @@ namespace upolynomial {
set(sz, p, p1);
normalize(p1);
TRACE("upolynomial",
TRACE(upolynomial,
scoped_numeral U(m());
root_upper_bound(p1.size(), p1.data(), U);
unsigned U_k = m().log2(U) + 1;
@ -2395,7 +2395,7 @@ namespace upolynomial {
scoped_mpbq curr_upper(bqm);
sturm_seq(sz, p, seq);
ss_frame_stack s(bqm);
TRACE("upolynomial", tout << "p: "; display(tout, sz, p); tout << "\nSturm seq:\n"; display(tout, seq); tout << "\n";);
TRACE(upolynomial, tout << "p: "; display(tout, sz, p); tout << "\nSturm seq:\n"; display(tout, seq); tout << "\n";);
unsigned lower_sv = sign_variations_at_minus_inf(seq);
unsigned zero_sv = sign_variations_at_zero(seq);
@ -2449,7 +2449,7 @@ namespace upolynomial {
SASSERT(lower_sv > upper_sv + 1);
bqm.add(curr_lower, curr_upper, mid);
bqm.div2(mid);
TRACE("upolynomial",
TRACE(upolynomial,
tout << "depth: " << s.size() << "\n";
tout << "lower_sv: " << lower_sv << "\n";
tout << "upper_sv: " << upper_sv << "\n";
@ -2495,7 +2495,7 @@ namespace upolynomial {
roots.push_back(mpbq(0));
scoped_numeral_vector nz_p(m());
remove_zero_roots(sz, p, nz_p);
TRACE("upolynomial", tout << "after removing zero root:\n"; display(tout, nz_p); tout << "\n";);
TRACE(upolynomial, tout << "after removing zero root:\n"; display(tout, nz_p); tout << "\n";);
SASSERT(!has_zero_roots(nz_p.size(), nz_p.data()));
sqf_nz_isolate_roots(nz_p.size(), nz_p.data(), bqm, roots, lowers, uppers);
}
@ -2506,10 +2506,10 @@ namespace upolynomial {
void manager::isolate_roots(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq_vector & roots, mpbq_vector & lowers, mpbq_vector & uppers) {
SASSERT(sz > 0);
TRACE("upolynomial", tout << "isolating roots of:\n"; display(tout, sz, p); tout << "\n";);
TRACE(upolynomial, tout << "isolating roots of:\n"; display(tout, sz, p); tout << "\n";);
scoped_numeral_vector sqf_p(m());
square_free(sz, p, sqf_p);
TRACE("upolynomial", tout << "square free part:\n"; display(tout, sqf_p); tout << "\n";);
TRACE(upolynomial, tout << "square free part:\n"; display(tout, sqf_p); tout << "\n";);
sqf_isolate_roots(sqf_p.size(), sqf_p.data(), bqm, roots, lowers, uppers);
}
@ -2596,7 +2596,7 @@ namespace upolynomial {
bool manager::isolating2refinable(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq & a, mpbq & b) {
int sign_a = eval_sign_at(sz, p, a);
int sign_b = eval_sign_at(sz, p, b);
TRACE("upolynomial", tout << "sign_a: " << sign_a << ", sign_b: " << sign_b << "\n";);
TRACE(upolynomial, tout << "sign_a: " << sign_a << ", sign_b: " << sign_b << "\n";);
if (sign_a != 0 && sign_b != 0) {
// CASE 1
SASSERT(sign_a == -sign_b); // p is square free
@ -2609,7 +2609,7 @@ namespace upolynomial {
bqm.add(a, b, new_a);
bqm.div2(new_a);
while (true) {
TRACE("upolynomial", tout << "CASE 2, a: " << bqm.to_string(a) << ", b: " << bqm.to_string(b) << ", new_a: " << bqm.to_string(new_a) << "\n";);
TRACE(upolynomial, tout << "CASE 2, a: " << bqm.to_string(a) << ", b: " << bqm.to_string(b) << ", new_a: " << bqm.to_string(new_a) << "\n";);
int sign_new_a = eval_sign_at(sz, p, new_a);
if (sign_new_a != sign_b) {
swap(new_a, a);
@ -2634,7 +2634,7 @@ namespace upolynomial {
bqm.add(a, b, new_b);
bqm.div2(new_b);
while (true) {
TRACE("upolynomial", tout << "CASE 3, a: " << bqm.to_string(a) << ", b: " << bqm.to_string(b) << ", new_b: " << bqm.to_string(new_b) << "\n";);
TRACE(upolynomial, tout << "CASE 3, a: " << bqm.to_string(a) << ", b: " << bqm.to_string(b) << ", new_b: " << bqm.to_string(new_b) << "\n";);
int sign_new_b = eval_sign_at(sz, p, new_b);
if (sign_new_b != sign_a) {
SASSERT(sign_new_b == 0 || // found the actual root
@ -2687,7 +2687,7 @@ namespace upolynomial {
bqm.div2(new_b2);
while (true) {
TRACE("upolynomial",
TRACE(upolynomial,
tout << "CASE 4\na1: " << bqm.to_string(a1) << ", b1: " << bqm.to_string(b1) << ", new_a1: " << bqm.to_string(new_a1) << "\n";
tout << "a2: " << bqm.to_string(a2) << ", b2: " << bqm.to_string(b2) << ", new_b2: " << bqm.to_string(new_b2) << "\n";);
int sign_new_a1 = eval_sign_at(sz, p, new_a1);
@ -2814,14 +2814,14 @@ namespace upolynomial {
SASSERT(!::is_zero(sign_a) && !::is_zero(sign_b));
SASSERT(sign_a == -sign_b);
bool found_d = false;
TRACE("convert_bug",
TRACE(convert_bug,
tout << "a: " << m().to_string(a.numerator()) << "/" << m().to_string(a.denominator()) << "\n";
tout << "b: " << m().to_string(b.numerator()) << "/" << m().to_string(b.denominator()) << "\n";
tout << "sign_a: " << sign_a << "\n";
tout << "sign_b: " << sign_b << "\n";);
scoped_mpbq lower(bqm), upper(bqm);
if (bqm.to_mpbq(a, lower)) {
TRACE("convert_bug", tout << "found c: " << lower << "\n";);
TRACE(convert_bug, tout << "found c: " << lower << "\n";);
// found c
swap(c, lower);
SASSERT(bqm.eq(c, a));
@ -2832,7 +2832,7 @@ namespace upolynomial {
bqm.mul2(upper);
if (m_manager.is_neg(a.numerator()))
::swap(lower, upper);
TRACE("convert_bug",
TRACE(convert_bug,
tout << "a: "; m().display(tout, a.numerator()); tout << "/"; m().display(tout, a.denominator()); tout << "\n";
tout << "lower: "; bqm.display(tout, lower); tout << ", upper: "; bqm.display(tout, upper); tout << "\n";);
SASSERT(bqm.lt(lower, a));
@ -2958,12 +2958,12 @@ namespace upolynomial {
void manager::factor_2_sqf_pp(numeral_vector & p, factors & r, unsigned k) {
SASSERT(p.size() == 3); // p has degree 2
TRACE("factor", tout << "factor square free (degree == 2):\n"; display(tout, p); tout << "\n";);
TRACE(factor, tout << "factor square free (degree == 2):\n"; display(tout, p); tout << "\n";);
numeral const & a = p[2];
numeral const & b = p[1];
numeral const & c = p[0];
TRACE("factor", tout << "a: " << m().to_string(a) << "\nb: " << m().to_string(b) << "\nc: " << m().to_string(c) << "\n";);
TRACE(factor, tout << "a: " << m().to_string(a) << "\nb: " << m().to_string(b) << "\nc: " << m().to_string(c) << "\n";);
// Create the discriminant: b^2 - 4*a*c
scoped_numeral b2(m());
scoped_numeral ac(m());
@ -2979,7 +2979,7 @@ namespace upolynomial {
r.push_back(p, k);
return;
}
TRACE("factor", tout << "disc_sqrt: " << m().to_string(disc_sqrt) << "\n";);
TRACE(factor, tout << "disc_sqrt: " << m().to_string(disc_sqrt) << "\n";);
// p = cont*(2*a*x + b - disc_sqrt)*(2*a*x + b + disc_sqrt)
scoped_numeral_vector f1(m());
scoped_numeral_vector f2(m());
@ -2992,7 +2992,7 @@ namespace upolynomial {
set_size(2, f2);
normalize(f1);
normalize(f2);
TRACE("factor", tout << "f1: "; display(tout, f1); tout << "\nf2: "; display(tout, f2); tout << "\n";);
TRACE(factor, tout << "f1: "; display(tout, f1); tout << "\nf2: "; display(tout, f2); tout << "\n";);
DEBUG_CODE({
scoped_numeral_vector f1f2(m());
mul(f1, f2, f1f2);
@ -3064,12 +3064,12 @@ namespace upolynomial {
else {
// B is of the form P_2 * P_3^2 * ... * P_k^{k-1}
VERIFY(exact_div(C, B, A));
TRACE("factor_bug", tout << "C: "; display(tout, C); tout << "\nB: "; display(tout, B); tout << "\nA: "; display(tout, A); tout << "\n";);
TRACE(factor_bug, tout << "C: "; display(tout, C); tout << "\nB: "; display(tout, B); tout << "\nA: "; display(tout, A); tout << "\n";);
// A is of the form P_1 * P_2 * ... * P_k
unsigned j = 1;
while (!is_const(A)) {
checkpoint();
TRACE("factor", tout << "factor_core main loop j: " << j << "\nA: "; display(tout, A); tout << "\nB: "; display(tout, B); tout << "\n";);
TRACE(factor, tout << "factor_core main loop j: " << j << "\nA: "; display(tout, A); tout << "\nB: "; display(tout, B); tout << "\n";);
// A is of the form P_j * P_{j+1} * P_{j+2} * ... * P_k
// B is of the form P_{j+1} * P_{j+2}^2 * ... * P_k^{k - j - 2}
gcd(A, B, D);
@ -3082,20 +3082,20 @@ namespace upolynomial {
result = false;
}
else {
TRACE("factor", tout << "const C: "; display(tout, C); tout << "\n";);
TRACE(factor, tout << "const C: "; display(tout, C); tout << "\n";);
SASSERT(C.size() == 1);
SASSERT(m().is_one(C[0]) || m().is_minus_one(C[0]));
if (m().is_minus_one(C[0]) && j % 2 == 1)
flip_sign(r);
}
TRACE("factor_bug", tout << "B: "; display(tout, B); tout << "\nD: "; display(tout, D); tout << "\n";);
TRACE(factor_bug, tout << "B: "; display(tout, B); tout << "\nD: "; display(tout, D); tout << "\n";);
VERIFY(exact_div(B, D, B));
// B is of the form P_{j+2} * ... * P_k^{k - j - 3}
A.swap(D);
// D is of the form P_{j+1} * P_{j+2} * ... * P_k
j++;
}
TRACE("factor_bug", tout << "A: "; display(tout, A); tout << "\n";);
TRACE(factor_bug, tout << "A: "; display(tout, A); tout << "\n";);
SASSERT(A.size() == 1 && m().is_one(A[0]));
}
return result;

118
src/math/polynomial/upolynomial_factorization.cpp
View file

@ -94,7 +94,7 @@ public:
m_row_pivot(m_size, -1) {
unsigned p = get_p_from_manager(m_zpm);
TRACE("polynomial::factorization::bughunt", tout << "polynomial::berlekamp_matrix("; m_upm.display(tout, f); tout << ", " << p << ")" << endl;);
TRACE(polynomial_factorization__bughunt, tout << "polynomial::berlekamp_matrix("; m_upm.display(tout, f); tout << ", " << p << ")" << endl;);
// the first row is always the vector [1, 0, ..., 0], since x^0 = 0 (modulo f)
m_matrix.push_back(1);
@ -139,7 +139,7 @@ public:
m_zpm.dec(get(i, i));
}
TRACE("polynomial::factorization::bughunt", tout << "polynomial::berlekamp_matrix():" << endl; display(tout); tout << endl;);
TRACE(polynomial_factorization__bughunt, tout << "polynomial::berlekamp_matrix():" << endl; display(tout); tout << endl;);
}
/**
@ -195,7 +195,7 @@ public:
}
}
TRACE("polynomial::factorization::bughunt", tout << "polynomial::diagonalize():" << endl; display(tout); tout << endl;);
TRACE(polynomial_factorization__bughunt, tout << "polynomial::diagonalize():" << endl; display(tout); tout << endl;);
return null_rank;
}
@ -253,7 +253,7 @@ void zp_square_free_factor(zp_manager & upm, numeral_vector const & f, zp_factor
zp_numeral_manager & nm = upm.m();
unsigned p = get_p_from_manager(upm.m());
TRACE("polynomial::factorization", tout << "polynomial::square_free_factor("; upm.display(tout, f); tout << ") over Z_" << p << endl;);
TRACE(polynomial_factorization, tout << "polynomial::square_free_factor("; upm.display(tout, f); tout << ") over Z_" << p << endl;);
scoped_numeral_vector div_tmp(nm);
@ -265,7 +265,7 @@ void zp_square_free_factor(zp_manager & upm, numeral_vector const & f, zp_factor
scoped_numeral constant(nm);
upm.mk_monic(T_0.size(), T_0.data(), constant);
sq_free_factors.set_constant(constant);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "Initial factors: " << sq_free_factors << endl;
tout << "R.<x> = GF(" << p << ")['x']" << endl;
tout << "T_0 = "; upm.display(tout, T_0); tout << endl;
@ -282,25 +282,25 @@ void zp_square_free_factor(zp_manager & upm, numeral_vector const & f, zp_factor
{
// [initialize e-loop] T = gcd(T_0, T_0'), V / T_0/T, k = 0
unsigned k = 0;
TRACE("polynomial::factorization::bughunt", tout << "k = 0" << endl;);
TRACE(polynomial_factorization__bughunt, tout << "k = 0" << endl;);
// T_0_d = T_0'
upm.derivative(T_0.size(), T_0.data(), T_0_d);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "T_0_d = T_0.derivative(x)" << endl;
tout << "T_0_d == "; upm.display(tout, T_0_d); tout << endl;
);
// T = gcd(T_0, T_0')
upm.gcd(T_0.size(), T_0.data(), T_0_d.size(), T_0_d.data(), T);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "T = T_0.gcd(T_0_d)" << endl;
tout << "T == "; upm.display(tout, T); tout << endl;
);
// V = T_0 / T
upm.div(T_0.size(), T_0.data(), T.size(), T.data(), V);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "V = T_0.quo_rem(T)[0]" << endl;
tout << "V == "; upm.display(tout, V); tout << endl;
);
@ -311,7 +311,7 @@ void zp_square_free_factor(zp_manager & upm, numeral_vector const & f, zp_factor
++ k;
// T = T/V
upm.div(T.size(), T.data(), V.size(), V.data(), T);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "T = T.quo_rem(V)[0]" << endl;
tout << "T == "; upm.display(tout, T); tout << endl;
);
@ -321,35 +321,35 @@ void zp_square_free_factor(zp_manager & upm, numeral_vector const & f, zp_factor
// W = gcd(T, V)
upm.gcd(T.size(), T.data(), V.size(), V.data(), W);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "W = T.gcd(V)" << endl;
upm.display(tout, W); tout << endl;
);
// A_ek = V/W
upm.div(V.size(), V.data(), W.size(), W.data(), A_ek);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "A_ek = V.quo_rem(W)[0]" << endl;
tout << "A_ek == "; upm.display(tout, A_ek); tout << endl;
);
// V = W
V.swap(W);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "V = W" << endl;
tout << "V == "; upm.display(tout, V); tout << endl;
);
// T = T/V
upm.div(T.size(), T.data(), V.size(), V.data(), T);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "T = T.quo_rem(V)[0]" << endl;
tout << "T == "; upm.display(tout, T); tout << endl;
);
// if not constant, we output it
if (A_ek.size() > 1) {
TRACE("polynomial::factorization::bughunt", tout << "Factor: ("; upm.display(tout, A_ek); tout << ")^" << e*k << endl;);
TRACE(polynomial_factorization__bughunt, tout << "Factor: ("; upm.display(tout, A_ek); tout << ")^" << e*k << endl;);
sq_free_factors.push_back(A_ek, e*k);
}
}
@ -361,18 +361,18 @@ void zp_square_free_factor(zp_manager & upm, numeral_vector const & f, zp_factor
T_0.push_back(numeral());
nm.set(T_0.back(), T[deg_p]);
}
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "T_0 = "; upm.display(tout, T_0); tout << endl;
);
}
TRACE("polynomial::factorization", tout << "polynomial::square_free_factor("; upm.display(tout, f); tout << ") => " << sq_free_factors << endl;);
TRACE(polynomial_factorization, tout << "polynomial::square_free_factor("; upm.display(tout, f); tout << ") => " << sq_free_factors << endl;);
}
bool zp_factor(zp_manager & upm, numeral_vector const & f, zp_factors & factors) {
TRACE("polynomial::factorization",
TRACE(polynomial_factorization,
unsigned p = get_p_from_manager(upm.m());
tout << "polynomial::factor("; upm.display(tout, f); tout << ") over Z_" << p << endl;);
@ -396,7 +396,7 @@ bool zp_factor(zp_manager & upm, numeral_vector const & f, zp_factors & factors)
// add the constant
factors.set_constant(sq_free_factors.get_constant());
TRACE("polynomial::factorization", tout << "polynomial::factor("; upm.display(tout, f); tout << ") => " << factors << endl;);
TRACE(polynomial_factorization, tout << "polynomial::factor("; upm.display(tout, f); tout << ") => " << factors << endl;);
return factors.total_factors() > 1;
}
@ -411,7 +411,7 @@ bool zp_factor_square_free_berlekamp(zp_manager & upm, numeral_vector const & f,
mpzzp_manager & zpm = upm.m();
unsigned p = get_p_from_manager(zpm);
TRACE("polynomial::factorization", tout << "upolynomial::factor_square_free_berlekamp("; upm.display(tout, f); tout << ", " << p << ")" << endl;);
TRACE(polynomial_factorization, tout << "upolynomial::factor_square_free_berlekamp("; upm.display(tout, f); tout << ", " << p << ")" << endl;);
SASSERT(zpm.is_one(f.back()));
// construct the berlekamp Q matrix to get the null space
@ -425,11 +425,11 @@ bool zp_factor_square_free_berlekamp(zp_manager & upm, numeral_vector const & f,
unsigned r = Q_I.diagonalize();
if (r == 1) {
// since r == 1 == number of factors, then f is irreducible
TRACE("polynomial::factorization", tout << "upolynomial::factor_square_free_berlekamp("; upm.display(tout, f); tout << ", " << p << ") => " << factors << endl;);
TRACE(polynomial_factorization, tout << "upolynomial::factor_square_free_berlekamp("; upm.display(tout, f); tout << ", " << p << ") => " << factors << endl;);
return false;
}
TRACE("polynomial::factorization::bughunt", tout << "upolynomial::factor_square_free_berlekamp(): computing factors, expecting " << r << endl;);
TRACE(polynomial_factorization__bughunt, tout << "upolynomial::factor_square_free_berlekamp(): computing factors, expecting " << r << endl;);
scoped_numeral_vector gcd(zpm);
scoped_numeral_vector div(zpm);
@ -440,7 +440,7 @@ bool zp_factor_square_free_berlekamp(zp_manager & upm, numeral_vector const & f,
scoped_numeral_vector v_k(zpm);
while (Q_I.next_null_space_vector(v_k)) {
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "null vector: ";
for(unsigned j = 0; j < d; ++ j) {
tout << zpm.to_string(v_k[j]) << " ";
@ -449,7 +449,7 @@ bool zp_factor_square_free_berlekamp(zp_manager & upm, numeral_vector const & f,
);
upm.trim(v_k);
// TRACE("polynomial::factorization", tout << "v_k = "; upm.display(tout, v_k); tout << endl;);
// TRACE(polynomial_factorization, tout << "v_k = "; upm.display(tout, v_k); tout << endl;);
unsigned current_factor_end = factors.distinct_factors();
for (unsigned current_factor_i = first_factor; current_factor_i < current_factor_end; ++ current_factor_i) {
@ -479,7 +479,7 @@ bool zp_factor_square_free_berlekamp(zp_manager & upm, numeral_vector const & f,
// get the divisor also (no need to normalize the div, both are monic)
upm.div(current_factor.size(), current_factor.data(), gcd.size(), gcd.data(), div);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "current_factor = "; upm.display(tout, current_factor); tout << endl;
tout << "gcd_norm = "; upm.display(tout, gcd); tout << endl;
tout << "div = "; upm.display(tout, div); tout << endl;
@ -495,7 +495,7 @@ bool zp_factor_square_free_berlekamp(zp_manager & upm, numeral_vector const & f,
// at the point where we have all the factors, we are done
if (factors.distinct_factors() - first_factor == r) {
TRACE("polynomial::factorization", tout << "polynomial::factor("; upm.display(tout, f); tout << ", " << p << ") => " << factors << " of degree " << factors.get_degree() << endl;);
TRACE(polynomial_factorization, tout << "polynomial::factor("; upm.display(tout, f); tout << ", " << p << ") => " << factors << " of degree " << factors.get_degree() << endl;);
return true;
}
}
@ -533,7 +533,7 @@ bool check_hansel_lift(z_manager & upm, numeral_vector const & C,
upm.sub(C.size(), C.data(), test1.size(), test1.data(), test1);
to_zp_manager(br_upm, test1);
if (!test1.empty()) {
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "sage: R.<x> = ZZ['x']" << endl;
tout << "sage: A = "; upm.display(tout, A); tout << endl;
tout << "sage: B = "; upm.display(tout, B); tout << endl;
@ -586,7 +586,7 @@ void hensel_lift(z_manager & upm, numeral const & a, numeral const & b, numeral
z_numeral_manager & nm = upm.zm();
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "polynomial::hensel_lift(";
tout << "a = " << nm.to_string(a) << ", ";
tout << "b = " << nm.to_string(b) << ", ";
@ -611,7 +611,7 @@ void hensel_lift(z_manager & upm, numeral const & a, numeral const & b, numeral
upm.sub(C.size(), C.data(), f.size(), f.data(), f);
upm.div(f, b);
to_zp_manager(r_upm, f);
TRACE("polynomial::factorization",
TRACE(polynomial_factorization,
tout << "f = "; upm.display(tout, f); tout << endl;
);
@ -634,33 +634,33 @@ void hensel_lift(z_manager & upm, numeral const & a, numeral const & b, numeral
// compute the S, T (compute in Z_r[x])
scoped_numeral_vector Vf(r_upm.m()), t(r_upm.m()), S(r_upm.m());
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "V == "; upm.display(tout, V); tout << endl;
);
r_upm.mul(V.size(), V.data(), f.size(), f.data(), Vf);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "Vf = V*f" << endl;
tout << "Vf == "; upm.display(tout, Vf); tout << endl;
);
r_upm.div_rem(Vf.size(), Vf.data(), A.size(), A.data(), t, S);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "[t, S] = Vf.quo_rem(A)" << endl;
tout << "t == "; upm.display(tout, t); tout << endl;
tout << "S == "; upm.display(tout, S); tout << endl;
);
scoped_numeral_vector T(r_upm.m()), tmp(r_upm.m());
r_upm.mul(U.size(), U.data(), f.size(), f.data(), T); // T = fU
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "T == U*f" << endl;
tout << "T == "; upm.display(tout, T); tout << endl;
);
r_upm.mul(B.size(), B.data(), t.size(), t.data(), tmp); // tmp = Bt
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "tmp = B*t" << endl;
tout << "tmp == "; upm.display(tout, tmp); tout << endl;
);
r_upm.add(T.size(), T.data(), tmp.size(), tmp.data(), T); // T = Uf + Bt
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "T = B*tmp" << endl;
tout << "T == "; upm.display(tout, T); tout << endl;
);
@ -685,7 +685,7 @@ bool check_quadratic_hensel(zp_manager & zpe_upm, numeral_vector const & U, nume
scoped_mpz_vector one(nm);
zpe_upm.add(tmp1.size(), tmp1.data(), tmp2.size(), tmp2.data(), one);
if (one.size() != 1 || !nm.is_one(one[0])) {
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "sage: R.<x> = Zmod(" << nm.to_string(zpe_upm.m().p()) << ")['x']" << endl;
tout << "sage: A = "; zpe_upm.display(tout, A); tout << endl;
tout << "sage: B = "; zpe_upm.display(tout, B); tout << endl;
@ -710,7 +710,7 @@ void hensel_lift_quadratic(z_manager& upm, numeral_vector const & C,
zp_manager & zpe_upm, numeral_vector & A, numeral_vector & B, unsigned e) {
z_numeral_manager & nm = upm.zm();
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "polynomial::hansel_lift_quadratic(";
tout << "A = "; upm.display(tout, A); tout << ", ";
tout << "B = "; upm.display(tout, B); tout << ", ";
@ -741,7 +741,7 @@ void hensel_lift_quadratic(z_manager& upm, numeral_vector const & C,
hensel_lift(upm, pe, pe, pe, U, A, V, B, C, A_lifted, B_lifted);
// now we have C = AB (mod b*r)
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "A_lifted = "; upm.display(tout, A_lifted); tout << endl;
tout << "B_lifted = "; upm.display(tout, B_lifted); tout << endl;
tout << "C = "; upm.display(tout, C); tout << endl;
@ -768,7 +768,7 @@ void hensel_lift_quadratic(z_manager& upm, numeral_vector const & C,
upm.sub(g.size(), g.data(), tmp1.size(), tmp1.data(), g); // g = 1 - UA - VB
upm.div(g, pe);
to_zp_manager(zpe_upm, g);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "g = (1 - A_lifted*U - B_lifted*V)/" << nm.to_string(pe) << endl;
tout << "g == "; upm.display(tout, g); tout << endl;
);
@ -823,7 +823,7 @@ bool check_hensel_lift(z_manager & upm, numeral_vector const & f, zp_factors con
to_zp_manager(zp_upm, f, f_zp);
zp_upm.mul(mult_zp, f_zp.back());
if (!upm.eq(mult_zp, f_zp)) {
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "f = "; upm.display(tout, f); tout << endl;
tout << "zp_fs = " << zp_fs << endl;
tout << "sage: R.<x> = Zmod(" << nm.to_string(p) << ")['x']" << endl;
@ -845,7 +845,7 @@ bool check_hensel_lift(z_manager & upm, numeral_vector const & f, zp_factors con
to_zp_manager(zpe_upm, f, f_zpe);
zpe_upm.mul(mult_zpe, f_zpe.back());
if (!upm.eq(mult_zpe, f_zpe)) {
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "f = "; upm.display(tout, f); tout << endl;
tout << "zpe_fs = " << zpe_fs << endl;
tout << "sage: R.<x> = Zmod(" << nm.to_string(pe) << ")['x']" << endl;
@ -883,7 +883,7 @@ void hensel_lift(z_manager & upm, numeral_vector const & f, zp_factors const & z
zp_manager & zpe_upm = zpe_fs.upm();
zpe_nm.set_zp(zp_nm.p());
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "polynomial::hansel_lift("; upm.display(tout, f); tout << ", " << zp_fs << ")" << endl;
);
@ -903,7 +903,7 @@ void hensel_lift(z_manager & upm, numeral_vector const & f, zp_factors const & z
// F_i = A (mod Z_p)
zp_upm.set(zp_fs[i].size(), zp_fs[i].data(), A);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "A = "; upm.display(tout, A); tout << endl;
);
@ -918,13 +918,13 @@ void hensel_lift(z_manager & upm, numeral_vector const & f, zp_factors const & z
zp_nm.set(lc, f.back());
zp_upm.mul(C, lc);
}
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "C = "; upm.display(tout, C); tout << endl;
);
// we take B to be what's left from C and A
zp_upm.div(C.size(), C.data(), A.size(), A.data(), B);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "B = "; upm.display(tout, B); tout << endl;
);
@ -932,7 +932,7 @@ void hensel_lift(z_manager & upm, numeral_vector const & f, zp_factors const & z
zpe_nm.set_zp(zp_nm.p());
hensel_lift_quadratic(upm, f_parts, zpe_upm, A, B, e);
CASSERT("polynomial::factorizatio::bughunt", check_individual_lift(zp_upm, zp_fs[i], zpe_upm, A));
TRACE("polynomial::factorization",
TRACE(polynomial_factorization,
tout << "lifted to " << nm.to_string(zpe_upm.m().p()) << endl;
tout << "A = "; upm.display(tout, A); tout << endl;
tout << "B = "; upm.display(tout, B); tout << endl;
@ -957,7 +957,7 @@ void hensel_lift(z_manager & upm, numeral_vector const & f, zp_factors const & z
zpe_upm.mul(B, lc_inv);
zpe_fs.push_back_swap(B, 1);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "polynomial::hansel_lift("; upm.display(tout, f); tout << ", " << zp_fs << ") => " << zpe_fs << endl;
);
@ -1017,7 +1017,7 @@ static unsigned mignotte_bound(z_manager & upm, numeral_vector const & f, numera
This method also assumes f is primitive.
*/
bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs, unsigned k, factor_params const & params) {
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "sage: f = "; upm.display(tout, f); tout << endl;
tout << "sage: if (not f.is_squarefree()): print 'Error, f is not square-free'" << endl;
tout << "sage: print 'Factoring :', f" << endl;
@ -1049,7 +1049,7 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
}
}
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "sage: f_pp = "; upm.display(tout, f_pp); tout << endl;
tout << "sage: if (not (f_pp * 1 == f): print 'Error, content computation wrong'" << endl;
);
@ -1072,7 +1072,7 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
prime_iterator prime_it;
scoped_numeral gcd_tmp(nm);
unsigned trials = 0;
TRACE("polynomial::factorization::bughunt", tout << "trials: " << params.m_p_trials << "\n";);
TRACE(polynomial_factorization__bughunt, tout << "trials: " << params.m_p_trials << "\n";);
while (trials <= params.m_p_trials) {
upm.checkpoint();
// construct prime to check
@ -1086,7 +1086,7 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
// we need gcd(lc(f_pp), p) = 1
nm.gcd(p, f_pp.back(), gcd_tmp);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "sage: if (not (gcd(" << nm.to_string(p) << ", " << nm.to_string(f_pp.back()) << ")) == " <<
nm.to_string(gcd_tmp) << "): print 'Error, wrong gcd'" << endl;
);
@ -1098,7 +1098,7 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
scoped_numeral_vector f_pp_zp(nm);
to_zp_manager(zp_upm, f_pp, f_pp_zp);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "sage: Rp.<x_p> = GF(" << nm.to_string(p) << ")['x_p']"; tout << endl;
tout << "sage: f_pp_zp = "; zp_upm.display(tout, f_pp_zp, "x_p"); tout << endl;
);
@ -1139,7 +1139,7 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
zp_fs.swap(current_fs);
nm.set(zp_fs_p, p);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "best zp factorization (Z_" << nm.to_string(zp_fs_p) << "): ";
tout << zp_fs << endl;
tout << "best degree set: "; degree_set.display(tout); tout << endl;
@ -1153,7 +1153,7 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
// make sure to set the zp_manager back to modulo zp_fs_p
zp_upm.set_zp(zp_fs_p);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "best zp factorization (Z_" << nm.to_string(zp_fs_p) << "): " << zp_fs << endl;
tout << "best degree set: "; degree_set.display(tout); tout << endl;
);
@ -1161,7 +1161,7 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
// get a bound on B for the factors of f_pp with degree less or equal to deg(f)/2
// and then choose e to be smallest such that p^e > 2*lc(f)*B, we use the mignotte
unsigned e = mignotte_bound(upm, f_pp, zp_fs_p);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "out p = " << nm.to_string(zp_fs_p) << ", and we'll work p^e for e = " << e << endl;
);
@ -1190,7 +1190,7 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
ufactorization_combination_iterator it(zpe_fs, degree_set);
scoped_numeral_vector trial_factor(nm), trial_factor_quo(nm);
scoped_numeral trial_factor_cont(nm);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "STARTING TRIAL DIVISION" << endl;
tout << "zpe_fs" << zpe_fs << endl;
tout << "degree_set = "; degree_set.display(tout); tout << endl;
@ -1237,14 +1237,14 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
// add the lc(f_pp) to the trial divisor
zpe_upm.mul(trial_factor, f_pp_lc);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "f_pp*lc(f_pp) = "; upm.display(tout, f_pp); tout << endl;
tout << "trial_factor = "; upm.display(tout, trial_factor); tout << endl;
);
bool true_factor = upm.exact_div(f_pp, trial_factor, trial_factor_quo);
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "trial_factor = "; upm.display(tout, trial_factor); tout << endl;
tout << "trial_factor_quo = "; upm.display(tout, trial_factor_quo); tout << endl;
tout << "result = " << (true_factor ? "true" : "false") << endl;
@ -1272,7 +1272,7 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
// don't remove this combination
remove = false;
}
TRACE("polynomial::factorization::bughunt",
TRACE(polynomial_factorization__bughunt,
tout << "factors = " << fs << endl;
tout << "f_pp*lc(f_pp) = "; upm.display(tout, f_pp); tout << endl;
tout << "lc(f_pp) = " << f_pp_lc << endl;