From 9f9c72cbd6a01adeb4bdfb864c93402808f3af20 Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Wed, 5 Sep 2007 21:53:56 +0200 Subject: [PATCH] separate computation of Bernoulli coefficients and Faulhaber polynomials --- bernoulli.c | 163 ++++++++++++++++++++++++++++++++++++------------------------ bernoulli.h | 28 ++++++++--- testlib.cc | 27 +++++----- 3 files changed, 134 insertions(+), 84 deletions(-) diff --git a/bernoulli.c b/bernoulli.c index df1e774..d7328d6 100644 --- a/bernoulli.c +++ b/bernoulli.c @@ -6,93 +6,126 @@ #define ALLOC(type) (type*)malloc(sizeof(type)) #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type)) -static struct bernoulli bernoulli; +static struct bernoulli_coef bernoulli_coef; +static struct poly_list faulhaber; -struct bernoulli *bernoulli_compute(int n) +struct bernoulli_coef *bernoulli_coef_compute(int n) { int i, j; - Value lcm, factor, tmp; + Value factor, tmp; - if (n < bernoulli.n) - return &bernoulli; + if (n < bernoulli_coef.n) + return &bernoulli_coef; - if (n >= bernoulli.size) { - int size = n + 5; + if (n >= bernoulli_coef.size) { + int size = 3*(n + 5)/2; Vector *b; - Vector **sum; b = Vector_Alloc(size); - if (bernoulli.n) { - Vector_Copy(bernoulli.b_num->p, b->p, bernoulli.n); - Vector_Free(bernoulli.b_num); + if (bernoulli_coef.n) { + Vector_Copy(bernoulli_coef.num->p, b->p, bernoulli_coef.n); + Vector_Free(bernoulli_coef.num); } - bernoulli.b_num = b; + bernoulli_coef.num = b; b = Vector_Alloc(size); - if (bernoulli.n) { - Vector_Copy(bernoulli.b_den->p, b->p, bernoulli.n); - Vector_Free(bernoulli.b_den); + if (bernoulli_coef.n) { + Vector_Copy(bernoulli_coef.den->p, b->p, bernoulli_coef.n); + Vector_Free(bernoulli_coef.den); } - bernoulli.b_den = b; - sum = ALLOCN(Vector *, size); - for (i = 0; i < bernoulli.n; ++i) - sum[i] = bernoulli.sum[i]; - free(bernoulli.sum); - bernoulli.sum = sum; - - bernoulli.size = size; + bernoulli_coef.den = b; + b = Vector_Alloc(size); + if (bernoulli_coef.n) { + Vector_Copy(bernoulli_coef.lcm->p, b->p, bernoulli_coef.n); + Vector_Free(bernoulli_coef.lcm); + } + bernoulli_coef.lcm = b; + + bernoulli_coef.size = size; } - value_init(lcm); value_init(factor); value_init(tmp); - value_set_si(lcm, 1); - for (i = 0; i < bernoulli.n; ++i) - value_lcm(lcm, bernoulli.b_den->p[i], &lcm); - for (i = bernoulli.n; i <= n; ++i) { + for (i = bernoulli_coef.n; i <= n; ++i) { if (i == 0) { - value_set_si(bernoulli.b_num->p[0], 1); - value_set_si(bernoulli.b_den->p[0], 1); + value_set_si(bernoulli_coef.num->p[0], 1); + value_set_si(bernoulli_coef.den->p[0], 1); + value_set_si(bernoulli_coef.lcm->p[0], 1); continue; } - value_set_si(bernoulli.b_num->p[i], 0); + value_set_si(bernoulli_coef.num->p[i], 0); value_set_si(factor, -(i+1)); for (j = i-1; j >= 0; --j) { mpz_mul_ui(factor, factor, j+1); mpz_divexact_ui(factor, factor, i+1-j); - value_division(tmp, lcm, bernoulli.b_den->p[j]); - value_multiply(tmp, tmp, bernoulli.b_num->p[j]); + value_division(tmp, bernoulli_coef.lcm->p[i-1], + bernoulli_coef.den->p[j]); + value_multiply(tmp, tmp, bernoulli_coef.num->p[j]); value_multiply(tmp, tmp, factor); - value_addto(bernoulli.b_num->p[i], bernoulli.b_num->p[i], tmp); + value_addto(bernoulli_coef.num->p[i], bernoulli_coef.num->p[i], tmp); } - mpz_mul_ui(bernoulli.b_den->p[i], lcm, i+1); - Gcd(bernoulli.b_num->p[i], bernoulli.b_den->p[i], &tmp); + mpz_mul_ui(bernoulli_coef.den->p[i], bernoulli_coef.lcm->p[i-1], i+1); + Gcd(bernoulli_coef.num->p[i], bernoulli_coef.den->p[i], &tmp); if (value_notone_p(tmp)) { - value_division(bernoulli.b_num->p[i], bernoulli.b_num->p[i], tmp); - value_division(bernoulli.b_den->p[i], bernoulli.b_den->p[i], tmp); + value_division(bernoulli_coef.num->p[i], + bernoulli_coef.num->p[i], tmp); + value_division(bernoulli_coef.den->p[i], + bernoulli_coef.den->p[i], tmp); } - value_lcm(lcm, bernoulli.b_den->p[i], &lcm); + value_lcm(bernoulli_coef.lcm->p[i-1], bernoulli_coef.den->p[i], + &bernoulli_coef.lcm->p[i]); } - for (i = bernoulli.n; i <= n; ++i) { - bernoulli.sum[i] = Vector_Alloc(i+3); - value_assign(bernoulli.sum[i]->p[i+1], lcm); - mpz_mul_ui(bernoulli.sum[i]->p[i+2], lcm, i+1); + bernoulli_coef.n = n+1; + value_clear(factor); + value_clear(tmp); + + return &bernoulli_coef; +} + +struct poly_list *faulhaber_compute(int n) +{ + int i, j; + Value factor; + struct bernoulli_coef *bc; + + if (n < faulhaber.n) + return &faulhaber; + + if (n >= faulhaber.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; + + faulhaber.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); value_set_si(factor, 1); for (j = 1; j <= i; ++j) { mpz_mul_ui(factor, factor, i+2-j); mpz_divexact_ui(factor, factor, j); - value_division(bernoulli.sum[i]->p[i+1-j], lcm, bernoulli.b_den->p[j]); - value_multiply(bernoulli.sum[i]->p[i+1-j], - bernoulli.sum[i]->p[i+1-j], bernoulli.b_num->p[j]); - value_multiply(bernoulli.sum[i]->p[i+1-j], - bernoulli.sum[i]->p[i+1-j], factor); + value_division(faulhaber.poly[i]->p[i+1-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); } - Vector_Normalize(bernoulli.sum[i]->p, i+3); + Vector_Normalize(faulhaber.poly[i]->p, i+3); } - bernoulli.n = n+1; - value_clear(lcm); value_clear(factor); - value_clear(tmp); + pl->n = n+1; - return &bernoulli; + return &faulhaber; } /* shift variables in polynomial one down */ @@ -117,7 +150,7 @@ static evalue *shifted_copy(evalue *src) return e; } -static evalue *power_sums(struct bernoulli *bernoulli, evalue *poly, +static evalue *power_sums(struct poly_list *faulhaber, evalue *poly, Vector *arg, Value denom, int neg, int alt_neg) { int i; @@ -125,7 +158,7 @@ static evalue *power_sums(struct bernoulli *bernoulli, evalue *poly, evalue *sum = evalue_zero(); for (i = 1; i < poly->x.p->size; ++i) { - evalue *term = evalue_polynomial(bernoulli->sum[i], base); + evalue *term = evalue_polynomial(faulhaber->poly[i], base); evalue *factor = shifted_copy(&poly->x.p->arr[i]); emul(factor, term); if (alt_neg && (i % 2)) @@ -213,10 +246,10 @@ static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data) } else { evalue *poly_u = NULL, *poly_l = NULL; Polyhedron *D; - struct bernoulli *bernoulli; + struct poly_list *faulhaber; assert(data->e->x.p->type == polynomial); assert(data->e->x.p->pos == 1); - bernoulli = bernoulli_compute(data->e->x.p->size-1); + faulhaber = faulhaber_compute(data->e->x.p->size-1); /* lower is the constraint * b i - l' >= 0 i >= l'/b = l * upper is the constraint @@ -241,10 +274,10 @@ static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data) Vector_Copy(upper+2, row->p, dim+1); value_oppose(tmp, upper[1]); value_addto(row->p[dim], row->p[dim], tmp); - poly_u = power_sums(bernoulli, data->e, row, tmp, 0, 0); + poly_u = power_sums(faulhaber, data->e, row, tmp, 0, 0); } Vector_Oppose(lower+2, row->p, dim+1); - extra = power_sums(bernoulli, data->e, row, lower[1], 1, 0); + extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0); eadd(poly_u, extra); eadd(linear, extra); @@ -267,11 +300,11 @@ static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data) if (!poly_l) { Vector_Copy(lower+2, row->p, dim+1); value_addto(row->p[dim], row->p[dim], lower[1]); - poly_l = power_sums(bernoulli, data->e, row, lower[1], 0, 1); + poly_l = power_sums(faulhaber, data->e, row, lower[1], 0, 1); } Vector_Oppose(upper+2, row->p, dim+1); value_oppose(tmp, upper[1]); - extra = power_sums(bernoulli, data->e, row, tmp, 1, 1); + extra = power_sums(faulhaber, data->e, row, tmp, 1, 1); eadd(poly_l, extra); eadd(linear, extra); @@ -293,13 +326,13 @@ static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data) if (!poly_l) { Vector_Copy(lower+2, row->p, dim+1); value_addto(row->p[dim], row->p[dim], lower[1]); - poly_l = power_sums(bernoulli, data->e, row, lower[1], 0, 1); + poly_l = power_sums(faulhaber, data->e, row, lower[1], 0, 1); } if (!poly_u) { Vector_Copy(upper+2, row->p, dim+1); value_oppose(tmp, upper[1]); value_addto(row->p[dim], row->p[dim], tmp); - poly_u = power_sums(bernoulli, data->e, row, tmp, 0, 0); + poly_u = power_sums(faulhaber, data->e, row, tmp, 0, 0); } data->s[data->ns].E = ALLOC(evalue); @@ -328,7 +361,7 @@ static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data) Vector_Copy(upper+2, row->p, dim+1); value_oppose(tmp, upper[1]); value_addto(row->p[dim], row->p[dim], tmp); - poly_u = power_sums(bernoulli, data->e, row, tmp, 0, 0); + poly_u = power_sums(faulhaber, data->e, row, tmp, 0, 0); } eadd(linear, poly_u); @@ -357,7 +390,7 @@ static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data) if (!poly_l) { Vector_Copy(lower+2, row->p, dim+1); value_addto(row->p[dim], row->p[dim], lower[1]); - poly_l = power_sums(bernoulli, data->e, row, lower[1], 0, 1); + poly_l = power_sums(faulhaber, data->e, row, lower[1], 0, 1); } eadd(linear, poly_l); diff --git a/bernoulli.h b/bernoulli.h index 9bf99ec..4df25d1 100644 --- a/bernoulli.h +++ b/bernoulli.h @@ -6,15 +6,31 @@ extern "C" { struct barvinok_options; -struct bernoulli { - Vector *b_num; - Vector *b_den; - Vector **sum; +struct bernoulli_coef { + Vector *num; + Vector *den; + Vector *lcm; /* lcm of this and previous denominators */ int size; - int n; + int n; /* The number of Bernoulli coefficients */ }; -struct bernoulli *bernoulli_compute(int n); +struct poly_list { + Vector **poly; + int size; + int n; /* The number of polynomials */ +}; + +/* + * Compute Bernoulli coefficients up to the nth. + * The returned structure will contain at least n+1 coefficients (0..n). + */ +struct bernoulli_coef *bernoulli_coef_compute(int n); + +/* + * Compute Faulhaber polynomials up to the nth. + * The returned structure will contain at least n+1 polynomials (0..n). + */ +struct poly_list *faulhaber_compute(int n); evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C, struct barvinok_options *options); diff --git a/testlib.cc b/testlib.cc index 72e6eb5..e115a1d 100644 --- a/testlib.cc +++ b/testlib.cc @@ -204,18 +204,19 @@ int test_todd(struct barvinok_options *options) int test_bernoulli(struct barvinok_options *options) { - struct bernoulli *bernoulli; - bernoulli = bernoulli_compute(2); - bernoulli = bernoulli_compute(4); - bernoulli = bernoulli_compute(8); - assert(value_cmp_si(bernoulli->b_num->p[6], 1) == 0); - assert(value_cmp_si(bernoulli->b_den->p[6], 42) == 0); - assert(value_cmp_si(bernoulli->sum[3]->p[0], 0) == 0); - assert(value_cmp_si(bernoulli->sum[3]->p[1], 0) == 0); - assert(value_cmp_si(bernoulli->sum[3]->p[2], 1) == 0); - assert(value_cmp_si(bernoulli->sum[3]->p[3], -2) == 0); - assert(value_cmp_si(bernoulli->sum[3]->p[4], 1) == 0); - assert(value_cmp_si(bernoulli->sum[3]->p[5], 4) == 0); + struct bernoulli_coef *bernoulli_coef; + struct poly_list *faulhaber; + bernoulli_coef = bernoulli_coef_compute(2); + faulhaber = faulhaber_compute(4); + bernoulli_coef = bernoulli_coef_compute(8); + assert(value_cmp_si(bernoulli_coef->num->p[6], 1) == 0); + assert(value_cmp_si(bernoulli_coef->den->p[6], 42) == 0); + assert(value_cmp_si(faulhaber->poly[3]->p[0], 0) == 0); + assert(value_cmp_si(faulhaber->poly[3]->p[1], 0) == 0); + assert(value_cmp_si(faulhaber->poly[3]->p[2], 1) == 0); + assert(value_cmp_si(faulhaber->poly[3]->p[3], -2) == 0); + assert(value_cmp_si(faulhaber->poly[3]->p[4], 1) == 0); + assert(value_cmp_si(faulhaber->poly[3]->p[5], 4) == 0); unsigned nvar, nparam; char **all_vars; @@ -223,7 +224,7 @@ int test_bernoulli(struct barvinok_options *options) base = evalue_read_from_str("(1 * n + 1)", NULL, &all_vars, &nvar, &nparam, options->MaxRays); - sum1 = evalue_polynomial(bernoulli->sum[3], base); + sum1 = evalue_polynomial(faulhaber->poly[3], base); Free_ParamNames(all_vars, nvar+nparam); sum2 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)", -- 2.11.4.GIT