From 57a1206dcc4bb7d218a072bbccc76ee9b350cf55 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Sun, 27 Oct 2019 09:45:08 +0200
Subject: [PATCH] A new polynomial class with "square root of rationals"
 coeffs.

Former-commit-id: 4f06dc080b45dc0e89f5899dce728fb340cd6ef8
---
 qpms/polynomials.c        | 158 +++++++------------------------------
 qpms/polynomials.h        |  66 ++++++++++++++++
 qpms/polynomials.template | 159 ++++++++++++++++++++++++++++++++++++++
 3 files changed, 252 insertions(+), 131 deletions(-)
 create mode 100644 qpms/polynomials.template

diff --git a/qpms/polynomials.c b/qpms/polynomials.c
index 9e87775..e8f16ec 100644
--- a/qpms/polynomials.c
+++ b/qpms/polynomials.c
@@ -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);
diff --git a/qpms/polynomials.h b/qpms/polynomials.h
index a99ef74..2f7adff 100644
--- a/qpms/polynomials.h
+++ b/qpms/polynomials.h
@@ -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;
diff --git a/qpms/polynomials.template b/qpms/polynomials.template
new file mode 100644
index 0000000..944f158
--- /dev/null
+++ b/qpms/polynomials.template
@@ -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