Compare commits

...

18 Commits

Author SHA1 Message Date
Marek Nečada 6cecad6d18 mpqs refactoring minor stuff
Former-commit-id: c17570e41047a309cc76bbe485309acfcac6e82a
2019-11-10 20:10:26 +02:00
Marek Nečada 345352dbd5 mpqs_t minor refactoring
Former-commit-id: 79a252c00ab7bf0f8d6c05d68be5b0adae51f9a6
2019-11-10 15:12:42 +02:00
Marek Nečada 1ea9d9cb47 mpqs_t multiplication and simple initialisation.
Former-commit-id: 2d2a5df78e1d0b47e300569c352e7758a0e64f94
2019-11-10 00:24:08 +02:00
Marek Nečada fa15871bd9 mpqs_t basics (addition, substraction) now compile (not yet tested)
Former-commit-id: 6bf0dd8553d514fe187ea5b4f6f8ff7a2752362e
2019-11-08 20:56:39 +02:00
Marek Nečada 34e52a42cf WIP mpqs_t
Former-commit-id: c0a46a1a7f0f00e1950a067f74440bd40e695db9
2019-11-06 00:22:20 +02:00
Marek Nečada 706073315b WIP mpzs_hh_t
Former-commit-id: 5df48ddb3c9aac26887c81c54d71b47173b9ec0d
2019-11-01 04:40:17 +02:00
Marek Nečada c0ec5389e0 WIP polynomials etc.
Former-commit-id: 546c139f9cce0686c1b3065403beaa769bd7b55b
2019-10-29 14:00:09 +02:00
Troy D. Hanson & al 6b5b773d0b Add uthash.h
Former-commit-id: 8bbe0ac9bdbc0e6b30c5b3c2af51b34f4f4dd4a3
2019-10-28 13:41:54 +02:00
Marek Nečada 57a1206dcc A new polynomial class with "square root of rationals" coeffs.
Former-commit-id: 4f06dc080b45dc0e89f5899dce728fb340cd6ef8
2019-10-27 09:45:08 +02:00
Marek Nečada 0c5a38371f qpq_t ointer mixing for _mul and _div, fix last merge.
Former-commit-id: 4c7a5de21c8c783ae4af4781b77f8b0d3babd6e3
2019-10-25 21:03:55 +03:00
Marek Nečada 5f77fd9386 Merge branch 'gmp_integration' of necada:~/repo/qpms into gmp_integration
Former-commit-id: 7d6118efbb06d1cee0eb3e89afc4aeb65665d023
2019-10-25 20:43:46 +03:00
Marek Nečada 8db97686e5 Polynomial division, bug fixes, enhancements.
- Addition and substraction fixed, supports result and operand mixing.


Former-commit-id: 0a53db4cd1fee70c6cf77a6ea2c65840de6d3ca1
2019-10-25 20:40:05 +03:00
Marek Nečada d6b3951dc2 "Legendroid" data structure
Former-commit-id: 34908e2933a0aae79c46d665c2abe20ab5504eed
2019-10-22 06:42:36 +03:00
Marek Nečada b4ac597771 qpq cython wrapper class, bug fixes
The new cython module qpms.cypolynomials requires gmpy2 version 2.1 or
higher.


Former-commit-id: 10198a7341437620f5a493d16a217bf6014e2759
2019-10-20 21:26:20 +03:00
Marek Nečada 25d1eced26 qpq_t element assignment functionality
Former-commit-id: 63edc3b1f27209af18932f6fd8f284c1859eb40d
2019-10-19 14:34:09 +03:00
Marek Nečada a4df0139b9 Polynomials with rational coeffs: basic functionality.
Untested.


Former-commit-id: d89daafeb73b498aab73c4a1c67823f3044499b2
2019-10-15 16:30:13 +03:00
Marek Nečada 3f5aaacbf2 First draft of essential polynomial arithmetic functionality.
Former-commit-id: f34d258aaffbd52c4e4d251f5b6f1ff3733e3a08
2019-10-13 18:08:15 +03:00
Marek Nečada 4ca003ec8b CMake scripts for finding GMP
Former-commit-id: 7ac8c4899598b44afa36300dd962a2173c07d1c1
2019-10-09 14:29:38 +03:00
11 changed files with 2320 additions and 4 deletions

View File

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

View File

@ -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.

View File

@ -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)

3
cmake-scripts/ORIGINS Normal file
View File

@ -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/

View File

@ -2,6 +2,7 @@
find_package(GSL 2.0 REQUIRED) find_package(GSL 2.0 REQUIRED)
find_package(BLAS REQUIRED) find_package(BLAS REQUIRED)
find_package(LAPACK REQUIRED) find_package(LAPACK REQUIRED)
find_package(GMP REQUIRED)
#includes #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 ewaldsf.c pointgroups.c latticegens.c
lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.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 bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
lll.c lll.c polynomials.c
) )
use_c99() use_c99()
@ -23,6 +24,7 @@ target_link_libraries (qpms
lapack lapack
blas blas
amos amos
gmp
) )
target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})

146
qpms/cypolynomials.pyx Normal file
View File

@ -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

546
qpms/polynomials.c Normal file
View File

@ -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);
}

179
qpms/polynomials.h Normal file
View File

@ -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

159
qpms/polynomials.template Normal file
View File

@ -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

1230
qpms/uthash.h Normal file

File diff suppressed because it is too large Load Diff

View File

@ -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 #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', setup(name='qpms',
version = "0.2.996", version = "0.2.996",
packages=['qpms'], packages=['qpms'],
# libraries = [('amos', {'sources': amos_sources} )], # libraries = [('amos', {'sources': amos_sources} )],
setup_requires=['cython>=0.28',], setup_requires=['cython>=0.28', 'gmpy2>=2.1b'],
install_requires=['cython>=0.28', install_requires=['cython>=0.28',
#'quaternion','spherical_functions', #'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'], #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}, cmdclass = {'build_ext': build_ext},
zip_safe=False zip_safe=False
) )