Rosetta Code

Faulhaber's formula: Difference between revisions

Page Discussion

Faulhaber's formula (view source)

Revision as of 13:06, 17 March 2024

51,843 bytes added ,  2 months ago
Added FreeBASIC
(Added FreeBASIC)
 
(36 intermediate revisions by 14 users not shown)
Line 1:
{{draft task|Mathematics}}
In mathematics, [[wp:Faulhaber's formula|Faulhaber's formula]], named after Johann Faulhaber, expresses the sum of the ''p''-th powers of the first ''n'' positive integers as a ''(p + 1)''th-degree polynomial function of n, the coefficients involving [[Bernoulli numbers]].
 
In mathematics,   Faulhaber's formula,   named after Johann Faulhaber,   expresses the sum of the ''p''-th powers of the first ''n'' positive integers as a ''(p + 1)''th-degree polynomial function of n,   the coefficients involving [[Bernoulli numbers]].
;Task
 
 
;Task:
Generate the first 10 closed-form expressions, starting with ''p = 0''.
 
 
;See also
;Related tasks:
* [https://en.wikipedia.org/wiki/Bernoulli_number Bernoulli numbers] and [https://en.wikipedia.org/wiki/Binomial_coefficient binomial coefficients] on Wikipedia.
:*   [[Bernoulli numbers]].
:*   [[Evaluate_binomial_coefficients|evaluate binomial coefficients]].
 
 
;See also:
:*   The Wikipedia entry:   [[wp:Faulhaber's formula|Faulhaber's formula]].
:*   The Wikipedia entry:   [https://en.wikipedia.org/wiki/Bernoulli_number Bernoulli numbers].
:*   The Wikipedia entry:   [https://en.wikipedia.org/wiki/Binomial_coefficient binomial coefficients].
<br><br>
 
=={{header|C}}==
{{trans|Modula-2}}
<langsyntaxhighlight lang="c">#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
Line 178 ⟶ 189:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>0 : n
Line 191 ⟶ 202:
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2</pre>
 
=={{header|C++ sharp|C#}}==
{{trans|D}}
Uses C++17
<lang cpp>#include <iostream>
#include <numeric>
#include <sstream>
#include <vector>
 
class Frac {
public:
Frac(long n, long d) {
if (d == 0) {
throw new std::runtime_error("d must not be zero");
}
 
long nn = n;
long dd = d;
if (nn == 0) {
dd = 1;
} else if (dd < 0) {
nn = -nn;
dd = -dd;
}
 
long g = abs(std::gcd(nn, dd));
if (g > 1) {
nn /= g;
dd /= g;
}
 
num = nn;
denom = dd;
}
 
Frac operator-() const {
return Frac(-num, denom);
}
 
Frac operator+(const Frac& rhs) const {
return Frac(num*rhs.denom + denom * rhs.num, rhs.denom*denom);
}
 
Frac operator-(const Frac& rhs) const {
return Frac(num*rhs.denom - denom * rhs.num, rhs.denom*denom);
}
 
Frac operator*(const Frac& rhs) const {
return Frac(num*rhs.num, denom*rhs.denom);
}
 
bool operator==(const Frac& rhs) const {
return num == rhs.num && denom == rhs.denom;
}
 
bool operator!=(const Frac& rhs) const {
return num != rhs.num || denom != rhs.denom;
}
 
bool operator<(const Frac& rhs) const {
if (denom == rhs.denom) {
return num < rhs.num;
}
return num * rhs.denom < rhs.num * denom;
}
 
friend std::ostream& operator<<(std::ostream&, const Frac&);
 
static Frac ZERO() {
return Frac(0, 1);
}
 
static Frac ONE() {
return Frac(1, 1);
}
 
private:
long num;
long denom;
};
 
std::ostream & operator<<(std::ostream & os, const Frac &f) {
if (f.num == 0 || f.denom == 1) {
return os << f.num;
}
 
std::stringstream ss;
ss << f.num << "/" << f.denom;
return os << ss.str();
}
 
Frac bernoulli(int n) {
if (n < 0) {
throw new std::runtime_error("n may not be negative or zero");
}
 
std::vector<Frac> a;
for (int m = 0; m <= n; m++) {
a.push_back(Frac(1, m + 1));
for (int j = m; j >= 1; j--) {
a[j - 1] = (a[j - 1] - a[j]) * Frac(j, 1);
}
}
 
// returns 'first' Bernoulli number
if (n != 1) return a[0];
return -a[0];
}
 
int binomial(int n, int k) {
if (n < 0 || k < 0 || n < k) {
throw new std::runtime_error("parameters are invalid");
}
if (n == 0 || k == 0) return 1;
 
int num = 1;
for (int i = k + 1; i <= n; i++) {
num *= i;
}
 
int denom = 1;
for (int i = 2; i <= n - k; i++) {
denom *= i;
}
 
return num / denom;
}
 
void faulhaber(int p) {
using namespace std;
cout << p << " : ";
 
auto q = Frac(1, p + 1);
int sign = -1;
for (int j = 0; j < p + 1; j++) {
sign *= -1;
auto coeff = q * Frac(sign, 1) * Frac(binomial(p + 1, j), 1) * bernoulli(j);
if (coeff == Frac::ZERO()) {
continue;
}
if (j == 0) {
if (coeff == -Frac::ONE()) {
cout << "-";
} else if (coeff != Frac::ONE()) {
cout << coeff;
}
} else {
if (coeff == Frac::ONE()) {
cout << " + ";
} else if (coeff == -Frac::ONE()) {
cout << " - ";
} else if (coeff < Frac::ZERO()) {
cout << " - " << -coeff;
} else {
cout << " + " << coeff;
}
}
int pwr = p + 1 - j;
if (pwr > 1) {
cout << "n^" << pwr;
} else {
cout << "n";
}
}
cout << endl;
}
 
int main() {
for (int 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#|C sharp}}==
{{trans|Java}}
<langsyntaxhighlight lang="csharp">using System;
 
namespace FaulhabersFormula {
Line 557 ⟶ 383:
}
}
}</langsyntaxhighlight>
{{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++}}==
{{trans|D}}
Uses C++17
<syntaxhighlight lang="cpp">#include <iostream>
#include <numeric>
#include <sstream>
#include <vector>
 
class Frac {
public:
Frac(long n, long d) {
if (d == 0) {
throw new std::runtime_error("d must not be zero");
}
 
long nn = n;
long dd = d;
if (nn == 0) {
dd = 1;
} else if (dd < 0) {
nn = -nn;
dd = -dd;
}
 
long g = abs(std::gcd(nn, dd));
if (g > 1) {
nn /= g;
dd /= g;
}
 
num = nn;
denom = dd;
}
 
Frac operator-() const {
return Frac(-num, denom);
}
 
Frac operator+(const Frac& rhs) const {
return Frac(num*rhs.denom + denom * rhs.num, rhs.denom*denom);
}
 
Frac operator-(const Frac& rhs) const {
return Frac(num*rhs.denom - denom * rhs.num, rhs.denom*denom);
}
 
Frac operator*(const Frac& rhs) const {
return Frac(num*rhs.num, denom*rhs.denom);
}
 
bool operator==(const Frac& rhs) const {
return num == rhs.num && denom == rhs.denom;
}
 
bool operator!=(const Frac& rhs) const {
return num != rhs.num || denom != rhs.denom;
}
 
bool operator<(const Frac& rhs) const {
if (denom == rhs.denom) {
return num < rhs.num;
}
return num * rhs.denom < rhs.num * denom;
}
 
friend std::ostream& operator<<(std::ostream&, const Frac&);
 
static Frac ZERO() {
return Frac(0, 1);
}
 
static Frac ONE() {
return Frac(1, 1);
}
 
private:
long num;
long denom;
};
 
std::ostream & operator<<(std::ostream & os, const Frac &f) {
if (f.num == 0 || f.denom == 1) {
return os << f.num;
}
 
std::stringstream ss;
ss << f.num << "/" << f.denom;
return os << ss.str();
}
 
Frac bernoulli(int n) {
if (n < 0) {
throw new std::runtime_error("n may not be negative or zero");
}
 
std::vector<Frac> a;
for (int m = 0; m <= n; m++) {
a.push_back(Frac(1, m + 1));
for (int j = m; j >= 1; j--) {
a[j - 1] = (a[j - 1] - a[j]) * Frac(j, 1);
}
}
 
// returns 'first' Bernoulli number
if (n != 1) return a[0];
return -a[0];
}
 
int binomial(int n, int k) {
if (n < 0 || k < 0 || n < k) {
throw new std::runtime_error("parameters are invalid");
}
if (n == 0 || k == 0) return 1;
 
int num = 1;
for (int i = k + 1; i <= n; i++) {
num *= i;
}
 
int denom = 1;
for (int i = 2; i <= n - k; i++) {
denom *= i;
}
 
return num / denom;
}
 
void faulhaber(int p) {
using namespace std;
cout << p << " : ";
 
auto q = Frac(1, p + 1);
int sign = -1;
for (int j = 0; j < p + 1; j++) {
sign *= -1;
auto coeff = q * Frac(sign, 1) * Frac(binomial(p + 1, j), 1) * bernoulli(j);
if (coeff == Frac::ZERO()) {
continue;
}
if (j == 0) {
if (coeff == -Frac::ONE()) {
cout << "-";
} else if (coeff != Frac::ONE()) {
cout << coeff;
}
} else {
if (coeff == Frac::ONE()) {
cout << " + ";
} else if (coeff == -Frac::ONE()) {
cout << " - ";
} else if (coeff < Frac::ZERO()) {
cout << " - " << -coeff;
} else {
cout << " + " << coeff;
}
}
int pwr = p + 1 - j;
if (pwr > 1) {
cout << "n^" << pwr;
} else {
cout << "n";
}
}
cout << endl;
}
 
int main() {
for (int i = 0; i < 10; i++) {
faulhaber(i);
}
 
return 0;
}</syntaxhighlight>
{{out}}
<pre>0 : n
Line 572 ⟶ 583:
=={{header|D}}==
{{trans|Kotlin}}
<langsyntaxhighlight Dlang="d">import std.algorithm : fold;
import std.exception : enforce;
import std.format : formattedWrite;
Line 721 ⟶ 732:
faulhaber(i);
}
}</langsyntaxhighlight>
{{out}}
<pre>0 : n
Line 735 ⟶ 746:
 
=={{header|EchoLisp}}==
<langsyntaxhighlight lang="scheme">
(lib 'math) ;; for bernoulli numbers
(string-delimiter "")
Line 753 ⟶ 764:
(define (Faulcomp n p)
(printf "Σ(1..%d) n^%d = %d" n p (/ (poly n (Faulhaber p)) (1+ p) )))
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 782 ⟶ 793:
</pre>
 
=={{header|Factor}}==
<syntaxhighlight lang="factor">USING: formatting kernel math math.combinatorics math.extras
math.functions regexp sequences ;
 
: faulhaber ( p -- seq )
1 + dup recip swap dup <iota>
[ [ nCk ] [ -1 swap ^ ] [ bernoulli ] tri * * * ] 2with map ;
 
: (poly>str) ( seq -- str )
reverse [ 1 + "%un^%d" sprintf ] map-index reverse " + " join ;
 
: clean-up ( str -- str' )
R/ n\^1\z/ "n" re-replace ! Change n^1 to n.
R/ 1n/ "n" re-replace ! Change 1n to n.
R/ \+ -/ "- " re-replace ! Change + - to - .
R/ [+-] 0n(\^\d+ )?/ "" re-replace ; ! Remove terms of zero.
 
: poly>str ( seq -- str ) (poly>str) clean-up ;
 
10 [ dup faulhaber poly>str "%d: %s\n" printf ] each-integer</syntaxhighlight>
{{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|FreeBASIC}}==
{{trans|C}}
<syntaxhighlight lang="vbnet">Type Fraction
As Integer num
As Integer den
End Type
 
Function Binomial(n As Integer, k As Integer) As Integer
If n < 0 Or k < 0 Then Print "Arguments must be non-negative integers": Exit Function
If n < k Then Print "The second argument cannot be more than the first.": Exit Function
If n = 0 Or k = 0 Then Return 1
Dim As Integer i, num, den
num = 1
For i = k + 1 To n
num *= i
Next i
den = 1
For i = 2 To n - k
den *= i
Next i
Return num / den
End Function
 
Function GCD(n As Integer, d As Integer) As Integer
Return Iif(d = 0, n, GCD(d, n Mod d))
End Function
 
Function makeFrac(n As Integer, d As Integer) As Fraction
Dim As Fraction result
Dim As Integer g
If d = 0 Then
result.num = 0
result.den = 0
Return result
End If
If n = 0 Then
d = 1
Elseif d < 0 Then
n = -n
d = -d
End If
g = Abs(gcd(n, d))
If g > 1 Then
n /= g
d /= g
End If
result.num = n
result.den = d
Return result
End Function
 
Function negateFrac(f As Fraction) As Fraction
Return makeFrac(-f.num, f.den)
End Function
 
Function subFrac(lhs As Fraction, rhs As Fraction) As Fraction
Return makeFrac(lhs.num * rhs.den - lhs.den * rhs.num, rhs.den * lhs.den)
End Function
 
Function multFrac(lhs As Fraction, rhs As Fraction) As Fraction
Return makeFrac(lhs.num * rhs.num, lhs.den * rhs.den)
End Function
 
Function equalFrac(lhs As Fraction, rhs As Fraction) As Integer
Return (lhs.num = rhs.num) And (lhs.den = rhs.den)
End Function
 
Function lessFrac(lhs As Fraction, rhs As Fraction) As Integer
Return (lhs.num * rhs.den) < (rhs.num * lhs.den)
End Function
 
Sub printFrac(f As Fraction)
Print Str(f.num);
If f.den <> 1 Then Print "/" & f.den;
End Sub
 
Function Bernoulli(n As Integer) As Fraction
If n < 0 Then Print "Argument must be non-negative": Exit Function
Dim As Fraction a(16)
Dim As Integer j, m
If (n < 0) Then
a(0).num = 0
a(0).den = 0
Return a(0)
End If
For m = 0 To n
a(m) = makeFrac(1, m + 1)
For j = m To 1 Step -1
a(j - 1) = multFrac(subFrac(a(j - 1), a(j)), makeFrac(j, 1))
Next j
Next m
If n <> 1 Then Return a(0)
Return negateFrac(a(0))
End Function
 
Sub Faulhaber(p As Integer)
Dim As Fraction coeff, q
Dim As Integer j, pwr, sign
Print p & " : ";
q = makeFrac(1, p + 1)
sign = -1
For j = 0 To p
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))) Then Continue For
If j = 0 Then
If Not equalFrac(coeff, makeFrac(1, 1)) Then
If equalFrac(coeff, makeFrac(-1, 1)) Then
Print "-";
Else
printFrac(coeff)
End If
End If
Else
If equalFrac(coeff, makeFrac(1, 1)) Then
Print " + ";
Elseif equalFrac(coeff, makeFrac(-1, 1)) Then
Print " - ";
Elseif lessFrac(makeFrac(0, 1), coeff) Then
Print " + ";
printFrac(coeff)
Else
Print " - ";
printFrac(negateFrac(coeff))
End If
End If
pwr = p + 1 - j
Print Iif(pwr > 1, "n^" & pwr, "n");
Next j
Print
End Sub
 
For i As Integer = 0 To 9
Faulhaber(i)
Next i
 
Sleep</syntaxhighlight>
{{out}}
<pre>Same as C entry.</pre>
 
=={{header|Fōrmulæ}}==
{{FormulaeEntry|page=https://formulae.org/?script=examples/Faulhaber}}
 
'''Solution'''
The following function creates the Faulhaber's coefficients up to a given number of rows, according to the [http://www.ww.ingeniousmathstat.org/sites/default/files/Torabi-Dashti-CMJ-2011.pdf paper] of of Mohammad Torabi Dashti:
 
(This is exactly the as than the task [[Faulhaber%27s_triangle#F%C5%8Drmul%C3%A6|Faulhaber's triangle]])
 
[[File:Fōrmulæ - Faulhaber 01.png]]
 
The following function creates the sum of the p-th powers of the first n positive integers as a (p + 1)th-degree polynomial function of n:
 
Notes. The -1 index means the last element (-2 is the penultimate element, and so on). So it retrieves the last row of the triangle. |x| is the cardinality (number of elements) of x.
 
(This is exactly the as than the task [[Faulhaber%27s_triangle#F%C5%8Drmul%C3%A6|Faulhaber's triangle]])
 
This function can be used for both symbolic or numeric computation of the polynomial:
 
[[File:Fōrmulæ - Faulhaber 06.png]]
 
To generate the first 10 closed-form expressions, starting with p = 0:
 
[[File:Fōrmulæ - Faulhaber 07.png]]
 
[[File:Fōrmulæ - Faulhaber 08.png]]
 
=={{header|GAP}}==
Straightforward implementation using GAP polynomials, and two different formulas: one based on Stirling numbers of the second kind (sum1, see Python implementation below in this page), and the usual Faulhaber formula (sum2). No optimization is made (one could compute Stirling numbers row by row, or the product in sum1 may be kept from one call to the other). Notice the Bernoulli term in the first formula is here only to correct the value of sum1(0), which is off by one because sum1 computes sums from 0 to n.
 
<langsyntaxhighlight lang="gap">n := X(Rationals, "n");
sum1 := p -> Sum([0 .. p], k -> Stirling2(p, k) * Product([0 .. k], j -> n + 1 - j) / (k + 1)) + 2 * Bernoulli(2 * p + 1);
sum2 := p -> Sum([0 .. p], j -> (-1)^j * Binomial(p + 1, j) * Bernoulli(j) * n^(p + 1 - j)) / (p + 1);
Line 804 ⟶ 1,026:
1/8*n^8+1/2*n^7+7/12*n^6-7/24*n^4+1/12*n^2
1/9*n^9+1/2*n^8+2/3*n^7-7/15*n^5+2/9*n^3-1/30*n
1/10*n^10+1/2*n^9+3/4*n^8-7/10*n^6+1/2*n^4-3/20*n^2</langsyntaxhighlight>
 
=={{header|Go}}==
<langsyntaxhighlight Golang="go">package main
 
import (
Line 866 ⟶ 1,088:
fmt.Println()
}
}</langsyntaxhighlight>
{{out}}
<pre>
Line 880 ⟶ 1,102:
9 : 1/10×n^10 -1/2 ×n^9 +3/4 ×n^8 -7/10×n^6 +1/2 ×n^4 -3/20×n^2
</pre>
 
=={{header|Groovy}}==
{{trans|Java}}
<syntaxhighlight lang="groovy">import java.util.stream.IntStream
 
class FaulhabersFormula {
private static long gcd(long a, long b) {
if (b == 0) {
return a
}
return gcd(b, a % b)
}
 
private static class Frac implements Comparable<Frac> {
private long num
private long denom
 
public static final Frac ZERO = new Frac(0, 1)
public static final Frac ONE = new Frac(1, 1)
 
Frac(long n, long d) {
if (d == 0) throw new IllegalArgumentException("d must not be zero")
long nn = n
long dd = d
if (nn == 0) {
dd = 1
} else if (dd < 0) {
nn = -nn
dd = -dd
}
long g = Math.abs(gcd(nn, dd))
if (g > 1) {
nn /= g
dd /= g
}
num = nn
denom = dd
}
 
Frac plus(Frac rhs) {
return new Frac(num * rhs.denom + denom * rhs.num, rhs.denom * denom)
}
 
Frac negative() {
return new Frac(-num, denom)
}
 
Frac minus(Frac rhs) {
return this + -rhs
}
 
Frac multiply(Frac rhs) {
return new Frac(this.num * rhs.num, this.denom * rhs.denom)
}
 
@Override
int compareTo(Frac o) {
double diff = toDouble() - o.toDouble()
return Double.compare(diff, 0.0)
}
 
@Override
boolean equals(Object obj) {
return null != obj && obj instanceof Frac && this == (Frac) obj
}
 
@Override
String toString() {
if (denom == 1) {
return Long.toString(num)
}
return String.format("%d/%d", num, denom)
}
 
private double toDouble() {
return (double) num / denom
}
}
 
private static Frac bernoulli(int n) {
if (n < 0) throw new IllegalArgumentException("n may not be negative or zero")
Frac[] a = new Frac[n + 1]
Arrays.fill(a, Frac.ZERO)
for (int m = 0; m <= n; ++m) {
a[m] = new Frac(1, m + 1)
for (int j = m; j >= 1; --j) {
a[j - 1] = (a[j - 1] - a[j]) * new Frac(j, 1)
}
}
// returns 'first' Bernoulli number
if (n != 1) return a[0]
return -a[0]
}
 
private static int binomial(int n, int k) {
if (n < 0 || k < 0 || n < k) throw new IllegalArgumentException()
if (n == 0 || k == 0) return 1
int num = IntStream.rangeClosed(k + 1, n).reduce(1, { a, b -> a * b })
int den = IntStream.rangeClosed(2, n - k).reduce(1, { acc, i -> acc * i })
return num / den
}
 
private static void faulhaber(int p) {
print("$p : ")
Frac q = new Frac(1, p + 1)
int sign = -1
for (int j = 0; j <= p; ++j) {
sign *= -1
Frac coeff = q * new Frac(sign, 1) * new Frac(binomial(p + 1, j), 1) * bernoulli(j)
if (Frac.ZERO == coeff) continue
if (j == 0) {
if (Frac.ONE != coeff) {
if (-Frac.ONE == coeff) {
print("-")
} else {
print(coeff)
}
}
} else {
if (Frac.ONE == coeff) {
print(" + ")
} else if (-Frac.ONE == coeff) {
print(" - ")
} else if (coeff > Frac.ZERO) {
print(" + $coeff")
} else {
print(" - ${-coeff}")
}
}
int pwr = p + 1 - j
if (pwr > 1) {
print("n^$pwr")
} else {
print("n")
}
}
println()
}
 
static void main(String[] args) {
for (int i = 0; i <= 9; ++i) {
faulhaber(i)
}
}
}</syntaxhighlight>
{{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|Haskell}}==
====Bernouilli polynomials====
<langsyntaxhighlight Haskelllang="haskell">import Data.Ratio ((%), numerator, denominator)
import Data.List (intercalate, transpose)
import ControlData.ArrowBifunctor ((&&&), (***)bimap)
import Data.Char (isSpace)
import Data.Monoid ((<>))
import Data.Bool (bool)
 
-- FAULHABER -------------------------------------------------------------------
------------------------- FAULHABER ------------------------
faulhaber :: [[Rational]]
faulhaber =
Line 899 ⟶ 1,278:
[]
[0 ..]
 
-- EXPRESSION STRINGS ----------------------------------------------------------
polynomials :: [[(String, String)]]
polynomials = fmap ((ratioPower =<<) . reverse . flip zip [1 ..]) faulhaber
 
---------------------------- TEST --------------------------
main :: IO ()
main = (putStrLn . unlines . expressionTable . take 10) polynomials
 
--------------------- EXPRESSION STRINGS -------------------
 
-- Rows of (Power string, Ratio string) tuples -> Printable lines
Line 914 ⟶ 1,300:
zipWith
(\(lw, rw) col ->
(fmap (bimap (justifyLeft lw ' ') *** (justifyLeft rw ' ')) col))
(colWidths cols)
cols)
 
-- Value pair -> String pair (lifted into list for use with >>=)
ratioPower :: (Rational, Integer) -> [(String, String)]
ratioPower (nd, j) =
let (num, den) = ((,) . numerator &&&<*> denominator) nd
sn
| num == 0 = []
Line 932 ⟶ 1,318:
| otherwise = intercalate "/" [show num, show den]
s = sr <> sn
in ifbool [(sn, sr)] [] (null s)
then []
else [(sn, sr)]
 
-- Rows of uneven length -> All rows padded to length of longest
fullTable :: [[(String, String)]] -> [[(String, String)]]
Line 941 ⟶ 1,325:
let lng = maximum $ length <$> xs
in (<>) <*> (flip replicate ([], []) . (-) lng . length) <$> xs
 
justifyLeft :: Int -> Char -> String -> String
justifyLeft n c s = take n (s <> replicate n c)
 
-- (Row index, Expression pairs) -> String joined by conjunctions
expressionRow :: (Int, [(String, String)]) -> String
Line 952 ⟶ 1,336:
, " -> "
, foldr
(\s a -> concat [s, bool " + " " " (blank a || head a == '-'), a])
concat[]
[ s
, if blank a || head a == '-'
then " "
else " + "
, a
])
""
(polyTerm <$> row)
]
 
-- (Power string, Ratio String) -> Combined string with possible '*'
polyTerm :: (String, String) -> String
Line 970 ⟶ 1,347:
| head r == '-' = concat ["- ", l, " * ", tail r]
| otherwise = intercalate " * " [l, r]
 
blank :: String -> Bool
blank = all isSpace
 
-- Maximum widths of power and ratio elements in each column
colWidths :: [[(String, String)]] -> [(Int, Int)]
Line 981 ⟶ 1,358:
(\(ls, rs) (lMax, rMax) -> (max (length ls) lMax, max (length rs) rMax))
(0, 0))
 
-- Length of string excluding any leading '-'
unsignedLength :: String -> Int
unsignedLength xs =
let l = length xs
in bool (bool l (l - 1) ('-' == head xs)) 0 (0 == l)</syntaxhighlight>
in case l of
0 -> 0
_ ->
case head xs of
'-' -> l - 1
_ -> l
 
-- TEST ------------------------------------------------------------------------
main :: IO ()
main = (putStrLn . unlines . expressionTable . take 10) polynomials</lang>
{{Out}}
<pre>0 -> n
Line 1,012 ⟶ 1,380:
Implementation:
 
<langsyntaxhighlight Jlang="j">Bsecond=:verb define"0
+/,(<:*(_1^[)*!*(y^~1+[)%1+])"0/~i.1x+y
)
Line 1,020 ⟶ 1,388:
Faul=:adverb define
(0,|.(%m+1x) * (_1x&^ * !&(m+1) * Bfirst) i.1+m)&p.
)</langsyntaxhighlight>
 
Task example:
 
<langsyntaxhighlight Jlang="j"> 0 Faul
0 1x&p.
1 Faul
Line 1,043 ⟶ 1,411:
0 _1r30 0 2r9 0 _7r15 0 2r3 1r2 1r9&p.
9 Faul
0 0 _3r20 0 1r2 0 _7r10 0 3r4 1r2 1r10&p.</langsyntaxhighlight>
 
Double checking our work:
 
<langsyntaxhighlight Jlang="j"> Fcheck=: dyad def'+/(1+i.y)^x'"0
9 Faul i.5
0 1 513 20196 282340
Line 1,055 ⟶ 1,423:
0 1 5 14 30
2 Fcheck i.5
0 1 5 14 30</langsyntaxhighlight>
 
=={{header|Java}}==
{{trans|Kotlin}}
{{works with|Java|8}}
<langsyntaxhighlight Javalang="java">import java.util.Arrays;
import java.util.stream.IntStream;
 
Line 1,201 ⟶ 1,569:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>0 : n
Line 1,213 ⟶ 1,581:
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|jq}}==
{{works with|jq}}
'''Works with gojq, the Go implementation of jq, and with fq'''
 
The following assumes the following rational arithmetic functions
as provided by the jq "rational" module at [[Arithmetic/Rational#jq]]:
 
* r/2 (constructor)
* requal/2 (equality)
* rlessthan/1
* rminus/1
* rminus/2 (subtraction)
* rmult/0 and rmult/2 (multiplication)
 
In addition, the pretty-printing function defined here (rpp) assumes
rationals are JSON objects of the form {n,d}
<syntaxhighlight lang=jq>
include "rational" {search: "."}; # see [[Arithmetic/Rational#jq]]:
 
# Preliminaries
# for gojq
def idivide($j):
. as $i
| ($i % $j) as $mod
| ($i - $mod) / $j ;
 
# use idivide for precision
def binomial(n; k):
if k > n / 2 then binomial(n; n-k)
else reduce range(1; k+1) as $i (1; . * (n - $i + 1) | idivide($i))
end;
 
# pretty print a Rational assumed to have the {n,d} form
def rpp:
if .n == 0 then "0"
elif .d == 1 then .n | tostring
else "\(.n)/\(.d)"
end;
 
# The following definition reflects the "modern view" that B(1) is 1 // 2
def bernoulli:
if type != "number" or . < 0 then "bernoulli must be given a non-negative number vs \(.)" | error
else . as $n
| reduce range(0; $n+1) as $i ([];
.[$i] = r(1; $i + 1)
| reduce range($i; 0; -1) as $j (.;
.[$j-1] = rmult($j; rminus(.[$j-1]; .[$j])) ) )
| .[0] # the modern view
end;
 
 
# The task
def faulhaber($p):
# The traditional view:
def bernouilli($n):
$n | bernouilli | if $n==1 then rminus else . end;
r(1; $p+1) as $q
| { sign: -1 }
| reduce range(0; 1+$p) as $j (.;
.sign *= -1
| r(binomial($p+1; $j); 1) as $b
| ([$q, .sign, $b, ($j|bernoulli)] | rmult) as $coeff
| if requal($coeff; r(0;1))|not
then .emit +=
(if $j == 0
then (if requal($coeff; r( 1;1)) then ""
elif requal($coeff; r(-1;1)) then "-"
else "\($coeff|rpp) " end)
else (if requal($coeff; r(1;1)) then " + "
elif requal($coeff; r(-1;1)) then " - "
elif r(0;1)|rlessthan($coeff) then " + \($coeff|rpp) "
else " - \($coeff|rminus|rpp) "
end)
end)
| ($p + 1 - $j) as $pwr
| .emit += (if 1 < $pwr then "n^\($pwr)" else "n" end)
else .
end )
| .emit ;
 
range(0;10) | "\(.) : \(faulhaber(.))"
</syntaxhighlight>
{{output}}
<pre>
0 : n
1 : 1/2 n^2 - 1/2 n
2 : 1/3 n^3 - 1/2 n^2 + 1/6 n
3 : 1/4 n^4 - 1/2 n^3 + 1/4 n^2
4 : 1/5 n^5 - 1/2 n^4 + 1/3 n^3 - 1/30 n
5 : 1/6 n^6 - 1/2 n^5 + 5/12 n^4 - 1/12 n^2
6 : 1/7 n^7 - 1/2 n^6 + 1/2 n^5 - 1/6 n^3 + 1/42 n
7 : 1/8 n^8 - 1/2 n^7 + 7/12 n^6 - 7/24 n^4 + 1/12 n^2
8 : 1/9 n^9 - 1/2 n^8 + 2/3 n^7 - 7/15 n^5 + 2/9 n^3 - 1/30 n
9 : 1/10 n^10 - 1/2 n^9 + 3/4 n^8 - 7/10 n^6 + 1/2 n^4 - 3/20 n^2
</pre>
 
=={{header|Julia}}==
Line 1,218 ⟶ 1,683:
 
'''Module''':
<langsyntaxhighlight lang="julia">module Faulhaber
 
function bernoulli(n::Integer)
Line 1,258 ⟶ 1,723:
end
 
end # module Faulhaber</langsyntaxhighlight>
 
'''Main''':
<syntaxhighlight lang ="julia">Faulhaber.formula.(1:10)</langsyntaxhighlight>
 
{{out}}
Line 1,277 ⟶ 1,742:
=={{header|Kotlin}}==
As Kotlin doesn't have support for rational numbers built in, a cut-down version of the Frac class in the Arithmetic/Rational task has been used in order to express the polynomial coefficients as fractions.
<langsyntaxhighlight lang="scala">// version 1.1.2
 
fun gcd(a: Long, b: Long): Long = if (b == 0L) a else gcd(b, a % b)
Line 1,395 ⟶ 1,860:
fun main(args: Array<String>) {
for (i in 0..9) faulhaber(i)
}</langsyntaxhighlight>
 
{{out}}
Line 1,410 ⟶ 1,875:
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2
</pre>
 
=={{header|Lua}}==
{{trans|C}}
<syntaxhighlight lang="lua">function binomial(n,k)
if n<0 or k<0 or n<k then return -1 end
if n==0 or k==0 then return 1 end
 
local num = 1
for i=k+1,n do
num = num * i
end
 
local denom = 1
for i=2,n-k do
denom = denom * i
end
 
return num / denom
end
 
function gcd(a,b)
while b ~= 0 do
local temp = a % b
a = b
b = temp
end
return a
end
 
function makeFrac(n,d)
local result = {}
 
if d==0 then
result.num = 0
result.denom = 0
return result
end
 
if n==0 then
d = 1
elseif d < 0 then
n = -n
d = -d
end
 
local g = math.abs(gcd(n, d))
if g>1 then
n = n / g
d = d / g
end
 
result.num = n
result.denom = d
return result
end
 
function negateFrac(f)
return makeFrac(-f.num, f.denom)
end
 
function subFrac(lhs, rhs)
return makeFrac(lhs.num * rhs.denom - lhs.denom * rhs.num, rhs.denom * lhs.denom)
end
 
function multFrac(lhs, rhs)
return makeFrac(lhs.num * rhs.num, lhs.denom * rhs.denom)
end
 
function equalFrac(lhs, rhs)
return (lhs.num == rhs.num) and (lhs.denom == rhs.denom)
end
 
function lessFrac(lhs, rhs)
return (lhs.num * rhs.denom) < (rhs.num * lhs.denom)
end
 
function printFrac(f)
io.write(f.num)
if f.denom ~= 1 then
io.write("/"..f.denom)
end
return nil
end
 
function bernoulli(n)
if n<0 then
return {num=0, denom=0}
end
 
local a = {}
for m=0,n do
a[m] = makeFrac(1, m+1)
for j=m,1,-1 do
a[j-1] = multFrac(subFrac(a[j-1], a[j]), makeFrac(j, 1))
end
end
 
if n~=1 then
return a[0]
end
return negateFrac(a[0])
end
 
function faulhaber(p)
io.write(p.." : ")
local q = makeFrac(1, p+1)
local sign = -1
for j=0,p do
sign = -1 * sign
local coeff = multFrac(multFrac(multFrac(q, makeFrac(sign, 1)), makeFrac(binomial(p + 1, j), 1)), bernoulli(j))
if not equalFrac(coeff, makeFrac(0, 1)) then
if j==0 then
if not equalFrac(coeff, makeFrac(1, 1)) then
if equalFrac(coeff, makeFrac(-1, 1)) then
io.write("-")
else
printFrac(coeff)
end
end
else
if equalFrac(coeff, makeFrac(1, 1)) then
io.write(" + ")
elseif equalFrac(coeff, makeFrac(-1, 1)) then
io.write(" - ")
elseif lessFrac(makeFrac(0, 1), coeff) then
io.write(" + ")
printFrac(coeff)
else
io.write(" - ")
printFrac(negateFrac(coeff))
end
end
 
local pwr = p + 1 - j
if pwr>1 then
io.write("n^"..pwr)
else
io.write("n")
end
end
end
print()
return nil
end
 
-- main
for i=0,9 do
faulhaber(i)
end</syntaxhighlight>
{{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|Mathematica}} / {{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">ClearAll[Faulhaber]
Faulhaber[n_, 0] := n
Faulhaber[n_, p_] := n^(p + 1)/(p + 1) + 1/2 n^p + Sum[BernoulliB[k]/k! p!/(p - k + 1)! n^(p - k + 1), {k, 2, p}]
Table[{p, Faulhaber[n, p]}, {p, 0, 9}] // Grid</syntaxhighlight>
{{out}}
<pre>0 n
1 n/2+n^2/2
2 n/6+n^2/2+n^3/3
3 n^2/4+n^3/2+n^4/4
4 -(n/30)+n^3/3+n^4/2+n^5/5
5 -(n^2/12)+(5 n^4)/12+n^5/2+n^6/6
6 n/42-n^3/6+n^5/2+n^6/2+n^7/7
7 n^2/12-(7 n^4)/24+(7 n^6)/12+n^7/2+n^8/8
8 -(n/30)+(2 n^3)/9-(7 n^5)/15+(2 n^7)/3+n^8/2+n^9/9
9 -((3 n^2)/20)+n^4/2-(7 n^6)/10+(3 n^8)/4+n^9/2+n^10/10</pre>
 
=={{header|Maxima}}==
<langsyntaxhighlight lang="maxima">sum1(p):=sum(stirling2(p,k)*pochhammer(n-k+1,k+1)/(k+1),k,0,p)$
sum2(p):=sum((-1)^j*binomial(p+1,j)*bern(j)*n^(p-j+1),j,0,p)/(p+1)$
 
Line 1,418 ⟶ 2,060:
[0,0,0,0,0,0,0,0,0,0]
 
for p from 0 thru 9 do print(expand(sum2(p)));</langsyntaxhighlight>
{{out}}
<pre>
Line 1,435 ⟶ 2,077:
=={{header|Modula-2}}==
{{trans|C#}}
<langsyntaxhighlight lang="modula2">MODULE Faulhaber;
FROM EXCEPTIONS IMPORT AllocateSource,ExceptionSource,GetMessage,RAISE;
FROM FormatString IMPORT FormatString;
Line 1,616 ⟶ 2,258:
END;
ReadChar
END Faulhaber.</langsyntaxhighlight>
{{out}}
<pre>0 : n
Line 1,628 ⟶ 2,270:
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|Nim}}==
{{trans|Kotlin}}
<syntaxhighlight lang="nim">import math, rationals
 
type
Fraction = Rational[int]
FaulhaberSequence = seq[Fraction]
 
const
Zero = 0 // 1
One = 1 // 1
MinusOne = -1 // 1
Powers = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"]
 
#---------------------------------------------------------------------------------------------------
 
func bernoulli(n: Natural): Fraction =
## Return nth Bernoulli coefficient.
 
var a = newSeq[Fraction](n + 1)
for m in 0..n:
a[m] = 1 // (m + 1)
for k in countdown(m, 1):
a[k - 1] = (a[k - 1] - a[k]) * k
result = if n != 1: a[0] else: -a[0]
 
#---------------------------------------------------------------------------------------------------
 
func faulhaber(n: Natural): FaulhaberSequence =
## Return nth Faulhaber sequence.
 
var a = 1 // (n + 1)
var sign = -1
for k in 0..n:
sign = -sign
result.add(a * sign * binom(n + 1, k) * bernoulli(k))
 
#---------------------------------------------------------------------------------------------------
 
func npower(k: Natural): string =
## Return the string representing "n" at power "k".
 
if k == 0: return ""
if k == 1: return "n"
var k = k
result = "n"
while k != 0:
result.insert(Powers[k mod 10], 1)
k = k div 10
 
#---------------------------------------------------------------------------------------------------
 
func `$`(fs: FaulhaberSequence): string =
## Return the string representing a Faulhaber sequence.
for i, coeff in fs:
 
# Process coefficient.
if coeff.num == 0: continue
if i == 0:
if coeff == MinusOne: result.add(" - ")
elif coeff != One: result.add($coeff)
else:
if coeff == One: result.add(" + ")
elif coeff == MinusOne: result.add(" - ")
elif coeff > Zero: result.add(" + " & $coeff)
else: result.add(" - " & $(-coeff))
 
# Process power of "n".
let pwr = fs.len - i
result.add(npower(pwr))
 
#———————————————————————————————————————————————————————————————————————————————————————————————————
 
for n in 0..9:
echo n, ": ", faulhaber(n)</syntaxhighlight>
 
{{out}}
<pre>0: n
1: 1/2n² + 1/2n
2: 1/3n³ + 1/2n² + 1/6n
3: 1/4n⁴ + 1/2n³ + 1/4n²
4: 1/5n⁵ + 1/2n⁴ + 1/3n³ - 1/30n
5: 1/6n⁶ + 1/2n⁵ + 5/12n⁴ - 1/12n²
6: 1/7n⁷ + 1/2n⁶ + 1/2n⁵ - 1/6n³ + 1/42n
7: 1/8n⁸ + 1/2n⁷ + 7/12n⁶ - 7/24n⁴ + 1/12n²
8: 1/9n⁹ + 1/2n⁸ + 2/3n⁷ - 7/15n⁵ + 2/9n³ - 1/30n
9: 1/10n¹⁰ + 1/2n⁹ + 3/4n⁸ - 7/10n⁶ + 1/2n⁴ - 3/20n²</pre>
 
=={{header|PARI/GP}}==
Line 1,637 ⟶ 2,368:
It's not worth using Bernoulli numbers in PARI/GP, because too much cleaning if you are avoiding "dirty" (but correct) result.<br>
Note: Find ssubstr() function here on RC.
<langsyntaxhighlight lang="parigp">
\\ Using "Faulhaber's" formula based on Bernoulli numbers. aev 2/7/17
\\ In str string replace all occurrences of the search string ssrch with the replacement string srepl. aev 3/8/16
Line 1,669 ⟶ 2,400:
{\\ Testing:
for(i=0,9, Faulhaber2(i))}
</syntaxhighlight>
</lang>
{{Output}}
<pre>
Line 1,687 ⟶ 2,418:
This version is using the sums of pth powers formula from [[wp:Bernoulli_polynomials| Bernoulli polynomials]].
It has small, simple and clear code, and produces instant result.
<langsyntaxhighlight lang="parigp">
\\ Using a formula based on Bernoulli polynomials. aev 2/5/17
Faulhaber1(m)={
Line 1,698 ⟶ 2,429:
{\\ Testing:
for(i=0,9, Faulhaber1(i))}
</syntaxhighlight>
</lang>
{{Output}}
<pre>
Line 1,713 ⟶ 2,444:
> ##
*** last result computed in 0 ms
</pre>
 
=={{header|Pascal}}==
A console program that runs under Lazarus or Delphi. Does not make use of Bernoulli numbers.
<syntaxhighlight lang="pascal">
program Faulhaber;
 
{$IFDEF FPC} // Lazarus
{$MODE Delphi} // ensure Lazarus accepts Delphi-style code
{$ASSERTIONS+} // by default, Lazarus does not compile 'Assert' statements
{$ELSE} // Delphi
{$APPTYPE CONSOLE}
{$ENDIF}
 
uses SysUtils;
 
type TRational = record
Num, Den : integer; // where Den > 0 and Num, Den are coprime
end;
 
const
ZERO : TRational = ( Num: 0; Den : 1);
HALF : TRational = ( Num: 1; Den : 2);
 
// Construct rational a/b, assuming b > 0.
function Rational( const a, b : integer) : TRational;
var
t, x, y : integer;
begin
if b <= 0 then raise SysUtils.Exception.Create( 'Denominator must be > 0');
// Find HCF of a and b (Euclid's algorithm) and cancel it out.
x := Abs(a);
y := b;
while y <> 0 do begin
t := x mod y;
x := y;
y := t;
end;
result.Num := a div x;
result.Den := b div x
end;
 
function Prod( r, s : TRational) : TRational; // result := r*s
begin
result := Rational( r.Num*s.Num, r.Den*s.Den);
end;
 
procedure DecRat( var r : TRational;
const s : TRational); // r := r - s
begin
r := Rational( r.Num*s.Den - s.Num*r.Den, r.Den * s.Den);
end;
 
// Write a term such as ' - (7/10)n^6' to the console.
procedure WriteTerm( coeff : TRational;
index : integer;
printPlus : boolean);
begin
if Coeff.Num = 0 then exit;
with coeff do begin
if Num < 0 then Write(' - ')
else if printPlus then Write(' + ');
// Put brackets round a fractional coefficient
if (Den > 1) then Write('(');
// If coefficient is 1, don't write it
if (Den > 1) or (Abs(Num) > 1) then Write( Abs(Num));
// Write denominator if it's not 1
if (Den > 1) then Write('/', Den, ')');
end;
Write('n');
if index > 1 then Write('^', index);
end;
 
{-------------------------------------------------------------------------------
Main routine. Calculation of Faulhaber polynomials
F_p(n) = 1^p + 2^p + ... + n^p, p = 0, 1, ..., p_max
}
var
p_max : integer;
c : array of array of TRational;
i, j, p : integer;
coeff_of_n : TRational;
begin
// User types program name, optionally followed by maximum power p (defaults to 9)
if ParamCount = 0 then p_max := 9
else p_max := SysUtils.StrToInt( ParamStr(1));
 
// c[p, i] is coefficient of n^i in the polynomial F_p(n).
// Initialize all coefficients to 0.
SetLength( c, p_max + 1, p_max + 2);
for i := 0 to p_max do
for j := 0 to p_max + 1 do
c[i, j] := ZERO;
 
c[0, 1] := Rational(1, 1); // F_0(n) = n, special case
for p := 1 to p_max do begin
// Initialize calculation of coefficient of n, needed if p is even.
// If p is odd, still calculate it as a check on the working (should be 0).
// Calculation uses the fact that F_p(1) = 1.
coeff_of_n := Rational(1, 1);
 
c[p, p+1] := Rational(1, p + 1);
DecRat( coeff_of_n, c[p, p + 1]);
c[p, p] := HALF;
DecRat( coeff_of_n, c[p, p]);
i := p - 1;
while (i >= 2) do begin
c[p, i] := Prod( Rational(p, i), c[p - 1, i - 1]);
DecRat( coeff_of_n, c[p, i]);
dec(i, 2);
end;
if i = 1 then // p is even
c[p, 1] := coeff_of_n // = the Bernoulli number B_p
else // p is odd
Assert( coeff_of_n.Num = 0); // just checking
end; // for p
 
// Print the result
for p := 0 to p_max do begin
Write( 'F_', p, '(n) = ');
for j := p + 1 downto 1 do WriteTerm( c[p, j], j, j <= p);
WriteLn;
end;
end.
</syntaxhighlight>
{{out}}
<pre>
F_0(n) = n
F_1(n) = (1/2)n^2 + (1/2)n
F_2(n) = (1/3)n^3 + (1/2)n^2 + (1/6)n
F_3(n) = (1/4)n^4 + (1/2)n^3 + (1/4)n^2
F_4(n) = (1/5)n^5 + (1/2)n^4 + (1/3)n^3 - (1/30)n
F_5(n) = (1/6)n^6 + (1/2)n^5 + (5/12)n^4 - (1/12)n^2
F_6(n) = (1/7)n^7 + (1/2)n^6 + (1/2)n^5 - (1/6)n^3 + (1/42)n
F_7(n) = (1/8)n^8 + (1/2)n^7 + (7/12)n^6 - (7/24)n^4 + (1/12)n^2
F_8(n) = (1/9)n^9 + (1/2)n^8 + (2/3)n^7 - (7/15)n^5 + (2/9)n^3 - (1/30)n
F_9(n) = (1/10)n^10 + (1/2)n^9 + (3/4)n^8 - (7/10)n^6 + (1/2)n^4 - (3/20)n^2
</pre>
 
=={{header|Perl}}==
<langsyntaxhighlight lang="perl">use 5.014;
use Math::Algebra::Symbols;
 
Line 1,758 ⟶ 2,626:
foreach my $i (0 .. 9) {
say "$i: ", faulhaber_s_formula($i);
}</langsyntaxhighlight>
{{out}}
<pre>
Line 1,771 ⟶ 2,639:
8: -1/30 n+2/9 n^3-7/15 n^5+2/3 n^7+1/2 n^8+1/9 n^9
9: -3/20 n^2+1/2 n^4-7/10 n^6+3/4 n^8+1/2 n^9+1/10 n^10
</pre>
 
=={{header|Perl 6}}==
{{works with|Rakudo|2018.04.01}}
<lang perl6>sub bernoulli_number($n) {
 
return 1/2 if $n == 1;
return 0/1 if $n % 2;
 
my @A;
for 0..$n -> $m {
@A[$m] = 1 / ($m + 1);
for $m, $m-1 ... 1 -> $j {
@A[$j - 1] = $j * (@A[$j - 1] - @A[$j]);
}
}
 
return @A[0];
}
 
sub binomial($n, $k) {
$k == 0 || $n == $k ?? 1 !! binomial($n-1, $k-1) + binomial($n-1, $k);
}
 
sub faulhaber_s_formula($p) {
 
my @formula = gather for 0..$p -> $j {
take '('
~ join('/', (binomial($p+1, $j) * bernoulli_number($j)).Rat.nude)
~ ")*n^{$p+1 - $j}";
}
 
my $formula = join(' + ', @formula.grep({!m{'(0/1)*'}}));
 
$formula .= subst(rx{ '(1/1)*' }, '', :g);
$formula .= subst(rx{ '^1'» }, '', :g);
 
"1/{$p+1} * ($formula)";
}
 
for 0..9 -> $p {
say "f($p) = ", faulhaber_s_formula($p);
}</lang>
{{out}}
<pre>
f(0) = 1/1 * (n)
f(1) = 1/2 * (n^2 + n)
f(2) = 1/3 * (n^3 + (3/2)*n^2 + (1/2)*n)
f(3) = 1/4 * (n^4 + (2/1)*n^3 + n^2)
f(4) = 1/5 * (n^5 + (5/2)*n^4 + (5/3)*n^3 + (-1/6)*n)
f(5) = 1/6 * (n^6 + (3/1)*n^5 + (5/2)*n^4 + (-1/2)*n^2)
f(6) = 1/7 * (n^7 + (7/2)*n^6 + (7/2)*n^5 + (-7/6)*n^3 + (1/6)*n)
f(7) = 1/8 * (n^8 + (4/1)*n^7 + (14/3)*n^6 + (-7/3)*n^4 + (2/3)*n^2)
f(8) = 1/9 * (n^9 + (9/2)*n^8 + (6/1)*n^7 + (-21/5)*n^5 + (2/1)*n^3 + (-3/10)*n)
f(9) = 1/10 * (n^10 + (5/1)*n^9 + (15/2)*n^8 + (-7/1)*n^6 + (5/1)*n^4 + (-3/2)*n^2)
</pre>
 
=={{header|Phix}}==
{{trans|C#}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<lang Phix>include builtins\pfrac.e -- (0.8.0+)
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
 
<span style="color: #008080;">include</span> <span style="color: #000000;">builtins</span><span style="color: #0000FF;">\</span><span style="color: #000000;">pfrac</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> <span style="color: #000080;font-style:italic;">-- (0.8.0+)</span>
function bernoulli(integer n)
sequence a = {}
<span style="color: #008080;">function</span> <span style="color: #000000;">bernoulli</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
for m=0 to n do
<span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
a = append(a,{1,m+1})
<span style="color: #008080;">for</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
for j=m to 1 by -1 do
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">})</span>
a[j] = frac_mul({j,1},frac_sub(a[j+1],a[j]))
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">m</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
end for
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">frac_mul</span><span style="color: #0000FF;">({</span><span style="color: #000000;">j</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span><span style="color: #000000;">frac_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]))</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
if n!=1 then return a[1] end if
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
return frac_uminus(a[1])
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">frac_uminus</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function binomial(integer n, k)
if n<0 or k<0 or n<k then ?9/0 end if
<span style="color: #008080;">function</span> <span style="color: #000000;">binomial</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">)</span>
if n=0 or k=0 then return 1 end if
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">k</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">n</span><span style="color: #0000FF;"><</span><span style="color: #000000;">k</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
integer num = 1,
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
denom = 1
<span style="color: #004080;">integer</span> <span style="color: #000000;">num</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>
for i=k+1 to n do
<span style="color: #000000;">denom</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
num *= i
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
end for
<span style="color: #000000;">num</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">i</span>
for i=2 to n-k do
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
denom *= i
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">k</span> <span style="color: #008080;">do</span>
end for
<span style="color: #000000;">denom</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">i</span>
return num / denom
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">num</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">denom</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
procedure faulhaber(integer p)
string res = sprintf("%d : ", p)
<span style="color: #008080;">procedure</span> <span style="color: #000000;">faulhaber</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
frac q = {1, p+1}
<span style="color: #004080;">string</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"%d : "</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
for j=0 to p do
<span style="color: #000000;">frac</span> <span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}</span>
frac bj = bernoulli(j)
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">p</span> <span style="color: #008080;">do</span>
if frac_ne(bj,frac_zero) then
<span style="color: #000000;">frac</span> <span style="color: #000000;">bj</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">bernoulli</span><span style="color: #0000FF;">(</span><span style="color: #000000;">j</span><span style="color: #0000FF;">)</span>
frac coeff = frac_mul({binomial(p+1,j),p+1},bj)
<span style="color: #008080;">if</span> <span style="color: #000000;">frac_ne</span><span style="color: #0000FF;">(</span><span style="color: #000000;">bj</span><span style="color: #0000FF;">,</span><span style="color: #000000;">frac_zero</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
string s = frac_sprint(coeff)
<span style="color: #000000;">frac</span> <span style="color: #000000;">coeff</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">frac_mul</span><span style="color: #0000FF;">({</span><span style="color: #000000;">binomial</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">),</span><span style="color: #000000;">p</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span><span style="color: #000000;">bj</span><span style="color: #0000FF;">)</span>
if j=0 then
<span style="color: #004080;">string</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">frac_sprint</span><span style="color: #0000FF;">(</span><span style="color: #000000;">coeff</span><span style="color: #0000FF;">)</span>
if s="1" then
<span style="color: #008080;">if</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
s = ""
<span style="color: #008080;">if</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">=</span><span style="color: #008000;">"1"</span> <span style="color: #008080;">then</span>
end if
<span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">""</span>
else
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
if s[1]='-' then
<span s[1..1] style= "color: - #008080;">else</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]=</span><span style="color: #008000;">'-'</span> <span style="color: #008080;">then</span>
else
<span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">" - "</span>
s[1..0] = " + "
end if<span style="color: #008080;">else</span>
<span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">0</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">" + "</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
res &= s&"n"
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
integer pwr = p+1-j
<span style="color: #000000;">res</span> <span style="color: #0000FF;">&=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">&</span><span style="color: #008000;">"n"</span>
if pwr>1 then
<span style="color: #004080;">integer</span> <span style="color: #000000;">pwr</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">j</span>
res &= sprintf("^%d", pwr)
<span style="color: #008080;">if</span> <span style="color: #000000;">pwr</span><span style="color: #0000FF;">></span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
end if
<span style="color: #000000;">res</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"^%d"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pwr</span><span style="color: #0000FF;">)</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
printf(1,"%s\n",{res})
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end procedure
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">res</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
for i=0 to 9 do
faulhaber(i)
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">9</span> <span style="color: #008080;">do</span>
end for</lang>
<span style="color: #000000;">faulhaber</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Line 1,904 ⟶ 2,720:
</pre>
 
== {{header|Python}} ==
 
The following implementation does not use [https://en.wikipedia.org/wiki/Bernoulli_number Bernoulli numbers], but [https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind Stirling numbers of the second kind], based on the relation: <math>m^n=\sum_{k=0}^n S_n^k (m)_k=\sum_{k=0}^n S_n^k k!{m\choose k}</math>.
Line 1,913 ⟶ 2,729:
 
 
'''Note:''' a number of the formulae above are '''invisible''' to the majority of browsers, including ''Chrome'', ''IE/Edge'', ''Safari'' and ''Opera''. They may (subject to the installation of necessary frontsfonts) be visible to some ''Firefox'' installations.
 
 
<langsyntaxhighlight lang="python">from fractions import Fraction
 
def nextu(a):
Line 1,978 ⟶ 2,794:
 
for i, p in enumerate(sumpol(10)):
print(i, ":", polstr(p))</langsyntaxhighlight>
 
{{out}}
Line 1,996 ⟶ 2,812:
Racket will simplify rational numbers; if this code simplifies the expressions too much for your tastes (e.g. you like <code>1/1 * (n)</code>) then tweak the simplify... clauses to taste.
 
<langsyntaxhighlight lang="racket">#lang racket/base
 
(require racket/match
Line 2,052 ⟶ 2,868:
(for ((p (in-range 0 (add1 9))))
(printf "f(~a) = ~a~%" p (expression->infix-string (faulhaber p))))
</syntaxhighlight>
</lang>
 
{{out}}
Line 2,068 ⟶ 2,884:
</pre>
 
=={{header|SidefRaku}}==
(formerly Perl 6)
<lang ruby>func faulhaber_s_formula(p) {
{{works with|Rakudo|2018.04.01}}
<syntaxhighlight lang="raku" line>sub bernoulli_number($n) {
 
varreturn formula1/2 =if gather$n == {1;
return 0/1 if $n {% |j|2;
 
take "(#{binomial(p+1, j) * j.bernfrac -> as_rat})*n^#{p+1 - j}"
my } << 0..p@A;
for 0..$n -> $m {
@A[$m] = 1 / ($m + 1);
for $m, $m-1 ... 1 -> $j {
@A[$j - 1] = $j * (@A[$j - 1] - @A[$j]);
}
}
 
return @A[0];
formula.grep! { !.contains('(0)*') }.join!(' + ')
}
 
sub binomial($n, $k) {
formula -= /\(1\)\*/g
$k == 0 || $n == $k ?? 1 !! binomial($n-1, $k-1) + binomial($n-1, $k);
formula -= /\^1\b/g
}
formula.gsub!(/\(([^+]*?)\)/, { _ })
 
sub faulhaber_s_formula($p) {
"1/#{p + 1} * (#{formula})"
 
my @formula = gather for 0..$p -> $j {
take '('
~ join('/', (binomial($p+1, $j) * bernoulli_number($j)).Rat.nude)
~ ")*n^{$p+1 - $j}";
}
 
my $formula = join(' + ', @formula.grep({!m{'(0/1)*'}}));
 
$formula .= subst(rx{ '(1/1)*' }, '', :g);
$formula .= subst(rx{ '^1'» }, '', :g);
 
"1/{$p+1} * ($formula)";
}
 
for 0..9 -> $p {
{ |p|
printf(say "%2d:f($p) = %s\n", p, faulhaber_s_formula($p));
}</syntaxhighlight>
} << ^10</lang>
{{out}}
<pre>
f(0:) = 1/1 * (n)
f(1:) = 1/2 * (n^2 + n)
f(2:) = 1/3 * (n^3 + (3/2)*n^2 + (1/2)*n)
f(3:) = 1/4 * (n^4 + (2/1)*n^3 + n^2)
f(4:) = 1/5 * (n^5 + (5/2)*n^4 + (5/3)*n^3 + (-1/6)*n)
f(5:) = 1/6 * (n^6 + (3/1)*n^5 + (5/2)*n^4 + (-1/2)*n^2)
f(6:) = 1/7 * (n^7 + (7/2)*n^6 + (7/2)*n^5 + (-7/6)*n^3 + (1/6)*n)
f(7:) = 1/8 * (n^8 + (4/1)*n^7 + (14/3)*n^6 + (-7/3)*n^4 + (2/3)*n^2)
f(8:) = 1/9 * (n^9 + (9/2)*n^8 + (6/1)*n^7 + (-21/5)*n^5 + (2/1)*n^3 + (-3/10)*n)
f(9:) = 1/10 * (n^10 + (5/1)*n^9 + (15/2)*n^8 + (-7/1)*n^6 + (5/1)*n^4 + (-3/2)*n^2)
</pre>
 
=={{header|RPL}}==
By not simplifying the formulas, we can have a much cleaner code:
{{works with|HP|49}}
<lang ruby>func faulhaber_s_formula(p) {
≪ → p
"1/#{p + 1} * (" + gather {
≪ { |j|0
p 0 '''FOR''' m
take "#{binomial(p+1, j) * j.bernfrac -> as_rat}*n^#{p+1 - j}"
} << 0.. p 1 + m 1 + COMB
}.join(' + ') + ")" p m - IBERNOULLI
'''IF''' LASTARG 1 == '''THEN''' NEG '''END'''
* EVAL + 'n' *
-1 '''STEP'''
p 1 + / EXPAN
≫ ≫ '<span style="color:blue>FAULH</span>' STO
 
≪ P <span style="color:blue>FAULH</span> ≫ 'P' 0 9 1 SEQ
{{out}}
<pre>
1: {'n' '1/2n^2+1/2n' '1/3n^3+1/2n^2+1/6n' '1/4n^4+1/2n^3+1/4n^2' '1/5n^5+1/2n^4+1/3n^3-1/30n' '1/6n^6+1/2n^5+5/12n^4-1/12n^2' '1/7n^7+1/2n^6+1/2n^5-1/6n^3+1/42n' '1/8n^8+1/2n^7+7/12n^6-7/24n^4+1/12n^2' '1/9n^9+1/2n^8+2/3n^7-7/15n^5+2/9n^3-1/30n' '1/10n^10+1/2n^9+3/4n^8-7/10n^6+1/2n^4-3/20n^2'}
</pre>
 
=={{header|Ruby}}==
{{trans|C}}
<syntaxhighlight lang="ruby">def binomial(n,k)
if n < 0 or k < 0 or n < k then
return -1
end
if n == 0 or k == 0 then
return 1
end
 
num = 1
for i in k+1 .. n do
num = num * i
end
 
denom = 1
for i in 2 .. n-k do
denom = denom * i
end
 
return num / denom
end
 
def bernoulli(n)
if n < 0 then
raise "n cannot be less than zero"
end
 
a = Array.new(16)
for m in 0 .. n do
a[m] = Rational(1, m + 1)
for j in m.downto(1) do
a[j-1] = (a[j-1] - a[j]) * Rational(j)
end
end
 
if n != 1 then
return a[0]
end
return -a[0]
end
 
def faulhaber(p)
print("%d : " % [p])
q = Rational(1, p + 1)
sign = -1
for j in 0 .. p do
sign = -1 * sign
coeff = q * Rational(sign) * Rational(binomial(p+1, j)) * bernoulli(j)
if coeff == 0 then
next
end
if j == 0 then
if coeff != 1 then
if coeff == -1 then
print "-"
else
print coeff
end
end
else
if coeff == 1 then
print " + "
elsif coeff == -1 then
print " - "
elsif 0 < coeff then
print " + "
print coeff
else
print " - "
print -coeff
end
end
pwr = p + 1 - j
if pwr > 1 then
print "n^%d" % [pwr]
else
print "n"
end
end
print "\n"
end
 
def main
for i in 0 .. 9 do
faulhaber(i)
end
end
 
main()</syntaxhighlight>
{{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|Scala}}==
{{trans|Java}}
<syntaxhighlight lang="scala">import scala.math.Ordering.Implicits.infixOrderingOps
 
abstract class Frac extends Comparable[Frac] {
val num: BigInt
val denom: BigInt
 
def unary_-(): Frac = {
Frac(-num, denom)
}
 
def +(rhs: Frac): Frac = {
Frac(
num * rhs.denom + rhs.num * denom,
denom * rhs.denom
)
}
 
def -(rhs: Frac): Frac = {
Frac(
num * rhs.denom - rhs.num * denom,
denom * rhs.denom
)
}
 
def *(rhs: Frac): Frac = {
Frac(num * rhs.num, denom * rhs.denom)
}
 
override def compareTo(rhs: Frac): Int = {
val ln = num * rhs.denom
val rn = rhs.num * denom
ln.compare(rn)
}
 
def canEqual(other: Any): Boolean = other.isInstanceOf[Frac]
 
override def equals(other: Any): Boolean = other match {
case that: Frac =>
(that canEqual this) &&
num == that.num &&
denom == that.denom
case _ => false
}
 
override def hashCode(): Int = {
val state = Seq(num, denom)
state.map(_.hashCode()).foldLeft(0)((a, b) => 31 * a + b)
}
 
override def toString: String = {
if (denom == 1) {
return s"$num"
}
s"$num/$denom"
}
}
 
object Frac {
{ |p|
val ZERO: Frac = Frac(0)
printf("%2d: %s\n", p, faulhaber_s_formula(p))
val ONE: Frac = Frac(1)
} << ^10</lang>
 
def apply(n: BigInt): Frac = new Frac {
val num: BigInt = n
val denom: BigInt = 1
}
 
def apply(n: BigInt, d: BigInt): Frac = {
if (d == 0) {
throw new IllegalArgumentException("Parameter d may not be zero.")
}
 
var nn = n
var dd = d
 
if (nn == 0) {
dd = 1
} else if (dd < 0) {
nn = -nn
dd = -dd
}
 
val g = nn.gcd(dd)
if (g > 0) {
nn /= g
dd /= g
}
 
new Frac {
val num: BigInt = nn
val denom: BigInt = dd
}
}
}
 
object Faulhaber {
def bernoulli(n: Int): Frac = {
if (n < 0) {
throw new IllegalArgumentException("n may not be negative or zero")
}
 
val a = Array.fill(n + 1)(Frac.ZERO)
for (m <- 0 to n) {
a(m) = Frac(1, m + 1)
for (j <- m to 1 by -1) {
a(j - 1) = (a(j - 1) - a(j)) * Frac(j)
}
}
 
// returns 'first' Bernoulli number
if (n != 1) {
return a(0)
}
-a(0)
}
 
def binomial(n: Int, k: Int): Int = {
if (n < 0 || k < 0 || n < k) {
throw new IllegalArgumentException()
}
if (n == 0 || k == 0) {
return 1
}
val num = (k + 1 to n).product
val den = (2 to n - k).product
num / den
}
 
def faulhaber(p: Int): Unit = {
print(s"$p : ")
val q = Frac(1, p + 1)
var sign = -1
for (j <- 0 to p) {
sign *= -1
val coeff = q * Frac(sign) * Frac(binomial(p + 1, j)) * bernoulli(j)
if (Frac.ZERO != coeff) {
if (j == 0) {
if (Frac.ONE != coeff) {
if (-Frac.ONE == coeff) {
print('-')
} else {
print(coeff)
}
}
} else {
if (Frac.ONE == coeff) {
print(" + ")
} else if (-Frac.ONE == coeff) {
print(" - ")
} else if (coeff > Frac.ZERO) {
print(s" + ${coeff}")
} else {
print(s" - ${-coeff}")
}
}
val pwr = p + 1 - j
if (pwr > 1) {
print(s"n^${pwr}")
} else {
print('n')
}
}
}
println()
}
 
def main(args: Array[String]): Unit = {
for (i <- 0 to 9) {
faulhaber(i)
}
}
}</syntaxhighlight>
{{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|Sidef}}==
<syntaxhighlight lang="ruby">func faulhaber_formula(p) {
(p+1).of { |j|
Poly(p - j + 1 => 1) * bernoulli(j) * binomial(p+1, j)
}.sum / (p+1)
}
 
for p in (^10) {
printf("%2d: %s\n", p, faulhaber_formula(p))
}</syntaxhighlight>
{{out}}
<pre>
0: 1/1 * (1*n^1)x
1: 1/2 * (1*nx^2 + 1/2*n^1)x
2: 1/3 * (1*nx^3 + 31/2*nx^2 + 1/26*n^1)x
3: 1/4 * (1*nx^4 + 1/2*nx^3 + 1/4*nx^2 + 0*n^1)
4: 1/5 * (1*nx^5 + 51/2*nx^4 + 51/3*nx^3 +- 0*n^2 + -1/630*n^1)x
5: 1/6 * (1*nx^6 + 31/2*nx^5 + 5/212*nx^4 +- 0*n^3 + -1/212*nx^2 + 0*n^1)
6: 1/7 * (1*nx^7 + 71/2*nx^6 + 71/2*nx^5 +- 0*n^4 + -71/6*nx^3 + 0*n^2 + 1/642*n^1)x
7: 1/8 * (1*nx^8 + 41/2*nx^7 + 147/312*nx^6 +- 0*n^5 + -7/324*nx^4 + 0*n^3 + 21/312*nx^2 + 0*n^1)
8: 1/9 * (1*nx^9 + 91/2*nx^8 + 62/3*nx^7 +- 0*n^6 + -217/515*nx^5 + 0*n^4 + 2/9*nx^3 +- 0*n^2 + -31/1030*n^1)x
9: 1/10 * (1*nx^10 + 51/2*nx^9 + 153/24*nx^8 +- 0*n^7 + -7/10*nx^6 + 01/2*n^5 + 5*nx^4 +- 0*n^3 + -3/220*nx^2 + 0*n^1)
</pre>
 
=={{header|Visual Basic .NET}}==
{{trans|C#}}
<syntaxhighlight lang="vbnet">Module Module1
Function Gcd(a As Long, b As Long)
If b = 0 Then
Return a
End If
Return Gcd(b, a Mod b)
End Function
 
Class Frac
ReadOnly num As Long
ReadOnly denom As Long
 
Public Shared ReadOnly ZERO As New Frac(0, 1)
Public Shared ReadOnly ONE As New Frac(1, 1)
 
Public Sub New(n As Long, d As Long)
If d = 0 Then Throw New ArgumentException("d must not be zero")
Dim nn = n
Dim dd = d
If nn = 0 Then
dd = 1
ElseIf dd < 0 Then
nn = -nn
dd = -dd
End If
Dim g = Math.Abs(Gcd(nn, dd))
If g > 1 Then
nn /= g
dd /= g
End If
num = nn
denom = dd
End Sub
 
Public Shared Operator -(self As Frac) As Frac
Return New Frac(-self.num, self.denom)
End Operator
 
Public Shared Operator +(lhs As Frac, rhs As Frac) As Frac
Return New Frac(lhs.num * rhs.denom + lhs.denom * rhs.num, rhs.denom * lhs.denom)
End Operator
 
Public Shared Operator -(lhs As Frac, rhs As Frac) As Frac
Return lhs + -rhs
End Operator
 
Public Shared Operator *(lhs As Frac, rhs As Frac) As Frac
Return New Frac(lhs.num * rhs.num, lhs.denom * rhs.denom)
End Operator
 
Public Shared Operator <(lhs As Frac, rhs As Frac) As Boolean
Dim x = lhs.num / lhs.denom
Dim y = rhs.num / rhs.denom
Return x < y
End Operator
 
Public Shared Operator >(lhs As Frac, rhs As Frac) As Boolean
Dim x = lhs.num / lhs.denom
Dim y = rhs.num / rhs.denom
Return x > y
End Operator
 
Public Shared Operator =(lhs As Frac, rhs As Frac) As Boolean
Return lhs.num = rhs.num AndAlso lhs.denom = rhs.denom
End Operator
 
Public Shared Operator <>(lhs As Frac, rhs As Frac) As Boolean
Return lhs.num <> rhs.num OrElse lhs.denom <> rhs.denom
End Operator
 
Public Overloads Function Equals(obj As Object) As Boolean
Dim frac = CType(obj, Frac)
Return Not IsNothing(frac) AndAlso num = frac.num AndAlso denom = frac.denom
End Function
 
Public Overloads Function GetHashCode() As Integer
Dim hashCode = 1317992671
hashCode = hashCode * -1521134295 + num.GetHashCode()
hashCode = hashCode * -1521134295 + denom.GetHashCode()
Return hashCode
End Function
 
Public Overloads Function ToString() As String
If denom = 1 Then Return num.ToString()
Return String.Format("{0}/{1}", num, denom)
End Function
End Class
 
Function Bernoulli(n As Integer) As Frac
If n < 0 Then Throw New ArgumentException("n may not be negative or zero")
Dim a(n + 1) As Frac
For m = 0 To n
a(m) = New Frac(1, m + 1)
For j = m To 1 Step -1
a(j - 1) = (a(j - 1) - a(j)) * New Frac(j, 1)
Next
Next
'returns the first Bernoulli number
If n <> 1 Then Return a(0)
Return -a(0)
End Function
 
Function Binomial(n As Integer, k As Integer) As Integer
If n < 0 OrElse k < 0 OrElse n < k Then
Throw New ArgumentException()
End If
If n = 0 OrElse k = 0 Then
Return 1
End If
Dim num = 1
For i = k + 1 To n
num *= i
Next
Dim denom = 1
For i = 2 To n - k
denom *= i
Next
Return num / denom
End Function
 
Sub Faulhaber(p As Integer)
Console.Write("{0} : ", p)
Dim q As New Frac(1, p + 1)
Dim sign = -1
For j = 0 To p
sign *= -1
Dim coeff = q * New Frac(sign, 1) * New Frac(Binomial(p + 1, j), 1) * Bernoulli(j)
If Frac.ZERO = coeff Then Continue For
If j = 0 Then
If Frac.ONE <> coeff Then
If -Frac.ONE = coeff Then
Console.Write("-")
Else
Console.Write(coeff.ToString())
End If
End If
Else
If Frac.ONE = coeff Then
Console.Write(" + ")
ElseIf -Frac.ONE = coeff Then
Console.Write(" - ")
ElseIf Frac.ZERO < coeff Then
Console.Write(" + {0}", coeff.ToString())
Else
Console.Write(" - {0}", (-coeff).ToString())
End If
End If
Dim pwr = p + 1 - j
If pwr > 1 Then
Console.Write("n^{0}", pwr)
Else
Console.Write("n")
End If
Next
Console.WriteLine()
End Sub
 
Sub Main()
For i = 0 To 9
Faulhaber(i)
Next
End Sub
End Module</syntaxhighlight>
{{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|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-math}}
{{libheader|Wren-rat}}
<syntaxhighlight lang="wren">import "./math" for Int
import "./rat" for Rat
 
var bernoulli = Fn.new { |n|
if (n < 0) Fiber.abort("Argument must be non-negative")
var a = List.filled(n+1, null)
for (m in 0..n) {
a[m] = Rat.new(1, m+1)
var j = m
while (j >= 1) {
a[j-1] = (a[j-1] - a[j]) * Rat.new(j, 1)
j = j - 1
}
}
return (n != 1) ? a[0] : -a[0] // 'first' Bernoulli number
}
 
var binomial = Fn.new { |n, k|
if (n < 0 || k < 0) Fiber.abort("Arguments must be non-negative integers")
if (n < k) Fiber.abort("The second argument cannot be more than the first.")
if (n == k) return 1
var prod = 1
var i = n - k + 1
while (i <= n) {
prod = prod * i
i = i + 1
}
return prod / Int.factorial(k)
}
 
var faulhaber = Fn.new { |p|
System.write("%(p) : ")
var q = Rat.new(1, p+1)
var sign = -1
for (j in 0..p) {
sign = sign * -1
var b = Rat.new(binomial.call(p+1, j), 1)
var coeff = q * Rat.new(sign, 1) * b * bernoulli.call(j)
if (coeff != Rat.zero) {
if (j == 0) {
System.write((coeff == Rat.one) ? "" : (coeff == Rat.minusOne) ? "-" : "%(coeff)")
} else {
System.write((coeff == Rat.one) ? " + " : (coeff == Rat.minusOne) ? " - " :
(coeff > Rat.zero) ? " + %(coeff)" : " - %(-coeff)")
}
var pwr = p + 1 - j
System.write((pwr > 1) ? "n^%(pwr)" : "n")
}
}
System.print()
}
 
for (i in 0..9) faulhaber.call(i)</syntaxhighlight>
 
{{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>
 
Line 2,132 ⟶ 3,522:
{{libheader|GMP}} GNU Multiple Precision Arithmetic Library
Uses code from the Bernoulli numbers task (copied here).
<langsyntaxhighlight lang="zkl">var [const] BN=Import("zklBigNum"); // libGMP (GNU MP Bignum Library)
 
fcn faulhaberFormula(p){ //-->(Rational,Rational...)
[p..0,-1].pump(List(),'wrap(k){ B(k)*BN(p+1).binomial(k) })
.apply('*(Rational(1,p+1)))
}</langsyntaxhighlight>
<langsyntaxhighlight lang="zkl">foreach p in (10){
println("F(%d) --> %s".fmt(p,polyRatString(faulhaberFormula(p))))
}</langsyntaxhighlight>
<langsyntaxhighlight lang="zkl">class Rational{ // Weenie Rational class, can handle BigInts
fcn init(_a,_b){ var a=_a, b=_b; normalize(); }
fcn toString{
Line 2,164 ⟶ 3,554:
}
fcn __opDiv(n){ self(a*n.b,b*n.a) } // Rat / Rat
}</langsyntaxhighlight>
<langsyntaxhighlight lang="zkl">fcn B(N){ // calculate Bernoulli(n) --> Rational
var A=List.createLong(100,0); // aka static aka not thread safe
foreach m in (N+1){
Line 2,181 ⟶ 3,571:
if(str[0]=="+") str[1,*]; // leave leading space
else String("-",str[2,*]);
}</langsyntaxhighlight>
{{out}}
<pre>
Jjuanhdez
2,122

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/263237...362129"