From a4df0139b90d0fd91497ecf3dba8167c327c9d6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 15 Oct 2019 16:30:13 +0300 Subject: [PATCH] Polynomials with rational coeffs: basic functionality. Untested. Former-commit-id: d89daafeb73b498aab73c4a1c67823f3044499b2 --- qpms/CMakeLists.txt | 2 +- qpms/polynomials.c | 136 ++++++++++++++++++++++++++++++++++++++++++++ qpms/polynomials.h | 29 +++++++++- 3 files changed, 163 insertions(+), 4 deletions(-) create mode 100644 qpms/polynomials.c diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index 856c41e..bf65f76 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -13,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() diff --git a/qpms/polynomials.c b/qpms/polynomials.c new file mode 100644 index 0000000..d5140c1 --- /dev/null +++ b/qpms/polynomials.c @@ -0,0 +1,136 @@ +#include "polynomials.h" +#include +#include "qpms_error.h" +#include + +#define MAX(x, y) (((x) > (y)) ? (x) : (y)) +#define MIN(x, y) (((x) <= (y)) ? (x) : (y)) + +// qpq_t internal consistency check +static inline void qpq_cc(const qpq_t *p) { + if (!p->coeffs) return; + QPMS_ENSURE(p->capacity >= p->order - p->offset + 1, "qpq_t inconsistency detected; have you initialised it properly?"); +} + +_Bool qpq_nonzero(const qpq_t *p) { + qpq_cc(p); + if (p->capacity <= 0) return false; + + for(int i = 0; i <= p->order - p->offset; ++i) + if (mpq_sgn(p->coeffs[i])) + return true; + return false; +} + +void qpq_init(qpq_t *p, int capacity) { + *p = QPQ_ZERO; + if (capacity <= 0) + return; + QPMS_CRASHING_MALLOC(p->coeffs, capacity * sizeof(mpq_t)); + for(int i = 0; i < capacity; ++i) + mpq_init(p->coeffs[i]); + p->capacity = capacity; +} + +void qpq_extend(qpq_t *p, int cap) { + if (cap > 0 && cap > p->capacity) { + QPMS_CRASHING_REALLOC(p->coeffs, sizeof(mpq_t) * cap); + for(int i = p->capacity; i < cap; ++i) + mpq_init(p->coeffs[i]); + p->capacity = cap; + } +} + +void qpq_clear(qpq_t *p) { + if (p->capacity > 0) { + for (int i = p->capacity; i >= 0; --i) + mpq_clear(p->coeffs[i]); + free(p->coeffs); + } + *p = QPQ_ZERO; +} + +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); + 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 - x->offset]); + else { + // 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(sum->coeffs[i - minoffset]); + mpq_init(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); + 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 - x->offset]); + mpq_neg(dif->coeffs[i - minoffset], dif->coeffs[i - minoffset]); + } else { + // Maybe not the best way to set it to zero. + mpq_clear(dif->coeffs[i - minoffset]); + mpq_init(dif->coeffs[i - minoffset]); + } + } + } + dif->offset = minoffset; + dif->order = maxorder; +} + +void qpq_mul(qpq_t *p, const qpq_t *x, const qpq_t *y) { + const int maxorder = x->order + y->order; + const int minoffset = x->offset + y->offset; + // Easiest way to set p to a zero polynomial... + qpq_clear(p); + qpq_init(p, maxorder - minoffset + 1); + for (int xi = x->offset; xi <= x->order; ++xi) + for (int yi = y->offset; yi <= y->order; ++yi) + mpq_mul(p->coeffs[xi + yi - minoffset], + x->coeffs[xi - x->offset], y->coeffs[yi - y->offset]); + p->order = maxorder; + p->offset = minoffset; +} + +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), 1); // 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; +} diff --git a/qpms/polynomials.h b/qpms/polynomials.h index 8101578..0dcb089 100644 --- a/qpms/polynomials.h +++ b/qpms/polynomials.h @@ -6,6 +6,7 @@ #include /// Polynomial with rational coeffs. +// TODO more docs about initialisation etc. typedef struct qpq_t { int order; int offset; @@ -13,8 +14,24 @@ typedef struct qpq_t { mpq_t *coeffs; } qpq_t; +const static qpq_t QPQ_ZERO = {-1, 0, 0, NULL}; + /// Initiasise the coefficients array in qpq_t. -void qpq_init(qpq_t *x, int maxorder); +/** 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 *x, 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 *x, int new_capacity); + +void qpq_set(qpq_t *copy, const qpq_t *orig); /// Deinitialise the coefficients array in qpq_t. void qpq_clear(qpq_t *x); @@ -26,7 +43,12 @@ 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); /// Polynomial multiplication. -void qpq_mul(qpq_t product, const qpq_t *multiplier, const qpq_t *multiplicand); +void qpq_mul(qpq_t *product, const qpq_t *multiplier, const qpq_t *multiplicand); + +/// Polynomial derivative. +void qpq_deriv(qpq_t *dPdx, const qpq_t *P); + +_Bool qpq_nonzero(const qpq_t *); /// Polynomial with double coeffs. @@ -56,7 +78,7 @@ void qpz_mul(qpz_t product, const qpz_t *multiplier, const qpz_t *multiplicand); 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; @@ -75,5 +97,6 @@ double qpms_legendre_table_eval(const qpms_legendre_table_t *, // TODO pre-calculate also the products?? +#endif #endif //QPMS_POLYNOMIALS_H