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

mpqs_t multiplication and simple initialisation. 

Former-commit-id: 2d2a5df78e1d0b47e300569c352e7758a0e64f94
This commit is contained in:
Marek Nečada 2019-11-10 00:24:08 +02:00
parent fa15871bd9
commit 1ea9d9cb47
2 changed files with 63 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
qpms/polynomials.c
View File

@ -93,6 +93,34 @@ void mpqs_nt_append(mpqs_t x, const struct _qp_mpzs_hashed *numelem) {
mpz_size(n->_2) * sizeof(mp_limb_t), n);
}
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);
}
void mpqs_nt_add(mpqs_t x, const struct _qp_mpzs_hashed *addend) {
struct _qp_mpzs_hashed *s;
HASH_FIND(hh, x->nt, mpz_limbs_read(addend->_2),
@ -204,6 +232,37 @@ void mpqs_sub(mpqs_t dif, const mpqs_t x, const mpqs_t y) {
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));
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->nt->hh.next;
nx != NULL; nx = nx->hh.next) {
for (const struct _qp_mpzs_hashed *ny = x->nt->hh.next;
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);
mpqs_nt_add(prod, addend);
}
}
mpzs_hh_clear(addend);
mpz_clear(mx); mpz_clear(my); mpz_clear(gcd);
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.

4
qpms/polynomials.h
View File

@ -100,6 +100,10 @@ typedef struct _qp_mpqs {
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);