1 /**
2 * Copyright 2015, SRI International.
3 *
4 * This file is part of LibPoly.
5 *
6 * LibPoly is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * LibPoly is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with LibPoly. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 #include "polypyPolynomial.h"
21
22 #include "polypyInteger.h"
23 #include "polypyVariable.h"
24 #include "polypyVariableOrder.h"
25 #include "polypyAssignment.h"
26 #include "polypyValue.h"
27 #include "polypyInterval.h"
28 #include "polypyFeasibilitySet.h"
29
30 #include "utils.h"
31 #include "variable_list.h"
32 #include "sign_condition.h"
33 #include "feasibility_set.h"
34 #include "value.h"
35 #include "polynomial_vector.h"
36
37 #include <structmember.h>
38
39 static lp_polynomial_context_t* default_ctx = 0;
40
Polynomial_get_default_context(void)41 const lp_polynomial_context_t* Polynomial_get_default_context(void) {
42 if (!default_ctx) {
43 default_ctx = lp_polynomial_context_new(0, Variable_get_default_db(), (lp_variable_order_t*) VariableOrder_get_default_order());
44 }
45 return default_ctx;
46 }
47
48 static void
49 Polynomial_dealloc(Polynomial* self);
50
51 static PyObject*
52 Polynomial_new(PyTypeObject *type, PyObject *args, PyObject *kwds);
53
54 static int
55 Polynomial_cmp(PyObject* self, PyObject* args);
56
57 static PyObject*
58 Polynomial_richcompare(PyObject* self, PyObject* args, int op);
59
60 static long
61 Polynomial_hash(PyObject* self);
62
63 static PyObject*
64 Polynomial_degree(PyObject* self);
65
66 static PyObject*
67 Polynomial_coefficients(PyObject* self);
68
69 static PyObject*
70 Polynomial_reductum(PyObject* self, PyObject* args);
71
72 static PyObject*
73 Polynomial_sgn(PyObject* self, PyObject* arguments);
74
75 static PyObject*
76 Polynomial_sgn_check(PyObject* self, PyObject* args);
77
78 static PyObject*
79 Polynomial_rem(PyObject* self, PyObject* args);
80
81 static PyObject*
82 Polynomial_prem(PyObject* self, PyObject* args);
83
84 static PyObject*
85 Polynomial_sprem(PyObject* self, PyObject* args);
86
87 static PyObject*
88 Polynomial_gcd(PyObject* self, PyObject* args);
89
90 static PyObject*
91 Polynomial_lcm(PyObject* self, PyObject* args);
92
93 static PyObject*
94 Polynomial_derivative(PyObject* self);
95
96 static PyObject*
97 Polynomial_extended_gcd(PyObject* self, PyObject* args);
98
99 static PyObject*
100 Polynomial_factor(PyObject* self);
101
102 static PyObject*
103 Polynomial_factor_square_free(PyObject* self);
104
105 static PyObject*
106 Polynomial_roots_count(PyObject* self, PyObject* args);
107
108 static PyObject*
109 Polynomial_roots_isolate(PyObject* self, PyObject* args);
110
111 static PyObject*
112 Polynomial_sturm_sequence(PyObject* self);
113
114 static PyObject*
115 Polynomial_str(PyObject* self);
116
117 static int
118 Polynomial_nonzero(PyObject* self);
119
120 static PyObject*
121 Polynomial_add(PyObject* self, PyObject* args);
122
123 static PyObject*
124 Polynomial_neg(PyObject* self);
125
126 static PyObject*
127 Polynomial_sub(PyObject* self, PyObject* args);
128
129 static PyObject*
130 Polynomial_mul(PyObject* self, PyObject* args);
131
132 static PyObject*
133 Polynomial_div(PyObject* self, PyObject* args);
134
135 static PyObject*
136 Polynomial_rem_operator(PyObject* self, PyObject* args);
137
138 static PyObject*
139 Polynomial_divmod(PyObject* self, PyObject* args);
140
141 static PyObject*
142 Polynomial_pow(PyObject* self, PyObject* args);
143
144 static PyObject*
145 Polynomial_resultant(PyObject* self, PyObject* args);
146
147 static PyObject*
148 Polynomial_psc(PyObject* self, PyObject* args);
149
150 static PyObject*
151 Polynomial_mgcd(PyObject* self, PyObject* args);
152
153 static PyObject*
154 Polynomial_evaluate(PyObject* self, PyObject* args);
155
156 static PyObject*
157 Polynomial_vars(PyObject* self);
158
159 static PyObject*
160 Polynomial_var(PyObject* self);
161
162 static PyObject*
163 Polynomial_pp(PyObject* self);
164
165 static PyObject*
166 Polynomial_cont(PyObject* self);
167
168 static PyObject*
169 Polynomial_pp_cont(PyObject* self);
170
171 static PyObject*
172 Polynomial_feasible_intervals(PyObject* self, PyObject* args);
173
174 static PyObject*
175 Polynomial_feasible_set(PyObject* self, PyObject* args);
176
177 PyMethodDef Polynomial_methods[] = {
178 {"degree", (PyCFunction)Polynomial_degree, METH_NOARGS, "Returns the degree of the polynomial in its top variable"},
179 {"coefficients", (PyCFunction)Polynomial_coefficients, METH_NOARGS, "Returns a dictionary from degrees to coefficients"},
180 {"reductum", (PyCFunction)Polynomial_reductum, METH_VARARGS, "Returns the reductum of the polynomial"},
181 {"sgn", (PyCFunction)Polynomial_sgn, METH_VARARGS, "Returns the sign of the polynomials in the given model"},
182 {"sgn_check", (PyCFunction)Polynomial_sgn_check, METH_VARARGS, "Returns true if the sign of the polynomail respects the sign condition."},
183 {"rem", (PyCFunction)Polynomial_rem, METH_VARARGS, "Returns the remainder of current and given polynomial"},
184 {"prem", (PyCFunction)Polynomial_prem, METH_VARARGS, "Returns the pseudo remainder of current and given polynomial"},
185 {"sprem", (PyCFunction)Polynomial_sprem, METH_VARARGS, "Returns the sparse pseudo remainder of current and given polynomial"},
186 {"gcd", (PyCFunction)Polynomial_gcd, METH_VARARGS, "Returns the gcd of current and given polynomial"},
187 {"lcm", (PyCFunction)Polynomial_lcm, METH_VARARGS, "Returns the lcm of current and given polynomial"},
188 {"extended_gcd", (PyCFunction)Polynomial_extended_gcd, METH_VARARGS, "Returns the extended gcd, i.e. (gcd, u, v), of current and given polynomial"},
189 {"factor", (PyCFunction)Polynomial_factor, METH_NOARGS, "Returns the factorization of the polynomial"},
190 {"factor_square_free", (PyCFunction)Polynomial_factor_square_free, METH_NOARGS, "Returns the square-free factorization of the polynomial"},
191 {"roots_count", (PyCFunction)Polynomial_roots_count, METH_VARARGS, "Returns the number of real roots in the given interval"},
192 {"roots_isolate", (PyCFunction)Polynomial_roots_isolate, METH_VARARGS, "Returns the list of real roots (has to be univariate modulo the assignment)"},
193 {"sturm_sequence", (PyCFunction)Polynomial_sturm_sequence, METH_NOARGS, "Returns the Sturm sequence"},
194 {"derivative", (PyCFunction)Polynomial_derivative, METH_NOARGS, "Returns the derivative of the polynomial"},
195 {"resultant", (PyCFunction)Polynomial_resultant, METH_VARARGS, "Returns the resultant of the current and given polynomial"},
196 {"psc", (PyCFunction)Polynomial_psc, METH_VARARGS, "Returns the principal subresultant coefficients of the current and given polynomial"},
197 {"mgcd", (PyCFunction)Polynomial_mgcd, METH_VARARGS, "Returns assumptions that the GCD of two polynomials is of same degree"},
198 {"factor_square_free", (PyCFunction)Polynomial_factor_square_free, METH_NOARGS, "Returns the square-free factorization of the polynomial"},
199 {"evaluate", (PyCFunction)Polynomial_evaluate, METH_VARARGS, "Returns the value of the polynomial in the given assignment (or null if it doesn't fully evaluate"},
200 {"vars", (PyCFunction)Polynomial_vars, METH_NOARGS, "Returns the list of variables in the polynomial"},
201 {"var", (PyCFunction)Polynomial_var, METH_NOARGS, "Returns the top variable of the polynomial"},
202 {"pp", (PyCFunction)Polynomial_pp, METH_NOARGS, "Returns the primitive part of the polynomial"},
203 {"cont", (PyCFunction)Polynomial_cont, METH_NOARGS, "Returns the content of the polynomial"},
204 {"pp_cont", (PyCFunction)Polynomial_pp_cont, METH_NOARGS, "Returns the tuple (pp, cont) of the polynomial"},
205 {"feasible_intervals", (PyCFunction)Polynomial_feasible_intervals, METH_VARARGS, "Returns feasible intervals (list) of the polynomial (has to be univariate modulo the assignment)"},
206 {"feasible_set", (PyCFunction)Polynomial_feasible_set, METH_VARARGS, "Returns feasible set of the polynomial (has to be univariate modulo the assignment)"},
207 {NULL} /* Sentinel */
208 };
209
210 PyNumberMethods Polynomial_NumberMethods = {
211 Polynomial_add, // binaryfunc nb_add;
212 Polynomial_sub, // binaryfunc nb_subtract;
213 Polynomial_mul, // binaryfunc nb_multiply;
214 Polynomial_div, // binaryfunc nb_divide;
215 Polynomial_rem_operator, // binaryfunc nb_remainder;
216 Polynomial_divmod, // binaryfunc nb_divmod;
217 (ternaryfunc)Polynomial_pow, // ternaryfunc nb_power;
218 Polynomial_neg, // unaryfunc nb_negative;
219 0, // unaryfunc nb_positive;
220 0, // unaryfunc nb_absolute;
221 Polynomial_nonzero, // inquiry nb_nonzero; /* Used by PyObject_IsTrue */
222 0, // unaryfunc nb_invert;
223 0, // binaryfunc nb_lshift;
224 0, // binaryfunc nb_rshift;
225 0, // binaryfunc nb_and;
226 0, // binaryfunc nb_xor;
227 0, // binaryfunc nb_or;
228 0, // coercion nb_coerce; /* Used by the coerce() function */
229 0, // unaryfunc nb_int;
230 0, // unaryfunc nb_long;
231 0, // unaryfunc nb_float;
232 0, // unaryfunc nb_oct;
233 0, // unaryfunc nb_hex;
234
235 /* Added in release 2.0 */
236 0, // binaryfunc nb_inplace_add;
237 0, // binaryfunc nb_inplace_subtract;
238 0, // binaryfunc nb_inplace_multiply;
239 0, // binaryfunc nb_inplace_divide;
240 0, // binaryfunc nb_inplace_remainder;
241 0, // ternaryfunc nb_inplace_power;
242 0, // binaryfunc nb_inplace_lshift;
243 0, // binaryfunc nb_inplace_rshift;
244 0, // binaryfunc nb_inplace_and;
245 0, // binaryfunc nb_inplace_xor;
246 0, // binaryfunc nb_inplace_or;
247
248 /* Added in release 2.2 */
249 0, // binaryfunc nb_floor_divide;
250 0, // binaryfunc nb_true_divide;
251 0, // binaryfunc nb_inplace_floor_divide;
252 0, // binaryfunc nb_inplace_true_divide;
253
254 /* Added in release 2.5 */
255 0 // unaryfunc nb_index;
256 };
257
258 PyTypeObject PolynomialType = {
259 PyObject_HEAD_INIT(NULL)
260 0, /*ob_size*/
261 "polypy.Polynomial", /*tp_name*/
262 sizeof(Polynomial), /*tp_basicsize*/
263 0, /*tp_itemsize*/
264 (destructor)Polynomial_dealloc, /*tp_dealloc*/
265 0, /*tp_print*/
266 0, /*tp_getattr*/
267 0, /*tp_setattr*/
268 Polynomial_cmp, /*tp_compare*/
269 Polynomial_str, /*tp_repr*/
270 &Polynomial_NumberMethods, /*tp_as_number*/
271 0, /*tp_as_sequence*/
272 0, /*tp_as_mapping*/
273 Polynomial_hash, /*tp_hash */
274 0, /*tp_call*/
275 Polynomial_str, /*tp_str*/
276 0, /*tp_getattro*/
277 0, /*tp_setattro*/
278 0, /*tp_as_buffer*/
279 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_CHECKTYPES, /*tp_flags*/
280 "Polynomial objects", /* tp_doc */
281 0, /* tp_traverse */
282 0, /* tp_clear */
283 Polynomial_richcompare, /* tp_richcompare */
284 0, /* tp_weaklistoffset */
285 0, /* tp_iter */
286 0, /* tp_iternext */
287 Polynomial_methods, /* tp_methods */
288 0, /* tp_members */
289 0, /* tp_getset */
290 0, /* tp_base */
291 0, /* tp_dict */
292 0, /* tp_descr_get */
293 0, /* tp_descr_set */
294 0, /* tp_dictoffset */
295 0, /* tp_init */
296 0, /* tp_alloc */
297 Polynomial_new, /* tp_new */
298 };
299
300 static void
Polynomial_dealloc(Polynomial * self)301 Polynomial_dealloc(Polynomial* self)
302 {
303 if (self->p) {
304 lp_polynomial_destruct(self->p);
305 free(self->p);
306 }
307 self->ob_type->tp_free((PyObject*)self);
308 }
309
310 PyObject*
Polynomial_create(lp_polynomial_t * p)311 Polynomial_create(lp_polynomial_t* p) {
312 Polynomial *self;
313 self = (Polynomial*)PolynomialType.tp_alloc(&PolynomialType, 0);
314 lp_polynomial_set_external(p);
315 self->p = p;
316 return (PyObject*) self;
317 }
318
319 static PyObject*
Polynomial_new(PyTypeObject * type,PyObject * args,PyObject * kwds)320 Polynomial_new(PyTypeObject *type, PyObject *args, PyObject *kwds) {
321 return Polynomial_create(0);
322 }
323
324 static PyObject*
Polynomial_richcompare(PyObject * self,PyObject * other,int op)325 Polynomial_richcompare(PyObject* self, PyObject* other, int op) {
326 PyObject *result = 0;
327
328 const lp_polynomial_context_t* ctx = 0;
329
330 // One of them is a polynomial
331 if (PyPolynomial_CHECK(self)) {
332 ctx = lp_polynomial_get_context(((Polynomial*) self)->p);
333 } else {
334 ctx = lp_polynomial_get_context(((Polynomial*) other)->p);
335 }
336
337 int dec_other = 0;
338 int dec_self = 0;
339
340 // Check arguments
341 if (!PyPolynomial_CHECK(self)) {
342 if (PyVariable_CHECK(self)) {
343 self = PyPolynomial_FromVariable(self, ctx);
344 dec_self = 1;
345 } else if (PyLong_or_Int_Check(self)) {
346 self = PyPolynomial_FromLong_or_Int(self, ctx);
347 dec_self = 1;
348 } else {
349 Py_INCREF(Py_NotImplemented);
350 return Py_NotImplemented;
351 }
352 }
353
354 // Check arguments
355 if (!PyPolynomial_CHECK(other)) {
356 if (PyVariable_CHECK(other)) {
357 other = PyPolynomial_FromVariable(other, ctx);
358 dec_other = 1;
359 } else if (PyLong_or_Int_Check(other)) {
360 other = PyPolynomial_FromLong_or_Int(other, ctx);
361 dec_other = 1;
362 } else {
363 Py_INCREF(Py_NotImplemented);
364 return Py_NotImplemented;
365 }
366 }
367
368 lp_polynomial_t* self_p = ((Polynomial*) self)->p;
369 lp_polynomial_t* other_p = ((Polynomial*) other)->p;
370
371 int cmp = lp_polynomial_cmp(self_p, other_p);
372
373 switch (op) {
374 case Py_LT:
375 result = cmp < 0 ? Py_True : Py_False;
376 break;
377 case Py_LE:
378 result = cmp <= 0 ? Py_True : Py_False;
379 break;
380 case Py_EQ:
381 result = cmp == 0 ? Py_True : Py_False;
382 break;
383 case Py_NE:
384 result = cmp != 0 ? Py_True : Py_False;
385 break;
386 case Py_GT:
387 result = cmp > 0 ? Py_True : Py_False;
388 break;
389 case Py_GE:
390 result = cmp >= 0 ? Py_True : Py_False;
391 break;
392 }
393
394 if (dec_self) {
395 Py_DECREF(self);
396 }
397 if (dec_other) {
398 Py_DECREF(other);
399 }
400
401 Py_INCREF(result);
402 return result;
403 }
404
405 static int
Polynomial_cmp(PyObject * self,PyObject * other)406 Polynomial_cmp(PyObject* self, PyObject* other) {
407
408 // Check arguments
409 if (!PyPolynomial_CHECK(self) || !PyPolynomial_CHECK(other)) {
410 // should return -1 and set an exception condition when an error occurred
411 return -1;
412 }
413 // Get arguments
414 Polynomial* p1 = (Polynomial*) self;
415 Polynomial* p2 = (Polynomial*) other;
416 // Compare
417 int cmp = lp_polynomial_cmp(p1->p, p2->p);
418 return cmp > 0 ? 1 : cmp < 0 ? -1 : 0;
419 }
420
421 static long
Polynomial_hash(PyObject * self)422 Polynomial_hash(PyObject* self) {
423 Polynomial* p = (Polynomial*) self;
424 long hash = lp_polynomial_hash(p->p);
425 if (hash == -1) {
426 // value -1 should not be returned as a normal return value
427 hash = 0;
428 }
429 return hash;
430 }
431
Polynomial_str(PyObject * self)432 static PyObject* Polynomial_str(PyObject* self) {
433 Polynomial* p = (Polynomial*) self;
434 if (p) {
435 char* p_str = lp_polynomial_to_string(p->p);
436 PyObject* str = PyString_FromString(p_str);
437 free(p_str);
438 return str;
439 } else {
440 Py_RETURN_NONE;
441 }
442 }
443
444 PyObject*
PyPolynomial_FromVariable(PyObject * variable,const lp_polynomial_context_t * ctx)445 PyPolynomial_FromVariable(PyObject* variable, const lp_polynomial_context_t* ctx) {
446 // The variable
447 lp_variable_t x = ((Variable*) variable)->x;
448
449 // The constant
450 lp_integer_t one;
451 lp_integer_construct_from_int(ctx->K, &one, 1);
452
453 // The x polynomial
454 lp_polynomial_t* p_x = lp_polynomial_alloc();
455 lp_polynomial_construct_simple(p_x, ctx, &one, x, 1);
456
457 // Remove temps
458 lp_integer_destruct(&one);
459
460 // Return the polynomial
461 PyObject* result = Polynomial_create(p_x);
462 Py_INCREF(result);
463 return result;
464 }
465
466 PyObject*
PyPolynomial_FromLong_or_Int(PyObject * number,const lp_polynomial_context_t * ctx)467 PyPolynomial_FromLong_or_Int(PyObject* number, const lp_polynomial_context_t* ctx) {
468 // The constants
469 lp_integer_t c;
470 PyLong_or_Int_to_integer(number, 0, &c);
471
472 // The c polynomial
473 lp_polynomial_t* p_c = lp_polynomial_alloc();
474 lp_polynomial_construct_simple(p_c, ctx, &c, 0, 0);
475
476 // Remove temps
477 lp_integer_destruct(&c);
478
479 // Return the polynomial
480 PyObject* result = Polynomial_create(p_c);
481 Py_INCREF(result);
482 return result;
483 }
484
485 static PyObject*
Polynomial_add(PyObject * self,PyObject * other)486 Polynomial_add(PyObject* self, PyObject* other) {
487
488 int dec_other = 0;
489
490 if (!PyPolynomial_CHECK(self)) {
491 return Polynomial_add(other, self);
492 }
493
494 Polynomial* p1 = (Polynomial*) self;
495 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
496
497 // Check argument
498 if (!PyPolynomial_CHECK(other)) {
499 if (PyVariable_CHECK(other)) {
500 other = PyPolynomial_FromVariable(other, p1_ctx);
501 dec_other = 1;
502 } else if (PyLong_or_Int_Check(other)) {
503 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
504 dec_other = 1;
505 } else {
506 Py_INCREF(Py_NotImplemented);
507 return Py_NotImplemented;
508 }
509 }
510
511 // Get arguments
512 Polynomial* p2 = (Polynomial*) other;
513 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
514 if (p1_ctx != p2_ctx) {
515 Py_INCREF(Py_NotImplemented);
516 return Py_NotImplemented;
517 }
518
519 // Add the polynomials
520 lp_polynomial_t* sum = lp_polynomial_new(p1_ctx);
521 lp_polynomial_add(sum, p1->p, p2->p);
522
523 if (dec_other) {
524 Py_DECREF(other);
525 }
526
527 // Return the result
528 return Polynomial_create(sum);
529 }
530
531 static PyObject*
Polynomial_neg(PyObject * self)532 Polynomial_neg(PyObject* self) {
533 Polynomial* p = (Polynomial*) self;
534 const lp_polynomial_context_t* p_ctx = lp_polynomial_get_context(p->p);
535 lp_polynomial_t* neg = lp_polynomial_new(p_ctx);
536 lp_polynomial_neg(neg, p->p);
537 return Polynomial_create(neg);
538 }
539
540 static PyObject*
Polynomial_sub(PyObject * self,PyObject * other)541 Polynomial_sub(PyObject* self, PyObject* other) {
542
543 int dec_other = 0;
544
545 if (!PyPolynomial_CHECK(self)) {
546 Polynomial* sub = (Polynomial*) Polynomial_sub(other, self);
547 lp_polynomial_neg(sub->p, sub->p);
548 return (PyObject*) sub;
549 }
550
551 Polynomial* p1 = (Polynomial*) self;
552 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
553
554 // Check argument
555 if (!PyPolynomial_CHECK(other)) {
556 if (PyVariable_CHECK(other)) {
557 other = PyPolynomial_FromVariable(other, p1_ctx);
558 dec_other = 1;
559 } else if (PyLong_or_Int_Check(other)) {
560 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
561 dec_other = 1;
562 } else {
563 Py_INCREF(Py_NotImplemented);
564 return Py_NotImplemented;
565 }
566 }
567
568 // Get arguments
569 Polynomial* p2 = (Polynomial*) other;
570 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
571 if (p1_ctx != p2_ctx) {
572 Py_INCREF(Py_NotImplemented);
573 return Py_NotImplemented;
574 }
575
576 // Subtract the polynomials
577 lp_polynomial_t* sub = lp_polynomial_new(p1_ctx);
578 lp_polynomial_sub(sub, p1->p, p2->p);
579
580 if (dec_other) {
581 Py_DECREF(other);
582 }
583
584 // Return the result
585 return Polynomial_create(sub);
586 }
587
588 static PyObject*
Polynomial_mul(PyObject * self,PyObject * other)589 Polynomial_mul(PyObject* self, PyObject* other) {
590
591 int dec_other = 0;
592
593 if (!PyPolynomial_CHECK(self)) {
594 return Polynomial_mul(other, self);
595 }
596
597 Polynomial* p1 = (Polynomial*) self;
598 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
599
600 // Check argument
601 if (!PyPolynomial_CHECK(other)) {
602 if (PyVariable_CHECK(other)) {
603 other = PyPolynomial_FromVariable(other, p1_ctx);
604 dec_other = 1;
605 } else if (PyLong_or_Int_Check(other)) {
606 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
607 dec_other = 1;
608 } else {
609 Py_INCREF(Py_NotImplemented);
610 return Py_NotImplemented;
611 }
612 }
613
614 // Get arguments
615 Polynomial* p2 = (Polynomial*) other;
616 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
617 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
618 Py_INCREF(Py_NotImplemented);
619 return Py_NotImplemented;
620 }
621
622 // Multiply the polynomials
623 lp_polynomial_t* mul = lp_polynomial_new(p1_ctx);
624 lp_polynomial_mul(mul, p1->p, p2->p);
625
626 if (dec_other) {
627 Py_DECREF(other);
628 }
629
630 // Return the result
631 return Polynomial_create(mul);
632 }
633
634
635 static PyObject*
Polynomial_pow(PyObject * self,PyObject * other)636 Polynomial_pow(PyObject* self, PyObject* other) {
637 // Check arguments
638 if (!PyPolynomial_CHECK(self) || !PyInt_Check(other)) {
639 Py_INCREF(Py_NotImplemented);
640 return Py_NotImplemented;
641 }
642 // Get arguments
643 Polynomial* p = (Polynomial*) self;
644 long n = PyInt_AsLong(other);
645 if (n < 0) {
646 Py_INCREF(Py_NotImplemented);
647 return Py_NotImplemented;
648 }
649 const lp_polynomial_context_t* p_ctx = lp_polynomial_get_context(p->p);
650 // Compute
651 lp_polynomial_t* pow = lp_polynomial_new(p_ctx);
652 lp_polynomial_pow(pow, p->p, n);
653 // Return the result
654 return Polynomial_create(pow);
655 }
656
657 static PyObject*
Polynomial_div(PyObject * self,PyObject * other)658 Polynomial_div(PyObject* self, PyObject* other) {
659
660 int dec_other = 0;
661
662 if (!PyPolynomial_CHECK(self)) {
663 Py_INCREF(Py_NotImplemented);
664 return Py_NotImplemented;
665 }
666
667 // self is always a polynomial
668 Polynomial* p1 = (Polynomial*) self;
669 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
670
671 // Make sure other is a polynomial
672 if (!PyPolynomial_CHECK(other)) {
673 if (PyVariable_CHECK(other)) {
674 other = PyPolynomial_FromVariable(other, p1_ctx);
675 dec_other = 1;
676 } else if (PyLong_or_Int_Check(other)) {
677 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
678 dec_other = 1;
679 } else {
680 Py_INCREF(Py_NotImplemented);
681 return Py_NotImplemented;
682 }
683 }
684
685 // other can be a variable or a number
686 Polynomial* p2 = (Polynomial*) other;
687 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
688 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
689 Py_INCREF(Py_NotImplemented);
690 return Py_NotImplemented;
691 }
692
693 // Multiply the polynomials
694 lp_polynomial_t* div = lp_polynomial_new(p1_ctx);
695 lp_polynomial_div(div, p1->p, p2->p);
696
697 if (dec_other) {
698 Py_DECREF(other);
699 }
700
701 // Return the result
702 return Polynomial_create(div);
703 }
704
705 static PyObject*
Polynomial_rem_operator(PyObject * self,PyObject * other)706 Polynomial_rem_operator(PyObject* self, PyObject* other) {
707 int dec_other = 0;
708
709 if (!PyPolynomial_CHECK(self)) {
710 Py_INCREF(Py_NotImplemented);
711 return Py_NotImplemented;
712 }
713
714 // self is always a polynomial
715 Polynomial* p1 = (Polynomial*) self;
716 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
717
718 // Make sure other is a polynomial
719 if (!PyPolynomial_CHECK(other)) {
720 if (PyVariable_CHECK(other)) {
721 other = PyPolynomial_FromVariable(other, p1_ctx);
722 dec_other = 1;
723 } else if (PyLong_or_Int_Check(other)) {
724 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
725 dec_other = 1;
726 } else {
727 Py_INCREF(Py_NotImplemented);
728 return Py_NotImplemented;
729 }
730 }
731
732 // other can be a variable or a number
733 Polynomial* p2 = (Polynomial*) other;
734 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
735 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
736 Py_INCREF(Py_NotImplemented);
737 return Py_NotImplemented;
738 }
739
740 // Multiply the polynomials
741 lp_polynomial_t* rem = lp_polynomial_new(p1_ctx);
742 lp_polynomial_rem(rem, p1->p, p2->p);
743
744 if (dec_other) {
745 Py_DECREF(other);
746 }
747
748 // Return the result
749 return Polynomial_create(rem);
750 }
751
752 static PyObject*
Polynomial_divmod(PyObject * self,PyObject * other)753 Polynomial_divmod(PyObject* self, PyObject* other) {
754
755 int dec_other = 0;
756
757 if (!PyPolynomial_CHECK(self)) {
758 Py_INCREF(Py_NotImplemented);
759 return Py_NotImplemented;
760 }
761
762 // self is always a polynomial
763 Polynomial* p1 = (Polynomial*) self;
764 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
765
766 // Make sure other is a polynomial
767 if (!PyPolynomial_CHECK(other)) {
768 if (PyVariable_CHECK(other)) {
769 other = PyPolynomial_FromVariable(other, p1_ctx);
770 dec_other = 1;
771 } else if (PyLong_or_Int_Check(other)) {
772 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
773 dec_other = 1;
774 } else {
775 Py_INCREF(Py_NotImplemented);
776 return Py_NotImplemented;
777 }
778 }
779
780 // other can be a variable or a number
781 Polynomial* p2 = (Polynomial*) other;
782 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
783 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
784 Py_INCREF(Py_NotImplemented);
785 return Py_NotImplemented;
786 }
787
788 // Multiply the polynomials
789 lp_polynomial_t* rem = lp_polynomial_new(p1_ctx);
790 lp_polynomial_t* div = lp_polynomial_new(p1_ctx);
791 lp_polynomial_divrem(div, rem, p1->p, p2->p);
792
793 if (dec_other) {
794 Py_DECREF(other);
795 }
796
797 // Return the result
798 PyObject* pair = PyTuple_New(2);
799 PyObject* divObj = Polynomial_create(div);
800 PyObject* remObj = Polynomial_create(rem);
801 Py_INCREF(divObj);
802 Py_INCREF(remObj);
803 PyTuple_SetItem(pair, 0, divObj);
804 PyTuple_SetItem(pair, 1, remObj);
805 return pair;
806 }
807
808 static int
Polynomial_nonzero(PyObject * self)809 Polynomial_nonzero(PyObject* self) {
810 // Get arguments
811 Polynomial* p = (Polynomial*) self;
812 // Return the result
813 return !lp_polynomial_is_zero(p->p);
814 }
815
816 enum rem_type {
817 REM_EXACT,
818 REM_PSEUDO,
819 REM_SPARSE_PSEUDO
820 };
821
822 static PyObject*
Polynomial_rem_general(PyObject * self,PyObject * args,enum rem_type type)823 Polynomial_rem_general(PyObject* self, PyObject* args, enum rem_type type) {
824 int dec_other = 0;
825
826 // self is always a polynomial
827 Polynomial* p1 = (Polynomial*) self;
828 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
829
830 if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
831 Py_INCREF(Py_NotImplemented);
832 return Py_NotImplemented;
833 }
834
835 PyObject* other = PyTuple_GetItem(args, 0);
836
837 // Make sure other is a polynomial
838 if (!PyPolynomial_CHECK(other)) {
839 if (PyVariable_CHECK(other)) {
840 other = PyPolynomial_FromVariable(other, p1_ctx);
841 dec_other = 1;
842 } else if (PyLong_or_Int_Check(other)) {
843 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
844 dec_other = 1;
845 } else {
846 Py_INCREF(Py_NotImplemented);
847 return Py_NotImplemented;
848 }
849 }
850
851 // other can be a variable or a number
852 Polynomial* p2 = (Polynomial*) other;
853 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
854 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
855 Py_INCREF(Py_NotImplemented);
856 return Py_NotImplemented;
857 }
858
859 // Multiply the polynomials
860 lp_polynomial_t* rem = lp_polynomial_new(p1_ctx);
861 switch (type) {
862 case REM_EXACT:
863 lp_polynomial_rem(rem, p1->p, p2->p);
864 break;
865 case REM_PSEUDO:
866 lp_polynomial_prem(rem, p1->p, p2->p);
867 break;
868 case REM_SPARSE_PSEUDO:
869 lp_polynomial_sprem(rem, p1->p, p2->p);
870 break;
871 }
872
873 if (dec_other) {
874 Py_DECREF(other);
875 }
876
877 // Return the result
878 return Polynomial_create(rem);
879 }
880
881 static PyObject*
Polynomial_rem(PyObject * self,PyObject * other)882 Polynomial_rem(PyObject* self, PyObject* other) {
883 return Polynomial_rem_general(self, other, REM_EXACT);
884 }
885
886 static PyObject*
Polynomial_prem(PyObject * self,PyObject * other)887 Polynomial_prem(PyObject* self, PyObject* other) {
888 return Polynomial_rem_general(self, other, REM_PSEUDO);
889 }
890
891 static PyObject*
Polynomial_sprem(PyObject * self,PyObject * other)892 Polynomial_sprem(PyObject* self, PyObject* other) {
893 return Polynomial_rem_general(self, other, REM_SPARSE_PSEUDO);
894 }
895
896 static PyObject*
Polynomial_gcd(PyObject * self,PyObject * args)897 Polynomial_gcd(PyObject* self, PyObject* args) {
898
899 int dec_other = 0;
900
901 // self is always a polynomial
902 Polynomial* p1 = (Polynomial*) self;
903 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
904
905 if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
906 Py_INCREF(Py_NotImplemented);
907 return Py_NotImplemented;
908 }
909
910 PyObject* other = PyTuple_GetItem(args, 0);
911
912 // Make sure other is a polynomial
913 if (!PyPolynomial_CHECK(other)) {
914 if (PyVariable_CHECK(other)) {
915 other = PyPolynomial_FromVariable(other, p1_ctx);
916 dec_other = 1;
917 } else if (PyLong_or_Int_Check(other)) {
918 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
919 dec_other = 1;
920 } else {
921 Py_INCREF(Py_NotImplemented);
922 return Py_NotImplemented;
923 }
924 }
925
926 // other can be a variable or a number
927 Polynomial* p2 = (Polynomial*) other;
928 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
929 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
930 Py_INCREF(Py_NotImplemented);
931 return Py_NotImplemented;
932 }
933
934 // Multiply the polynomials
935 lp_polynomial_t* gcd = lp_polynomial_new(p1_ctx);
936 lp_polynomial_gcd(gcd, p1->p, p2->p);
937
938 if (dec_other) {
939 Py_DECREF(other);
940 }
941
942 // Return the result
943 return Polynomial_create(gcd);
944 }
945
946 static PyObject*
Polynomial_lcm(PyObject * self,PyObject * args)947 Polynomial_lcm(PyObject* self, PyObject* args) {
948
949 int dec_other = 0;
950
951 // self is always a polynomial
952 Polynomial* p1 = (Polynomial*) self;
953 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
954
955 if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
956 Py_INCREF(Py_NotImplemented);
957 return Py_NotImplemented;
958 }
959
960 PyObject* other = PyTuple_GetItem(args, 0);
961
962 // Make sure other is a polynomial
963 if (!PyPolynomial_CHECK(other)) {
964 if (PyVariable_CHECK(other)) {
965 other = PyPolynomial_FromVariable(other, p1_ctx);
966 dec_other = 1;
967 } else if (PyLong_or_Int_Check(other)) {
968 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
969 dec_other = 1;
970 } else {
971 Py_INCREF(Py_NotImplemented);
972 return Py_NotImplemented;
973 }
974 }
975
976 // other can be a variable or a number
977 Polynomial* p2 = (Polynomial*) other;
978 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
979 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
980 Py_INCREF(Py_NotImplemented);
981 return Py_NotImplemented;
982 }
983
984 // Multiply the polynomials
985 lp_polynomial_t* lcm = lp_polynomial_new(p1_ctx);
986 lp_polynomial_lcm(lcm, p1->p, p2->p);
987
988 if (dec_other) {
989 Py_DECREF(other);
990 }
991
992 // Return the result
993 return Polynomial_create(lcm);
994 }
995
996 static PyObject*
Polynomial_psc(PyObject * self,PyObject * args)997 Polynomial_psc(PyObject* self, PyObject* args) {
998
999 // self is always a polynomial
1000 Polynomial* p1 = (Polynomial*) self;
1001 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
1002
1003 if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
1004 Py_INCREF(Py_NotImplemented);
1005 return Py_NotImplemented;
1006 }
1007
1008 PyObject* other = PyTuple_GetItem(args, 0);
1009
1010 int dec_other = 0;
1011
1012 // Make sure other is a polynomial
1013 if (!PyPolynomial_CHECK(other)) {
1014 if (PyVariable_CHECK(other)) {
1015 other = PyPolynomial_FromVariable(other, p1_ctx);
1016 dec_other = 1;
1017 } else if (PyLong_or_Int_Check(other)) {
1018 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
1019 dec_other = 1;
1020 } else {
1021 Py_INCREF(Py_NotImplemented);
1022 return Py_NotImplemented;
1023 }
1024 }
1025
1026 // Othe polynomial
1027 Polynomial* p2 = (Polynomial*) other;
1028 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
1029 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
1030 Py_INCREF(Py_NotImplemented);
1031 return Py_NotImplemented;
1032 }
1033
1034 // Check the arguments (must be same top variable)
1035 if (lp_polynomial_is_constant(p1->p) || lp_polynomial_is_constant(p2->p)) {
1036 Py_INCREF(Py_NotImplemented);
1037 return Py_NotImplemented;
1038 }
1039
1040 if (lp_polynomial_top_variable(p1->p) != lp_polynomial_top_variable(p2->p)) {
1041 Py_INCREF(Py_NotImplemented);
1042 return Py_NotImplemented;
1043 }
1044
1045 // Allocate the polynomials for the sequence
1046 size_t p1_deg = lp_polynomial_degree(p1->p);
1047 size_t p2_deg = lp_polynomial_degree(p2->p);
1048 int size = p1_deg > p2_deg ? p2_deg + 1 : p1_deg + 1;
1049
1050 lp_polynomial_t** psc = malloc(sizeof(lp_polynomial_t*)*size);
1051 int i;
1052 for (i = 0; i < size; ++ i) {
1053 psc[i] = lp_polynomial_new(p1_ctx);
1054 }
1055
1056 // Compute the psc
1057 lp_polynomial_psc(psc, p1->p, p2->p);
1058
1059 // Copy the polynomials into a list
1060 PyObject* list = PyList_New(size);
1061 for (i = 0; i < size; ++i) {
1062 PyObject* p = Polynomial_create(psc[i]);
1063 PyList_SetItem(list, i, p);
1064 }
1065
1066 if (dec_other) {
1067 Py_DECREF(other);
1068 }
1069
1070 // Return the result
1071 return list;
1072 }
1073
1074 static PyObject*
Polynomial_mgcd(PyObject * self,PyObject * args)1075 Polynomial_mgcd(PyObject* self, PyObject* args) {
1076
1077 // self is always a polynomial
1078 Polynomial* p1 = (Polynomial*) self;
1079 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
1080
1081 if (!PyTuple_Check(args) || PyTuple_Size(args) != 2) {
1082 PyErr_SetString(PyExc_RuntimeError, "mgcd(): Need two arguments.");
1083 return NULL;
1084 }
1085
1086 // Assignment
1087 PyObject* py_assignment = PyTuple_GetItem(args, 1);
1088 if (!PyAssignment_CHECK(py_assignment)) {
1089 PyErr_SetString(PyExc_RuntimeError, "mgcd(): Second argument should be an assignment.");
1090 return NULL;
1091 }
1092 const lp_assignment_t* assignment = ((Assignment*) py_assignment)->assignment;
1093
1094 // Other polynomial
1095 PyObject* other = PyTuple_GetItem(args, 0);
1096
1097 int dec_other = 0;
1098
1099 // Make sure other is a polynomial
1100 if (!PyPolynomial_CHECK(other)) {
1101 if (PyVariable_CHECK(other)) {
1102 other = PyPolynomial_FromVariable(other, p1_ctx);
1103 dec_other = 1;
1104 } else {
1105 PyErr_SetString(PyExc_RuntimeError, "mgcd(): First argument should be a polynomial.");
1106 return NULL;
1107 }
1108 }
1109
1110 // Othe polynomial
1111 Polynomial* p2 = (Polynomial*) other;
1112 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
1113 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
1114 PyErr_SetString(PyExc_RuntimeError, "mgcd(): Polynomials should be over the same context.");
1115 return NULL;
1116 }
1117
1118 // Check the arguments (must be same top variable)
1119 if (lp_polynomial_is_constant(p1->p) || lp_polynomial_is_constant(p2->p)) {
1120 PyErr_SetString(PyExc_RuntimeError, "mgcd(): Polynomials should be over the same top variables.");
1121 return NULL;
1122 }
1123
1124 if (lp_polynomial_top_variable(p1->p) != lp_polynomial_top_variable(p2->p)) {
1125 PyErr_SetString(PyExc_RuntimeError, "mgcd(): Polynomials should be over the same top variables.");
1126 return NULL;
1127 }
1128
1129 // Compute the gcd
1130 lp_polynomial_vector_t* mgcd = lp_polynomial_mgcd(p1->p, p2->p, assignment);
1131
1132 // Copy the polynomials into a list
1133 size_t size = lp_polynomial_vector_size(mgcd);
1134 PyObject* list = PyList_New(size);
1135 size_t i;
1136 for (i = 0; i < size; ++i) {
1137 lp_polynomial_t* mgcd_i = lp_polynomial_vector_at(mgcd, i);
1138 PyObject* p = Polynomial_create(mgcd_i);
1139 PyList_SetItem(list, i, p);
1140 }
1141 lp_polynomial_vector_delete(mgcd);
1142
1143 if (dec_other) {
1144 Py_DECREF(other);
1145 }
1146
1147 // Return the result
1148 return list;
1149 }
1150
1151 static PyObject*
Polynomial_resultant(PyObject * self,PyObject * args)1152 Polynomial_resultant(PyObject* self, PyObject* args) {
1153
1154 // self is always a polynomial
1155 Polynomial* p1 = (Polynomial*) self;
1156 const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
1157
1158 if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
1159 Py_INCREF(Py_NotImplemented);
1160 return Py_NotImplemented;
1161 }
1162
1163 PyObject* other = PyTuple_GetItem(args, 0);
1164
1165 int dec_other = 0;
1166
1167 // Make sure other is a polynomial
1168 if (!PyPolynomial_CHECK(other)) {
1169 if (PyVariable_CHECK(other)) {
1170 other = PyPolynomial_FromVariable(other, p1_ctx);
1171 dec_other = 1;
1172 } else if (PyLong_or_Int_Check(other)) {
1173 other = PyPolynomial_FromLong_or_Int(other, p1_ctx);
1174 dec_other = 1;
1175 } else {
1176 Py_INCREF(Py_NotImplemented);
1177 return Py_NotImplemented;
1178 }
1179 }
1180
1181 // Othe polynomial
1182 Polynomial* p2 = (Polynomial*) other;
1183 const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
1184 if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
1185 Py_INCREF(Py_NotImplemented);
1186 return Py_NotImplemented;
1187 }
1188
1189 // Check the arguments (must be same top variable)
1190 if (lp_polynomial_is_constant(p1->p) || lp_polynomial_is_constant(p2->p)) {
1191 Py_INCREF(Py_NotImplemented);
1192 return Py_NotImplemented;
1193 }
1194
1195 if (lp_polynomial_top_variable(p1->p) != lp_polynomial_top_variable(p2->p)) {
1196 Py_INCREF(Py_NotImplemented);
1197 return Py_NotImplemented;
1198 }
1199
1200 // Allocate the resultant
1201 lp_polynomial_t* resultant = lp_polynomial_new(p1_ctx);
1202
1203 // Compute the psc
1204 lp_polynomial_resultant(resultant, p1->p, p2->p);
1205
1206 if (dec_other) {
1207 Py_DECREF(other);
1208 }
1209
1210 // Return the result
1211 return Polynomial_create(resultant);
1212 }
1213
1214
1215 static PyObject*
Polynomial_extended_gcd(PyObject * self,PyObject * args)1216 Polynomial_extended_gcd(PyObject* self, PyObject* args) {
1217 return 0;
1218 }
1219
1220 static PyObject*
Polynomial_factor(PyObject * self)1221 Polynomial_factor(PyObject* self) {
1222 return 0;
1223 }
1224
1225 // Creates a python list from the factors, taking over the polynomials
factors_to_PyList(lp_polynomial_t ** factors,size_t * multiplicities,size_t size)1226 PyObject* factors_to_PyList(lp_polynomial_t** factors, size_t* multiplicities, size_t size) {
1227 // Construct the result
1228 PyObject* factors_list = PyList_New(size);
1229
1230 // Copy the constant
1231 // Copy over the factors
1232 int i;
1233 for (i = 0; i < size; ++ i) {
1234 PyObject* p_i = Polynomial_create(factors[i]);
1235 Py_INCREF(p_i);
1236 PyObject* d = PyInt_FromSize_t(multiplicities[i]);
1237 PyObject* pair = PyTuple_New(2);
1238 PyTuple_SetItem(pair, 0, p_i);
1239 PyTuple_SetItem(pair, 1, d);
1240 PyList_SetItem(factors_list, i, pair);
1241 }
1242
1243 // Return the list
1244 return factors_list;
1245 }
1246
1247 static PyObject*
Polynomial_factor_square_free(PyObject * self)1248 Polynomial_factor_square_free(PyObject* self) {
1249 // Get arguments
1250 Polynomial* p = (Polynomial*) self;
1251 // Factor
1252 lp_polynomial_t** factors = 0;
1253 size_t* multiplicities = 0;
1254 size_t factors_size = 0;
1255 lp_polynomial_factor_square_free(p->p, &factors, &multiplicities, &factors_size);
1256 // Create the list
1257 PyObject* factors_list = factors_to_PyList(factors, multiplicities, factors_size);
1258 // Get rid of the factors (not the polynomials)
1259 free(factors);
1260 free(multiplicities);
1261 // Return the list
1262 return factors_list;
1263 }
1264
1265 static PyObject*
Polynomial_roots_count(PyObject * self,PyObject * args)1266 Polynomial_roots_count(PyObject* self, PyObject* args) {
1267 return 0;
1268 }
1269
1270 static PyObject*
Polynomial_roots_isolate(PyObject * self,PyObject * args)1271 Polynomial_roots_isolate(PyObject* self, PyObject* args) {
1272
1273 int i;
1274
1275 if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
1276 Py_INCREF(Py_NotImplemented);
1277 return Py_NotImplemented;
1278 }
1279
1280 PyObject* assignment_obj = PyTuple_GetItem(args, 0);
1281
1282 if (!PyAssignment_CHECK(assignment_obj)) {
1283 Py_INCREF(Py_NotImplemented);
1284 return Py_NotImplemented;
1285 }
1286
1287 lp_polynomial_t* p = ((Polynomial*) self)->p;
1288 lp_assignment_t* assignment = ((Assignment*) assignment_obj)->assignment;
1289
1290 // Check that the top variable is the only unassigned
1291 if (!lp_polynomial_is_univariate_m(p, assignment)) {
1292 PyErr_SetString(PyExc_RuntimeError, "roots_count(): Polynomial must be univariate modulo the assignment.");
1293 return NULL;
1294 }
1295
1296 // Get the degree of the polynomial and allocate the values
1297 lp_value_t* roots = malloc(sizeof(lp_value_t)*lp_polynomial_degree(p));
1298 size_t roots_size = 0;
1299
1300 // Get the roots
1301 lp_polynomial_roots_isolate(p, assignment, roots, &roots_size);
1302
1303 // Generate a list of roots
1304 PyObject* list = PyList_New(roots_size);
1305
1306 for (i = 0; i < roots_size; ++ i) {
1307 PyObject* c = PyValue_create(roots + i);
1308 PyList_SetItem(list, i, c);
1309 }
1310
1311 // Get rid of the temporaries
1312 for (i = 0; i < roots_size; ++ i) {
1313 lp_value_destruct(roots + i);
1314 }
1315 free(roots);
1316
1317 // Return the list
1318 return list;
1319 }
1320
1321 static PyObject*
Polynomial_derivative(PyObject * self)1322 Polynomial_derivative(PyObject* self) {
1323 lp_polynomial_t* p = ((Polynomial*) self)->p;
1324 lp_polynomial_t* p_derivative = lp_polynomial_new(lp_polynomial_get_context(p));
1325 lp_polynomial_derivative(p_derivative, p);
1326 return Polynomial_create(p_derivative);
1327 }
1328
1329 static PyObject*
Polynomial_sturm_sequence(PyObject * self)1330 Polynomial_sturm_sequence(PyObject* self) {
1331 return 0;
1332 }
1333
1334 static PyObject*
Polynomial_degree(PyObject * self)1335 Polynomial_degree(PyObject* self) {
1336 Polynomial* p = (Polynomial*) self;
1337 return PyInt_FromLong(lp_polynomial_degree(p->p));
1338 }
1339
1340 static PyObject*
Polynomial_coefficients(PyObject * self)1341 Polynomial_coefficients(PyObject* self) {
1342 int i;
1343
1344 lp_polynomial_t* p = ((Polynomial*) self)->p;
1345 size_t size = lp_polynomial_degree(p) + 1;
1346
1347 // Get the coefficients
1348 const lp_polynomial_context_t* ctx = lp_polynomial_get_context(p);
1349
1350 // Copy the polynomials into a list
1351 PyObject* list = PyList_New(size);
1352 for (i = 0; i < size; ++i) {
1353 lp_polynomial_t* c_p = lp_polynomial_new(ctx);
1354 lp_polynomial_get_coefficient(c_p, p, i);
1355 PyObject* c = Polynomial_create(c_p);
1356 PyList_SetItem(list, i, c);
1357 }
1358
1359 return list;
1360 }
1361
1362 static PyObject*
Polynomial_vars(PyObject * self)1363 Polynomial_vars(PyObject* self) {
1364
1365 lp_polynomial_t* p = ((Polynomial*) self)->p;
1366
1367 lp_variable_list_t p_vars;
1368 lp_variable_list_construct(&p_vars);
1369
1370 // Get the variables
1371 lp_polynomial_get_variables(p, &p_vars);
1372
1373 // Copy the polynomials into a list
1374 PyObject* list = PyList_New(p_vars.list_size);
1375 size_t i;
1376 for (i = 0; i < p_vars.list_size; ++i) {
1377 PyObject* c = PyVariable_create(p_vars.list[i]);
1378 PyList_SetItem(list, i, c);
1379 }
1380
1381 lp_variable_list_destruct(&p_vars);
1382
1383 return list;
1384 }
1385
1386 static PyObject*
Polynomial_var(PyObject * self)1387 Polynomial_var(PyObject* self) {
1388 lp_polynomial_t* p = ((Polynomial*) self)->p;
1389 return PyVariable_create(lp_polynomial_top_variable(p));
1390 }
1391
1392
1393 static PyObject*
Polynomial_reductum(PyObject * self,PyObject * args)1394 Polynomial_reductum(PyObject* self, PyObject* args) {
1395 lp_polynomial_t* p = ((Polynomial*) self)->p;
1396 const lp_polynomial_context_t* ctx = lp_polynomial_get_context(p);
1397
1398 if (!PyTuple_Check(args) || PyTuple_Size(args) > 1) {
1399 Py_INCREF(Py_NotImplemented);
1400 return Py_NotImplemented;
1401 }
1402
1403 lp_assignment_t* assignment = 0;
1404
1405 if (PyTuple_Size(args) == 1) {
1406 PyObject* assignment_obj = PyTuple_GetItem(args, 0);
1407 if (!PyAssignment_CHECK(assignment_obj)) {
1408 Py_INCREF(Py_NotImplemented);
1409 return Py_NotImplemented;
1410 } else {
1411 assignment = ((Assignment*) assignment_obj)->assignment;
1412 }
1413 }
1414
1415 lp_polynomial_t* result = lp_polynomial_new(ctx);
1416 if (assignment) {
1417 lp_polynomial_reductum_m(result, p, assignment);
1418 } else {
1419 lp_polynomial_reductum(result, p);
1420 }
1421
1422 return Polynomial_create(result);
1423 }
1424
1425 static PyObject*
Polynomial_sgn(PyObject * self,PyObject * args)1426 Polynomial_sgn(PyObject* self, PyObject* args) {
1427
1428 if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
1429 Py_INCREF(Py_NotImplemented);
1430 return Py_NotImplemented;
1431 }
1432
1433 PyObject* assignment_obj = PyTuple_GetItem(args, 0);
1434
1435 if (!PyAssignment_CHECK(assignment_obj)) {
1436 Py_INCREF(Py_NotImplemented);
1437 return Py_NotImplemented;
1438 }
1439
1440 lp_polynomial_t* p = ((Polynomial*) self)->p;
1441 lp_assignment_t* assignment = ((Assignment*) assignment_obj)->assignment;
1442
1443 int sgn = lp_polynomial_sgn(p, assignment);
1444
1445 return PyInt_FromLong(sgn);
1446 }
1447
1448
1449 static PyObject*
Polynomial_evaluate(PyObject * self,PyObject * args)1450 Polynomial_evaluate(PyObject* self, PyObject* args) {
1451
1452 if (!PyTuple_Check(args) || PyTuple_Size(args) != 1) {
1453 Py_INCREF(Py_NotImplemented);
1454 return Py_NotImplemented;
1455 }
1456
1457 PyObject* assignment_obj = PyTuple_GetItem(args, 0);
1458
1459 if (!PyAssignment_CHECK(assignment_obj)) {
1460 Py_INCREF(Py_NotImplemented);
1461 return Py_NotImplemented;
1462 }
1463
1464 lp_polynomial_t* p = ((Polynomial*) self)->p;
1465 lp_assignment_t* assignment = ((Assignment*) assignment_obj)->assignment;
1466
1467 lp_value_t* value = lp_polynomial_evaluate(p, assignment);
1468 PyObject* value_obj = PyValue_create(value);
1469 lp_value_delete(value);
1470
1471 return value_obj;
1472 }
1473
1474 static PyObject*
Polynomial_feasible_intervals(PyObject * self,PyObject * args)1475 Polynomial_feasible_intervals(PyObject* self, PyObject* args) {
1476
1477 if (!PyTuple_Check(args) || PyTuple_Size(args) != 2) {
1478 PyErr_SetString(PyExc_RuntimeError, "feasible_intervals(): Needs two arguments, an assignment and a sign condition.");
1479 return NULL;
1480 }
1481
1482 PyObject* assignment_obj = PyTuple_GetItem(args, 0);
1483 if (!PyAssignment_CHECK(assignment_obj)) {
1484 PyErr_SetString(PyExc_RuntimeError, "feasible_intervals(): First argument must be an assignment.");
1485 return NULL;
1486 }
1487
1488 PyObject* sgn_condition_obj = PyTuple_GetItem(args, 1);
1489 if (!PyInt_Check(sgn_condition_obj)) {
1490 PyErr_SetString(PyExc_RuntimeError, "feasible_intervals(): Second argument must be a sign-condition.");
1491 return NULL;
1492 }
1493
1494 // Get the arguments
1495 lp_polynomial_t* p = ((Polynomial*) self)->p;
1496 lp_assignment_t* assignment = ((Assignment*) assignment_obj)->assignment;
1497 lp_sign_condition_t sgn_condition = PyInt_AsLong(sgn_condition_obj);
1498
1499 // Check if all but the top variable are unassigned
1500 if (!lp_polynomial_is_univariate_m(p, assignment)) {
1501 PyErr_SetString(PyExc_RuntimeError, "feasible_intervals(): Polynomial must be univariate modulo the assignment.");
1502 return NULL;
1503 }
1504
1505 // Get the feasible intervals
1506 lp_feasibility_set_t* feasible = lp_polynomial_constraint_get_feasible_set(p, sgn_condition, 0, assignment);
1507
1508 // The list where we return the arguments
1509 PyObject* list = PyList_New(feasible->size);
1510 // Copy over to the list
1511 size_t i;
1512 for (i = 0; i < feasible->size; ++i) {
1513 PyObject* p = PyInterval_create(feasible->intervals + i);
1514 PyList_SetItem(list, i, p);
1515 }
1516 // Remove temp
1517 lp_feasibility_set_delete(feasible);
1518
1519 // Return the list
1520 return list;
1521 }
1522
1523 static PyObject*
Polynomial_feasible_set(PyObject * self,PyObject * args)1524 Polynomial_feasible_set(PyObject* self, PyObject* args) {
1525
1526 if (!PyTuple_Check(args) || PyTuple_Size(args) != 2) {
1527 Py_INCREF(Py_NotImplemented);
1528 return Py_NotImplemented;
1529 }
1530
1531 PyObject* assignment_obj = PyTuple_GetItem(args, 0);
1532 if (!PyAssignment_CHECK(assignment_obj)) {
1533 Py_INCREF(Py_NotImplemented);
1534 return Py_NotImplemented;
1535 }
1536
1537 PyObject* sgn_condition_obj = PyTuple_GetItem(args, 1);
1538 if (!PyInt_Check(sgn_condition_obj)) {
1539 Py_INCREF(Py_NotImplemented);
1540 return Py_NotImplemented;
1541 }
1542
1543 // Get the arguments
1544 lp_polynomial_t* p = ((Polynomial*) self)->p;
1545 lp_assignment_t* assignment = ((Assignment*) assignment_obj)->assignment;
1546 lp_sign_condition_t sgn_condition = PyInt_AsLong(sgn_condition_obj);
1547
1548 // Check if all but the top variable are unassigned
1549 if (!lp_polynomial_is_univariate_m(p, assignment)) {
1550 PyErr_SetString(PyExc_RuntimeError, "feasible_set(): Polynomial must be univariate modulo the assignment.");
1551 return NULL;
1552 }
1553
1554 // Get the feasible intervals
1555 lp_feasibility_set_t* feasible = lp_polynomial_constraint_get_feasible_set(p, sgn_condition, 0, assignment);
1556
1557 // Return the list
1558 return PyFeasibilitySet_create(feasible);
1559 }
1560
1561 static PyObject*
Polynomial_sgn_check(PyObject * self,PyObject * args)1562 Polynomial_sgn_check(PyObject* self, PyObject* args) {
1563
1564 if (!PyTuple_Check(args) || PyTuple_Size(args) != 2) {
1565 Py_INCREF(Py_NotImplemented);
1566 return Py_NotImplemented;
1567 }
1568
1569 PyObject* assignment_obj = PyTuple_GetItem(args, 0);
1570 if (!PyAssignment_CHECK(assignment_obj)) {
1571 Py_INCREF(Py_NotImplemented);
1572 return Py_NotImplemented;
1573 }
1574
1575 PyObject* sgn_condition_obj = PyTuple_GetItem(args, 1);
1576 if (!PyInt_Check(sgn_condition_obj)) {
1577 Py_INCREF(Py_NotImplemented);
1578 return Py_NotImplemented;
1579 }
1580
1581 // Get the arguments
1582 lp_polynomial_t* p = ((Polynomial*) self)->p;
1583 lp_assignment_t* assignment = ((Assignment*) assignment_obj)->assignment;
1584 lp_sign_condition_t sgn_condition = PyInt_AsLong(sgn_condition_obj);
1585
1586 // Check if all but the top variable are unassigned
1587 if (!lp_polynomial_is_assigned(p, assignment)) {
1588 PyErr_SetString(PyExc_RuntimeError, "sgn_check(): All polynomial variables should be assigned by the given assignment.");
1589 return NULL;
1590 }
1591
1592 // Check the sign
1593 int sgn = lp_polynomial_sgn(p, assignment);
1594 if (lp_sign_condition_consistent(sgn_condition, sgn)) {
1595 Py_RETURN_TRUE;
1596 } else {
1597 Py_RETURN_FALSE;
1598 }
1599 }
1600
1601 static PyObject*
Polynomial_pp(PyObject * self)1602 Polynomial_pp(PyObject* self) {
1603 lp_polynomial_t* p = ((Polynomial*) self)->p;
1604 const lp_polynomial_context_t* p_ctx = lp_polynomial_get_context(p);
1605 lp_polynomial_t* pp = lp_polynomial_new(p_ctx);
1606 lp_polynomial_pp(pp, p);
1607 PyObject* pp_py = Polynomial_create(pp);
1608 return pp_py;
1609 }
1610
1611 static PyObject*
Polynomial_cont(PyObject * self)1612 Polynomial_cont(PyObject* self) {
1613 lp_polynomial_t* p = ((Polynomial*) self)->p;
1614 const lp_polynomial_context_t* p_ctx = lp_polynomial_get_context(p);
1615 lp_polynomial_t* cont = lp_polynomial_new(p_ctx);
1616 lp_polynomial_cont(cont, p);
1617 PyObject* cont_py = Polynomial_create(cont);
1618 return cont_py;
1619 }
1620
1621 static PyObject*
Polynomial_pp_cont(PyObject * self)1622 Polynomial_pp_cont(PyObject* self) {
1623 lp_polynomial_t* p = ((Polynomial*) self)->p;
1624 const lp_polynomial_context_t* p_ctx = lp_polynomial_get_context(p);
1625 lp_polynomial_t* pp = lp_polynomial_new(p_ctx);
1626 lp_polynomial_t* cont = lp_polynomial_new(p_ctx);
1627 lp_polynomial_pp_cont(pp, cont, p);
1628 PyObject* pp_py = Polynomial_create(pp);
1629 PyObject* cont_py = Polynomial_create(cont);
1630 PyObject* tuple = PyTuple_New(2);
1631 PyTuple_SetItem(tuple, 0, pp_py);
1632 PyTuple_SetItem(tuple, 1, cont_py);
1633 return tuple;
1634 }
1635