symbolic_polynomials/
functions.rs

1use traits::*;
2use monomial::Monomial;
3use polynomial::Polynomial;
4use composite::Composite;
5use ::std::collections::HashMap;
6
7/// Returns a polynomial representing 1 * x^1 + 0,
8/// where 'x' is a variable uniquely identifiable by the provided `id`.
9pub fn variable<I, C, P>(id: I) -> Polynomial<I, C, P>
10    where I: Id, C: Coefficient, P: Power {
11    Polynomial {
12        monomials: vec![Monomial {
13                            coefficient: C::one(),
14                            powers: vec![(Composite::Variable(id), P::one())],
15                        }],
16    }
17}
18
19/// Computes a symbolic `max` between two polynomials.
20pub fn max<I, C, P>(left: &Polynomial<I, C, P>, right: &Polynomial<I, C, P>) -> Polynomial<I, C, P>
21    where I: Id, C: Coefficient, P: Power {
22    if left.is_constant() && right.is_constant() {
23        let v1 = left.eval(&HashMap::default()).ok().unwrap();
24        let v2 = right.eval(&HashMap::default()).ok().unwrap();
25        Polynomial::from(if v1 > v2 { v1 } else { v2 })
26    } else {
27        Polynomial {
28            monomials: vec![Monomial {
29                                coefficient: C::one(),
30                                powers: vec![(Composite::Max(::std::rc::Rc::new(left.clone()),
31                                                             ::std::rc::Rc::new(right.clone())),
32                                              P::one())],
33                            }],
34        }
35    }
36}
37
38/// Computes a symbolic `min` between two polynomials.
39pub fn min<I, C, P>(left: &Polynomial<I, C, P>, right: &Polynomial<I, C, P>) -> Polynomial<I, C, P>
40    where I: Id, C: Coefficient, P: Power {
41    if left.is_constant() && right.is_constant() {
42        let v1 = left.eval(&HashMap::default()).ok().unwrap();
43        let v2 = right.eval(&HashMap::default()).ok().unwrap();
44        Polynomial::from(if v1 < v2 { v1 } else { v2 })
45    } else {
46        Polynomial {
47            monomials: vec![Monomial {
48                                coefficient: C::one(),
49                                powers: vec![(Composite::Min(::std::rc::Rc::new(left.clone()),
50                                                             ::std::rc::Rc::new(right.clone())),
51                                              P::one())],
52                            }],
53        }
54    }
55}
56
57/// Computes a symbolic `ceil` between two polynomials.
58pub fn ceil<I, C, P>(left: &Polynomial<I, C, P>, right: &Polynomial<I, C, P>) -> Polynomial<I, C, P>
59    where I: Id, C: Coefficient, P: Power {
60    if left.is_constant() && right.is_constant() {
61        let v1 = left.eval(&HashMap::default()).ok().unwrap();
62        let v2 = right.eval(&HashMap::default()).ok().unwrap();
63        let (d, rem) = v1.div_rem(&v2);
64        if rem == C::zero() {
65            Polynomial::from(d)
66        } else {
67            Polynomial::from(d + C::one())
68        }
69    } else {
70        let (result, reminder) = left.div_rem(right);
71        if reminder.monomials.is_empty() {
72            result
73        } else {
74            Polynomial {
75                monomials: vec![Monomial {
76                                    coefficient: C::one(),
77                                    powers: vec![(Composite::Ceil(::std::rc::Rc::new(left.clone()),
78                                                                  ::std::rc::Rc::new(right.clone())),
79                                                  P::one())],
80                                }],
81            }
82        }
83    }
84}
85
86/// Computes a symbolic `floor` between two polynomials.
87pub fn floor<I, C, P>(left: &Polynomial<I, C, P>, right: &Polynomial<I, C, P>) -> Polynomial<I, C, P>
88    where I: Id, C: Coefficient, P: Power {
89    if left.is_constant() && right.is_constant() {
90        let v1 = left.eval(&HashMap::default()).ok().unwrap();
91        let v2 = right.eval(&HashMap::default()).ok().unwrap();
92        Polynomial::from(C::div_floor(&v1, &v2))
93    } else {
94        let (result, reminder) = left.div_rem(right);
95        if reminder.monomials.is_empty() {
96            result
97        } else {
98            Polynomial {
99                monomials: vec![Monomial {
100                                    coefficient: C::one(),
101                                    powers: vec![(Composite::Floor(::std::rc::Rc::new(left.clone()),
102                                                                   ::std::rc::Rc::new(right.clone())),
103                                                  P::one())],
104                                }],
105            }
106        }
107    }
108}
109
110/// Reduces the monomial, given the variable assignments provided.
111pub fn reduce_monomial<I, C, P>(monomial: &Monomial<I, C, P>, values: &HashMap<I, C>) -> Monomial<I, C, P>
112    where I: Id, C: Coefficient, P: Power {
113    if monomial.is_constant() {
114        monomial.clone()
115    } else {
116        let mut result = Monomial::<I, C, P> {
117            coefficient: monomial.coefficient.clone(),
118            powers: Vec::new(),
119        };
120        for &(ref c, ref p) in &monomial.powers {
121            match *c {
122                Composite::Variable(ref id) => {
123                    match values.get(id) {
124                        Some(value) => {
125                            result.coefficient *= ::num::pow(value.clone(), p.to_usize().unwrap());
126                        }
127                        None => {
128                            result *= &Monomial::<I, C, P> {
129                                coefficient: C::one(),
130                                powers: vec![(c.clone(), p.clone())],
131                            };
132                        }
133                    }
134                }
135                Composite::Max(ref left, ref right) => {
136                    let mut reduced_left = ::std::rc::Rc::new(reduce(&*left, &*values));
137                    let mut reduced_right = ::std::rc::Rc::new(reduce(&*right, &*values));
138                    if reduced_left.eq(left) {
139                        reduced_left = left.clone();
140                    }
141                    if reduced_right.eq(right) {
142                        reduced_right = right.clone();
143                    }
144                    let c = Composite::Max(reduced_left.clone(), reduced_right.clone());
145                    if reduced_left.is_constant() && reduced_right.is_constant() {
146                        result.coefficient *= c.eval(values).unwrap();
147                    } else {
148                        result *= &Monomial::<I, C, P> {
149                            coefficient: C::one(),
150                            powers: vec![(c, p.clone())],
151                        };
152                    }
153                }
154                Composite::Min(ref left, ref right) => {
155                    let mut reduced_left = ::std::rc::Rc::new(reduce(&*left, &*values));
156                    let mut reduced_right = ::std::rc::Rc::new(reduce(&*right, &*values));
157                    if reduced_left.eq(left) {
158                        reduced_left = left.clone();
159                    }
160                    if reduced_right.eq(right) {
161                        reduced_right = right.clone();
162                    }
163                    let c = Composite::Min(reduced_left.clone(), reduced_right.clone());
164                    if reduced_left.is_constant() && reduced_right.is_constant() {
165                        result.coefficient *= c.eval(values).unwrap();
166                    } else {
167                        result *= &Monomial::<I, C, P> {
168                            coefficient: C::one(),
169                            powers: vec![(c, p.clone())],
170                        };
171                    }
172                }
173                Composite::Ceil(ref left, ref right) => {
174                    let mut reduced_left = ::std::rc::Rc::new(reduce(&*left, &*values));
175                    let mut reduced_right = ::std::rc::Rc::new(reduce(&*right, &*values));
176                    if reduced_left.eq(left) {
177                        reduced_left = left.clone();
178                    }
179                    if reduced_right.eq(right) {
180                        reduced_right = right.clone();
181                    }
182                    let c = Composite::Ceil(reduced_left.clone(), reduced_right.clone());
183                    if reduced_left.is_constant() && reduced_right.is_constant() {
184                        result.coefficient *= c.eval(values).unwrap();
185                    } else {
186                        result *= &Monomial::<I, C, P> {
187                            coefficient: C::one(),
188                            powers: vec![(c, p.clone())],
189                        };
190                    }
191                }
192                Composite::Floor(ref left, ref right) => {
193                    let mut reduced_left = ::std::rc::Rc::new(reduce(&*left, &*values));
194                    let mut reduced_right = ::std::rc::Rc::new(reduce(&*right, &*values));
195                    if reduced_left.eq(left) {
196                        reduced_left = left.clone();
197                    }
198                    if reduced_right.eq(right) {
199                        reduced_right = right.clone();
200                    }
201                    let c = Composite::Floor(reduced_left.clone(), reduced_right.clone());
202                    if reduced_left.is_constant() && reduced_right.is_constant() {
203                        result.coefficient *= c.eval(values).unwrap();
204                    } else {
205                        result *= &Monomial::<I, C, P> {
206                            coefficient: C::one(),
207                            powers: vec![(c, p.clone())],
208                        };
209                    }
210                }
211            }
212        }
213        result
214    }
215}
216
217/// Reduces the polynomial, given the variable assignments provided.
218pub fn reduce<I, C, P>(polynomial: &Polynomial<I, C, P>, values: &HashMap<I, C>) -> Polynomial<I, C, P>
219    where I: Id, C: Coefficient, P: Power {
220    let mut result = Polynomial::<I, C, P> { monomials: Vec::new() };
221    for m in &polynomial.monomials {
222        result += &reduce_monomial(m, values);
223    }
224    result
225}
226
227/// Automatically deduces the individual variable assignments based on the
228/// system of equations specified by the mapping of `Polynomial` to a constant value.
229pub fn deduce_values<I, C, P>(original_values: &Vec<(Polynomial<I, C, P>, C)>) -> Result<HashMap<I, C>, String>
230    where I: Id, C: Coefficient, P: Power {
231    let mut implicit_values = original_values.clone();
232    let mut indexes: Vec<usize> = (0..original_values.len()).collect();
233    let mut values: HashMap<I, C> = HashMap::new();
234    let mut i = 0;
235    while i < implicit_values.len() {
236        let (to_remove, to_reduce) = {
237            let (ref p, ref c) = implicit_values[i];
238            if p.is_constant() {
239                let value = p.eval(&HashMap::new()).unwrap();
240                if value != *c {
241                    return Err(format!("Value deduction failed for {} = {}, as it was deduced to {}.",
242                                       original_values[indexes[i]].0,
243                                       c,
244                                       value));
245                } else {
246                    (true, false)
247                }
248            } else if (p.monomials.len() == 1 || (p.monomials.len() == 2 && p.monomials[1].is_constant())) &&
249                      p.monomials[0].powers.len() == 1 {
250                if let Composite::Variable(ref id) = p.monomials[0].powers[0].0 {
251                    // The polynomial is in the form a * x^n + b, so we can deduce value of 'x'
252                    let b = p.monomials.get(1).map_or(C::zero(), |m| m.eval(&values).unwrap());
253                    let a = &p.monomials[0].coefficient;
254                    let n = &p.monomials[0].powers[0].1;
255                    let inferred = match nth_root((c.clone() - b.clone()) / a.clone(), n.clone()) {
256                        Some(val) => val,
257                        None => {
258                            return Err(format!("Could not find integer solution to {} * {}^{} + {} = {}.",
259                                               a,
260                                               id,
261                                               n,
262                                               b,
263                                               c))
264                        }
265                    };
266                    values.insert(id.clone(), inferred);
267                    (true, true)
268                } else {
269                    (false, false)
270                }
271            } else {
272                (false, false)
273            }
274        };
275        if to_remove {
276            implicit_values.remove(i);
277            indexes.remove(i);
278            if implicit_values.is_empty() {
279                return Ok(values);
280            }
281        }
282        if to_reduce {
283            for &mut (ref mut p, _) in &mut implicit_values {
284                *p = reduce(p, &values);
285            }
286            i = 0;
287        } else if !to_remove {
288            i += 1;
289        }
290    }
291    Err("Could not deduce all variables.".into())
292}
293
294fn nth_root<C, P>(value: C, n: P) -> Option<C>
295    where C: Coefficient, P: Power {
296    let result = if value < C::zero() {
297        C::from_f64(-(-value.to_f64().unwrap()).powf(n.to_f64().unwrap().recip())).unwrap()
298    } else {
299        C::from_f64(value.to_f64().unwrap().powf(n.to_f64().unwrap().recip())).unwrap()
300    };
301    if ::num::pow(result.clone(), n.to_usize().unwrap()) == value {
302        Some(result)
303    } else {
304        None
305    }
306}