From b4ac597771eafef5a245137b8bdf1ec029a31d4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 20 Oct 2019 21:06:46 +0300 Subject: [PATCH] qpq cython wrapper class, bug fixes The new cython module qpms.cypolynomials requires gmpy2 version 2.1 or higher. Former-commit-id: 10198a7341437620f5a493d16a217bf6014e2759 --- qpms/CMakeLists.txt | 1 + qpms/cypolynomials.pxd | 23 ------- qpms/cypolynomials.pyx | 141 +++++++++++++++++++++++++++++++++++++++++ qpms/polynomials.c | 4 +- setup.py | 11 +++- 5 files changed, 152 insertions(+), 28 deletions(-) delete mode 100644 qpms/cypolynomials.pxd create mode 100644 qpms/cypolynomials.pyx diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index bf65f76..825bd97 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -24,6 +24,7 @@ target_link_libraries (qpms lapack blas amos + gmp ) target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/qpms/cypolynomials.pxd b/qpms/cypolynomials.pxd deleted file mode 100644 index 4d84409..0000000 --- a/qpms/cypolynomials.pxd +++ /dev/null @@ -1,23 +0,0 @@ -cdef extern from "": - cdef struct __mpq_struct: - pass - ctypedef __mpq_struct mpq_t[1] - - cdef struct __mpz_struct: - pass - ctypedef __mpz_struct mpz_t[1] - - double mpz_get_d(const mpz_t op) - - -cdef extern from "polynomials.h": - cdef struct qpq_t: - int order - int offset - int capacity - mpq_t *coeffs - - void qpq_init(qpq_t *x, int capacity) - void qpq_extend(qpq_t *x, int new_capacity) - void qpq_set(qpq_t *copy, const qpq_t *orig) - void qpq_clear(qpq_t diff --git a/qpms/cypolynomials.pyx b/qpms/cypolynomials.pyx new file mode 100644 index 0000000..f5687ad --- /dev/null +++ b/qpms/cypolynomials.pyx @@ -0,0 +1,141 @@ +from libc.limits cimport INT_MIN, INT_MAX, UINT_MAX, LONG_MIN, LONG_MAX, ULONG_MAX +from gmpy2 cimport * # requires gmpy2>=2.1.0 or later (pip has currently 2.0.something) + +cdef extern from "gmp.h": + # cdef struct __mpq_struct: + # pass + # ctypedef __mpq_struct mpq_t[1] + # cdef struct __mpz_struct: + # pass + # ctypedef __mpz_struct mpz_t[1] + double mpz_get_d(const mpz_t op) + void mpq_init(mpq_t q) + void mpq_clear(mpq_t q) + int mpq_sgn(mpq_t q) + +cdef extern from "polynomials.h": + cdef struct qpq_t: + int order + int offset + int capacity + mpq_t *coeffs + + void qpq_init(qpq_t *p, int capacity) + void qpq_extend(qpq_t *p, int new_capacity) + void qpq_set(qpq_t *copy, const qpq_t *orig) + void qpq_set_elem(qpq_t *p, int exponent, const mpq_t coeff) + void qpq_set_elem_si(qpq_t *p, int exponent, long numerator, unsigned long denominator) + void qpq_get_elem(mpq_t coeff, const qpq_t *p, int exponent) + int qpq_get_elem_si(long *numerator, unsigned long *denominator, const qpq_t *p, int exponent) + void qpq_clear(qpq_t *p) + void qpq_add(qpq_t *sum, const qpq_t *addend1, const qpq_t *addend2) + void qpq_sub(qpq_t *difference, const qpq_t *minuend, const qpq_t *substrahend) + void qpq_mul(qpq_t *product, const qpq_t *multiplier, const qpq_t *multiplicand) + void qpq_deriv(qpq_t *dPdx, const qpq_t *P) + bint qpq_nonzero(const qpq_t *) + + +import_gmpy2() # needed to initialise the C-API + +cdef class qpq: + """ Polynomials with rational coefficients """ + cdef qpq_t p + + def __cinit__(self, *args, **kwargs): + qpq_init(&self.p, 0) + + def __init__(self, *args, **kwargs): + cdef int offset, order + cdef dict thedict + cdef mpq coeff + if len(args) > 0 and isinstance(args[0], dict) and len(args[0]) > 0: + thedict = args[0] + keys = thedict.keys() + for key in keys: + if not isinstance(key, int) or key < 0 or key > INT_MAX: + raise TypeError("Coefficient dictionary keys must be non-negative integers.") + offset = min(keys) + order = max(keys) + qpq_extend(&self.p, order - offset + 1) + self.p.order = order + self.p.offset = offset + for key, val in thedict.items(): + if MPQ_Check(val): + coeff = val + else: + coeff = mpq(val) + qpq_set_elem(&self.p, key, MPQ(coeff)) + return + + def __dealloc__(self): + qpq_clear(&self.p) + + def __getitem__(self, key): + # Only one coefficient a time supported right now + cdef mpq q + cdef mpq_t cq + if isinstance(key, int): + if key >= 0 and key <= INT_MAX: + """ Can we alternatively do it like: + q = mpq() # or GMPy_MPQ_New(NULL) ? + qpq_get_elem(MPQ(q), &self.p, key) + return q + ? """ + mpq_init(cq) + qpq_get_elem(cq, &self.p, key) + q = GMPy_MPQ_From_mpq(cq) + mpq_clear(cq) + return q + else: raise IndexError("Only non-negative int exponents (indices) allowed.") + else: raise TypeError("Only integer exponents (indices) allowed.") + + def __setitem__(self, key, value): + # Only one coefficient a time supported right now + if not MPQ_Check(value): + value = mpq(value) + if (isinstance(key, int)): + if key >= 0 and key <= INT_MAX: + qpq_set_elem(&self.p, key, MPQ(value)) + else: raise IndexError("Only non-negative int exponents (indices) allowed.") + else: raise TypeError("Only integer exponents (indices) allowed.") + + def __add__(qpq self, qpq other): + cdef qpq result = qpq() + qpq_add(&result.p, &self.p, &other.p) + return result + def __sub__(qpq self, qpq other): + cdef qpq result = qpq() + qpq_sub(&result.p, &self.p, &other.p) + return result + def __mul__(qpq self, qpq other): + cdef qpq result = qpq() + qpq_mul(&result.p, &self.p, &other.p) + return result + def derivative(self): + cdef qpq result = qpq() + qpq_deriv(&result.p, &self.p) + return result + + @property + def order(self): + return self.p.order + @property + def offset(self): + return self.p.offset + @property + def capacity(self): + return self.p.capacity + + def coeffdict(self): + """ Returns a dictionary of all non-zero coefficients """ + cdef int exponent, i + cdef dict result = dict() + cdef mpq coeff + if not qpq_nonzero(&self.p): + return result + for exponent in range(self.p.offset, self.p.order + 1): + i = exponent - self.p.offset + if mpq_sgn(self.p.coeffs[i]): + result[i] = GMPy_MPQ_From_mpq(self.p.coeffs[i]) + return result + diff --git a/qpms/polynomials.c b/qpms/polynomials.c index 0d05587..eebf3ea 100644 --- a/qpms/polynomials.c +++ b/qpms/polynomials.c @@ -116,7 +116,7 @@ int qpq_get_elem_si(long *num, unsigned long *den, const qpq_t *p, const int exp void qpq_clear(qpq_t *p) { if (p->capacity > 0) { - for (int i = p->capacity; i >= 0; --i) + for (int i = p->capacity - 1; i >= 0; --i) mpq_clear(p->coeffs[i]); free(p->coeffs); } @@ -195,7 +195,7 @@ void qpq_deriv(qpq_t *dp, const qpq_t *p) { mpq_t qi; mpq_init(qi); // qi is now 0 / 1 for(int i = p->offset + !p->offset; i <= p->order; ++i) { - mpz_set_si(mpq_numref(qi), 1); // qi is now i / 1 + mpz_set_si(mpq_numref(qi), i); // qi is now i / 1 mpq_mul(dp->coeffs[i-1], qi, p->coeffs[i]); } mpq_clear(qi); diff --git a/setup.py b/setup.py index 4ae09d2..9146933 100755 --- a/setup.py +++ b/setup.py @@ -135,16 +135,21 @@ qpms_c = Extension('qpms.qpms_c', #extra_objects = ['amos/libamos.a'], # FIXME apparently, I would like to eliminate the need to cmake/make first ) +cypolynomials = Extension('qpms.cypolynomials', + sources = ['qpms/cypolynomials.pyx'], + libraries=['qpms', 'gsl', 'lapacke', 'blas', 'gslcblas', 'pthread', 'gmp'] + ) + setup(name='qpms', version = "0.2.996", packages=['qpms'], # libraries = [('amos', {'sources': amos_sources} )], - setup_requires=['cython>=0.28',], + setup_requires=['cython>=0.28', 'gmpy2>=2.1b'], install_requires=['cython>=0.28', #'quaternion','spherical_functions', - 'scipy>=0.18.0', 'sympy>=1.2'], + 'scipy>=0.18.0', 'sympy>=1.2', 'gmpy2>=2.1b'], #dependency_links=['https://github.com/moble/quaternion/archive/v2.0.tar.gz','https://github.com/moble/spherical_functions/archive/master.zip'], - ext_modules=cythonize([qpms_c, cywaves, cytranslations, cytmatrices, cycommon, cyquaternions, cybspec, cymaterials], include_path=['qpms', 'amos'], gdb_debug=True), + ext_modules=cythonize([qpms_c, cywaves, cytranslations, cytmatrices, cycommon, cyquaternions, cybspec, cymaterials, cypolynomials], include_path=['qpms', 'amos'], gdb_debug=True), cmdclass = {'build_ext': build_ext}, zip_safe=False )