From 1ea9d9cb470b36c07214562ea2c5f75de991027f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Sun, 10 Nov 2019 00:24:08 +0200
Subject: [PATCH] mpqs_t multiplication and simple initialisation.

Former-commit-id: 2d2a5df78e1d0b47e300569c352e7758a0e64f94
---
 qpms/polynomials.c | 59 ++++++++++++++++++++++++++++++++++++++++++++++
 qpms/polynomials.h |  4 ++++
 2 files changed, 63 insertions(+)

diff --git a/qpms/polynomials.c b/qpms/polynomials.c
index b29fdd1..7e16bda 100644
--- a/qpms/polynomials.c
+++ b/qpms/polynomials.c
@@ -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.
diff --git a/qpms/polynomials.h b/qpms/polynomials.h
index f1220d8..a0d57fc 100644
--- a/qpms/polynomials.h
+++ b/qpms/polynomials.h
@@ -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);