1 /*
2  * This file is part of the Yices SMT Solver.
3  * Copyright (C) 2017 SRI International.
4  *
5  * Yices is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * Yices is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with Yices.  If not, see <http://www.gnu.org/licenses/>.
17  */
18 
19 /*
20  * Power products: x_1^d_1 ... x_n^d_n
21  * - each x_i is a 32bit integer
22  * - each exponent d_i is a positive integer
23  *
24  * There's a limit on the total degree (YICES_MAX_DEGREE) defined in 'yices_limits.h'.
25  * If power products p1 and p2 have degree less than MAX_DEGREE
26  * then degree(p1) + degree(p2) can be computed without overflow on 32bit numbers.
27  * Power products (and polynomials) of degree >= YICES_MAX_DEGREE are not supported.
28  */
29 
30 #ifndef __POWER_PRODUCTS_H
31 #define __POWER_PRODUCTS_H
32 
33 #include <stdint.h>
34 #include <stdbool.h>
35 #include <stddef.h>
36 
37 #include "utils/tagged_pointers.h"
38 
39 /*
40  * Product:
41  * - len = number of components
42  * - prod = array of components, each component is a pair <var, exp>
43  * - degree = sum of the exponents
44  * - each variable is represented as a 32bit signed integer
45  * When normalized:
46  * - the variables are in increasing order
47  * - there are no duplicate variables; each var occurs in a single pair
48  * - the exponents are non-zero
49  */
50 // pair <variable, exponent>
51 typedef struct {
52   int32_t var;
53   uint32_t exp;
54 } varexp_t;
55 
56 // power product
57 typedef struct {
58   uint32_t len;
59   uint32_t degree;
60   varexp_t prod[0]; // real size is len
61 } pprod_t;
62 
63 
64 /*
65  * Buffer for intermediate computations.
66  */
67 typedef struct {
68   uint32_t size;    // size of the array prod
69   uint32_t len;     // elements of prod currently used
70   varexp_t *prod;
71 } pp_buffer_t;
72 
73 
74 /*
75  * Maximal number of pairs <variable, exponent> in a product
76  * - this is to ensure we can compute len * sizeof(varexp_t)
77  *   on 32bits without overflow
78  */
79 #define PPROD_MAX_LENGTH ((UINT32_MAX-sizeof(pprod_t))/sizeof(varexp_t))
80 
81 
82 /*
83  * Maximal variable index: variables should be between 0 and MAX_PROD_VAR
84  */
85 #define MAX_PPROD_VAR (INT32_MAX/2)
86 
87 
88 
89 
90 
91 /*
92  * POINTER TAGGING
93  */
94 
95 /*
96  * Power product are part of the polynomial representation. But in many
97  * cases, the polynomials are linear expressions (i.e., each power product
98  * is either the empty product or a single variable with exponent 1). To
99  * compactly encode the common case, we use tagged pointers and the
100  * following conventions.
101  * - the empty power product is the NULL pointer
102  * - a single variable x is packed in a pprod_t pointer with tag bit set to 1
103  *   (this gives a compact representation for x^1)
104  * - otherwise, the product p is a pointer to an actual pprod_t object
105  *
106  * We also use a special end-marker in polynomial representations.
107  * The end marker is larger than any other power product in the deglex ordering.
108  */
109 
110 /*
111  * Empty product
112  */
113 #define empty_pp ((pprod_t *) NULL)
114 
115 /*
116  * End marker: all bits are one
117  */
118 #define end_pp ((pprod_t *) ~((uintptr_t) 0))
119 
120 /*
121  * Variable x encoded as a power product
122  */
var_pp(int32_t x)123 static inline pprod_t *var_pp(int32_t x) {
124   assert(0 <= x && x <= MAX_PPROD_VAR);
125   return (pprod_t *) tag_i32(x);
126 }
127 
128 /*
129  * Check whether p is empty or a variable
130 ' */
pp_is_empty(pprod_t * p)131 static inline bool pp_is_empty(pprod_t *p) {
132   return p == empty_pp;
133 }
134 
pp_is_var(pprod_t * p)135 static inline bool pp_is_var(pprod_t *p) {
136   return has_int_tag(p);
137 }
138 
139 /*
140  * Get the variable index of p if pp_is_var(p) holds
141  */
var_of_pp(pprod_t * p)142 static inline int32_t var_of_pp(pprod_t *p) {
143   assert(pp_is_var(p));
144   return untag_i32(p);
145 }
146 
147 
148 
149 
150 /*
151  * PRODUCT BUFFERS
152  */
153 
154 /*
155  * Initialize a buffer of initial size n:
156  * - initial value = empty product
157  */
158 extern void init_pp_buffer(pp_buffer_t *b, uint32_t n);
159 extern void delete_pp_buffer(pp_buffer_t *b);
160 
161 /*
162  * Operations:
163  * - reset to empty product
164  * - set to or multiply by v
165  * - set to or multiply v^d
166  * - set to or multiply by v[0] * ... * v[n-1]
167  * - set to or multiply by v[0]^d[0] * ... * v[n-1]^d[n-1]
168  * all modify buffer b.
169  *
170  * The set_xxx operations do not normalize b.
171  * The mul_xxx operations do normalize b.
172  */
173 extern void pp_buffer_reset(pp_buffer_t *b);
174 
175 extern void pp_buffer_set_var(pp_buffer_t *b, int32_t v);
176 extern void pp_buffer_set_varexp(pp_buffer_t *b, int32_t v, uint32_t d);
177 extern void pp_buffer_set_vars(pp_buffer_t *b, uint32_t n, int32_t *v);
178 extern void pp_buffer_set_varexps(pp_buffer_t *b, uint32_t n, int32_t *v, uint32_t *d);
179 
180 extern void pp_buffer_mul_var(pp_buffer_t *b, int32_t v);
181 extern void pp_buffer_mul_varexp(pp_buffer_t *b, int32_t v, uint32_t d);
182 extern void pp_buffer_mul_vars(pp_buffer_t *b, uint32_t n, int32_t *v);
183 extern void pp_buffer_mul_varexps(pp_buffer_t *b, uint32_t n, int32_t *v, uint32_t *d);
184 
185 
186 /*
187  * Lower-level operation: add v or v^d as a factor of b
188  * - this does not normalize b
189  */
190 extern void pp_buffer_push_var(pp_buffer_t *b, int32_t v);
191 extern void pp_buffer_push_varexp(pp_buffer_t *n, int32_t v, uint32_t d);
192 
193 
194 /*
195  * Raise b to power d and normalize the result
196  */
197 extern void pp_buffer_exponentiate(pp_buffer_t *b, uint32_t d);
198 
199 
200 /*
201  * Assign or multiply b by a power-product p
202  * - p must follow the pointer tagging conventions
203  */
204 extern void pp_buffer_set_pprod(pp_buffer_t *b, pprod_t *p);
205 extern void pp_buffer_mul_pprod(pp_buffer_t *b, pprod_t *p);
206 
207 /*
208  * Copy src into b or multiply b by a
209  */
210 extern void pp_buffer_copy(pp_buffer_t *b, pp_buffer_t *src);
211 extern void pp_buffer_mul_buffer(pp_buffer_t *b, pp_buffer_t *a);
212 
213 
214 /*
215  * Normalize: remove x^0 and replace x^n * x^m by x^(n+m)
216  * also sort in increasing variable order
217  */
218 extern void pp_buffer_normalize(pp_buffer_t *b);
219 
220 /*
221  * Check whether b contains a trivial product
222  * (either empty or a single variable with exponent 1)
223  * - b must be normalized
224  */
pp_buffer_is_trivial(pp_buffer_t * b)225 static inline bool pp_buffer_is_trivial(pp_buffer_t *b) {
226   return (b->len == 0) || (b->len == 1 && b->prod[0].exp == 1);
227 }
228 
229 /*
230  * Degree computation
231  */
232 extern uint32_t pp_buffer_degree(pp_buffer_t *b);
233 
234 /*
235  * Check whether the degree is less than YICES_MAX_DEGREE.
236  */
237 extern bool pp_buffer_below_max_degree(pp_buffer_t *b);
238 
239 /*
240  * Degree of a variable x in product p
241  * (0 if x does not occur in p).
242  */
243 extern uint32_t pp_buffer_var_degree(pp_buffer_t *b, int32_t x);
244 
245 /*
246  * Check whether b1 and b2 are equal
247  * - both must be normalized
248  */
249 extern bool pp_buffer_equal(pp_buffer_t *b1, pp_buffer_t *b2);
250 
251 /*
252  * Convert b's content to a power-product object.
253  * - b must be normalized.
254  * - if b is empty or is of the form x^1 then the result is
255  *   the appropriate tagged pointer.
256  */
257 extern pprod_t *pp_buffer_getprod(pp_buffer_t *b);
258 
259 
260 /*
261  * Compute the common factors of b and b1. Store the result in b.
262  * - on entry: b and b1 must be normalized
263  * - on exit: b is modified to store the common factors and it's normalized
264  */
265 extern void pp_buffer_gcd(pp_buffer_t *b, pp_buffer_t *b1);
266 
267 /*
268  * Divide b by b1. Store the result in b.
269  * - on entry b and b1 must be normalized
270  * - on exit: b stores the division.
271  * - this should be used only if b1 is a divisor of b but we don't check.
272  */
273 extern void pp_buffer_divide(pp_buffer_t *b, pp_buffer_t *b1);
274 
275 /*
276  * Divide b by x. Store the result in b.
277  * - on entry, b must be normalized.
278  * - on exit, b store the division and is normalized.
279  */
280 extern void pp_buffer_divide_by_var(pp_buffer_t *b, int32_t x);
281 
282 
283 /*
284  * POWER PRODUCT OPERATIONS
285  */
286 
287 /*
288  * All the operations assume p follows the tagged pointer
289  * conventions.
290  * - p = NULL denotes the empty product
291  * - a tagged pointer x denotes the product (x^1)
292  * - otherwise, p is a pointer to an actual pprod_t structure
293  *   that must be normalized
294  */
295 
296 /*
297  * Degree of p
298  * - p must be distinct from end_pp
299  */
300 extern uint32_t pprod_degree(pprod_t *p);
301 
302 
303 /*
304  * Check whether p's degree is below MAX_DEGREE
305  * - p must be distinct from end_pp
306  */
307 extern bool pprod_below_max_degree(pprod_t *p);
308 
309 
310 /*
311  * Degree of variable x in product p
312  * - returns 0 if x does not occur in p
313  * - p must be distinct from end_pp
314  */
315 extern uint32_t pprod_var_degree(pprod_t *p, int32_t x);
316 
317 
318 /*
319  * Lexicographic ordering:
320  * - lex_cmp(p1, p2) < 0 means p1 < p2
321  * - lex_cmp(p1, p2) = 0 means p1 = p2
322  * - lex_cmp(p1, p2) > 0 means p2 < p1
323  *
324  * The order is defined by
325  *   x_1^d_1 ... x_n^d_n < x_1^c_1 ... x_n^c_n
326  * if, for some index i, we have
327  *   d_1 = c_1 ... d_{i-1} = c_{i-1} and d_i > c_i.
328  *
329  * p1 and p2 must be distinct from end_pp
330  */
331 extern int32_t pprod_lex_cmp(pprod_t *p1, pprod_t *p2);
332 
333 
334 /*
335  * Degree then lex ordering. This ordering is defined by
336  *  p1 < p2 if degree(p1) < degree(p2)
337  *          or degree(p1) = degree(p2) and (p1 < p2 in lex ordering)
338  *
339  * That's the ordering we use for normalizing polynomials.
340  * p1 or p2 (or both) may be equal to the end marker end_pp.
341  *
342  * There are three cases:
343  * - if p1 = end_pp then
344  *     pprod_precedes(p1, p2) is false for any p2
345  * - if p1 != end_pp and p2 = end_pp then
346  *     pprod_precedes(p1, p2) is true
347  * - otherwise,
348  *     pprod_precedes(p1, p2) is true
349  *        iff p1 < p2 in the deglex ordering
350  */
351 extern bool pprod_precedes(pprod_t *p1, pprod_t *p2);
352 
353 
354 /*
355  * Check whether p1 and p2 are equal
356  * - p1 and p2 must be distinct from end_pp
357  */
358 extern bool pprod_equal(pprod_t *p1, pprod_t *p2);
359 
360 
361 /*
362  * Divisibility test: returns true if p1 is a divisor of p2
363  * - p1 and p2 must be distinct from end_pp
364  */
365 extern bool pprod_divides(pprod_t *p1, pprod_t *p2);
366 
367 
368 /*
369  * Same thing but also computes the quotient (p2/p1) in b
370  * - p1 and p2 must be distinct from end_pp
371  */
372 extern bool pprod_divisor(pp_buffer_t *b, pprod_t *p1, pprod_t *p2);
373 
374 
375 /*
376  * Free a power-product
377  * - p must be the result of pp_buffer_getprod or make_pprod.
378  * - do nothing if p is empty or a variable.
379  */
380 extern void free_pprod(pprod_t *p);
381 
382 
383 
384 
385 
386 /*
387  * LOW-LEVEL OPERATIONS
388  */
389 
390 /*
391  * Check whether a and b (both of length n) are equal
392  * - both a and b must be normalized
393  */
394 extern bool varexp_array_equal(varexp_t *a, varexp_t *b, uint32_t n);
395 
396 /*
397  * Build a pproduct by making a copy of a
398  * - a must be normalized
399  * - n must be the length of a
400  * - n must be less than PPROD_MAX_LENGTH
401  */
402 extern pprod_t *make_pprod(varexp_t *a, uint32_t n);
403 
404 
405 
406 
407 #endif /* __POWER_PRODUCTS_H */
408