From 34e52a42cf2c48b492c898724c5e9a7bb3d53b2a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 6 Nov 2019 00:22:20 +0200 Subject: [PATCH] WIP mpqs_t Former-commit-id: c0a46a1a7f0f00e1950a067f74440bd40e695db9 --- qpms/polynomials.c | 148 ++++++++++++++++++++++++++++++++++++++++----- qpms/polynomials.h | 50 ++------------- 2 files changed, 139 insertions(+), 59 deletions(-) diff --git a/qpms/polynomials.c b/qpms/polynomials.c index 37c8a21..fc60f46 100644 --- a/qpms/polynomials.c +++ b/qpms/polynomials.c @@ -7,6 +7,7 @@ #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); @@ -23,39 +24,156 @@ void mpzs_set(mpzs_t x, const mpzs_t y) { } // 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 -void mpqs_clear_num(mpqs_t x) { - TODO; +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); +} + +//===== mpqs_t ===== + void mpqs_init(mpqs_t x) { - x->sz = 0; + //x->sz = 0; x->nt = 0; - mpz_init_set_ui(x->den, 1); + mpq_init(x->f); + mpq_set_ui(x->f, 1, 1); +} + +void mpqs_clear_num(mpqs_t x) { + struct _qp_mpzs_hashed *current, *tmp; + HASH_ITER(hh, *(x->nt), current, tmp) { + HASH_DEL(x->nt, current); + free(current); + //x->sz--; + } } void mpqs_clear(mpqs_t x) { mpqs_clear_num(mpqs_t x); - mpz_clear(x->den); - x->sz = -1 + mpq_clear(x->f); + //x->sz = 0; x->nt = 0; } -static void mpqs_num_add_elem(mpqs_t q, mpzs_t n) { - mpzs_t *n_owned; - QPMS_CRASHING_MALLOC(n_owned, sizeof(mpzs_t)); - mpzs_init(*n_owned); - mpzs_set(*n_owned, n); - void *added = - tsearch(n_owned, &(t->nt), mpzs_cmp2); - QPMS_ENSURE(added, "Failed to add numerator element. Memory error?"); - QPMS_ENSURE(added != n_owned, "FIXME another numerator elements with the same square root found."); // TODO how to handle this? +void mpqs_nt_append(mpqs_t x, const mpzs_hh_t numelem) { + mpzs_hh_t *n; + QPMS_CRASHING_MALLOC(n, sizeof(mpzs_hh_t)); + mpzs_hh_init(*n); + mpzs_hh_set(*n, numelem); + HASH_ADD_KEYPTR(hh, x->nt, mpz_limbs_read(*n->_2), + mpz_size(*n->_2) * sizeof(mp_limb_t), n); } +void mpqs_nt_add(mpqs_t x, const mpzs_hh_t addend) { + mpzs_hh_t *s; + HASH_FIND(hh, x->nt, mpz_limbs_read(addend->_2), + mpz_size(addend->_2), s); + if (!s) mpqs_nt_append(addend, x); // if not found + else { + mpz_add(*s->_1, *s->_1, addend); + if(!mpz_sgn(*s->_1)) { // If zero, annihilate + HASH_DEL(x->nt, s); + mpzs_hh_clear(*s); + free(*s); + } + } +} +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_nt_gcd(mpz_t gcd, const mpqs_t x) { + if(x->nt) { + mpz_set(gcd, *(x->nt)->_1); + for(mpzs_hh_t *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); + 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); + + mpzs_hh_t hh_cur; + + // 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(y->f)); + + // Left operand numerator factor + mpz_t tmp; mpz_init(tmp); + mpz_set(tmp, mpq_denref(y->f)); + mpz_divexact(tmp, 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(mpzs_hh_t *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_set(tmp, mpq_denref(x->f)); + mpz_divexact(tmp, den_gcd); + mpz_mul(tmp, tmp, mpq_numref(y->f)); + + mpzs_hh_t addend; mpzs_hh_init(addend); + for(const mpzs_hh_t *n = y->nt; n != NULL; n = n->hh.next) { + mpz_mul(addend->_1, tmp, *n->_1); + mpz_set(addend->_2, *n->_2); + mpqs_nt_add(y, addend); + } + mpzs_hh_clear(addend); + mpz_clear(tmp); + mpz_clear(den_gcd); + + mpqs_canonicalise(sum); + mpqs_set(sum_final, sum); + mpqs_clear(sum); +} diff --git a/qpms/polynomials.h b/qpms/polynomials.h index 383b1c4..c53a585 100644 --- a/qpms/polynomials.h +++ b/qpms/polynomials.h @@ -4,6 +4,7 @@ */ #ifndef QPMS_POLYNOMIALS_H #include +#include "uthash.h" /// Polynomial with rational coeffs. // TODO more docs about initialisation etc. @@ -72,47 +73,6 @@ void qpq_deriv(qpq_t *dPdx, const qpq_t *P); _Bool qpq_nonzero(const qpq_t *); -/// 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); #if 0 /// Type representing a number of form \f$ a \sqrt{b}; a \in \ints, b \in \nats \f$. @@ -137,21 +97,23 @@ typedef struct _qp_mpzs_hashed { 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); /// 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. + //int sz; ///< Used size of the numtree. mpzs_hh_t *nt; ///< List of numerator components.. - mpz_t den; ///< Denominator. Always positive. + 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_add(mpqs_t sum, const mpqs_t addend1, const mpqs_t addend2); 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(mpqs_t quotient, const mpqs_t dividend, const mpqs_t divisor); +//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.