1 /*++
2 Copyright (c) 2017 Microsoft Corporation
3
4 Author:
5 Lev Nachmanson (levnach)
6 Nikolaj Bjorner (nbjorner)
7
8 --*/
9 #include "math/lp/nla_tangent_lemmas.h"
10 #include "math/lp/nla_core.h"
11
12 namespace nla {
13
14 class tangent_imp {
15 point m_a;
16 point m_b;
17 point m_xy;
18 rational m_correct_v;
19 // "below" means that the incorrect value is less than the correct one, that is m_v < m_correct_v
20 bool m_below;
21 rational m_v; // the monomial value
22 lpvar m_j; // the monic variable
23 const monic& m_m;
24 const factor& m_x;
25 const factor& m_y;
26 lpvar m_jx;
27 lpvar m_jy;
28 tangents& m_tang;
29 bool m_is_mon;
30
31 public:
tangent_imp(point xy,const rational & v,const monic & m,const factorization & f,tangents & tang)32 tangent_imp(point xy,
33 const rational& v,
34 const monic& m,
35 const factorization& f,
36 tangents& tang) : m_xy(xy),
37 m_correct_v(xy.x * xy.y),
38 m_below(v < m_correct_v),
39 m_v(v),
40 m_j(m.var()),
41 m_m(m),
42 m_x(f[0]),
43 m_y(f[1]),
44 m_jx(m_x.var()),
45 m_jy(m_y.var()),
46 m_tang(tang),
47 m_is_mon(f.is_mon()) {
48 SASSERT(f.size() == 2);
49 }
50
operator ()()51 void operator()() {
52 get_points();
53 TRACE("nla_solver", print_tangent_domain(tout << "tang domain = ") << std::endl;);
54 generate_line1();
55 generate_line2();
56 generate_plane(m_a);
57 generate_plane(m_b);
58 }
59
60 private:
61
c()62 core & c() { return m_tang.c(); }
63
explain(new_lemma & lemma)64 void explain(new_lemma& lemma) {
65 if (!m_is_mon) {
66 lemma &= m_m;
67 lemma &= m_x;
68 lemma &= m_y;
69 }
70 }
71
generate_plane(const point & pl)72 void generate_plane(const point & pl) {
73 new_lemma lemma(c(), "generate tangent plane");
74 c().negate_relation(lemma, m_jx, m_x.rat_sign()*pl.x);
75 c().negate_relation(lemma, m_jy, m_y.rat_sign()*pl.y);
76 #if Z3DEBUG
77 SASSERT(c().val(m_x) == m_xy.x && c().val(m_y) == m_xy.y);
78 int mult_sign = nla::rat_sign(pl.x - m_xy.x)*nla::rat_sign(pl.y - m_xy.y);
79 SASSERT((mult_sign == 1) == m_below);
80 // If "mult_sign is 1" then (a - x)(b-y) > 0 and ab - bx - ay + xy > 0
81 // or -ab + bx + ay < xy or -ay - bx + xy > -ab
82 // val(j) stands for xy. So, finally we have -ay - bx + j > - ab
83 #endif
84
85 lp::lar_term t;
86 t.add_monomial(- m_y.rat_sign()*pl.x, m_jy);
87 t.add_monomial(- m_x.rat_sign()*pl.y, m_jx);
88 t.add_var(m_j);
89 lemma |= ineq(t, m_below? llc::GT : llc::LT, - pl.x*pl.y);
90 explain(lemma);
91 }
92
generate_line1()93 void generate_line1() {
94 new_lemma lemma(c(), "tangent line 1");
95 // Should be v = val(m_x)*val(m_y), and val(factor) = factor.rat_sign()*var(factor.var())
96 lemma |= ineq(m_jx, llc::NE, c().val(m_jx));
97 lemma |= ineq(lp::lar_term(m_j, - m_y.rat_sign() * m_xy.x, m_jy), llc::EQ, 0);
98 explain(lemma);
99 }
100
generate_line2()101 void generate_line2() {
102 new_lemma lemma(c(), "tangent line 2");
103 lemma |= ineq(m_jy, llc::NE, c().val(m_jy));
104 lemma |= ineq(lp::lar_term(m_j, - m_x.rat_sign() * m_xy.y, m_jx), llc::EQ, 0);
105 explain(lemma);
106 }
107
108 // Get two planes tangent to surface z = xy, one at point a, and another at point b, creating a cut
get_initial_points()109 void get_initial_points() {
110 const rational& x = m_xy.x;
111 const rational& y = m_xy.y;
112 bool all_ints = m_v.is_int() && x.is_int() && y.is_int();
113 rational delta = rational(1);
114 if (!all_ints )
115 delta = std::min(delta, abs(m_correct_v - m_v));
116 TRACE("nla_solver", tout << "delta = " << delta << "\n";);
117 if (!m_below){
118 m_a = point(x - delta, y + delta);
119 m_b = point(x + delta, y - delta);
120 }
121 else {
122 // denote x = xy.x and y = xy.y, and vx, vy - the values of x and y.
123 // we have val(xy) < vx*y + vy*x - vx*vy = pl(x, y);
124 // The plane with delta (1, 1) is (vx + 1)y + (vy + 1)x - (vx + 1)(vy + 1) =
125 // vx*y + vy*x - vx*vy + y + x - xv*vy - vx - vy - 1 = pl(x, y) - 1
126 // For integers the last expression is greater than or equal to val(xy) when x = vx and y = vy.
127 // If x <= vx+1 and y <= vy+1 then (vx+1-x)*(vy+1-y) > 0, that creates a cut
128 // - (vx + 1)y - (vy + 1)x + xy > - (vx+1)*(vx+1).
129 // If all_ints is false then we use the fact that
130 // tang_plane() will not change more than on delta*delta
131 m_a = point(x - delta, y - delta);
132 m_b = point(x + delta, y + delta);
133 }
134 }
135
push_point(point & a)136 void push_point(point & a) {
137 SASSERT(plane_is_correct_cut(a));
138 int steps = 10;
139 point del = a - m_xy;
140 while (steps-- && !c().done()) {
141 del *= rational(2);
142 point na = m_xy + del;
143 TRACE("nla_solver_tp", tout << "del = " << del << std::endl;);
144 if (!plane_is_correct_cut(na)) {
145 TRACE("nla_solver_tp", tout << "exit\n";);
146 return;
147 }
148 a = na;
149 }
150 }
151
tang_plane(const point & a) const152 rational tang_plane(const point& a) const {
153 return a.x * m_xy.y + a.y * m_xy.x - a.x * a.y;
154 }
155
get_points()156 void get_points() {
157 get_initial_points();
158 TRACE("nla_solver", tout << "xy = " << m_xy << ", correct val = " << m_correct_v;
159 print_tangent_domain(tout << "\ntang points:") << std::endl;);
160 push_point(m_a);
161 push_point(m_b);
162 TRACE("nla_solver",
163 tout << "pushed a = " << m_a << std::endl
164 << "pushed b = " << m_b << std::endl
165 << "tang_plane(a) = " << tang_plane(m_a) << " , val = " << m_a << ", "
166 << "tang_plane(b) = " << tang_plane(m_b) << " , val = " << m_b << std::endl;);
167 }
168
print_tangent_domain(std::ostream & out)169 std::ostream& print_tangent_domain(std::ostream& out) {
170 return out << "(" << m_a << ", " << m_b << ")";
171 }
172
plane_is_correct_cut(const point & plane) const173 bool plane_is_correct_cut(const point& plane) const {
174 TRACE("nla_solver", tout << "plane = " << plane << "\n";
175 tout << "tang_plane() = " << tang_plane(plane) << ", v = " << m_v << ", correct_v = " << m_correct_v << "\n";);
176 SASSERT((m_below && m_v < m_correct_v) ||
177 ((!m_below) && m_v > m_correct_v));
178 rational sign = rational(m_below ? 1 : -1);
179 rational px = tang_plane(plane);
180 return ((m_correct_v - px)*sign).is_pos() && !((px - m_v)*sign).is_neg();
181 }
182 };
183
tangents(core * c)184 tangents::tangents(core * c) : common(c) {}
185
tangent_lemma()186 void tangents::tangent_lemma() {
187 factorization bf(nullptr);
188 const monic* m = nullptr;
189 if (c().m_nla_settings.run_tangents() && c().find_bfc_to_refine(m, bf)) {
190 lpvar j = m->var();
191 tangent_imp tangent(point(val(bf[0]), val(bf[1])), c().val(j), *m, bf, *this);
192 tangent();
193 }
194 }
195
196
197 }
198