QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

A new polynomial class with "square root of rationals" coeffs. 

Former-commit-id: 4f06dc080b45dc0e89f5899dce728fb340cd6ef8
This commit is contained in:
Marek Nečada 2019-10-27 09:45:08 +02:00
parent 0c5a38371f
commit 57a1206dcc
3 changed files with 252 additions and 131 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

158
qpms/polynomials.c
View File

@ -8,13 +8,6 @@
#define MIN(x, y) (((x) <= (y)) ? (x) : (y))
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);
}
// 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.
@ -23,109 +16,43 @@ static inline void mpq_zero(mpq_t q) {
mpq_init(q);
}
// TODO to the template
static inline void qpq_zero(qpq_t *q) {
qpq_clear(q);
}
// 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?");
#define TEMPLATE_QPQ
#include "polynomials.template"
#undef TEMPLATE_QPQ
// 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);
}
_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;
static inline void qpqs_zero(qpqs_t *q) {
qpqs_clear(q);
}
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;
#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);
}
void qpq_extend(qpq_t *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(mpq_t) * cap);
for(int i = p->capacity; i < cap; ++i)
mpq_init(p->coeffs[i]);
p->capacity = cap;
}
}
void qpq_shrink(qpq_t *p) {
if (p->capacity > 0) {
for(int n = p->capacity - 1; n > p->order - p->offset; --n)
mpq_clear(p->coeffs[n]);
p->capacity = p->order - p->offset + 1;
if (p->capacity > 0)
QPMS_CRASHING_REALLOC(p->coeffs, p->capacity * sizeof(mpq_t));
}
}
void qpq_canonicalise(qpq_t *p) {
qpq_cc(p);
// Lower the order if necessary (here one can get -1 if the polynomial is 0)
while (p->order >= p->offset && !mpq_sgn(p->coeffs[p->order - p->offset])) --p->order;
if (p->order < p->offset)
qpq_zero(p);
else { // remove the lowest-order coefficients which are in fact zero.
int i = 0;
while (i <= p->order - p->offset && !mpq_sgn(p->coeffs[i])) ++i;
if (i > 0)
for (int j = 0; j <= p->order - p->offset - i; ++j)
mpq_swap(p->coeffs[j], p->coeffs[j+i]);
p->offset += i;
// canonicalise the fractions
for (i = 0; i <= p->order - p->offset; ++i)
mpq_canonicalize(p->coeffs[i]);
}
}
void qpq_set(qpq_t *p, const qpq_t *orig) {
if(p != orig) {
const int order = orig->order, offset = orig->offset;
qpq_extend(p, order - offset + 1);
for (int i = orig->offset; i <= order; ++i)
mpq_set(p->coeffs[i - offset], orig->coeffs[i - offset]);
p->offset = offset;
p->order = order;
}
}
void qpq_set_elem(qpq_t *p, const int exponent, const mpq_t coeff) {
QPMS_ENSURE(exponent >= 0, "Exponent must be non-negative, got %d", exponent);
int offset = p->offset, order = p->order;
if(mpq_sgn(coeff) == 0 && (exponent < offset || exponent > order))
return; // exponent out of range, but zero needs not to be assigned explicitly
if(exponent < p->offset) {
qpq_extend(p, p->order - exponent + 1);
offset = exponent;
for(int i = order - offset; i >= p->offset - offset; --i)
mpq_swap(p->coeffs[i], p->coeffs[i - (p->offset - offset)]);
for(int i = p->offset - offset - 1; i > 0; --i)
mpq_zero(p->coeffs[i]);
p->offset = offset;
}
if(exponent > order) {
qpq_extend(p, exponent - p->offset + 1);
for(int i = p->order - p->offset + 1; i <= exponent - p->offset; ++i)
mpq_zero(p->coeffs[i]);
p->order = exponent;
}
mpq_set(p->coeffs[exponent - p->offset], coeff);
return;
}
void qpq_set_elem_si(qpq_t *p, const int exponent, const long int num, const unsigned long int den) {
mpq_t q;
@ -135,13 +62,6 @@ void qpq_set_elem_si(qpq_t *p, const int exponent, const long int num, const uns
mpq_clear(q);
}
void qpq_get_elem(mpq_t coeff, const qpq_t *p, const int exponent) {
if (exponent < p->offset || exponent > p->order)
mpq_zero(coeff);
else
mpq_set(coeff, p->coeffs[exponent-p->offset]);
}
int qpq_get_elem_si(long *num, unsigned long *den, const qpq_t *p, const int exponent) {
mpq_t q;
mpq_init(q);
@ -157,30 +77,6 @@ int qpq_get_elem_si(long *num, unsigned long *den, const qpq_t *p, const int exp
return retval;
}
void qpq_clear(qpq_t *p) {
if (p->capacity > 0) {
for (int i = p->capacity - 1; i >= 0; --i)
mpq_clear(p->coeffs[i]);
free(p->coeffs);
}
*p = QPQ_ZERO;
}
// Auxillary function for lowering the offset
void qpq_lower_offset(qpq_t *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) {
qpq_extend(p, p->order - p->offset + dec + 1);
for(int i = p->order; i >= p->offset; --i)
mpq_swap(p->coeffs[i+dec - p->offset], p->coeffs[i - p->offset]);
p->offset -= dec;
for(int i = dec - 1; i >= 0; --i)
mpq_set_si(p->coeffs[i], 0, 1);
}
}
void qpq_add(qpq_t *sum, const qpq_t *x, const qpq_t *y) {
const int maxorder = MAX(x->order, y->order);

66
qpms/polynomials.h
View File

@ -71,6 +71,72 @@ void qpq_deriv(qpq_t *dPdx, const qpq_t *P);
/// Tests whether a polynomial is non-zero.
_Bool qpq_nonzero(const qpq_t *);
/// Square root of a rational number.
typedef struct mpqs_struct_t {
mpq_t _2;
} mpqs_t[1];
static inline void mpqs_init(mpqs_t q) {mpq_init(q->_2);}
static inline void mpqs_clear(mpqs_t q) {mpq_clear(q->_2);}
static inline void mpqs_mul(mpqs_t product, const mpqs_t multiplier, const mpqs_t multiplicand) {
mpq_mul(product->_2, multiplier->_2, multiplicand->_2);
}
static inline void mpqs_div(mpqs_t quotient, const mpqs_t dividend, const mpqs_t divisor) {
mpq_div(quotient->_2, dividend->_2, divisor->_2);
}
static inline void mpqs_set(mpqs_t copy, const mpqs_t orig) { mpq_set(copy->_2, orig->_2); }
static inline void mpqs_set_mpq(mpqs_t qs, const mpq_t q) {
mpq_mul(qs->_2, q, q);
}
static inline int mpqs_sgn(const mpqs_t q) { return mpq_sgn(q->_2); }
static inline void mpqs_canonicalize(mpqs_t q) { mpq_canonicalize(q->_2); }
static inline void mpqs_swap(mpqs_t a, mpqs_t b) { mpq_swap(a->_2, b->_2); }
/// Polynomial with rational square root coeffs.
// TODO more docs about initialisation etc.
typedef struct qpqs_t {
int order;
int offset;
int capacity;
mpqs_t *coeffs;
} qpqs_t;
const static qpqs_t QPQS_ZERO = {-1, 0, 0, NULL};
/// Initiasise the coefficients array in qpqs_t.
/** Do not use on qpqs_t that has already been initialised
* (and not recently cleared),
* otherwise you can get a memory leak.
*/
void qpqs_init(qpqs_t *p, int capacity);
/// Extend capacity of a qpqs_t instance.
/** If the requested new_capacity is larger than the qpqs_t's
* capacity, the latter is extended to match new_capacity.
* Otherwise, nothing happend (this function does _not_ trim
* the capacity).
*/
void qpqs_extend(qpqs_t *p, int new_capacity);
/// Shrinks the capacity to the minimum that can store the current polynomial.
void qpqs_shrink(qpqs_t *p);
/// Canonicalises the coefficients and (re)sets the correct degree.
void qpqs_canonicalise(qpqs_t *p);
void qpqs_set(qpqs_t *copy, const qpqs_t *orig);
void qpqs_set_elem(qpqs_t *p, int exponent, const mpqs_t coeff);
void qpqs_set_elem_si(qpqs_t *p, int exponent, long numerator, unsigned long denominator);
void qpqs_get_elem(mpqs_t coeff, const qpqs_t *p, int exponent);
/** \returns zero if the result fits into long / unsigned long; non-zero otherwise. */
int qpqs_get_elem_si(long *numerator, unsigned long *denominator, const qpqs_t *p, int exponent);
/// Deinitialise the coefficients array in qpqs_t.
void qpqs_clear(qpqs_t *p);
/// 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;

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