mpqs refactoring minor stuff

Former-commit-id: c17570e41047a309cc76bbe485309acfcac6e82a

CMakeLists.txt

@@ -1,4 +1,5 @@
 cmake_minimum_required(VERSION 3.0.2)
+set (CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};${CMAKE_CURRENT_SOURCE_DIR}/cmake-scripts")
 include(CMakeAddFortranSubdirectory)
 include(version.cmake)
 include(GNUInstallDirs)

cmake-scripts/COPYING-CMAKE-SCRIPTS

@@ -0,0 +1,22 @@
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. The name of the author may not be used to endorse or promote products 
+   derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

cmake-scripts/FindGMP.cmake

@@ -0,0 +1,23 @@
+# Try to find the GMP librairies
+#  GMP_FOUND - system has GMP lib
+#  GMP_INCLUDE_DIR - the GMP include directory
+#  GMP_LIBRARIES - Libraries needed to use GMP
+
+# Copyright (c) 2006, Laurent Montel, <montel@kde.org>
+#
+# Redistribution and use is allowed according to the terms of the BSD license.
+# For details see the accompanying COPYING-CMAKE-SCRIPTS file.
+
+
+if (GMP_INCLUDE_DIR AND GMP_LIBRARIES)
+  # Already in cache, be silent
+  set(GMP_FIND_QUIETLY TRUE)
+endif (GMP_INCLUDE_DIR AND GMP_LIBRARIES)
+
+find_path(GMP_INCLUDE_DIR NAMES gmp.h )
+find_library(GMP_LIBRARIES NAMES gmp libgmp)
+
+include(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(GMP DEFAULT_MSG GMP_INCLUDE_DIR GMP_LIBRARIES)
+
+mark_as_advanced(GMP_INCLUDE_DIR GMP_LIBRARIES)

cmake-scripts/ORIGINS

@@ -0,0 +1,3 @@
+The FindGMP module comes from KDE's Extra CMake Modules v5.62.0, see:
+- https://api.kde.org/ecm/
+- https://cgit.kde.org/extra-cmake-modules.git/ 

qpms/CMakeLists.txt

@@ -2,6 +2,7 @@
 find_package(GSL 2.0 REQUIRED)
 find_package(BLAS REQUIRED)
 find_package(LAPACK REQUIRED)
+find_package(GMP REQUIRED)
 
 
 #includes
@@ -12,7 +13,7 @@ add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c e
 	ewaldsf.c pointgroups.c latticegens.c
 	lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
 	bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
-	lll.c
+	lll.c polynomials.c
 	)
 use_c99()
 
@@ -23,6 +24,7 @@ target_link_libraries (qpms
 	lapack
 	blas
 	amos
+	gmp
 	)
 target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 

qpms/cypolynomials.pyx

@@ -0,0 +1,146 @@
+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_div(qpq_t *quotient, qpq_t *remainder, const qpq_t *dividend, const qpq_t *divisor)
+    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 __divmod__(qpq self, qpq other):
+        cdef qpq q = qpq(), r = qpq()
+        qpq_div(&q.p, &r.p, &self.p, &other.p)
+        return (q, r)
+    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[exponent] = GMPy_MPQ_From_mpq(self.p.coeffs[i])
+        return result
+

qpms/polynomials.c

@@ -0,0 +1,546 @@
+#include "polynomials.h"
+#include <stdlib.h>
+#include "qpms_error.h"
+#include <stdbool.h>
+#include <stdio.h>
+
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#define MIN(x, y) (((x) <= (y)) ? (x) : (y))
+
+#if 0
+void mpzs_init(mpzs_t x) {
+  mpz_init(x->_1);
+  mpz_init(x->_2);
+}
+
+void mpzs_clear(mpzs_t x) {
+  mpz_clear(x->_1);
+  mpz_clear(x->_2);
+}
+
+void mpzs_set(mpzs_t x, const mpzs_t y) {
+  mpz_set(x->_1, y->_1);
+  mpz_set(x->_2, y->_2);
+}
+
+// Compares the square-rooted part of mpzs_t, so we can use it as a key in a search tree.
+// Probably deprecated, since we now use hashes.
+static int mpzs_cmp2(void *op1, void *op2) {
+  mpz_cmpabs(((mpzs_t) op1)->_2, ((mpzs_t) op2)->_2);
+}
+#endif
+
+/// Type representing a number of form \f$ a \sqrt{b}; a \in \ints, b \in \nats \f$.
+/** This includes UT hash handle structure */
+typedef struct _qp_mpzs_hashed {
+        mpz_t _1; ///< The integer factor \f$ a \f$.
+        mpz_t _2; ///< The square-rooted factor \f$ b \f$. Always positive.
+        UT_hash_handle hh;
+} mpzs_hh_t[1];
+
+void mpzs_hh_init(mpzs_hh_t x);
+void mpzs_hh_clear(mpzs_hh_t x);
+void mpzs_hh_set(mpzs_hh_t x, const mpzs_hh_t y);
+// void mpzs_hh_set_ui(mpzs_hh_t x, long integer_factor, unsigned long sqrt_part);
+
+
+void mpzs_hh_init(mpzs_hh_t x) {
+  mpz_init(x->_1);
+  mpz_init(x->_2);
+}
+
+void mpzs_hh_clear(mpzs_hh_t x) {
+  mpz_clear(x->_1);
+  mpz_clear(x->_2);
+}
+
+void mpzs_hh_set(mpzs_hh_t x, const mpzs_hh_t y) {
+  mpz_set(x->_1, y->_1);
+  mpz_set(x->_2, y->_2);
+}
+
+static inline void mpzs_hash_clear(struct _qp_mpzs_hashed **hash) {
+  struct _qp_mpzs_hashed *current, *tmp;
+  HASH_ITER(hh, *hash, current, tmp) {
+    HASH_DEL(*hash, current);
+    free(current);
+    //x->sz--;
+  }
+}
+
+/* Append an mpzs element to a mpzs sum that must not contain a matching key
+ * in advance.
+ */
+static inline void mpzs_hash_append(struct _qp_mpzs_hashed **hash, 
+    const struct _qp_mpzs_hashed *newelem) {
+  struct _qp_mpzs_hashed *n;
+  QPMS_CRASHING_MALLOC(n, sizeof(mpzs_hh_t));
+  mpzs_hh_init(n);
+  mpzs_hh_set(n, newelem);
+  HASH_ADD_KEYPTR(hh, *hash, mpz_limbs_read(n->_2),
+      mpz_size(n->_2) * sizeof(mp_limb_t), n);
+}
+
+// Arithmetically add an mpzs element to a mpzs sum/hash.
+static inline void mpzs_hash_addelem(struct _qp_mpzs_hashed **hash,
+    const struct _qp_mpzs_hashed *addend) {
+  struct _qp_mpzs_hashed *s;
+  HASH_FIND(hh, *hash, mpz_limbs_read(addend->_2),
+      mpz_size(addend->_2), s);
+  if (!s) mpzs_hash_append(hash, addend); // if not found
+  else {
+    mpz_add(s->_1, s->_1, addend->_1);
+    if(!mpz_sgn(s->_1)) { // If zero, annihilate
+      HASH_DEL(*hash, s);
+      mpzs_hh_clear(s);
+      free(s);
+    }
+  }
+}
+
+// Multiplies two mpzs hashes, adding them to *prod; 
+// *prod must differ from x and y (this is not checked)
+static inline void  mpzs_hash_mul(struct _qp_mpzs_hashed **prod, 
+    const struct _qp_mpzs_hashed *x, const struct _qp_mpzs_hashed *y) {
+  mpz_t gcd; mpz_init(gcd); // common denominator of the sqrt to factor out
+  mpz_t mx, my; mpz_init(mx); mpz_init(my);
+  mpzs_hh_t addend; mpzs_hh_init(addend);
+
+  for (const struct _qp_mpzs_hashed *nx = x; 
+      nx != NULL; nx = nx->hh.next) {
+    for (const struct _qp_mpzs_hashed *ny = y; 
+        ny != NULL; ny = ny->hh.next) {
+      mpz_gcd(gcd, nx->_2, ny->_2);
+      mpz_divexact(mx, nx->_2, gcd);
+      mpz_divexact(my, ny->_2, gcd);
+      mpz_mul(addend->_2, mx, my);
+      mpz_mul(addend->_1, ny->_1, nx->_1);
+      mpz_mul(addend->_1, addend->_1, gcd);
+      mpzs_hash_addelem(prod, addend);
+    }
+  }
+  mpzs_hh_clear(addend);
+  mpz_clear(mx); mpz_clear(my); mpz_clear(gcd);
+}
+
+
+//===== mpqs_t =====
+
+void mpqs_init(mpqs_t x) {
+  //x->sz = 0;
+  x->nt = 0;
+  mpq_init(x->f);
+  mpq_set_ui(x->f, 1, 1);
+}
+
+void mpqs_clear_num(struct _qp_mpqs *x) {
+  mpzs_hash_clear(&(x->nt));
+}
+
+void mpqs_clear(mpqs_t x) {
+  mpqs_clear_num(x);
+  mpq_clear(x->f);
+  //x->sz = 0;
+  x->nt = 0;
+}
+
+// Append an element to mpqs_t's numerator hash without any checks.
+void mpqs_nt_append(mpqs_t x, const struct _qp_mpzs_hashed *numelem) {
+  mpzs_hash_append(&(x->nt), numelem);
+}
+
+void mpqs_set_z(mpqs_t x, const mpz_t num1, const mpz_t num2,
+    const mpz_t den) {
+  mpqs_clear_num(x);
+  if (mpz_sgn(num1) == 0 || mpz_sgn(num2) == 0) return;
+  mpz_abs(mpq_numref(x->f), num1);
+  mpz_abs(mpq_denref(x->f), den);
+  mpq_canonicalize(x->f);
+  mpzs_hh_t tmp; mpzs_hh_init(tmp);
+  mpz_set_si(tmp->_1, mpz_sgn(num1) * mpz_sgn(den));
+  mpz_set(tmp->_2, num2);
+  mpqs_nt_append(x, tmp);
+  mpzs_hh_clear(tmp);
+}
+
+void mpqs_set_si(mpqs_t x, long num1, unsigned long num2,
+    unsigned long den) {
+  mpqs_clear_num(x);
+  if(num1 == 0 || num2 == 0) return;
+  mpz_set_ui(mpq_numref(x->f), labs(num1));
+  mpz_set_ui(mpq_denref(x->f), den);
+  mpq_canonicalize(x->f);
+  mpzs_hh_t tmp; mpzs_hh_init(tmp);
+  mpz_set_si(tmp->_1, (num1 < 0) ? -1 : 1);
+  mpz_set_ui(tmp->_2, num2);
+  mpqs_nt_append(x, tmp);
+  mpzs_hh_clear(tmp);
+}
+
+// Arithmetically adds an element to mpqs's numerator hash.
+void mpqs_nt_addelem(mpqs_t x, const struct _qp_mpzs_hashed *addend) {
+  mpzs_hash_addelem(&(x->nt), addend);
+}
+
+void mpqs_init_set(mpqs_t dest, const mpqs_t src) {
+  mpqs_init(dest);
+  mpq_set(dest->f, src->f);
+  struct _qp_mpzs_hashed *numitem, *tmp;
+  HASH_ITER(hh, src->nt, numitem, tmp) {
+    mpqs_nt_append(dest, numitem);
+  }
+}
+
+void mpqs_set(mpqs_t x, const mpqs_t y) {
+  mpqs_clear(x);
+  mpqs_init_set(x, y);
+}
+
+void mpqs_neg_inplace(mpqs_t x) {
+  for(struct _qp_mpzs_hashed *n = (x->nt)->hh.next; n != NULL; n = n->hh.next) {
+    mpz_neg(n->_1, n->_1);
+  }
+}
+
+void mpqs_neg(mpqs_t negated_operand, const mpqs_t operand) {
+  mpqs_set(negated_operand, operand);
+  mpqs_neg_inplace(negated_operand);
+}
+
+void mpqs_nt_gcd(mpz_t gcd, const mpqs_t x) {
+  if(x->nt) {
+    mpz_set(gcd, x->nt->_1);
+    for(struct _qp_mpzs_hashed *n = (x->nt)->hh.next; n != NULL; n = n->hh.next) {
+      mpz_gcd(gcd, gcd, n->_1);
+    }
+  } else 
+    mpz_set_ui(gcd, 1);
+}
+
+void mpqs_canonicalise(mpqs_t x) {
+  mpz_t gcd; mpz_init(gcd);
+  mpqs_nt_gcd(gcd, x);
+  mpz_mul(mpq_numref(x->f), mpq_numref(x->f), gcd);
+  for(struct _qp_mpzs_hashed *n = (x->nt); n != NULL; n = n->hh.next) {
+    mpz_divexact(n->_1, n->_1, gcd);
+  }
+  mpz_clear(gcd);
+  mpq_canonicalize(x->f);
+}
+
+void mpqs_add(mpqs_t sum_final, const mpqs_t x, const mpqs_t y) {
+  mpqs_t sum; mpqs_init(sum);
+  mpqs_set(sum, x);
+
+  // Denominators gcd
+  mpz_t den_gcd; mpz_init(den_gcd);
+  mpz_gcd(den_gcd, mpq_denref(x->f), mpq_denref(y->f));
+  // Common denominator
+  mpz_divexact(mpq_denref(sum->f), mpq_denref(sum->f), den_gcd);
+  mpz_mul(mpq_denref(sum->f), mpq_denref(sum->f), mpq_denref(y->f));
+  
+  // Left operand numerator factor
+  mpz_t tmp; mpz_init(tmp);
+  mpz_divexact(tmp, mpq_denref(y->f), den_gcd);
+  mpz_mul(tmp, tmp, mpq_numref(x->f));
+  /* Distribute the factor to numerator elements; this approach might be
+   * be suboptimal if x is a complicated expression
+   * and y is quite simple. Maybe optimise later 
+   */
+  for(struct _qp_mpzs_hashed *n = sum->nt; n != NULL; n = n->hh.next)
+    mpz_mul(n->_1, n->_1, tmp);
+  mpz_set_ui(mpq_numref(sum->f), 1);
+  
+  // At this point, sum is equal to x but in a form prepared to
+  // take summands from y.
+
+  // Right operand numerator factor
+  mpz_divexact(tmp, mpq_denref(x->f), den_gcd);
+  mpz_mul(tmp, tmp, mpq_numref(y->f));
+
+  mpzs_hh_t addend; mpzs_hh_init(addend);
+  for(const struct _qp_mpzs_hashed *n = y->nt; n != NULL; n = n->hh.next) {
+    mpz_mul(addend->_1, tmp, n->_1);
+    mpz_set(addend->_2, n->_2);
+    mpqs_nt_addelem(sum, addend);
+  }
+  mpzs_hh_clear(addend);
+  mpz_clear(tmp);
+  mpz_clear(den_gcd);
+  
+  mpqs_canonicalise(sum);
+  mpqs_set(sum_final, sum);
+  mpqs_clear(sum);
+}
+
+void mpqs_sub(mpqs_t dif, const mpqs_t x, const mpqs_t y) {
+  // This is probably not optimal
+  mpqs_t tmp; mpqs_init(tmp);
+  mpqs_neg(tmp, y);
+  mpqs_add(dif, x, tmp);
+  mpqs_clear(tmp);
+}
+
+
+void mpqs_mul(mpqs_t product_final, const mpqs_t x, const mpqs_t y) {
+  mpqs_t prod; mpqs_init(prod);
+
+  mpz_mul(mpq_numref(prod->f), mpq_numref(x->f), mpq_numref(y->f));
+  mpz_mul(mpq_denref(prod->f), mpq_denref(x->f), mpq_denref(y->f));
+
+  mpzs_hash_mul(&(prod->nt), x->nt, y->nt);
+
+  mpqs_canonicalise(prod);
+  mpqs_set(product_final, prod);
+  mpqs_clear(prod);
+}
+
+
+// Auxillary function to set a mpq_t to 0/1
+static inline void mpq_zero(mpq_t q) {
+// Maybe not the best way to set it to zero.
+// Alternatively, we could use just mpz_set_si(mpq_numref(sum->coeffs[i - minoffset]), 0);
+  mpq_clear(q);
+  mpq_init(q);
+}
+
+
+// TODO to the template
+
+static inline void qpq_zero(qpq_t *q) {
+  qpq_clear(q);
+}
+
+
+#define TEMPLATE_QPQ
+#include "polynomials.template"
+#undef TEMPLATE_QPQ
+
+#if 0
+// Used in the template but perhaps not too useful to put in the header.
+static inline void mpqs_set_si(mpqs_t q, int num, unsigned den) {mpq_set_si(q->_2, num, den) ;}
+
+// TODO Put into the template
+static inline void mpqs_zero(mpqs_t q) {
+  mpqs_clear(q);
+  mpqs_init(q);
+}
+
+static inline void qpqs_zero(qpqs_t *q) {
+  qpqs_clear(q);
+}
+
+#define TEMPLATE_QPQS
+#include "polynomials.template"
+#undef TEMPLATE_QPQS
+
+
+static void qpq_dbgprint(const qpq_t *p) {
+  for(int n = p->order; n >= p->offset; --n)
+    if(mpq_sgn(p->coeffs[n - p->offset]))
+      gmp_printf("%+Qdx**%d  ", p->coeffs[n - p->offset], n);
+  gmp_printf("[%d, %d, %d]\n", p->capacity, p->order, p->offset);
+}
+#endif
+
+void qpq_set_elem_si(qpq_t *p, const int exponent, const long int num, const unsigned long int den) {
+  mpq_t q;
+  mpq_init(q);
+  mpq_set_si(q, num, den);
+  qpq_set_elem(p, exponent, q);
+  mpq_clear(q);
+}
+
+int qpq_get_elem_si(long *num, unsigned long *den, const qpq_t *p, const int exponent) {
+  mpq_t q;
+  mpq_init(q);
+  qpq_get_elem(q, p, exponent);
+  int retval = 0;
+  *num = mpz_get_si(mpq_numref(q));
+  if (!mpz_fits_slong_p(mpq_numref(q)))
+    retval += 1;
+  *den = mpz_get_ui(mpq_denref(q));
+  if (!mpz_fits_ulong_p(mpq_denref(q)))
+    retval += 2;
+  mpq_clear(q);
+  return retval;
+}
+
+
+void qpq_add(qpq_t *sum, const qpq_t *x, const qpq_t *y) {
+  const int maxorder = MAX(x->order, y->order);
+  const int minoffset = MIN(x->offset, y->offset);
+  qpq_extend(sum, maxorder - minoffset + 1);
+  /* if sum is actually some of the summands and that summand has higher
+   * offset, we have to lower the offset.
+   */
+  if ((sum == x || sum == y) && sum->offset > minoffset) 
+    qpq_lower_offset(sum, sum->offset - minoffset);
+
+  for (int i = minoffset; i <= maxorder; ++i) {
+    if (i - x->offset >= 0 && i <= x->order) {
+      if (i - y->offset >= 0 && i <= y->order) 
+        mpq_add(sum->coeffs[i - minoffset], 
+        x->coeffs[i - x->offset], y->coeffs[i - y->offset]);
+      else
+        mpq_set(sum->coeffs[i - minoffset], x->coeffs[i - x->offset]);
+    } else {
+      if (i - y->offset >= 0 && i <= y->order)
+        mpq_set(sum->coeffs[i - minoffset], y->coeffs[i - y->offset]);
+      else {
+        mpq_zero(sum->coeffs[i - minoffset]);
+      }
+    }
+  }
+  sum->offset = minoffset;
+  sum->order = maxorder;
+}
+      
+void qpq_sub(qpq_t *dif, const qpq_t *x, const qpq_t *y) {
+  const int maxorder = MAX(x->order, y->order);
+  const int minoffset = MIN(x->offset, y->offset);
+  qpq_extend(dif, maxorder - minoffset + 1);
+  /* if dif is actually some of the summands and that summand has higher
+   * offset, we have to lower the offset.
+   */
+  if ((dif == x || dif == y) && dif->offset > minoffset) 
+    qpq_lower_offset(dif, dif->offset - minoffset);
+
+  for (int i = minoffset; i <= maxorder; ++i) {
+    if (i - x->offset >= 0 && i <= x->order) {
+      if (i - y->offset >= 0 && i <= y->order) 
+        mpq_sub(dif->coeffs[i - minoffset], 
+        x->coeffs[i - x->offset], y->coeffs[i - y->offset]);
+      else
+        mpq_set(dif->coeffs[i - minoffset], x->coeffs[i - x->offset]);
+    } else {
+      if (i - y->offset >= 0 && i <= y->order) {
+        mpq_set(dif->coeffs[i - minoffset], y->coeffs[i - y->offset]);
+        mpq_neg(dif->coeffs[i - minoffset], dif->coeffs[i - minoffset]);
+      } else {
+        mpq_zero(dif->coeffs[i - minoffset]);
+      }
+    }
+  }
+  dif->offset = minoffset;
+  dif->order = maxorder;
+}
+
+void qpq_mul(qpq_t *p_orig, const qpq_t *x, const qpq_t *y) {
+
+  const int maxorder = x->order + y->order;
+  const int minoffset = x->offset + y->offset;
+  
+  qpq_t p[1];
+  qpq_init(p, maxorder - minoffset + 1);
+  mpq_t tmp; mpq_init(tmp);
+  for (int xi = x->offset; xi <= x->order; ++xi)
+    for (int yi = y->offset; yi <= y->order; ++yi) {
+      mpq_mul(tmp, x->coeffs[xi - x->offset], y->coeffs[yi - y->offset]);
+      mpq_add(p->coeffs[xi + yi - minoffset], p->coeffs[xi + yi - minoffset], tmp);
+    }
+  mpq_clear(tmp);
+  p->order = maxorder;
+  p->offset = minoffset;
+  qpq_set(p_orig, p);
+  qpq_clear(p);
+} 
+
+void qpq_div(qpq_t *q_orig, qpq_t *r_orig, const qpq_t *dend, const qpq_t *dor) {
+  QPMS_ENSURE(q_orig != r_orig, 
+      "Quotient and remainder must be _different_ instances of qpq_t.");
+
+
+  // Split the divisor into "head" and "tail"
+  qpq_t dor_tail[1]; qpq_init(dor_tail, dor->order - dor->offset); // divisor tail
+  qpq_set(dor_tail, dor);
+  qpq_canonicalise(dor_tail);
+  QPMS_ENSURE(qpq_nonzero(dor_tail), "The divisor must be non-zero");
+
+  const int dor_order = dor_tail->order;
+  mpq_t dor_head; mpq_init(dor_head);
+  mpq_set(dor_head, dor_tail->coeffs[dor_order - dor_tail->offset]);
+  dor_tail->order--;
+
+  qpq_t q[1];  qpq_init(q, dend->order - dor_order + 1);
+  q->order = dend->order - dor_order;
+  q->offset = 0;
+
+  // Assign the dividend to r but with extended capacity (with zero offset)
+  qpq_t r[1]; qpq_init(r, dend->order + 1);
+  r->offset = 0;
+  r->order = dend->order;
+  for (int n = dend->offset; n <= dend->order; ++n) 
+    mpq_set(r->coeffs[n], dend->coeffs[n - dend->offset]);
+
+  qpq_t f[1]; qpq_init(f, 1);
+  qpq_t ftail[1]; qpq_init(ftail, dor_tail->order - dor_tail->offset + 1);
+  for(; r->order >= dor_order; --r->order) {
+    // Compute the current order (r->order - dor_order) of q
+    mpq_t * const hicoeff = &(q->coeffs[r->order - dor_order]);
+    mpq_div(*hicoeff, r->coeffs[r->order], dor_head);
+    mpq_canonicalize(*hicoeff);
+    
+    // Update the remainder
+    f->offset = f->order = r->order - dor_order;
+    mpq_set(f->coeffs[0], *hicoeff);
+    qpq_mul(ftail, f, dor_tail);
+    qpq_sub(r, r, ftail);
+  }
+  qpq_clear(ftail);
+  qpq_clear(f);
+  mpq_clear(dor_head);
+  qpq_clear(dor_tail);
+
+  qpq_canonicalise(r);
+  qpq_canonicalise(q);
+  qpq_set(q_orig, q);
+  qpq_set(r_orig, r);
+  qpq_clear(r); qpq_clear(f);
+}
+
+
+void qpq_deriv(qpq_t *dp, const qpq_t *p) {
+  if (p->order <= 0) { // p is constant, dp is zero.
+    qpq_clear(dp);
+    return;
+  }
+
+  qpq_extend(dp, p->order - p->offset + (p->offset > 0));
+
+  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), i); // qi is now i / 1
+    mpq_mul(dp->coeffs[i-1], qi, p->coeffs[i]);
+  }
+  mpq_clear(qi);
+  dp->order = p->order - 1;
+  dp->offset = p->offset - 1 + !p->offset;
+}
+
+
+void qpq_legendroid_init(qpq_legendroid_t *p) {
+  qpq_init(&p->p, 0);
+  p->f = 0;
+}
+
+void qpq_legendroid_clear(qpq_legendroid_t *p) {
+  qpq_clear(&p->p);
+  p->f = 0;
+}
+
+void qpq_legendroid_mul(qpq_legendroid_t *p, const qpq_legendroid_t *a, const qpq_legendroid_t *b) {
+  qpq_mul(&p->p, &a->p, &b->p);
+  if (a->f && b->f) {
+    // TODO make somehow a static representation of this constant polynomial
+    qpq_t ff;
+    qpq_init(&ff, 3);
+    qpq_set_elem_si(&ff, 0, -1, 1);
+    qpq_set_elem_si(&ff, 2, 1, 1);
+    qpq_mul(&p->p, &p->p, &ff);
+    qpq_clear(&ff);
+  }
+  p->f = !(a->f) != !(b->f);
+}
+

qpms/polynomials.h

@@ -0,0 +1,179 @@
+/** \file polynomials.h
+ * \brief Basic operations with polynomials.
+ *
+ */
+#ifndef QPMS_POLYNOMIALS_H
+#include <gmp.h>
+#include "uthash.h"
+
+/// Polynomial with rational coeffs.
+// TODO more docs about initialisation etc.
+typedef struct qpq_t {
+	int order;
+	int offset;
+	int capacity;
+	mpq_t *coeffs;
+} qpq_t;
+
+const static qpq_t QPQ_ZERO = {-1, 0, 0, NULL};
+
+/// Initiasise the coefficients array in qpq_t.
+/** Do not use on qpq_t that has already been initialised
+ * (and not recently cleared),
+ * otherwise you can get a memory leak.
+ */
+void qpq_init(qpq_t *p, int capacity);
+
+/// Extend capacity of a qpq_t instance.
+/** If the requested new_capacity is larger than the qpq_t's
+ * capacity, the latter is extended to match new_capacity.
+ * Otherwise, nothing happend (this function does _not_ trim
+ * the capacity).
+ */
+void qpq_extend(qpq_t *p, int new_capacity);
+
+/// Shrinks the capacity to the minimum that can store the current polynomial.
+void qpq_shrink(qpq_t *p);
+
+/// Canonicalises the coefficients and (re)sets the correct degree.
+void qpq_canonicalise(qpq_t *p);
+
+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);
+/** \returns zero if the result fits into long / unsigned long; non-zero otherwise. */
+int qpq_get_elem_si(long *numerator, unsigned long *denominator, const qpq_t *p, int exponent);
+
+/// Deinitialise the coefficients array in qpq_t.
+void qpq_clear(qpq_t *p);
+
+/// Polynomial addition.
+/** Supports operand and result pointer mixing. */
+void qpq_add(qpq_t *sum, const qpq_t *addend1, const qpq_t *addend2);
+
+/// Polynomial substraction.
+/** Supports operand and result pointer mixing. */
+void qpq_sub(qpq_t *difference, const qpq_t *minuend, const qpq_t *substrahend);
+
+/// Polynomial multiplication.
+/** Supports operand and result pointer mixing. */
+void qpq_mul(qpq_t *product, const qpq_t *multiplier, const qpq_t *multiplicand);
+
+/// Polynomial division with remainder.
+/** Supports operand and result pointer mixing. */
+void qpq_div(qpq_t *quotient, qpq_t *remainder, const qpq_t *dividend, const qpq_t *divisor);
+
+/// Polynomial derivative.
+/** Supports operand and result pointer mixing. */
+void qpq_deriv(qpq_t *dPdx, const qpq_t *P);
+
+/// Tests whether a polynomial is non-zero.
+_Bool qpq_nonzero(const qpq_t *);
+
+
+
+#if 0
+/// Type representing a number of form \f$ a \sqrt{b}; a \in \ints, b \in \nats \f$.
+typedef struct _qp_mpzs {
+	mpz_t _1; ///< The integer factor \f$ a \f$.
+	mpz_t _2; ///< The square-rooted factor \f$ b \f$. Always positive.
+} mpzs_t[1];
+
+void mpzs_init(mpzs_t x);
+void mpzs_clear(mpzs_t x);
+void mpzs_set(mpzs_t x, const mpzs_t y);
+#endif
+
+
+struct _qp_mpzs_hashed;
+
+/// Sum of square roots of rational numbers.
+/// Represented as \f$ \sum_s a_i \sqrt{b_i} / d \f$.
+typedef struct _qp_mpqs {
+	//int sz; ///< Used size of the numtree.
+	struct _qp_mpzs_hashed *nt; ///< List of numerator components..
+	mpq_t f; ///< Common rational factor. Always positive.
+} mpqs_t[1];
+
+void mpqs_init(mpqs_t x);
+void mpqs_clear(mpqs_t x);
+void mpqs_set(mpqs_t x, const mpqs_t y);
+void mpqs_set_z(mpqs_t x, const mpz_t numerator_int, 
+		const mpz_t numerator_sqrt, const mpz_t denominator);
+void mpqs_set_si(mpqs_t x, long numerator_int, 
+		unsigned long numerator_sqrt, unsigned long denominator);
+void mpqs_add(mpqs_t sum, const mpqs_t addend1, const mpqs_t addend2);
+void mpqs_neg(mpqs_t negated_operand, const mpqs_t operand);
+void mpqs_sub(mpqs_t difference, const mpqs_t minuend, const mpqs_t substrahend);
+void mpqs_mul(mpqs_t product, const mpqs_t multiplier, const mpqs_t multiplicand);
+//void mpqs_div_zs(mpqs_t quotient, const mpqs_t dividend, const mpzs_hh_t divisor);
+void mpqs_clear_num(mpqs_t x); ///< Sets the numerator to zero, clearing the numerator tree.
+
+
+/// A type representing a polynomial with rational coefficients times an optional factor \f$ \sqrt{1-x^2} \f$.
+typedef struct qpq_legendroid_t {
+	qpq_t p;
+	_Bool f;
+} qpq_legendroid_t;
+
+void qpq_legendroid_init(qpq_legendroid_t *p);
+void qpq_legendroid_clear(qpq_legendroid_t *p);
+
+/// Polynomial multiplication.
+void qpq_legendroid_mul(qpq_legendroid_t *product, const qpq_legendroid_t *multiplier, const qpq_legendroid_t *multiplicand);
+
+/// Polynomial derivative.
+void qpq_legendroid_deriv(qpq_legendroid_t *dP_dx, const qpq_legendroid_t *P);
+
+
+/// Polynomial with double coeffs.
+typedef struct qpz_t {
+	int order;
+	int offset;
+	int capacity;
+	double *coeffs;
+} qpz_t;
+
+/// Initiasise the coefficients array in qpz_t.
+void qpz_init(qpz_t *p, int maxorder);
+
+/// Deinitialise the coefficients array in qpz_t.
+void qpz_clear(qpz_t *p);
+
+/// Polynomial addition.
+void qpz_add(qpz_t *sum, const qpz_t *addend1, const qpz_t *addend2);
+
+/// Polynomial substraction.
+void qpz_sub(qpz_t *difference, const qpz_t *minuend, const qpz_t *substrahend);
+
+/// Polynomial multiplication.
+void qpz_mul(qpz_t product, const qpz_t *multiplier, const qpz_t *multiplicand);
+
+/// Convert rational coefficient polynomial to double coefficient polynomial
+void qpz_from_qpq(qpz_t *target, const qpq_t *src);
+
+
+#if 0 // This will go elsewhere
+/// Table with pre-calculated Ferrers function coeffs.
+typedef struct qpms_legendre_table_t {
+	qpms_l_t lMax;
+	// TODO
+} qpms_legendre_table_t;
+
+/// Constructor for qpms_legendre_table_t.
+qpms_legendre_table_t *qpms_legendre_table_init(qpms_l_t lMax);
+
+/// Destructor for qpms_legendre_table_t.
+void qpms_legendre_table_free(qpms_legendre_table_t *);
+
+/// Evaluates a Ferrers function.
+double qpms_legendre_table_eval(const qpms_legendre_table_t *,
+		qpms_l_t l, qpms_m_t m, double x);
+
+
+// TODO pre-calculate also the products??
+#endif
+
+#endif //QPMS_POLYNOMIALS_H

qpms/polynomials.template

@@ -0,0 +1,159 @@
+// C preprocessor sorcery inspired by gsl/specfunc/legendre_source.c
+
+
+
+#if defined(TEMPLATE_QPQ)
+#define POLYTYPE qpq_t
+#define COEFTYPE mpq_t
+#define POLYFUNC(x) qpq_ ## x
+#define COEFFUNC(x) mpq_ ## x
+#define POLYTYPE_ZERO QPQ_ZERO
+
+#elif defined(TEMPLATE_QPQS)
+#define POLYTYPE qpqs_t
+#define COEFTYPE mpqs_t
+#define POLYFUNC(x) qpqs_ ## x
+#define COEFFUNC(x) mpqs_ ## x
+#define POLYTYPE_ZERO QPQS_ZERO
+
+#endif
+
+
+// POLYTYPE internal consistency check
+static inline void POLYFUNC(cc)(const POLYTYPE *p) {
+  if (!p->coeffs) return;
+  QPMS_ENSURE(p->capacity >= p->order - p->offset + 1, "POLYTYPE inconsistency detected; have you initialised it properly?");
+}
+
+_Bool POLYFUNC(nonzero)(const POLYTYPE *p) {
+  POLYFUNC(cc)(p);
+  if (p->capacity == 0) return false;
+
+  for(int i = 0; i <= p->order - p->offset; ++i)
+    if (COEFFUNC(sgn)(p->coeffs[i]))
+      return true;
+  return false;
+}
+
+void POLYFUNC(init)(POLYTYPE *p, int capacity) {
+  *p = POLYTYPE_ZERO;
+  if (capacity <= 0)
+    return;
+  QPMS_CRASHING_MALLOC(p->coeffs, capacity * sizeof(COEFTYPE));
+  for(int i = 0; i < capacity; ++i)
+    COEFFUNC(init)(p->coeffs[i]);
+  p->capacity = capacity;
+}
+
+void POLYFUNC(extend)(POLYTYPE *p, int cap) {
+  QPMS_ENSURE(p->capacity >= 0, "Got polynomial with negative capacity (%d). Is this a manually allocated one?", p->capacity);
+  if (cap > 0 && cap > p->capacity) {
+    QPMS_CRASHING_REALLOC(p->coeffs, sizeof(COEFTYPE) * cap);
+    for(int i = p->capacity; i < cap; ++i)
+      COEFFUNC(init)(p->coeffs[i]);
+    p->capacity = cap;
+  }
+}
+
+void POLYFUNC(shrink)(POLYTYPE *p) {
+  if (p->capacity > 0) {
+      for(int n = p->capacity - 1; n > p->order - p->offset; --n)
+        COEFFUNC(clear)(p->coeffs[n]);
+      p->capacity = p->order - p->offset + 1;
+      if (p->capacity > 0)
+        QPMS_CRASHING_REALLOC(p->coeffs, p->capacity * sizeof(COEFTYPE));
+  }
+}
+
+
+void POLYFUNC(canonicalise)(POLYTYPE *p) {
+  POLYFUNC(cc)(p);
+  // Lower the order if necessary (here one can get -1 if the polynomial is 0)
+  while (p->order >= p->offset && !COEFFUNC(sgn)(p->coeffs[p->order - p->offset])) --p->order;
+  if (p->order < p->offset)
+    POLYFUNC(zero)(p);
+  else { // remove the lowest-order coefficients which are in fact zero.
+    int i = 0;
+    while (i <= p->order - p->offset && !COEFFUNC(sgn)(p->coeffs[i])) ++i;
+    if (i > 0)
+      for (int j = 0; j <= p->order - p->offset - i; ++j)
+        COEFFUNC(swap)(p->coeffs[j], p->coeffs[j+i]);
+    p->offset += i;
+    // canonicalise the fractions
+    for (i = 0; i <= p->order - p->offset; ++i)
+      COEFFUNC(canonicalize)(p->coeffs[i]);
+  }
+}
+
+void POLYFUNC(set)(POLYTYPE *p, const POLYTYPE *orig) {
+  if(p != orig) {
+    const int order = orig->order, offset = orig->offset;
+    POLYFUNC(extend)(p, order - offset + 1);
+    for (int i = orig->offset; i <= order; ++i)
+      COEFFUNC(set)(p->coeffs[i - offset], orig->coeffs[i - offset]);
+    p->offset = offset;
+    p->order = order;
+  }
+}
+
+void POLYFUNC(set_elem)(POLYTYPE *p, const int exponent, const COEFTYPE coeff)  {
+  QPMS_ENSURE(exponent >= 0, "Exponent must be non-negative, got %d", exponent);
+  int offset = p->offset, order = p->order;
+  if(COEFFUNC(sgn)(coeff) == 0 && (exponent < offset || exponent > order))
+    return; // exponent out of range, but zero needs not to be assigned explicitly
+  if(exponent < p->offset) {
+    POLYFUNC(extend)(p, p->order - exponent + 1);
+    offset = exponent;
+    for(int i = order - offset; i >= p->offset - offset; --i)
+      COEFFUNC(swap)(p->coeffs[i], p->coeffs[i - (p->offset - offset)]);
+    for(int i = p->offset - offset - 1; i > 0; --i)
+      COEFFUNC(zero)(p->coeffs[i]);
+    p->offset = offset;
+  }
+  if(exponent > order) {
+    POLYFUNC(extend)(p, exponent - p->offset + 1);
+    for(int i = p->order - p->offset + 1; i <= exponent - p->offset; ++i)
+      COEFFUNC(zero)(p->coeffs[i]);
+    p->order = exponent;
+  }
+  COEFFUNC(set)(p->coeffs[exponent - p->offset], coeff);
+  return;
+}
+
+
+void POLYFUNC(get_elem)(COEFTYPE coeff, const POLYTYPE *p, const int exponent) {
+  if (exponent < p->offset || exponent > p->order)
+    COEFFUNC(zero)(coeff);
+  else
+    COEFFUNC(set)(coeff, p->coeffs[exponent-p->offset]);
+}
+
+
+void POLYFUNC(clear)(POLYTYPE *p) {
+  if (p->capacity > 0) {
+    for (int i = p->capacity - 1; i >= 0; --i)
+      COEFFUNC(clear)(p->coeffs[i]);
+    free(p->coeffs);
+  }
+  *p = POLYTYPE_ZERO;
+}
+
+// Auxillary function for lowering the offset
+void POLYFUNC(lower_offset)(POLYTYPE *p, int dec) {
+  QPMS_ENSURE(dec >= 0, "Offset decrease must be positive, is %d!", dec);
+  QPMS_ENSURE(dec <= p->offset, "Offset can't be pushed below 0, (offset=%d, decr.=%d!)",
+      p->offset,  dec);
+  if(dec > 0) {
+    POLYFUNC(extend)(p, p->order - p->offset + dec + 1);
+    for(int i = p->order; i >= p->offset; --i)
+      COEFFUNC(swap)(p->coeffs[i+dec - p->offset], p->coeffs[i - p->offset]);
+    p->offset -= dec;
+    for(int i = dec - 1; i >= 0; --i)
+      COEFFUNC(set_si)(p->coeffs[i], 0, 1);
+  }
+}
+#undef POLYTYPE
+#undef COEFTYPE
+#undef POLYFUNC
+#undef COEFFUNC
+#undef POLYTYPE_ZERO

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
         )