Polynomials with rational coeffs: basic functionality.

Untested.


Former-commit-id: d89daafeb73b498aab73c4a1c67823f3044499b2

qpms/CMakeLists.txt

@@ -13,7 +13,7 @@ add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c e
 	ewaldsf.c pointgroups.c latticegens.c
 	lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
 	bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
-	lll.c
+	lll.c polynomials.c
 	)
 use_c99()
 

qpms/polynomials.c

@@ -0,0 +1,136 @@
+#include "polynomials.h"
+#include <stdlib.h>
+#include "qpms_error.h"
+#include <stdbool.h>
+
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#define MIN(x, y) (((x) <= (y)) ? (x) : (y))
+
+// 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?");
+}
+
+_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;
+}
+
+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;
+}
+
+void qpq_extend(qpq_t *p, int cap) {
+  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_clear(qpq_t *p) {
+  if (p->capacity > 0) {
+    for (int i = p->capacity; i >= 0; --i)
+      mpq_clear(p->coeffs[i]);
+    free(p->coeffs);
+  }
+  *p = QPQ_ZERO;
+}
+
+void qpq_add(qpq_t *sum, const qpq_t *x, const qpq_t *y) {
+  const int maxorder = MAX(x->order, y->order);
+  const int minoffset = MIN(x->offset, y->offset);
+  qpq_extend(sum, maxorder - minoffset + 1);
+  for (int i = minoffset; i <= maxorder; ++i) {
+    if (i - x->offset >= 0 && i <= x->order) {
+      if (i - y->offset >= 0 && i <= y->order) 
+        mpq_add(sum->coeffs[i - minoffset], 
+        x->coeffs[i - x->offset], y->coeffs[i - y->offset]);
+      else
+        mpq_set(sum->coeffs[i - minoffset], x->coeffs[i - x->offset]);
+    } else {
+      if (i - y->offset >= 0 && i <= y->order)
+        mpq_set(sum->coeffs[i - minoffset], y->coeffs[i - x->offset]);
+      else {
+        // Maybe not the best way to set it to zero.
+        // Alternatively, we could use just mpz_set_si(mpq_numref(sum->coeffs[i - minoffset]), 0);
+        mpq_clear(sum->coeffs[i - minoffset]);
+        mpq_init(sum->coeffs[i - minoffset]);
+      }
+    }
+  }
+  sum->offset = minoffset;
+  sum->order = maxorder;
+}
+      
+void qpq_sub(qpq_t *dif, const qpq_t *x, const qpq_t *y) {
+  const int maxorder = MAX(x->order, y->order);
+  const int minoffset = MIN(x->offset, y->offset);
+  qpq_extend(dif, maxorder - minoffset + 1);
+  for (int i = minoffset; i <= maxorder; ++i) {
+    if (i - x->offset >= 0 && i <= x->order) {
+      if (i - y->offset >= 0 && i <= y->order) 
+        mpq_sub(dif->coeffs[i - minoffset], 
+        x->coeffs[i - x->offset], y->coeffs[i - y->offset]);
+      else
+        mpq_set(dif->coeffs[i - minoffset], x->coeffs[i - x->offset]);
+    } else {
+      if (i - y->offset >= 0 && i <= y->order) {
+        mpq_set(dif->coeffs[i - minoffset], y->coeffs[i - x->offset]);
+        mpq_neg(dif->coeffs[i - minoffset], dif->coeffs[i - minoffset]);
+      } else {
+        // Maybe not the best way to set it to zero.
+        mpq_clear(dif->coeffs[i - minoffset]);
+        mpq_init(dif->coeffs[i - minoffset]);
+      }
+    }
+  }
+  dif->offset = minoffset;
+  dif->order = maxorder;
+}
+
+void qpq_mul(qpq_t *p, const qpq_t *x, const qpq_t *y) {
+  const int maxorder = x->order + y->order;
+  const int minoffset = x->offset + y->offset;
+  // Easiest way to set p to a zero polynomial...
+  qpq_clear(p);
+  qpq_init(p, maxorder - minoffset + 1);
+  for (int xi = x->offset; xi <= x->order; ++xi)
+    for (int yi = y->offset; yi <= y->order; ++yi)
+      mpq_mul(p->coeffs[xi + yi - minoffset],
+          x->coeffs[xi - x->offset], y->coeffs[yi - y->offset]);
+  p->order = maxorder;
+  p->offset = minoffset;
+}
+
+void qpq_deriv(qpq_t *dp, const qpq_t *p) {
+  if (p->order <= 0) { // p is constant, dp is zero.
+    qpq_clear(dp);
+    return;
+  }
+
+  qpq_extend(dp, p->order - p->offset + (p->offset > 0));
+
+  mpq_t qi;
+  mpq_init(qi); // qi is now 0 / 1
+  for(int i = p->offset + !p->offset; i <= p->order; ++i) {
+    mpz_set_si(mpq_numref(qi), 1); // qi is now i / 1
+    mpq_mul(dp->coeffs[i-1], qi, p->coeffs[i]);
+  }
+  mpq_clear(qi);
+  dp->order = p->order - 1;
+  dp->offset = p->offset - 1 + !p->offset;
+}

qpms/polynomials.h

@@ -6,6 +6,7 @@
 #include <gmp.h>
 
 /// Polynomial with rational coeffs.
+// TODO more docs about initialisation etc.
 typedef struct qpq_t {
 	int order;
 	int offset;
@@ -13,8 +14,24 @@ typedef struct qpq_t {
 	mpq_t *coeffs;
 } qpq_t;
 
+const static qpq_t QPQ_ZERO = {-1, 0, 0, NULL};
+
 /// Initiasise the coefficients array in qpq_t.
-void qpq_init(qpq_t *x, int maxorder);
+/** Do not use on qpq_t that has already been initialised
+ * (and not recently cleared),
+ * otherwise you can get a memory leak.
+ */
+void qpq_init(qpq_t *x, int capacity);
+
+/// Extend capacity of a qpq_t instance.
+/** If the requested new_capacity is larger than the qpq_t's
+ * capacity, the latter is extended to match new_capacity.
+ * Otherwise, nothing happend (this function does _not_ trim
+ * the capacity).
+ */
+void qpq_extend(qpq_t *x, int new_capacity);
+
+void qpq_set(qpq_t *copy, const qpq_t *orig);
 
 /// Deinitialise the coefficients array in qpq_t.
 void qpq_clear(qpq_t *x);
@@ -26,7 +43,12 @@ void qpq_add(qpq_t *sum, const qpq_t *addend1, const qpq_t *addend2);
 void qpq_sub(qpq_t *difference, const qpq_t *minuend, const qpq_t *substrahend);
 
 /// Polynomial multiplication.
-void qpq_mul(qpq_t product, const qpq_t *multiplier, const qpq_t *multiplicand);
+void qpq_mul(qpq_t *product, const qpq_t *multiplier, const qpq_t *multiplicand);
+
+/// Polynomial derivative.
+void qpq_deriv(qpq_t *dPdx, const qpq_t *P);
+
+_Bool qpq_nonzero(const qpq_t *);
 
 
 /// Polynomial with double coeffs.
@@ -56,7 +78,7 @@ void qpz_mul(qpz_t product, const qpz_t *multiplier, const qpz_t *multiplicand);
 void qpz_from_qpq(qpz_t *target, const qpq_t *src);
 
 
-
+#if 0 // This will go elsewhere
 /// Table with pre-calculated Ferrers function coeffs.
 typedef struct qpms_legendre_table_t {
 	qpms_l_t lMax;
@@ -75,5 +97,6 @@ double qpms_legendre_table_eval(const qpms_legendre_table_t *,
 
 
 // TODO pre-calculate also the products??
+#endif
 
 #endif //QPMS_POLYNOMIALS_H