Faulhaber's formula: Difference between revisions

Added solution using C
(Added solution using C)
Line 7:
;See also
* [https://en.wikipedia.org/wiki/Bernoulli_number Bernoulli numbers] and [https://en.wikipedia.org/wiki/Binomial_coefficient binomial coefficients] on Wikipedia.
 
=={{header|C}}==
{{trans|Modula-2}}
<lang c>#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
 
int binomial(int n, int k) {
int num, denom, i;
 
if (n < 0 || k < 0 || n < k) return -1;
if (n == 0 || k == 0) return 1;
 
num = 1;
for (i = k + 1; i <= n; ++i) {
num = num * i;
}
 
denom = 1;
for (i = 2; i <= n - k; ++i) {
denom *= i;
}
 
return num / denom;
}
 
int gcd(int a, int b) {
int temp;
while (b != 0) {
temp = a % b;
a = b;
b = temp;
}
return a;
}
 
typedef struct tFrac {
int num, denom;
} Frac;
 
Frac makeFrac(int n, int d) {
Frac result;
int g;
 
if (d == 0) {
result.num = 0;
result.denom = 0;
return result;
}
 
if (n == 0) {
d = 1;
} else if (d < 0) {
n = -n;
d = -d;
}
 
g = abs(gcd(n, d));
if (g > 1) {
n = n / g;
d = d / g;
}
 
result.num = n;
result.denom = d;
return result;
}
 
Frac negateFrac(Frac f) {
return makeFrac(-f.num, f.denom);
}
 
Frac subFrac(Frac lhs, Frac rhs) {
return makeFrac(lhs.num * rhs.denom - lhs.denom * rhs.num, rhs.denom * lhs.denom);
}
 
Frac multFrac(Frac lhs, Frac rhs) {
return makeFrac(lhs.num * rhs.num, lhs.denom * rhs.denom);
}
 
bool equalFrac(Frac lhs, Frac rhs) {
return (lhs.num == rhs.num) && (lhs.denom == rhs.denom);
}
 
bool lessFrac(Frac lhs, Frac rhs) {
return (lhs.num * rhs.denom) < (rhs.num * lhs.denom);
}
 
void printFrac(Frac f) {
printf("%d", f.num);
if (f.denom != 1) {
printf("/%d", f.denom);
}
}
 
Frac bernoulli(int n) {
Frac a[16];
int j, m;
 
if (n < 0) {
a[0].num = 0;
a[0].denom = 0;
return a[0];
}
 
for (m = 0; m <= n; ++m) {
a[m] = makeFrac(1, m + 1);
for (j = m; j >= 1; --j) {
a[j - 1] = multFrac(subFrac(a[j - 1], a[j]), makeFrac(j, 1));
}
}
 
if (n != 1) {
return a[0];
}
 
return negateFrac(a[0]);
}
 
void faulhaber(int p) {
Frac coeff, q;
int j, pwr, sign;
 
printf("%d : ", p);
q = makeFrac(1, p + 1);
sign = -1;
for (j = 0; j <= p; ++j) {
sign = -1 * sign;
coeff = multFrac(multFrac(multFrac(q, makeFrac(sign, 1)), makeFrac(binomial(p + 1, j), 1)), bernoulli(j));
if (equalFrac(coeff, makeFrac(0, 1))) {
continue;
}
if (j == 0) {
if (!equalFrac(coeff, makeFrac(1, 1))) {
if (equalFrac(coeff, makeFrac(-1, 1))) {
printf("-");
} else {
printFrac(coeff);
}
}
} else {
if (equalFrac(coeff, makeFrac(1, 1))) {
printf(" + ");
} else if (equalFrac(coeff, makeFrac(-1, 1))) {
printf(" - ");
} else if (lessFrac(makeFrac(0, 1), coeff)) {
printf(" + ");
printFrac(coeff);
} else {
printf(" - ");
printFrac(negateFrac(coeff));
}
}
pwr = p + 1 - j;
if (pwr > 1) {
printf("n^%d", pwr);
} else {
printf("n");
}
}
printf("\n");
}
 
int main() {
int i;
 
for (i = 0; i < 10; ++i) {
faulhaber(i);
}
 
return 0;
}</lang>
{{out}}
<pre>0 : n
1 : 1/2n^2 + 1/2n
2 : 1/3n^3 + 1/2n^2 + 1/6n
3 : 1/4n^4 + 1/2n^3 + 1/4n^2
4 : 1/5n^5 + 1/2n^4 + 1/3n^3 - 1/30n
5 : 1/6n^6 + 1/2n^5 + 5/12n^4 - 1/12n^2
6 : 1/7n^7 + 1/2n^6 + 1/2n^5 - 1/6n^3 + 1/42n
7 : 1/8n^8 + 1/2n^7 + 7/12n^6 - 7/24n^4 + 1/12n^2
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2</pre>
 
=={{header|C++}}==
1,452

edits