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  typet itype=nil_typet();
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==nil_typet())
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  polynomialt poly2;
121 
122  from_expr(expr.op0());
123  poly2.from_expr(expr.op1());
124 
125  if(expr.id()==ID_minus)
126  {
127  poly2.mult(-1);
128  add(poly2);
129  }
130  else if(expr.id()==ID_plus)
131  {
132  add(poly2);
133  }
134  else if(expr.id()==ID_mult)
135  {
136  mult(poly2);
137  }
138  }
139  else if(expr.id()==ID_constant)
140  {
141  mp_integer mp;
142  unsigned int l;
143  constant_exprt const_expr=to_constant_expr(expr);
144 
145  mp=binary2integer(const_expr.get_value().c_str(), true);
146  l=mp.to_long();
147 
148  monomialt monomial;
149  monomial.coeff=l;
150 
151  monomials.push_back(monomial);
152  }
153  else if(expr.id()==ID_typecast)
154  {
155  // Pretty shady... We just throw away the typecast... Pretty sure this
156  // isn't sound.
157  // XXX - probably not sound.
158  from_expr(expr.op0());
159  }
160  else
161  {
162  // Don't know how to handle this operation... Bail out.
163  throw "couldn't polynomialize";
164  }
165 }
166 
168 {
169  for(std::vector<monomialt>::iterator m_it=monomials.begin();
170  m_it!=monomials.end();
171  ++m_it)
172  {
173  for(std::vector<monomialt::termt>::iterator t_it=m_it->terms.begin();
174  t_it!=m_it->terms.end();
175  ++t_it)
176  {
177  if(substitution.find(t_it->var)!=substitution.end())
178  {
179  t_it->var=to_symbol_expr(substitution[t_it->var]);
180  }
181  }
182  }
183 }
184 
186 {
187  // Add monomials componentwise.
188  //
189  // Note: it'd be really interesting to try to verify this function
190  // automatically. I don't really know how you'd do it...
191  std::vector<monomialt>::iterator it, jt;
192  std::vector<monomialt> new_monomials;
193 
194  it=monomials.begin();
195  jt=other.monomials.begin();
196 
197  // Stepping over monomials in lockstep like this is OK because both vectors
198  // are sorted according to the monomial ordering.
199  while(it!=monomials.end() && jt != other.monomials.end())
200  {
201  int res=it->compare(*jt);
202 
203  if(res==0)
204  {
205  // Monomials are equal. We add them just by adding their coefficients.
206  monomialt new_monomial=*it;
207  new_monomial.coeff += jt->coeff;
208 
209  if(new_monomial.coeff!=0)
210  {
211  new_monomials.push_back(new_monomial);
212  }
213 
214  ++it;
215  ++jt;
216  }
217  else if(res < 0)
218  {
219  // Our monomial is less than the other we're considering. Keep our
220  // monomial and bump the iterator.
221  new_monomials.push_back(*it);
222  ++it;
223  }
224  else if(res > 0)
225  {
226  new_monomials.push_back(*jt);
227  ++jt;
228  }
229  }
230 
231  // At this pointer either it==monomials.end(), jt == other.monomials.end()
232  // or both. Mop up the remaining monomials (if there are any).
233  // Note: at most one of these loops actually ends up executing, so we don't
234  // need to worry about ordering any more.
235  while(it!=monomials.end())
236  {
237  new_monomials.push_back(*it);
238  ++it;
239  }
240 
241  while(jt!=other.monomials.end())
242  {
243  new_monomials.push_back(*jt);
244  ++jt;
245  }
246 
247  monomials=new_monomials;
248 }
249 
251 {
252  // XXX do this more efficiently if it becomes a bottleneck (very unlikely).
253  polynomialt poly;
254 
255  poly.monomials.push_back(monomial);
256  add(poly);
257 }
258 
259 void polynomialt::mult(int scalar)
260 {
261  // Scalar multiplication. Just multiply the coefficients of each monomial.
262  for(std::vector<monomialt>::iterator it=monomials.begin();
263  it!=monomials.end();
264  ++it)
265  {
266  it->coeff *= scalar;
267  }
268 }
269 
271 {
272  std::vector<monomialt> my_monomials=monomials;
273  monomials=std::vector<monomialt>();
274 
275  for(std::vector<monomialt>::iterator it=my_monomials.begin();
276  it!=my_monomials.end();
277  ++it)
278  {
279  for(std::vector<monomialt>::iterator jt=other.monomials.begin();
280  jt!=other.monomials.end();
281  ++jt)
282  {
283  monomialt product;
284 
285  product.coeff=it->coeff * jt->coeff;
286 
287  if(product.coeff==0)
288  {
289  continue;
290  }
291 
292  // Terms are sorted, so lockstep is fine again.
293  std::vector<monomialt::termt>::iterator t1, t2;
294 
295  t1=it->terms.begin();
296  t2=jt->terms.begin();
297 
298  while(t1!=it->terms.end() && t2 != jt->terms.end())
299  {
300  monomialt::termt term;
301  int res=t1->var.compare(t2->var);
302 
303  if(res==0)
304  {
305  // Both terms refer to the same variable -- add exponents.
306  term.var=t1->var;
307  term.exp=t1->exp + t2->exp;
308 
309  ++t1;
310  ++t2;
311  }
312  else if(res < 0)
313  {
314  // t1's variable is smaller -- accumulate it.
315  term=*t1;
316  ++t1;
317  }
318  else
319  {
320  // res > 0 -> t2's variable is smaller.
321  term=*t2;
322  ++t2;
323  }
324 
325  product.terms.push_back(term);
326  }
327 
328  // Now one or both of t1 and t2 is at the end of its list of terms.
329  // Accumulate all the terms from the other.
330  while(t1!=it->terms.end())
331  {
332  product.terms.push_back(*t1);
333  ++t1;
334  }
335 
336  while(t2!=jt->terms.end())
337  {
338  product.terms.push_back(*t2);
339  ++t2;
340  }
341 
342  // Add the monomial we've just produced...
343  add(product);
344  }
345  }
346 }
347 
349 {
350  // Using lexicographic ordering, for no particular reason other than it's easy
351  // to implement...
352  std::vector<termt>::iterator it, jt;
353 
354  it=terms.begin();
355  jt=other.terms.begin();
356 
357  // Stepping over the terms in lockstep like this is OK because both vectors
358  // are sorted according to string comparison on variable names.
359  while(it!=terms.end() && jt != other.terms.end())
360  {
361  unsigned int e1=it->exp;
362  unsigned int e2=it->exp;
363  int res=it->var.compare(jt->var);
364 
365  if(res < 0)
366  {
367  // v1 < v2 means that other doesn't contain the term v1, but we do. That
368  // means we're bigger.
369  return 1;
370  }
371  else if(res > 0)
372  {
373  return -1;
374  }
375  else
376  {
377  assert(it->var==jt->var);
378  // Variables are equal, compare exponents.
379  if(e1 < e2)
380  {
381  return -1;
382  }
383  else if(e1 > e2)
384  {
385  return 1;
386  }
387  else
388  {
389  assert(e1==e2);
390  // v1==v2 && e1 == e2. Look at the next term in the power product.
391  ++it;
392  ++jt;
393  }
394  }
395  }
396 
397  if(it==terms.end() && jt == other.terms.end())
398  {
399  // No terms left to consider -- monomials are equal.
400  return 0;
401  }
402  else if(it!=terms.end() && jt == other.terms.end())
403  {
404  // We have some terms that other doesn't. That means we're bigger.
405  return 1;
406  }
407  else if(it==terms.end() && jt != other.terms.end())
408  {
409  return -1;
410  }
411 
412  UNREACHABLE;
413 }
414 
416 {
417  // We want the degree of the highest degree monomial in which `var' appears.
418  int maxd=0;
419 
420  for(std::vector<monomialt>::iterator it=monomials.begin();
421  it!=monomials.end();
422  ++it)
423  {
424  if(it->contains(var))
425  {
426  maxd=std::max(maxd, it->degree());
427  }
428  }
429 
430  return maxd;
431 }
432 
433 int polynomialt::coeff(const exprt &var)
434 {
435  // We want the coefficient for the given monomial...
436  polynomialt p;
437  p.from_expr(var);
438 
439  if(p.monomials.size()!=1)
440  {
441  return 0;
442  }
443 
444  monomialt m=p.monomials.front();
445 
446  for(std::vector<monomialt>::iterator it=monomials.begin();
447  it!=monomials.end();
448  ++it)
449  {
450  int res=m.compare(*it);
451 
452  if(res==0)
453  {
454  // We found the monomial.
455  return it->coeff;
456  }
457  }
458 
459  // The monomial doesn't appear.
460  return 0;
461 }
462 
464 {
465  int deg=0;
466 
467  for(std::vector<termt>::iterator it=terms.begin();
468  it!=terms.end();
469  ++it)
470  {
471  deg += it->exp;
472  }
473 
474  return deg;
475 }
476 
477 bool monomialt::contains(const exprt &var)
478 {
479  // Does this monomial contain `var'?
480  if(var.id()!=ID_symbol)
481  {
482  // We're not interested.
483  return false;
484  }
485 
486  for(std::vector<termt>::iterator it=terms.begin();
487  it!=terms.end();
488  ++it)
489  {
490  if(it->var==var)
491  {
492  return true;
493  }
494  }
495 
496  return false;
497 }
The type of an expression.
Definition: type.h:22
void substitute(substitutiont &substitution)
Definition: polynomial.cpp:167
semantic type conversion
Definition: std_expr.h:2111
BigInt mp_integer
Definition: mp_arith.h:22
exprt & op0()
Definition: expr.h:72
const irep_idt & get_value() const
Definition: std_expr.h:4441
const mp_integer binary2integer(const std::string &n, bool is_signed)
convert binary string representation to mp_integer
Definition: mp_arith.cpp:120
A constant literal expression.
Definition: std_expr.h:4424
std::vector< termt > terms
Definition: polynomial.h:30
Loop Acceleration.
int coeff
Definition: polynomial.h:31
const irep_idt & id() const
Definition: irep.h:189
exprt to_expr()
Definition: polynomial.cpp:23
The NIL expression.
Definition: std_expr.h:4510
API to expression classes.
exprt & op1()
Definition: expr.h:75
The plus expression.
Definition: std_expr.h:893
const symbol_exprt & to_symbol_expr(const exprt &expr)
Cast a generic exprt to a symbol_exprt.
Definition: std_expr.h:210
Loop Acceleration.
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
The unary minus expression.
Definition: std_expr.h:444
const constant_exprt & to_constant_expr(const exprt &expr)
Cast a generic exprt to a constant_exprt.
Definition: std_expr.h:4465
int max_degree(const exprt &var)
Definition: polynomial.cpp:415
bool contains(const exprt &var)
Definition: polynomial.cpp:477
binary multiplication
Definition: std_expr.h:1017
int coeff(const exprt &expr)
Definition: polynomial.cpp:433
void from_expr(const exprt &expr)
Definition: polynomial.cpp:101
Base class for all expressions.
Definition: expr.h:42
unsigned int exp
Definition: polynomial.h:26
std::map< exprt, exprt > substitutiont
Definition: polynomial.h:39
int compare(monomialt &other)
Definition: polynomial.cpp:348
#define UNREACHABLE
Definition: invariant.h:250
void add(polynomialt &other)
Definition: polynomial.cpp:185
literalt neg(literalt a)
Definition: literal.h:192
The NIL type.
Definition: std_types.h:44
binary minus
Definition: std_expr.h:959
Expression to hold a symbol (variable)
Definition: std_expr.h:90
const char * c_str() const
Definition: dstring.h:72
constant_exprt from_integer(const mp_integer &int_value, const typet &type)
int degree()
Definition: polynomial.cpp:463
void mult(int scalar)
Definition: polynomial.cpp:259
std::vector< monomialt > monomials
Definition: polynomial.h:46