-
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathinequation.cpp
170 lines (147 loc) · 7.81 KB
/
inequation.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
#include "matchit.h"
#include "internal.h"
#include "mathiu/inequation.h"
#include "mathiu/diff.h"
#include "mathiu/solve.h"
#include "mathiu/polynomial.h"
#include "mathiu/setOp.h"
#include <numeric>
namespace mathiu
{
namespace impl
{
ExprPtr solveInequationImpl(ExprPtr const& ex, RelationalKind relKind, ExprPtr const& var, ExprPtr const& domain)
{
#if DEBUG
std::cout << "solveInequationImpl: " << toString(ex) << ",\t" << toString(relKind) << ",\t" << toString(var) << ",\t" << toString(domain) << std::endl;
#endif // DEBUG
using namespace matchit;
const auto freeOfVar = meet([&](auto&& e) { return freeOf(e, var); });
Id<PieceWise> iPieceWise;
auto result = match(ex)(
pattern | freeOfVar = [&]{
throw std::logic_error("Cannot solve inequation!");
return relational(relKind, ex, 0_i);
},
pattern | some(as<PieceWise>(iPieceWise)) = [&]
{
return std::accumulate((*iPieceWise).begin(), (*iPieceWise).end(), false_, [&](auto&& result, auto&& e)
{
return union_(result, solveInequationImpl(e.first, relKind, var, solveInequation(e.second, var, domain)));
});
},
pattern | _ = [&]
{
auto coeffList = coefficientList(ex, var);
Id<ExprPtr> c, b, a;
return match(std::get<List>(*coeffList))(
pattern | ds(c, b, a.at(some(isRational))) = [&]
{
auto a_ = evald(*a);
auto const rKind = a_ > 0 ? relKind : invertRelational(relKind);
auto _2a = 2_i * (*a);
auto b2_4ac = ((*b)^2_i) - 4_i * (*a) * (*c);
auto sqrtB2_4ac = sqrt(b2_4ac);
auto result = match(b2_4ac)(
pattern | asDouble | when([&]
{ return evald(b2_4ac) < 0; }) = [&]
{
switch (rKind)
{
case RelationalKind::kLESS:
return false_;
case RelationalKind::kLESS_EQUAL:
return false_;
case RelationalKind::kGREATER:
return true_;
case RelationalKind::kGREATER_EQUAL:
return true_;
case RelationalKind::kEQUAL:
throw std::logic_error{"Unreachable!"};
return false_;
}
throw std::logic_error{"Unreachable!"};
return false_;
},
pattern | _ = [&]
{
auto left = (-(*b) - sqrtB2_4ac) / _2a;
auto right = (-(*b) + sqrtB2_4ac) / _2a;
switch (rKind)
{
case RelationalKind::kLESS:
return interval(left, false, right, false);
case RelationalKind::kLESS_EQUAL:
return interval(left, true, right, true);
case RelationalKind::kGREATER:
return makeSharedExprPtr(SetOp{Union{{interval(-infinity, false, left, false),
interval(right, false, infinity, false)}}});
case RelationalKind::kGREATER_EQUAL:
return makeSharedExprPtr(SetOp{Union{{interval(-infinity, false, left, true),
interval(right, true, infinity, false)}}});
case RelationalKind::kEQUAL:
throw std::logic_error{"Unreachable!"};
}
throw std::logic_error{"Unreachable!"};
return false_;
});
return intersect(result, domain);
},
pattern | ds(c, b.at(some(isRational))) = [&]
{
auto const root = -(*c) / (*b);
auto b_ = evald(*b);
auto const rKind = b_ > 0 ? relKind : invertRelational(relKind);
auto result = [&]
{
switch (rKind)
{
case RelationalKind::kLESS:
return interval(-infinity, false, root, false);
case RelationalKind::kLESS_EQUAL:
return interval(-infinity, false, root, true);
case RelationalKind::kGREATER:
return interval(root, false, infinity, false);
case RelationalKind::kGREATER_EQUAL:
return interval(root, true, infinity, false);
case RelationalKind::kEQUAL:
throw std::logic_error{"Unreachable!"};
}
throw std::logic_error{"Unreachable!"};
return false_;
}();
return intersect(result, domain);
},
pattern | _ = [&]
{
throw std::runtime_error{"No match in solve!"};
return false_;
});
throw std::runtime_error{"Not implemented yet!"};
return false_;
});
#if DEBUG
std::cout << "solveInequationImpl: " << toString(ex) << ",\t" << toString(relKind) << ",\t" << toString(var) << ",\t" << toString(domain) << ",\tresult: " << toString(result) << std::endl;
#endif // DEBUG
return result;
}
ExprPtr solveInequation(ExprPtr const& ex, ExprPtr const& var, ExprPtr const& domain)
{
#if DEBUG
std::cout << "solveInequation: " << toString(ex) << ",\t" << toString(var) << ",\t" << toString(domain) << std::endl;
#endif // DEBUG
auto const var_ = std::get<Symbol>(*var);
using namespace matchit;
Id<RelationalKind> iRelKind;
Id<ExprPtr> iE1, iE2;
auto result = match(ex)(
pattern | or_(true_, false_) = expr(domain),
pattern | some(as<Relational>(ds(iRelKind, iE1, iE2))) = [&]
{ return solveInequationImpl(expand(*iE1 - *iE2), *iRelKind, var, domain); });
#if DEBUG
std::cout << "solveInequation: " << toString(ex) << ",\t" << toString(var) << ",\t" << toString(domain) << ",\tresult: " << toString(result) << std::endl;
#endif // DEBUG
return result;
}
} // namespace impl
} // namespace mathiu