Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixed divrem bugs, added pseudo divrem #72

Merged
merged 15 commits into from
Nov 3, 2023
Merged
4 changes: 2 additions & 2 deletions include/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
# Configure the version
configure_file(version.h.in ${CMAKE_CURRENT_SOURCE_DIR}/version.h)

# Find all the base headers and add them to install terget
# Find all the base headers and add them to install target
file(GLOB headers *.h)
install(FILES ${headers} DESTINATION include/poly)

# Find all the cxx headers and add them to install terget
# Find all the cxx headers and add them to install target
file(GLOB headers_cxx polyxx/*.h)
install(FILES ${headers_cxx} DESTINATION include/poly/polyxx)
15 changes: 9 additions & 6 deletions include/polynomial.h
Original file line number Diff line number Diff line change
Expand Up @@ -88,9 +88,6 @@ void lp_polynomial_swap(lp_polynomial_t* A1, lp_polynomial_t* A2);
/** Assign the polynomial a given polynomial. */
void lp_polynomial_assign(lp_polynomial_t* A, const lp_polynomial_t* from);

/** Returns the context of the polynomial */
const lp_polynomial_context_t* lp_polynomial_context(const lp_polynomial_t* A);

/** Returns the degree of the polynomial (in it's top variable) */
size_t lp_polynomial_degree(const lp_polynomial_t* A);

Expand Down Expand Up @@ -214,14 +211,14 @@ void lp_polynomial_add_mul(lp_polynomial_t* S, const lp_polynomial_t* A1, const
void lp_polynomial_sub_mul(lp_polynomial_t* S, const lp_polynomial_t* A1, const lp_polynomial_t* A2);

/**
* Reduce the polynomial A in Z[y,x] using B in Z[y,x] so that
* Reduce the polynomial A in Z[x,y] using B in Z[x,y] so that
*
* P*A = Q*B + R
*
* and
*
* P in Z[y]
* Q, R in Z[y,x]
* P in Z[x]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This conflicts with the comment of the coefficient_reduce method. Are you sure about this?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I forgot to change this at coefficient_reduce as well.

* Q, R in Z[x,y]
*
* with
*
Expand All @@ -239,9 +236,15 @@ void lp_polynomial_rem(lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_p
/** Compute a*A1 = D*A2 + R (pseudo remainder). */
void lp_polynomial_prem(lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);

/** Compute a*A1 = D*A2 + R (pseudo remainder). */
void lp_polynomial_pdivrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);

/** Compute a*A1 = D*A2 + R (sparse pseudo remainder). */
void lp_polynomial_sprem(lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);

/** Compute a*A1 = D*A2 + R (sparse pseudo remainder). */
void lp_polynomial_spdivrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);

/** Compute A1 = D*A2 + R (assumes that exact division). */
void lp_polynomial_divrem(lp_polynomial_t* D, lp_polynomial_t* R, const lp_polynomial_t* A1, const lp_polynomial_t* A2);

Expand Down
64 changes: 54 additions & 10 deletions python/polypyPolynomial2.c
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,12 @@ Polynomial_prem(PyObject* self, PyObject* args);
static PyObject*
Polynomial_sprem(PyObject* self, PyObject* args);

static PyObject*
Polynomial_pdivrem(PyObject* self, PyObject* args);

static PyObject*
Polynomial_spdivrem(PyObject* self, PyObject* args);

static PyObject*
Polynomial_gcd(PyObject* self, PyObject* args);

Expand Down Expand Up @@ -183,6 +189,8 @@ PyMethodDef Polynomial_methods[] = {
{"rem", (PyCFunction)Polynomial_rem, METH_VARARGS, "Returns the remainder of current and given polynomial"},
{"prem", (PyCFunction)Polynomial_prem, METH_VARARGS, "Returns the pseudo remainder of current and given polynomial"},
{"sprem", (PyCFunction)Polynomial_sprem, METH_VARARGS, "Returns the sparse pseudo remainder of current and given polynomial"},
{"pdivrem", (PyCFunction)Polynomial_pdivrem, METH_VARARGS, "Returns the pseudo quotient and pseudo remainder of current and given polynomial"},
{"spdivprem", (PyCFunction)Polynomial_spdivrem, METH_VARARGS, "Returns the sparse pseudo quotient and sparse pseudo remainder of current and given polynomial"},
{"gcd", (PyCFunction)Polynomial_gcd, METH_VARARGS, "Returns the gcd of current and given polynomial"},
{"lcm", (PyCFunction)Polynomial_lcm, METH_VARARGS, "Returns the lcm of current and given polynomial"},
{"extended_gcd", (PyCFunction)Polynomial_extended_gcd, METH_VARARGS, "Returns the extended gcd, i.e. (gcd, u, v), of current and given polynomial"},
Expand Down Expand Up @@ -749,8 +757,14 @@ Polynomial_rem_operator(PyObject* self, PyObject* other) {
return Polynomial_create(rem);
}

enum rem_type {
REM_EXACT,
REM_PSEUDO,
REM_SPARSE_PSEUDO
};

static PyObject*
Polynomial_divmod(PyObject* self, PyObject* other) {
Polynomial_divmod_general(PyObject* self, PyObject* args, enum rem_type type) {

int dec_other = 0;

Expand All @@ -763,7 +777,18 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
Polynomial* p1 = (Polynomial*) self;
const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);

// Make sure other is a polynomial
// check that there is only one other
PyObject *other;
if (PyTuple_Check(args)) {
if (PyTuple_Size(args) != 1) {
Py_INCREF(Py_NotImplemented);
return Py_NotImplemented;
}
other = PyTuple_GetItem(args, 0);
} else {
other = args;
}
// make sure other is a polynomial
if (!PyPolynomial_CHECK(other)) {
if (PyVariable_CHECK(other)) {
other = PyPolynomial_FromVariable(other, p1_ctx);
Expand All @@ -785,10 +810,20 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
return Py_NotImplemented;
}

// Multiply the polynomials
// Divide the polynomials
lp_polynomial_t* rem = lp_polynomial_new(p1_ctx);
lp_polynomial_t* div = lp_polynomial_new(p1_ctx);
lp_polynomial_divrem(div, rem, p1->p, p2->p);
switch (type) {
case REM_EXACT:
lp_polynomial_divrem(div, rem, p1->p, p2->p);
break;
case REM_PSEUDO:
lp_polynomial_pdivrem(div, rem, p1->p, p2->p);
break;
case REM_SPARSE_PSEUDO:
lp_polynomial_spdivrem(div, rem, p1->p, p2->p);
break;
}

if (dec_other) {
Py_DECREF(other);
Expand All @@ -805,6 +840,21 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
return pair;
}

static PyObject*
Polynomial_divmod(PyObject* self, PyObject* other) {
return Polynomial_divmod_general(self, other, REM_EXACT);
}

static PyObject*
Polynomial_pdivrem(PyObject* self, PyObject* other) {
return Polynomial_divmod_general(self, other, REM_PSEUDO);
}

static PyObject*
Polynomial_spdivrem(PyObject* self, PyObject* other) {
return Polynomial_divmod_general(self, other, REM_SPARSE_PSEUDO);
}

static int
Polynomial_nonzero(PyObject* self) {
// Get arguments
Expand All @@ -813,12 +863,6 @@ Polynomial_nonzero(PyObject* self) {
return !lp_polynomial_is_zero(p->p);
}

enum rem_type {
REM_EXACT,
REM_PSEUDO,
REM_SPARSE_PSEUDO
};

static PyObject*
Polynomial_rem_general(PyObject* self, PyObject* args, enum rem_type type) {
int dec_other = 0;
Expand Down
87 changes: 68 additions & 19 deletions python/polypyPolynomial3.c
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,12 @@ Polynomial_prem(PyObject* self, PyObject* args);
static PyObject*
Polynomial_sprem(PyObject* self, PyObject* args);

static PyObject*
Polynomial_pdivrem(PyObject* self, PyObject* args);

static PyObject*
Polynomial_spdivrem(PyObject* self, PyObject* args);

static PyObject*
Polynomial_gcd(PyObject* self, PyObject* args);

Expand Down Expand Up @@ -180,6 +186,8 @@ PyMethodDef Polynomial_methods[] = {
{"rem", (PyCFunction)Polynomial_rem, METH_VARARGS, "Returns the remainder of current and given polynomial"},
{"prem", (PyCFunction)Polynomial_prem, METH_VARARGS, "Returns the pseudo remainder of current and given polynomial"},
{"sprem", (PyCFunction)Polynomial_sprem, METH_VARARGS, "Returns the sparse pseudo remainder of current and given polynomial"},
{"pdivrem", (PyCFunction)Polynomial_pdivrem, METH_VARARGS, "Returns the pseudo quotient and pseudo remainder of current and given polynomial"},
{"spdivprem", (PyCFunction)Polynomial_spdivrem, METH_VARARGS, "Returns the sparse pseudo quotient and sparse pseudo remainder of current and given polynomial"},
{"gcd", (PyCFunction)Polynomial_gcd, METH_VARARGS, "Returns the gcd of current and given polynomial"},
{"lcm", (PyCFunction)Polynomial_lcm, METH_VARARGS, "Returns the lcm of current and given polynomial"},
{"extended_gcd", (PyCFunction)Polynomial_extended_gcd, METH_VARARGS, "Returns the extended gcd, i.e. (gcd, u, v), of current and given polynomial"},
Expand Down Expand Up @@ -730,8 +738,14 @@ Polynomial_rem_operator(PyObject* self, PyObject* other) {
return Polynomial_create(rem);
}

enum rem_type {
REM_EXACT,
REM_PSEUDO,
REM_SPARSE_PSEUDO
};

static PyObject*
Polynomial_divmod(PyObject* self, PyObject* other) {
Polynomial_divmod_general(PyObject* self, PyObject* args, enum rem_type type) {

int dec_other = 0;

Expand All @@ -741,10 +755,21 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
}

// self is always a polynomial
Polynomial* p1 = (Polynomial*) self;
const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
Polynomial *p1 = (Polynomial *) self;
const lp_polynomial_context_t *p1_ctx = lp_polynomial_get_context(p1->p);

// Make sure other is a polynomial
// check that there is only one other
PyObject *other;
if (PyTuple_Check(args)) {
if (PyTuple_Size(args) != 1) {
Py_INCREF(Py_NotImplemented);
return Py_NotImplemented;
}
other = PyTuple_GetItem(args, 0);
} else {
other = args;
}
// make sure other is a polynomial
if (!PyPolynomial_CHECK(other)) {
if (PyVariable_CHECK(other)) {
other = PyPolynomial_FromVariable(other, p1_ctx);
Expand All @@ -759,33 +784,58 @@ Polynomial_divmod(PyObject* self, PyObject* other) {
}

// other can be a variable or a number
Polynomial* p2 = (Polynomial*) other;
const lp_polynomial_context_t* p2_ctx = lp_polynomial_get_context(p2->p);
Polynomial *p2 = (Polynomial *) other;
const lp_polynomial_context_t *p2_ctx = lp_polynomial_get_context(p2->p);
if (!lp_polynomial_context_equal(p1_ctx, p2_ctx)) {
Py_INCREF(Py_NotImplemented);
return Py_NotImplemented;
}

// Multiply the polynomials
lp_polynomial_t* rem = lp_polynomial_new(p1_ctx);
lp_polynomial_t* div = lp_polynomial_new(p1_ctx);
lp_polynomial_divrem(div, rem, p1->p, p2->p);
// Divide the polynomials
lp_polynomial_t *rem = lp_polynomial_new(p1_ctx);
lp_polynomial_t *div = lp_polynomial_new(p1_ctx);
switch (type) {
case REM_EXACT:
lp_polynomial_divrem(div, rem, p1->p, p2->p);
break;
case REM_PSEUDO:
lp_polynomial_pdivrem(div, rem, p1->p, p2->p);
break;
case REM_SPARSE_PSEUDO:
lp_polynomial_spdivrem(div, rem, p1->p, p2->p);
break;
}

if (dec_other) {
Py_DECREF(other);
}

// Return the result
PyObject* pair = PyTuple_New(2);
PyObject* divObj = Polynomial_create(div);
PyObject* remObj = Polynomial_create(rem);
PyObject *pair = PyTuple_New(2);
PyObject *divObj = Polynomial_create(div);
PyObject *remObj = Polynomial_create(rem);
Py_INCREF(divObj);
Py_INCREF(remObj);
PyTuple_SetItem(pair, 0, divObj);
PyTuple_SetItem(pair, 1, remObj);
return pair;
}

static PyObject*
Polynomial_divmod(PyObject* self, PyObject* other) {
return Polynomial_divmod_general(self, other, REM_EXACT);
}

static PyObject*
Polynomial_pdivrem(PyObject* self, PyObject* other) {
return Polynomial_divmod_general(self, other, REM_PSEUDO);
}

static PyObject*
Polynomial_spdivrem(PyObject* self, PyObject* other) {
return Polynomial_divmod_general(self, other, REM_SPARSE_PSEUDO);
}

static int
Polynomial_nonzero(PyObject* self) {
// Get arguments
Expand All @@ -794,16 +844,15 @@ Polynomial_nonzero(PyObject* self) {
return !lp_polynomial_is_zero(p->p);
}

enum rem_type {
REM_EXACT,
REM_PSEUDO,
REM_SPARSE_PSEUDO
};

static PyObject*
Polynomial_rem_general(PyObject* self, PyObject* args, enum rem_type type) {
int dec_other = 0;

if (!PyPolynomial_CHECK(self)) {
Py_INCREF(Py_NotImplemented);
return Py_NotImplemented;
}

// self is always a polynomial
Polynomial* p1 = (Polynomial*) self;
const lp_polynomial_context_t* p1_ctx = lp_polynomial_get_context(p1->p);
Expand Down
4 changes: 2 additions & 2 deletions src/number/value.c
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ int lp_value_sgn(const lp_value_t* v) {
int lp_value_cmp(const lp_value_t* v1, const lp_value_t* v2) {

if (trace_is_enabled("value::cmp")) {
tracef("lp_value_cmp()\n")
tracef("lp_value_cmp()\n");
tracef("v1 = "); lp_value_print(v1, trace_out); tracef("\n");
tracef("v2 = "); lp_value_print(v2, trace_out); tracef("\n");
}
Expand Down Expand Up @@ -612,7 +612,7 @@ char* lp_value_to_string(const lp_value_t* v) {
void lp_value_get_value_between(const lp_value_t* a, int a_strict, const lp_value_t* b, int b_strict, lp_value_t* v) {

if (trace_is_enabled("value::get_value_between")) {
tracef("lp_value_get_value_between()\n")
tracef("lp_value_get_value_between()\n");
tracef("a = "); lp_value_print(a, trace_out); tracef(", a_strict = %d\n", a_strict);
tracef("b = "); lp_value_print(b, trace_out); tracef(", b_strict = %d\n", b_strict);
}
Expand Down
Loading