From c0f2f7eae7e590b00ccd55b506e1f1499022d7ec Mon Sep 17 00:00:00 2001
From: Sven Verdoolaege <skimo@kotnet.org>
Date: Fri, 31 Aug 2007 18:16:59 +0200
Subject: [PATCH] support computation of Bernoulli polynomials

---
 bernoulli.c | 63 ++++++++++++++++++++++++++++++++++++++++---------------------
 bernoulli.h |  6 ++++++
 testlib.cc  | 12 +++++++++++-
 3 files changed, 59 insertions(+), 22 deletions(-)

diff --git a/bernoulli.c b/bernoulli.c
index d7328d6..5beadef 100644
--- a/bernoulli.c
+++ b/bernoulli.c
@@ -7,6 +7,7 @@
 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
 
 static struct bernoulli_coef bernoulli_coef;
+static struct poly_list bernoulli;
 static struct poly_list faulhaber;
 
 struct bernoulli_coef *bernoulli_coef_compute(int n)
@@ -80,52 +81,72 @@ struct bernoulli_coef *bernoulli_coef_compute(int n)
     return &bernoulli_coef;
 }
 
-struct poly_list *faulhaber_compute(int n)
+/*
+ * Compute either Bernoulli B_n or Faulhaber F_n polynomials.
+ *
+ * B_n =         sum_{k=0}^n {  n  \choose k } b_k x^{n-k}
+ * F_n = 1/(n+1) sum_{k=0}^n { n+1 \choose k } b_k x^{n+1-k}
+ */
+static struct poly_list *bernoulli_faulhaber_compute(int n, struct poly_list *pl,
+						     int faulhaber)
 {
     int i, j;
     Value factor;
     struct bernoulli_coef *bc;
 
-    if (n < faulhaber.n)
-	return &faulhaber;
+    if (n < pl->n)
+	return pl;
 
-    if (n >= faulhaber.size) {
+    if (n >= pl->size) {
 	int size = 3*(n + 5)/2;
 	Vector **poly;
 
 	poly = ALLOCN(Vector *, size);
-	for (i = 0; i < faulhaber.n; ++i)
-	    poly[i] = faulhaber.poly[i];
-	free(faulhaber.poly);
-	faulhaber.poly = poly;
+	for (i = 0; i < pl->n; ++i)
+	    poly[i] = pl->poly[i];
+	free(pl->poly);
+	pl->poly = poly;
 
-	faulhaber.size = size;
+	pl->size = size;
     }
 
     bc = bernoulli_coef_compute(n);
 
     value_init(factor);
-    for (i = faulhaber.n; i <= n; ++i) {
-	faulhaber.poly[i] = Vector_Alloc(i+3);
-	value_assign(faulhaber.poly[i]->p[i+1], bc->lcm->p[i]);
-	mpz_mul_ui(faulhaber.poly[i]->p[i+2], bc->lcm->p[i], i+1);
+    for (i = pl->n; i <= n; ++i) {
+	pl->poly[i] = Vector_Alloc(i+faulhaber+2);
+	value_assign(pl->poly[i]->p[i+faulhaber], bc->lcm->p[i]);
+	if (faulhaber)
+	    mpz_mul_ui(pl->poly[i]->p[i+2], bc->lcm->p[i], i+1);
+	else
+	    value_assign(pl->poly[i]->p[i+1], bc->lcm->p[i]);
 	value_set_si(factor, 1);
 	for (j = 1; j <= i; ++j) {
-	    mpz_mul_ui(factor, factor, i+2-j);
+	    mpz_mul_ui(factor, factor, i+faulhaber+1-j);
 	    mpz_divexact_ui(factor, factor, j);
-	    value_division(faulhaber.poly[i]->p[i+1-j],
+	    value_division(pl->poly[i]->p[i+faulhaber-j],
 			   bc->lcm->p[i], bc->den->p[j]);
-	    value_multiply(faulhaber.poly[i]->p[i+1-j],
-			   faulhaber.poly[i]->p[i+1-j], bc->num->p[j]);
-	    value_multiply(faulhaber.poly[i]->p[i+1-j],
-			   faulhaber.poly[i]->p[i+1-j], factor);
+	    value_multiply(pl->poly[i]->p[i+faulhaber-j],
+			   pl->poly[i]->p[i+faulhaber-j], bc->num->p[j]);
+	    value_multiply(pl->poly[i]->p[i+faulhaber-j],
+			   pl->poly[i]->p[i+faulhaber-j], factor);
 	}
-	Vector_Normalize(faulhaber.poly[i]->p, i+3);
+	Vector_Normalize(pl->poly[i]->p, i+faulhaber+2);
     }
     value_clear(factor);
     pl->n = n+1;
 
-    return &faulhaber;
+    return pl;
+}
+
+struct poly_list *bernoulli_compute(int n)
+{
+    return bernoulli_faulhaber_compute(n, &bernoulli, 0);
+}
+
+struct poly_list *faulhaber_compute(int n)
+{
+    return bernoulli_faulhaber_compute(n, &faulhaber, 1);
 }
 
 /* shift variables in polynomial one down */
diff --git a/bernoulli.h b/bernoulli.h
index 4df25d1..6ff0d3f 100644
--- a/bernoulli.h
+++ b/bernoulli.h
@@ -32,6 +32,12 @@ struct bernoulli_coef *bernoulli_coef_compute(int n);
  */
 struct poly_list *faulhaber_compute(int n);
 
+/*
+ * Compute Bernoulli polynomials up to the nth.
+ * The returned structure will contain at least n+1 polynomials (0..n).
+ */
+struct poly_list *bernoulli_compute(int n);
+
 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
 			struct barvinok_options *options);
 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
diff --git a/testlib.cc b/testlib.cc
index e115a1d..58c44e1 100644
--- a/testlib.cc
+++ b/testlib.cc
@@ -205,7 +205,7 @@ int test_todd(struct barvinok_options *options)
 int test_bernoulli(struct barvinok_options *options)
 {
     struct bernoulli_coef *bernoulli_coef;
-    struct poly_list *faulhaber;
+    struct poly_list *faulhaber, *bernoulli;
     bernoulli_coef = bernoulli_coef_compute(2);
     faulhaber = faulhaber_compute(4);
     bernoulli_coef = bernoulli_coef_compute(8);
@@ -218,6 +218,16 @@ int test_bernoulli(struct barvinok_options *options)
     assert(value_cmp_si(faulhaber->poly[3]->p[4], 1) == 0);
     assert(value_cmp_si(faulhaber->poly[3]->p[5], 4) == 0);
 
+    bernoulli = bernoulli_compute(6);
+    assert(value_cmp_si(bernoulli->poly[6]->p[0], 1) == 0);
+    assert(value_cmp_si(bernoulli->poly[6]->p[1], 0) == 0);
+    assert(value_cmp_si(bernoulli->poly[6]->p[2], -21) == 0);
+    assert(value_cmp_si(bernoulli->poly[6]->p[3], 0) == 0);
+    assert(value_cmp_si(bernoulli->poly[6]->p[4], 105) == 0);
+    assert(value_cmp_si(bernoulli->poly[6]->p[5], -126) == 0);
+    assert(value_cmp_si(bernoulli->poly[6]->p[6], 42) == 0);
+    assert(value_cmp_si(bernoulli->poly[6]->p[7], 42) == 0);
+
     unsigned nvar, nparam;
     char **all_vars;
     evalue *base, *sum1, *sum2;
-- 
2.11.4.GIT