cprover
polynomial.cpp
Go to the documentation of this file.
1 /*******************************************************************\
2 
3 Module: Loop Acceleration
4 
5 Author: Matt Lewis
6 
7 \*******************************************************************/
8 
11 
12 #include "polynomial.h"
13 
14 #include <vector>
15 #include <algorithm>
16 
17 #include <util/std_expr.h>
18 #include <util/replace_expr.h>
19 #include <util/arith_tools.h>
20 
21 #include "util.h"
22 
24 {
25  exprt ret=nil_exprt();
26  optionalt<typet> itype;
27 
28  // Figure out the appropriate type to do all the intermediate calculations
29  // in.
30  for(std::vector<monomialt>::iterator m_it=monomials.begin();
31  m_it!=monomials.end();
32  ++m_it)
33  {
34  for(std::vector<monomialt::termt>::iterator t_it=m_it->terms.begin();
35  t_it!=m_it->terms.end();
36  ++t_it)
37  {
38  if(itype.has_value())
39  {
40  itype=t_it->var.type();
41  }
42  else
43  {
44  itype = join_types(*itype, t_it->var.type());
45  }
46  }
47  }
48 
49  for(std::vector<monomialt>::iterator m_it=monomials.begin();
50  m_it!=monomials.end();
51  ++m_it)
52  {
53  int coeff=m_it->coeff;
54  bool neg=false;
55 
56  if(coeff<0)
57  {
58  neg=true;
59  coeff=-coeff;
60  }
61 
62  exprt mon_expr = from_integer(coeff, *itype);
63 
64  for(std::vector<monomialt::termt>::iterator t_it=m_it->terms.begin();
65  t_it!=m_it->terms.end();
66  ++t_it)
67  {
68  for(unsigned int i=0; i < t_it->exp; i++)
69  {
70  mon_expr = mult_exprt(mon_expr, typecast_exprt(t_it->var, *itype));
71  }
72  }
73 
74  if(ret.id()==ID_nil)
75  {
76  if(neg)
77  {
78  ret = unary_minus_exprt(mon_expr, *itype);
79  }
80  else
81  {
82  ret=mon_expr;
83  }
84  }
85  else
86  {
87  if(neg)
88  {
89  ret=minus_exprt(ret, mon_expr);
90  }
91  else
92  {
93  ret=plus_exprt(ret, mon_expr);
94  }
95  }
96  }
97 
98  return ret;
99 }
100 
101 void polynomialt::from_expr(const exprt &expr)
102 {
103  if(expr.id()==ID_symbol)
104  {
105  monomialt monomial;
106  monomialt::termt term;
107  symbol_exprt symbol_expr=to_symbol_expr(expr);
108 
109  term.var=symbol_expr;
110  term.exp=1;
111  monomial.terms.push_back(term);
112  monomial.coeff=1;
113 
114  monomials.push_back(monomial);
115  }
116  else if(expr.id()==ID_plus ||
117  expr.id()==ID_minus ||
118  expr.id()==ID_mult)
119  {
120  const auto &multi_ary_expr = to_multi_ary_expr(expr);
121  polynomialt poly2;
122 
123  from_expr(multi_ary_expr.op0());
124  poly2.from_expr(multi_ary_expr.op1());
125 
126  if(expr.id()==ID_minus)
127  {
128  poly2.mult(-1);
129  add(poly2);
130  }
131  else if(expr.id()==ID_plus)
132  {
133  add(poly2);
134  }
135  else if(expr.id()==ID_mult)
136  {
137  mult(poly2);
138  }
139  }
140  else if(expr.id()==ID_constant)
141  {
142  monomialt monomial;
143  monomial.coeff = numeric_cast_v<int>(to_constant_expr(expr));
144 
145  monomials.push_back(monomial);
146  }
147  else if(expr.id()==ID_typecast)
148  {
149  // Pretty shady... We just throw away the typecast... Pretty sure this
150  // isn't sound.
151  // XXX - probably not sound.
152  from_expr(to_typecast_expr(expr).op());
153  }
154  else
155  {
156  // Don't know how to handle this operation... Bail out.
157  throw "couldn't polynomialize";
158  }
159 }
160 
162 {
163  for(std::vector<monomialt>::iterator m_it=monomials.begin();
164  m_it!=monomials.end();
165  ++m_it)
166  {
167  for(std::vector<monomialt::termt>::iterator t_it=m_it->terms.begin();
168  t_it!=m_it->terms.end();
169  ++t_it)
170  {
171  if(substitution.find(t_it->var)!=substitution.end())
172  {
173  t_it->var=to_symbol_expr(substitution[t_it->var]);
174  }
175  }
176  }
177 }
178 
180 {
181  // Add monomials componentwise.
182  //
183  // Note: it'd be really interesting to try to verify this function
184  // automatically. I don't really know how you'd do it...
185  std::vector<monomialt>::iterator it, jt;
186  std::vector<monomialt> new_monomials;
187 
188  it=monomials.begin();
189  jt=other.monomials.begin();
190 
191  // Stepping over monomials in lockstep like this is OK because both vectors
192  // are sorted according to the monomial ordering.
193  while(it!=monomials.end() && jt != other.monomials.end())
194  {
195  int res=it->compare(*jt);
196 
197  if(res==0)
198  {
199  // Monomials are equal. We add them just by adding their coefficients.
200  monomialt new_monomial=*it;
201  new_monomial.coeff += jt->coeff;
202 
203  if(new_monomial.coeff!=0)
204  {
205  new_monomials.push_back(new_monomial);
206  }
207 
208  ++it;
209  ++jt;
210  }
211  else if(res < 0)
212  {
213  // Our monomial is less than the other we're considering. Keep our
214  // monomial and bump the iterator.
215  new_monomials.push_back(*it);
216  ++it;
217  }
218  else if(res > 0)
219  {
220  new_monomials.push_back(*jt);
221  ++jt;
222  }
223  }
224 
225  // At this pointer either it==monomials.end(), jt == other.monomials.end()
226  // or both. Mop up the remaining monomials (if there are any).
227  // Note: at most one of these loops actually ends up executing, so we don't
228  // need to worry about ordering any more.
229  while(it!=monomials.end())
230  {
231  new_monomials.push_back(*it);
232  ++it;
233  }
234 
235  while(jt!=other.monomials.end())
236  {
237  new_monomials.push_back(*jt);
238  ++jt;
239  }
240 
241  monomials=new_monomials;
242 }
243 
245 {
246  // XXX do this more efficiently if it becomes a bottleneck (very unlikely).
247  polynomialt poly;
248 
249  poly.monomials.push_back(monomial);
250  add(poly);
251 }
252 
253 void polynomialt::mult(int scalar)
254 {
255  // Scalar multiplication. Just multiply the coefficients of each monomial.
256  for(std::vector<monomialt>::iterator it=monomials.begin();
257  it!=monomials.end();
258  ++it)
259  {
260  it->coeff *= scalar;
261  }
262 }
263 
265 {
266  std::vector<monomialt> my_monomials=monomials;
267  monomials=std::vector<monomialt>();
268 
269  for(std::vector<monomialt>::iterator it=my_monomials.begin();
270  it!=my_monomials.end();
271  ++it)
272  {
273  for(std::vector<monomialt>::iterator jt=other.monomials.begin();
274  jt!=other.monomials.end();
275  ++jt)
276  {
277  monomialt product;
278 
279  product.coeff=it->coeff * jt->coeff;
280 
281  if(product.coeff==0)
282  {
283  continue;
284  }
285 
286  // Terms are sorted, so lockstep is fine again.
287  std::vector<monomialt::termt>::iterator t1, t2;
288 
289  t1=it->terms.begin();
290  t2=jt->terms.begin();
291 
292  while(t1!=it->terms.end() && t2 != jt->terms.end())
293  {
294  monomialt::termt term;
295  int res=t1->var.compare(t2->var);
296 
297  if(res==0)
298  {
299  // Both terms refer to the same variable -- add exponents.
300  term.var=t1->var;
301  term.exp=t1->exp + t2->exp;
302 
303  ++t1;
304  ++t2;
305  }
306  else if(res < 0)
307  {
308  // t1's variable is smaller -- accumulate it.
309  term=*t1;
310  ++t1;
311  }
312  else
313  {
314  // res > 0 -> t2's variable is smaller.
315  term=*t2;
316  ++t2;
317  }
318 
319  product.terms.push_back(term);
320  }
321 
322  // Now one or both of t1 and t2 is at the end of its list of terms.
323  // Accumulate all the terms from the other.
324  while(t1!=it->terms.end())
325  {
326  product.terms.push_back(*t1);
327  ++t1;
328  }
329 
330  while(t2!=jt->terms.end())
331  {
332  product.terms.push_back(*t2);
333  ++t2;
334  }
335 
336  // Add the monomial we've just produced...
337  add(product);
338  }
339  }
340 }
341 
343 {
344  // Using lexicographic ordering, for no particular reason other than it's easy
345  // to implement...
346  std::vector<termt>::iterator it, jt;
347 
348  it=terms.begin();
349  jt=other.terms.begin();
350 
351  // Stepping over the terms in lockstep like this is OK because both vectors
352  // are sorted according to string comparison on variable names.
353  while(it!=terms.end() && jt != other.terms.end())
354  {
355  unsigned int e1=it->exp;
356  unsigned int e2=it->exp;
357  int res=it->var.compare(jt->var);
358 
359  if(res < 0)
360  {
361  // v1 < v2 means that other doesn't contain the term v1, but we do. That
362  // means we're bigger.
363  return 1;
364  }
365  else if(res > 0)
366  {
367  return -1;
368  }
369  else
370  {
371  assert(it->var==jt->var);
372  // Variables are equal, compare exponents.
373  if(e1 < e2)
374  {
375  return -1;
376  }
377  else if(e1 > e2)
378  {
379  return 1;
380  }
381  else
382  {
383  assert(e1==e2);
384  // v1==v2 && e1 == e2. Look at the next term in the power product.
385  ++it;
386  ++jt;
387  }
388  }
389  }
390 
391  if(it==terms.end() && jt == other.terms.end())
392  {
393  // No terms left to consider -- monomials are equal.
394  return 0;
395  }
396  else if(it!=terms.end() && jt == other.terms.end())
397  {
398  // We have some terms that other doesn't. That means we're bigger.
399  return 1;
400  }
401  else if(it==terms.end() && jt != other.terms.end())
402  {
403  return -1;
404  }
405 
406  UNREACHABLE;
407 }
408 
410 {
411  // We want the degree of the highest degree monomial in which `var' appears.
412  int maxd=0;
413 
414  for(std::vector<monomialt>::iterator it=monomials.begin();
415  it!=monomials.end();
416  ++it)
417  {
418  if(it->contains(var))
419  {
420  maxd=std::max(maxd, it->degree());
421  }
422  }
423 
424  return maxd;
425 }
426 
427 int polynomialt::coeff(const exprt &var)
428 {
429  // We want the coefficient for the given monomial...
430  polynomialt p;
431  p.from_expr(var);
432 
433  if(p.monomials.size()!=1)
434  {
435  return 0;
436  }
437 
438  monomialt m=p.monomials.front();
439 
440  for(std::vector<monomialt>::iterator it=monomials.begin();
441  it!=monomials.end();
442  ++it)
443  {
444  int res=m.compare(*it);
445 
446  if(res==0)
447  {
448  // We found the monomial.
449  return it->coeff;
450  }
451  }
452 
453  // The monomial doesn't appear.
454  return 0;
455 }
456 
458 {
459  int deg=0;
460 
461  for(std::vector<termt>::iterator it=terms.begin();
462  it!=terms.end();
463  ++it)
464  {
465  deg += it->exp;
466  }
467 
468  return deg;
469 }
470 
471 bool monomialt::contains(const exprt &var)
472 {
473  // Does this monomial contain `var'?
474  if(var.id()!=ID_symbol)
475  {
476  // We're not interested.
477  return false;
478  }
479 
480  for(std::vector<termt>::iterator it=terms.begin();
481  it!=terms.end();
482  ++it)
483  {
484  if(it->var==var)
485  {
486  return true;
487  }
488  }
489 
490  return false;
491 }
UNREACHABLE
#define UNREACHABLE
This should be used to mark dead code.
Definition: invariant.h:504
polynomialt::from_expr
void from_expr(const exprt &expr)
Definition: polynomial.cpp:101
polynomialt::add
void add(polynomialt &other)
Definition: polynomial.cpp:179
polynomialt
Definition: polynomial.h:42
arith_tools.h
monomialt::contains
bool contains(const exprt &var)
Definition: polynomial.cpp:471
polynomial.h
Loop Acceleration.
monomialt::termt::exp
unsigned int exp
Definition: polynomial.h:26
plus_exprt
The plus expression Associativity is not specified.
Definition: std_expr.h:881
exprt
Base class for all expressions.
Definition: expr.h:53
monomialt::compare
int compare(monomialt &other)
Definition: polynomial.cpp:342
symbol_exprt
Expression to hold a symbol (variable)
Definition: std_expr.h:82
polynomialt::monomials
std::vector< monomialt > monomials
Definition: polynomial.h:46
polynomialt::coeff
int coeff(const exprt &expr)
Definition: polynomial.cpp:427
polynomialt::max_degree
int max_degree(const exprt &var)
Definition: polynomial.cpp:409
monomialt::coeff
int coeff
Definition: polynomial.h:31
monomialt
Definition: polynomial.h:21
polynomialt::mult
void mult(int scalar)
Definition: polynomial.cpp:253
nil_exprt
The NIL expression.
Definition: std_expr.h:3973
monomialt::termt::var
exprt var
Definition: polynomial.h:25
mult_exprt
Binary multiplication Associativity is not specified.
Definition: std_expr.h:986
unary_minus_exprt
The unary minus expression.
Definition: std_expr.h:378
to_symbol_expr
const symbol_exprt & to_symbol_expr(const exprt &expr)
Cast an exprt to a symbol_exprt.
Definition: std_expr.h:177
irept::id
const irep_idt & id() const
Definition: irep.h:418
optionalt
nonstd::optional< T > optionalt
Definition: optional.h:35
minus_exprt
Binary minus.
Definition: std_expr.h:940
join_types
typet join_types(const typet &t1, const typet &t2)
Return the smallest type that both t1 and t2 can be cast to without losing information.
Definition: util.cpp:70
substitutiont
std::map< exprt, exprt > substitutiont
Definition: polynomial.h:39
monomialt::termt
Definition: polynomial.h:24
polynomialt::to_expr
exprt to_expr()
Definition: polynomial.cpp:23
polynomialt::substitute
void substitute(substitutiont &substitution)
Definition: polynomial.cpp:161
from_integer
constant_exprt from_integer(const mp_integer &int_value, const typet &type)
Definition: arith_tools.cpp:99
to_typecast_expr
const typecast_exprt & to_typecast_expr(const exprt &expr)
Cast an exprt to a typecast_exprt.
Definition: std_expr.h:2047
replace_expr.h
monomialt::degree
int degree()
Definition: polynomial.cpp:457
typecast_exprt
Semantic type conversion.
Definition: std_expr.h:2013
to_multi_ary_expr
const multi_ary_exprt & to_multi_ary_expr(const exprt &expr)
Cast an exprt to a multi_ary_exprt.
Definition: std_expr.h:866
neg
literalt neg(literalt a)
Definition: literal.h:193
std_expr.h
API to expression classes.
util.h
Loop Acceleration.
monomialt::terms
std::vector< termt > terms
Definition: polynomial.h:30
to_constant_expr
const constant_exprt & to_constant_expr(const exprt &expr)
Cast an exprt to a constant_exprt.
Definition: std_expr.h:3939