Rosetta Code

Faulhaber's triangle: Difference between revisions

Page Discussion

Faulhaber's triangle (view source)

Revision as of 11:46, 4 December 2023

57,634 bytes added ,  5 months ago
m
→‎{{header|Wren}}: Minor tidy
(Added Fōrmulæ)
m (→‎{{header|Wren}}: Minor tidy)
 
(34 intermediate revisions by 18 users not shown)
Line 1:
{{draft task|Mathematics}}
Named after [https://en.wikipedia.org/wiki/Johann_Faulhaber Johann Faulhaber], the rows of Faulhaber's triangle are the coefficients of polynomials that represent sums of integer powers, which are extracted from Faulhaber's formula:
 
Line 36:
* [http://www.ww.ingeniousmathstat.org/sites/default/files/Torabi-Dashti-CMJ-2011.pdf Faulhaber's triangle (PDF)]
<br>
 
=={{header|ALGOL 68}}==
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
Using code from the Algol 68 samples for the [[Arithmetic/Rational]] and [[Bernoulli numbers]] tasks and the Algol W sample for the [[Evaluate binomial coefficients]] task.<br>
Note that in the Bernoulli numbers task, the Algol 68 sample returns -1/2 for B(1) - this is modified here so B(1) is 1/2.<br>
Assumes LONG LONG INT is long enough to calculate the 17th power sum, the default precision of LONG LONG INT in ALGOL 68G is large enough.
<syntaxhighlight lang="algol68">
BEGIN # show some rows of Faulhaber's triangle #
 
# utility operators #
OP LENGTH = ( STRING a )INT: ( UPB a - LWB a ) + 1;
PRIO PAD = 9;
OP PAD = ( INT width, STRING v )STRING: # left blank pad v to width #
IF LENGTH v >= width THEN v ELSE ( " " * ( width - LENGTH v ) ) + v FI;
 
MODE INTEGER = LONG LONG INT; # mode for FRAC numberator & denominator #
OP TOINTEGER = ( INT n )INTEGER: n; # force widening n to INTEGER #
 
# Code from the Arithmetic/Rational task #
MODE FRAC = STRUCT( INTEGER num #erator#, den #ominator#);
 
PROC gcd = (INTEGER a, b) INTEGER: # greatest common divisor #
(a = 0 | b |: b = 0 | a |: ABS a > ABS b | gcd(b, a MOD b) | gcd(a, b MOD a));
PROC lcm = (INTEGER a, b)INTEGER: # least common multiple #
a OVER gcd(a, b) * b;
PRIO // = 9; # higher then the ** operator #
OP // = (INTEGER num, den)FRAC: ( # initialise and normalise #
INTEGER common = gcd(num, den);
IF den < 0 THEN
( -num OVER common, -den OVER common)
ELSE
( num OVER common, den OVER common)
FI
);
 
OP + = (FRAC a, b)FRAC: (
INTEGER common = lcm(den OF a, den OF b);
FRAC result := ( common OVER den OF a * num OF a + common OVER den OF b * num OF b, common );
num OF result//den OF result
);
OP - = (FRAC a, b)FRAC: a + -b,
* = (FRAC a, b)FRAC: (
INTEGER num = num OF a * num OF b,
den = den OF a * den OF b;
INTEGER common = gcd(num, den);
(num OVER common) // (den OVER common)
);
OP - = (FRAC frac)FRAC: (-num OF frac, den OF frac);
 
# end code from the Arithmetic/Rational task #
 
# alternative // operator for standard size INT values #
OP // = (INT num, den)FRAC: TOINTEGER num // TOINTEGER den;
# returns a * b #
OP * = ( INT a, FRAC b )FRAC: ( num OF b * a ) // den OF b;
OP * = ( INTEGER a, FRAC b )FRAC: ( num OF b * a ) // den OF b;
# sets a to a + b and returns a #
OP +:= = ( REF FRAC a, FRAC b )FRAC: a := a + b;
# sets a to - a and returns a #
OP -=: = ( REF FRAC a )FRAC: BEGIN num OF a := - num OF a; a END;
 
# returns the nth Bernoulli number, n must be >= 0 #
OP BERNOULLI = ( INT n )FRAC:
IF n < 0
THEN # n is out of range # 0 // 1
ELSE # n is valid #
[ 0 : n ]FRAC a;
FOR m FROM 0 TO n DO
a[ m ] := 1 // ( m + 1 );
FOR j FROM m BY -1 TO 1 DO
a[ j - 1 ] := j * ( a[ j - 1 ] - a[ j ] )
OD
OD;
IF n = 1 THEN - a[ 0 ] ELSE a[ 0 ] FI
FI # BERNOULLI # ;
 
# returns n! / k! #
PROC factorial over factorial = ( INT n, k )INTEGER:
IF k > n THEN 0
ELIF k = n THEN 1
ELSE # k < n #
INTEGER f := 1;
FOR i FROM k + 1 TO n DO f *:= i OD;
f
FI # factorial over Factorial # ;
 
# returns n! #
PROC factorial = ( INT n )INTEGER:
BEGIN
INTEGER f := 1;
FOR i FROM 2 TO n DO f *:= i OD;
f
END # factorial # ;
 
# returns the binomial coefficient of (n k) #
PROC binomial coefficient = ( INT n, k )INTEGER:
IF n - k > k
THEN factorial over factorial( n, n - k ) OVER factorial( k )
ELSE factorial over factorial( n, k ) OVER factorial( n - k )
FI # binomial coefficient # ;
 
# returns a string representation of a #
OP TOSTRING = ( FRAC a )STRING:
whole( num OF a, 0 ) + IF den OF a = 1 THEN "" ELSE "/" + whole( den OF a, 0 ) FI;
 
# returns the pth row of Faulhaber's triangle #
OP FAULHABER = ( INT p )[]FRAC:
BEGIN
FRAC q := -1 // ( p + 1 );
[ 0 : p ]FRAC coeffs;
FOR j FROM 0 TO p DO
coeffs[ p - j ] := binomial coefficient( p + 1, j ) * BERNOULLI j * -=: q
OD;
coeffs
END # faulhaber # ;
 
FOR i FROM 0 TO 9 DO # show the triabngle's first 10 rows #
[]FRAC frow = FAULHABER i;
FOR j FROM LWB frow TO UPB frow DO
print( ( " ", 6 PAD TOSTRING frow[ j ] ) )
OD;
print( ( newline ) )
OD;
BEGIN # compute the sum of k^17 for k = 1 to 1000 using triangle row 18 #
[]FRAC frow = FAULHABER 17;
FRAC sum := 0 // 1;
INTEGER kn := 1;
FOR j FROM LWB frow TO UPB frow DO
VOID( sum +:= ( kn *:= 1000 ) * frow[ j ] )
OD;
print( ( TOSTRING sum, newline ) )
END
END
</syntaxhighlight>
{{out}}
<pre>
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000
</pre>
 
=={{header|C}}==
{{trans|C++}}
<langsyntaxhighlight lang="c">#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
Line 196 ⟶ 349:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre> 1
Line 209 ⟶ 362:
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10</pre>
 
=={{header|C++ sharp|C#}}==
{{trans|C#}}
Uses C++ 17
<lang cpp>#include <exception>
#include <iomanip>
#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);
}
 
friend std::ostream& operator<<(std::ostream&, const Frac&);
 
static Frac ZERO() {
return Frac(0, 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;
}
 
std::vector<Frac> faulhaberTraingle(int p) {
std::vector<Frac> coeffs;
 
for (int i = 0; i < p + 1; i++) {
coeffs.push_back(Frac::ZERO());
}
 
Frac q{ 1, p + 1 };
int sign = -1;
for (int j = 0; j <= p; j++) {
sign *= -1;
coeffs[p - j] = q * Frac(sign, 1) * Frac(binomial(p + 1, j), 1) * bernoulli(j);
}
 
return coeffs;
}
 
int main() {
using namespace std;
 
for (int i = 0; i < 10; i++) {
vector<Frac> coeffs = faulhaberTraingle(i);
for (auto it = coeffs.begin(); it != coeffs.end(); it++) {
cout << right << setw(5) << *it << " ";
}
cout << endl;
}
 
return 0;
}</lang>
{{out}}
<pre> 1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10</pre>
 
=={{header|C#|C sharp}}==
{{trans|Java}}
<langsyntaxhighlight lang="csharp">using System;
 
namespace FaulhabersTriangle {
Line 517 ⟶ 518:
}
}
}</langsyntaxhighlight>
{{out}}
<pre> 1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10</pre>
 
=={{header|C++}}==
{{trans|C#}}
Uses C++ 17
<syntaxhighlight lang="cpp">#include <exception>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <sstream>
#include <vector>
class Frac {
public:
 
Frac() : num(0), denom(1) {}
 
Frac(int n, int d) {
if (d == 0) {
throw std::runtime_error("d must not be zero");
}
 
int sign_of_d = d < 0 ? -1 : 1;
int g = std::gcd(n, d);
 
num = sign_of_d * n / g;
denom = sign_of_d * d / g;
}
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);
}
 
Frac operator*(int rhs) const {
return Frac(num * rhs, denom);
}
 
friend std::ostream& operator<<(std::ostream&, const Frac&);
private:
int num;
int 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 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]) * j;
}
}
// 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 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;
}
std::vector<Frac> faulhaberTraingle(int p) {
std::vector<Frac> coeffs(p + 1);
Frac q{ 1, p + 1 };
int sign = -1;
for (int j = 0; j <= p; j++) {
sign *= -1;
coeffs[p - j] = q * sign * binomial(p + 1, j) * bernoulli(j);
}
return coeffs;
}
int main() {
for (int i = 0; i < 10; i++) {
std::vector<Frac> coeffs = faulhaberTraingle(i);
for (auto frac : coeffs) {
std::cout << std::right << std::setw(5) << frac << " ";
}
std::cout << std::endl;
}
return 0;
}</syntaxhighlight>
{{out}}
<pre> 1
Line 532 ⟶ 671:
=={{header|D}}==
{{trans|Kotlin}}
<langsyntaxhighlight Dlang="d">import std.algorithm : fold;
import std.conv : to;
import std.exception : enforce;
Line 665 ⟶ 804:
}
writeln;
}</langsyntaxhighlight>
{{out}}
<pre> 1
Line 680 ⟶ 819:
=={{header|F_Sharp|F#}}==
===The Function===
<langsyntaxhighlight lang="fsharp">
// Generate Faulhaber's Triangle. Nigel Galloway: May 8th., 2018
let Faulhaber=let fN n = (1N - List.sum n)::n
Line 687 ⟶ 826:
yield! Faul t (b+1N)}
Faul [] 0N
</syntaxhighlight>
</lang>
===The Task===
<langsyntaxhighlight lang="fsharp">
Faulhaber |> Seq.take 10 |> Seq.iter (printfn "%A")
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 707 ⟶ 846:
 
=={{header|Factor}}==
<langsyntaxhighlight lang="factor">USING: kernel math math.combinatorics math.extras math.functions
math.ranges prettyprint sequences ;
 
Line 714 ⟶ 853:
[ [ nCk ] [ -1 swap ^ ] [ bernoulli ] tri * * * ] 2with map ;
 
10 [ faulhaber . ] each-integer</langsyntaxhighlight>
{{out}}
<pre>
Line 728 ⟶ 867:
{ 0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 }
</pre>
 
=={{header|Fōrmulæ}}==
 
In [https://wiki.formulae.org/Faulhaber this] page you can see the solution of this task.
 
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text ([http://wiki.formulae.org/Editing_F%C5%8Drmul%C3%A6_expressions more info]). Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation &mdash;i.e. XML, JSON&mdash; they are intended for transportation effects more than visualization and edition.
 
The option to show Fōrmulæ programs and their results is showing images. Unfortunately images cannot be uploaded in Rosetta Code.
 
=={{header|FreeBASIC}}==
{{libheader|GMP}}
<langsyntaxhighlight lang="freebasic">' version 12-08-2017
' compile with: fbc -s console
' uses GMP
Line 826 ⟶ 957:
Print : Print "hit any key to end program"
Sleep
End</langsyntaxhighlight>
{{out}}
<pre>The first 10 rows
Line 841 ⟶ 972:
 
56056972216555580111030077961944183400198333273050000</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 same as the task [[Faulhaber%27s_formula#F%C5%8Drmul%C3%A6|Faulhaber's formula]])
 
[[File:Fōrmulæ - Faulhaber 01.png]]
 
'''Excecise 1.''' To show the first 11 rows (the first is the 0 row) of Faulhaber's triangle:
 
[[File:Fōrmulæ - Faulhaber 02.png]]
 
[[File:Fōrmulæ - Faulhaber 03.png]]
 
In order to show the previous result as a triangle:
 
[[File:Fōrmulæ - Faulhaber 04.png]]
 
[[File:Fōrmulæ - Faulhaber 05.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 same as the task [[Faulhaber%27s_formula#F%C5%8Drmul%C3%A6|Faulhaber's formula]])
 
This function can be used for both symbolic or numeric computation of the polynomial:
 
[[File:Fōrmulæ - Faulhaber 06.png]]
 
'''Excecise 2.''' Using the 18th row of Faulhaber's triangle, compute the sum <math>\sum_{k=1}^{1000} k^{17}</math>
 
[[File:Fōrmulæ - Faulhaber 09.png]]
 
[[File:Fōrmulæ - Faulhaber 10.png]]
 
Verification:
 
[[File:Fōrmulæ - Faulhaber 11.png]]
 
[[File:Fōrmulæ - Faulhaber 10.png]]
 
=={{header|Go}}==
{{trans|Kotlin}}
Except that there is no need to roll our own Frac type when we can use the big.Rat type from the Go standard library.
<langsyntaxhighlight lang="go">package main
 
import (
Line 926 ⟶ 1,103:
}
fmt.Println(sum.RatString())
}</langsyntaxhighlight>
 
{{out}}
Line 943 ⟶ 1,120:
56056972216555580111030077961944183400198333273050000
</pre>
 
=={{header|Groovy}}==
{{trans|Java}}
<syntaxhighlight lang="groovy">import java.math.MathContext
import java.util.stream.LongStream
 
class FaulhabersTriangle {
private static final MathContext MC = new MathContext(256)
 
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)
 
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)
}
 
double toDouble() {
return (double) num / denom
}
 
BigDecimal toBigDecimal() {
return BigDecimal.valueOf(num).divide(BigDecimal.valueOf(denom), MC)
}
}
 
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 long binomial(int n, int k) {
if (n < 0 || k < 0 || n < k) throw new IllegalArgumentException()
if (n == 0 || k == 0) return 1
long num = LongStream.rangeClosed(k + 1, n).reduce(1, { a, b -> a * b })
long den = LongStream.rangeClosed(2, n - k).reduce(1, { acc, i -> acc * i })
return num / den
}
 
private static Frac[] faulhaberTriangle(int p) {
Frac[] coeffs = new Frac[p + 1]
Arrays.fill(coeffs, Frac.ZERO)
Frac q = new Frac(1, p + 1)
int sign = -1
for (int j = 0; j <= p; ++j) {
sign *= -1
coeffs[p - j] = q * new Frac(sign, 1) * new Frac(binomial(p + 1, j), 1) * bernoulli(j)
}
return coeffs
}
 
static void main(String[] args) {
for (int i = 0; i <= 9; ++i) {
Frac[] coeffs = faulhaberTriangle(i)
for (Frac coeff : coeffs) {
printf("%5s ", coeff)
}
println()
}
println()
// get coeffs for (k + 1)th row
int k = 17
Frac[] cc = faulhaberTriangle(k)
int n = 1000
BigDecimal nn = BigDecimal.valueOf(n)
BigDecimal np = BigDecimal.ONE
BigDecimal sum = BigDecimal.ZERO
for (Frac c : cc) {
np = np * nn
sum = sum.add(np * c.toBigDecimal())
}
println(sum.toBigInteger())
}
}</syntaxhighlight>
{{out}}
<pre> 1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 </pre>
 
=={{header|Haskell}}==
{{works with|GHC|78.106.34}}
<langsyntaxhighlight lang="haskell">import Data.Ratio (Ratio, numeratordenominator, denominatornumerator, (%))
 
import Control.Arrow ((&&&))
------------------------ FAULHABER -----------------------
 
faulhaber :: Int -> Rational -> Rational
faulhaber p n =
sum $
zipWith ((*) . (n ^)) [1 ..] (faulhaberTriangle !! p)
 
 
-- FAULHABER -------------------------------------------------------------------
-- Infinite list of rows of Faulhaber's triangle
faulhaberTriangle :: [[Rational]]
faulhaberTriangle =
tail $
scanl
( \rs n ->
let xs = zipWith ((*) . (n %)) [2 ..] rs
in 1 - sum xs : xs)
[] )
[0 .. []
[0 ..]
 
-- p -> n -> Sum of the p-th powers of the first n positive integers
faulhaber :: Int -> Rational -> Rational
faulhaber p n = sum (zipWith ((*) . (n ^)) [1 ..] (faulhaberTriangle !! p))
 
-- DISPLAY -------------------------------------------- TEST -------------------------
main :: IO ()
-- (Max numerator+denominator widths) -> Column width -> Filler -> Ratio -> String
main = do
justifyRatio :: (Int, Int) -> Int -> Char -> Rational -> String
let triangle = take 10 faulhaberTriangle
widths = maxWidths triangle
mapM_
putStrLn
[ unlines
( (justifyRatio widths 8 ' ' =<<)
<$> triangle
),
(show . numerator) (faulhaber 17 1000)
]
 
------------------------- DISPLAY ------------------------
 
justifyRatio ::
(Int, Int) -> Int -> Char -> Rational -> String
justifyRatio (wn, wd) n c nd =
go $
let [num, den] = [numerator, denominator] <*> [nd]
[numerator, denominator] <*> [nd]
w = max n (wn + wd + 2) -- Minimum column width, or more if specified.
where
in if 1 == den
-- Minimum column thenwidth, centeror wmore cif (show num)specified.
w = max else letn (q,wn r)+ =wd quotRem+ (w - 12) 2
go [num, in concatden]
| 1 == den = center [ justifyRight qw c (show num)
| otherwise , "/"=
, justifyLeftlet (q +, r) c= quotRem (showw den- 1) 2
in ]concat
[ justifyRight q c (show num),
"/",
justifyLeft (q + r) c (show den)
]
 
center, justifyLeft, justifyRight :: Int -> Chara -> String[a] -> String[a]
justifyLeft n c s = take n (s <> replicate n c)
center n c s =
let (q, r) = quotRem (n - length s) 2
in concat [replicate q c, s, replicate (q + r) c]
 
justifyRight :: Int -> a -> [a] -> [a]
justifyLeft n c s = take n (s ++ replicate n c)
justifyRight n c = (drop . length) <*> (replicate n c <>)
 
center :: Int -> a -> [a] -> [a]
justifyRight n c s = drop (length s) (replicate n c ++ s)
center n c s =
let (q, r) = quotRem (n - length s) 2
pad = replicate q c
in concat [pad, s, pad, replicate r c]
 
-- List of Ratios -> (Max numerator width, Max denominator width)
maxWidths :: [[Rational]] -> (Int, Int)
maxWidths xss =
let widest f xs = maximum $ fmap (length . show . f) xs
in ((,) . widest numerator &&&<*> widest denominator) $ concat xss
concat xss</syntaxhighlight>
 
-- TEST ------------------------------------------------------------------------
main :: IO ()
main = do
let triangle = take 10 faulhaberTriangle
widths = maxWidths triangle
mapM_
putStrLn
[ unlines ((justifyRatio widths 8 ' ' =<<) <$> triangle)
, (show . numerator) (faulhaber 17 1000)
]</lang>
{{Out}}
<pre> 1
Line 1,018 ⟶ 1,363:
 
56056972216555580111030077961944183400198333273050000</pre>
 
=={{header|J}}==
<pre>
faulhaberTriangle=: ([: %. [: x: (1 _2 (p.) 2 | +/~) * >:/~ * (!~/~ >:))@:i.
 
faulhaberTriangle 10
1 0 0 0 0 0 0 0 0 0
1r2 1r2 0 0 0 0 0 0 0 0
1r6 1r2 1r3 0 0 0 0 0 0 0
0 1r4 1r2 1r4 0 0 0 0 0 0
_1r30 0 1r3 1r2 1r5 0 0 0 0 0
0 _1r12 0 5r12 1r2 1r6 0 0 0 0
1r42 0 _1r6 0 1r2 1r2 1r7 0 0 0
0 1r12 0 _7r24 0 7r12 1r2 1r8 0 0
_1r30 0 2r9 0 _7r15 0 2r3 1r2 1r9 0
0 _3r20 0 1r2 0 _7r10 0 3r4 1r2 1r10
 
NB.--------------------------------------------------------------
NB. matrix inverse (%.) of the extended precision (x:) product of
 
(!~/~ >:)@:i. 5 NB. Pascal's triangle
1 1 0 0 0
1 2 1 0 0
1 3 3 1 0
1 4 6 4 1
1 5 10 10 5
 
NB. second task in progress
(>:/~)@: i. 5 NB. lower triangle
1 0 0 0 0
1 1 0 0 0
1 1 1 0 0
1 1 1 1 0
1 1 1 1 1
(1 _2 p. (2 | +/~))@:i. 5 NB. 1 + (-2) * (parity table)
1 _1 1 _1 1
_1 1 _1 1 _1
1 _1 1 _1 1
_1 1 _1 1 _1
1 _1 1 _1 1
</pre>
 
=={{header|Java}}==
{{trans|Kotlin}}
{{works with|Java|8}}
<langsyntaxhighlight Javalang="java">import java.math.BigDecimal;
import java.math.MathContext;
import java.util.Arrays;
Line 1,163 ⟶ 1,551:
System.out.println(sum.toBigInteger());
}
}</langsyntaxhighlight>
{{out}}
<pre> 1
Line 1,184 ⟶ 1,572:
 
(Further progress would entail implementing some hand-crafted representation of arbitrary precision integers – perhaps a bit beyond the intended scope of this task, and good enough motivation to use a different language)
<langsyntaxhighlight JavaScriptlang="javascript">(() => {
 
// Order of Faulhaber's triangle -> rows of Faulhaber's triangle
Line 1,462 ⟶ 1,850:
]
);
})();</langsyntaxhighlight>
{{Out}}
<pre> 1
Line 1,483 ⟶ 1,871:
faulhaber(4, 1000)
200500333333300</pre>
 
=={{header|jq}}==
'''Works with [[jq]]''' (*)
'''Works with gojq, the Go implementation of jq'''
 
The solution presented here requires a "rational arithmetic" package, such
as the "Rational" module at [[Arithmetic/Rational#jq]].
As a reminder, the jq directive for including such a package appears as the first line of
the program below.
 
Note also that the function `bernoulli` is defined here in a way that ensures B(1) is `1 // 2`.
(*) The C implementation of jq does not have sufficient numeric precision for the "extra credit" task.
<syntaxhighlight lang="jq">include "Rational";
 
# Preliminaries
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
 
# 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;
 
# Here we conform to 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;
 
# Input: a non-negative integer, $p
# Output: an array of Rationals corresponding to the
# Faulhaber coefficients for row ($p + 1) (counting the first row as row 1).
def faulhabercoeffs:
def altBernoulli: # adjust B(1) for this task
bernoulli as $b
| if . == 1 then rmult(-1; $b) else $b end;
. as $p
| r(1; $p + 1) as $q
| { coeffs: [], sign: -1 }
| reduce range(0; $p+1) as $j (.;
.sign *= -1
| binomial($p + 1; $j) as $b
| .coeffs[$p - $j] = ([ .sign, $q, $b, ($j|altBernoulli) ] | rmult))
| .coeffs
;
 
# Calculate the sum for ($k|faulhabercoeffs)
def faulhabersum($n; $k):
($k|faulhabercoeffs) as $coe
| reduce range(0;$k+1) as $i ({sum: 0, power: 1};
.power *= $n
| .sum = radd(.sum; rmult(.power; $coe[$i]))
)
| .sum;
 
# 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;
def testfaulhaber:
(range(0; 10) as $i
| ($i | faulhabercoeffs | map(rpp | lpad(6)) | join(" "))),
"\nfaulhabersum(1000; 17):",
(faulhabersum(1000; 17) | rpp) ;
testfaulhaber</syntaxhighlight>
{{out}}
<pre>
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
 
faulhabersum(1000; 17):
56056972216555580111030077961944183400198333273050000
</pre>
 
 
=={{header|Julia}}==
<langsyntaxhighlight lang="julia">function bernoulli(n)
A = Vector{Rational{BigInt}}(undef, n + 1)
for i in 0:n
Line 1,522 ⟶ 2,009:
 
testfaulhaber()
</langsyntaxhighlight>{{out}}
<pre>
1
Line 1,537 ⟶ 2,024:
56056972216555580111030077961944183400198333273050000
</pre>
 
 
=={{header|Kotlin}}==
Uses appropriately modified code from the Faulhaber's Formula task:
<langsyntaxhighlight lang="scala">// version 1.1.2
 
import java.math.BigDecimal
Line 1,661 ⟶ 2,147:
}
println(sum.toBigInteger())
}</langsyntaxhighlight>
 
{{out}}
Line 1,681 ⟶ 2,167:
=={{header|Lua}}==
{{trans|C}}
<langsyntaxhighlight 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
Line 1,803 ⟶ 2,289:
for i=0,9 do
faulhaber(i)
end</langsyntaxhighlight>
{{out}}
<pre> 1
Line 1,815 ⟶ 2,301:
-1/30 -0 2/9 -0 -7/15 -0 2/3 1/2 1/9
-0 -3/20 -0 1/2 -0 -7/10 -0 3/4 1/2 1/10</pre>
 
 
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">ClearAll[Faulhaber]
bernoulliB[1] := 1/2
bernoulliB[n_] := BernoulliB[n]
Faulhaber[n_, p_] := 1/(p + 1) Sum[Binomial[p + 1, j] bernoulliB[j] n^(p + 1 - j), {j, 0, p}]
Table[Rest@CoefficientList[Faulhaber[n, t], n], {t, 0, 9}] // Grid
Faulhaber[1000, 17]</syntaxhighlight>
{{out}}
<pre>1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-(1/30) 0 1/3 1/2 1/5
0 -(1/12) 0 5/12 1/2 1/6
1/42 0 -(1/6) 0 1/2 1/2 1/7
0 1/12 0 -(7/24) 0 7/12 1/2 1/8
-(1/30) 0 2/9 0 -(7/15) 0 2/3 1/2 1/9
0 -(3/20) 0 1/2 0 -(7/10) 0 3/4 1/2 1/10
56056972216555580111030077961944183400198333273050000</pre>
 
=={{header|Maxima|}}==
<syntaxhighlight lang="maxima">
faulhaber_fraction(n, k) :=
if n = 0 and k = 1 then 1
else if k >= 2 and k <= n + 1 then (n/k) * faulhaber_fraction(n-1, k-1)
else if k = 1 then 1 - sum(faulhaber_fraction(n, i), i, 2, n+1)
else 0$
 
faulhaber_row(n):=makelist(faulhaber_fraction(n,k),k,1,n+1)$
/* Example */
triangle_faulhaber_first_ten_rows:block(makelist(faulhaber_row(i),i,0,9),table_form(%%));
</syntaxhighlight>
[[File:Faulhaber.png|thumb|center]]
 
=={{header|Nim}}==
{{libheader|bignum}}
For the task, we could use the standard module “rationals” but for the extra task we need big numbers (and big rationals). We use third party module “bignum” for this purpose.
<syntaxhighlight lang="nim">import algorithm, math, strutils
import bignum
 
type FaulhaberSequence = seq[Rat]
 
#---------------------------------------------------------------------------------------------------
 
func bernoulli(n: Natural): Rat =
## Return nth Bernoulli coefficient.
 
var a = newSeq[Rat](n + 1)
for m in 0..n:
a[m] = newRat(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 (high degree first).
 
var a = newRat(1, n + 1)
var sign = -1
for k in 0..n:
sign = -sign
result.add(a * sign * binom(n + 1, k) * bernoulli(k))
 
#---------------------------------------------------------------------------------------------------
 
proc display(fs: FaulhaberSequence) =
## Return the string representing a Faulhaber sequence.
 
var str = ""
for i, coeff in reversed(fs):
str.addSep(" ", 0)
str.add(($coeff).align(6))
echo str
 
#---------------------------------------------------------------------------------------------------
 
func evaluate(fs: FaulhaberSequence; n: int): Rat =
## Evaluate the polynomial associated to a sequence for value "n".
 
result = newRat(0)
for coeff in fs:
result = result * n + coeff
result *= n
 
#———————————————————————————————————————————————————————————————————————————————————————————————————
 
for n in 0..9:
display(faulhaber(n))
 
echo ""
let fs18 = faulhaber(17) # 18th row.
echo fs18.evaluate(1000)</syntaxhighlight>
 
{{out}}
<pre> 1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
 
56056972216555580111030077961944183400198333273050000</pre>
 
=={{header|Pascal}}==
{{libheader|IntXLib4Pascal}}
A console application in Free Pascal, created with the Lazarus IDE.
 
Row numbering below is 0-based, so row r has r+1 elements. Rather than use a rational number type, the program scales up row r by (r+1)!, which means that all the entries are integers.
<syntaxhighlight lang="pascal">
program FaulhaberTriangle;
uses uIntX, uEnums, // units in the library IntXLib4Pascal
SysUtils;
 
// Convert a rational num/den to a string, right-justified in the given width.
// Before converting, remove any common factor of num and den.
// For this application we can assume den > 0.
function RationalToString( num, den : TIntX;
minWidth : integer) : string;
var
num1, den1, divisor : TIntX;
w : integer;
begin
divisor := TIntX.GCD( num, den);
// TIntx.Divide requires the caller to specifiy the division mode
num1 := TIntx.Divide( num, divisor, uEnums.dmClassic);
den1 := TIntx.Divide( den, divisor, uEnums.dmClassic);
result := num1.ToString;
if not den1.IsOne then result := result + '/' + den1.ToString;
w := minWidth - Length( result);
if (w > 0) then result := StringOfChar(' ', w) + result;
end;
 
// Main routine
const
r_MAX = 17;
var
g : array [1..r_MAX + 1] of TIntX;
r, s, k : integer;
r_1_fac, sum, k_intx : TIntX;
begin
// Calculate rows 0..17 of Faulhaner's triangle, and show rows 0..9.
// For a given r, the subarray g[1..(r+1)] contains (r + 1)! times row r.
r_1_fac := 1; // (r + 1)!
g[1] := 1;
for r := 0 to r_MAX do begin
r_1_fac := r_1_fac * (r+1);
sum := 0;
for s := r downto 1 do begin
g[s + 1] := r*(r+1)*g[s] div (s+1);
sum := sum + g[s + 1];
end;
g[1] := r_1_fac - sum; // the scaled row must sum to (r + 1)!
if (r <= 9) then begin
for s := 1 to r + 1 do Write( RationalToString( g[s], r_1_fac, 7));
WriteLn;
end;
end;
 
// Use row 17 to sum 17th powers from 1 to 1000
sum := 0;
for s := r_MAX + 1 downto 1 do sum := (sum + g[s]) * 1000;
sum := TIntx.Divide( sum, r_1_fac, uEnums.dmClassic);
WriteLn;
WriteLn( 'Sum by Faulhaber = ' + sum.ToString);
 
// Check by direct calculation
sum := 0;
for k := 1 to 1000 do begin
k_intx := k;
sum := sum + TIntX.Pow( k_intx, r_MAX);
end;
WriteLn( 'by direct calc. = ' + sum.ToString);
end.
</syntaxhighlight>
{{out}}
<pre>
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
 
Sum by Faulhaber = 56056972216555580111030077961944183400198333273050000
by direct calc. = 56056972216555580111030077961944183400198333273050000
</pre>
 
 
=={{header|Perl}}==
{{libheader|ntheory}}
<langsyntaxhighlight lang="perl">use 5.010;
use List::Util qw(sum);
use Math::BigRat try => 'GMP';
Line 1,841 ⟶ 2,526:
my $n = Math::BigInt->new(1000);
my @r = faulhaber_triangle($p+1);
say "\n", sum(map { $r[$_] * $n**($_ + 1) } 0 .. $#r);</langsyntaxhighlight>
{{out}}
<pre>
Line 1,857 ⟶ 2,542:
56056972216555580111030077961944183400198333273050000
</pre>
 
=={{header|Perl 6}}==
{{works with|Rakudo|2017.05}}
{{trans|Sidef}}
 
<lang perl6># Helper subs
 
sub infix:<reduce> (\prev, \this) { this.key => this.key * (this.value - prev.value) }
 
sub next-bernoulli ( (:key($pm), :value(@pa)) ) {
$pm + 1 => [ map *.value, [\reduce] ($pm + 2 ... 1) Z=> 1 / ($pm + 2), |@pa ]
}
 
constant bernoulli = (0 => [1.FatRat], &next-bernoulli ... *).map: { .value[*-1] };
 
sub binomial (Int $n, Int $p) { combinations($n, $p).elems }
 
sub asRat (FatRat $r) { $r ?? $r.denominator == 1 ?? $r.numerator !! $r.nude.join('/') !! 0 }
 
 
# The task
sub faulhaber_triangle ($p) { map { binomial($p + 1, $_) * bernoulli[$_] / ($p + 1) }, ($p ... 0) }
 
# First 10 rows of Faulhaber's triangle:
say faulhaber_triangle($_)».&asRat.fmt('%5s') for ^10;
say '';
 
# Extra credit:
my $p = 17;
my $n = 1000;
say sum faulhaber_triangle($p).kv.map: { $^value * $n**($^key + 1) }</lang>
{{out}}
<pre> 1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
 
56056972216555580111030077961944183400198333273050000</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>
atom 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;">atom</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>
function faulhaber_triangle(integer p, bool asString=true)
sequence coeffs = repeat(frac_zero,p+1)
<span style="color: #008080;">function</span> <span style="color: #000000;">faulhaber_triangle</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">bool</span> <span style="color: #000000;">asString</span><span style="color: #0000FF;">=</span><span style="color: #004600;">true</span><span style="color: #0000FF;">)</span>
for j=0 to p do
<span style="color: #004080;">sequence</span> <span style="color: #000000;">coeffs</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">frac_zero</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 coeff = frac_mul({binomial(p+1,j),p+1},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>
coeffs[p-j+1] = iff(asString?sprintf("%5s",{frac_sprint(coeff)}):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;">bernoulli</span><span style="color: #0000FF;">(</span><span style="color: #000000;">j</span><span style="color: #0000FF;">))</span>
end for
<span style="color: #000000;">coeffs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">p</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: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">asString</span><span style="color: #0000FF;">?</span><span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"%5s"</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><span style="color: #000000;">coeff</span><span style="color: #0000FF;">)</span>
return coeffs
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">coeffs</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
for i=0 to 9 do
printf(1,"%s\n",{join(faulhaber_triangle(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
<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: #7060A8;">join</span><span style="color: #0000FF;">(</span><span style="color: #000000;">faulhaber_triangle</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">),</span><span style="color: #008000;">" "</span><span style="color: #0000FF;">)})</span>
puts(1,"\n")
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
 
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
sequence row18 = faulhaber_triangle(17,false)
frac res = frac_zero
<span style="color: #008080;">if</span> <span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()!=</span><span style="color: #004600;">JS</span> <span style="color: #008080;">then</span>
atom t1 = time()+1
<span style="color: #004080;">sequence</span> <span style="color: #000000;">row18</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">faulhaber_triangle</span><span style="color: #0000FF;">(</span><span style="color: #000000;">17</span><span style="color: #0000FF;">,</span><span style="color: #004600;">false</span><span style="color: #0000FF;">)</span>
integer lim = 1000
<span style="color: #000000;">frac</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">frac_zero</span>
for k=1 to lim do
<span style="color: #004080;">atom</span> <span style="color: #000000;">t1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()+</span><span style="color: #000000;">1</span>
bigatom nn = BA_ONE
<span style="color: #004080;">integer</span> <span style="color: #000000;">lim</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1000</span>
for i=1 to length(row18) do
<span style="color: #008080;">for</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;">lim</span> <span style="color: #008080;">do</span>
res = frac_add(res,frac_mul(row18[i],{nn,1}))
<span style="color: #000000;">bigatom</span> <span style="color: #000000;">nn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">BA_ONE</span>
nn = ba_mul(nn,lim)
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">row18</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
end for
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">frac_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">frac_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">row18</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],{</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}))</span>
if time()>t1 then printf(1,"calculating, k=%d...\r",k) t1 = time()+1 end if
<span style="color: #000000;">nn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ba_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lim</span><span style="color: #0000FF;">)</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
printf(1,"%s \n",{frac_sprint(res)})</lang>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()></span><span style="color: #000000;">t1</span> <span style="color: #008080;">then</span> <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;">"calculating, k=%d...\r"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> <span style="color: #000000;">t1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()+</span><span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<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;">frac_sprint</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;">if</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Line 1,973 ⟶ 2,619:
 
56056972216555580111030077961944183400198333273050000
</pre>
Note the extra credit takes about 90s, so I disabled it under pwa/p2js.
 
=={{header|Prolog}}==
<syntaxhighlight lang="prolog">
ft_rows(Lz) :-
lazy_list(ft_row, [], Lz).
 
ft_row([], R1, R1) :- R1 = [1].
ft_row(R0, R2, R2) :-
length(R0, P),
Jmax is 1 + P, numlist(2, Jmax, Qs),
maplist(term(P), Qs, R0, R1),
sum_list(R1, S), Bk is 1 - S, % Bk is Bernoulli number
R2 = [Bk | R1].
 
term(P, Q, R, S) :- S is R * (P rdiv Q).
 
show(N) :-
ft_rows(Rs),
length(Rows, N), prefix(Rows, Rs),
forall(
member(R, Rows),
(format(string(S), "~w", [R]),
re_replace(" rdiv "/g, "/", S, T),
re_replace(","/g, ", ", T, U),
write(U), nl)).
 
sum(N, K, S) :- % sum I=1,N (I ** K)
ft_rows(Rows), drop(K, Rows, [Coefs|_]),
reverse([0|Coefs], Poly),
foldl(horner(N), Poly, 0, S).
 
horner(N, A, S0, S1) :-
S1 is N*S0 + A.
 
drop(N, Lz1, Lz2) :-
append(Pfx, Lz2, Lz1), length(Pfx, N), !.
</syntaxhighlight>
{{Out}}
<pre>
?- show(10).
[1]
[1/2, 1/2]
[1/6, 1/2, 1/3]
[0, 1/4, 1/2, 1/4]
[-1/30, 0, 1/3, 1/2, 1/5]
[0, -1/12, 0, 5/12, 1/2, 1/6]
[1/42, 0, -1/6, 0, 1/2, 1/2, 1/7]
[0, 1/12, 0, -7/24, 0, 7/12, 1/2, 1/8]
[-1/30, 0, 2/9, 0, -7/15, 0, 2/3, 1/2, 1/9]
[0, -3/20, 0, 1/2, 0, -7/10, 0, 3/4, 1/2, 1/10]
true.
 
?- sum(1000, 17, S).
S = 56056972216555580111030077961944183400198333273050000.
</pre>
 
Line 1,978 ⟶ 2,680:
{{trans|Haskell}}
{{Works with|Python|3.7}}
<langsyntaxhighlight lang="python">'''Faulhaber's triangle'''
 
from itertools import (accumulate, chain, count, islice, starmap)
from fractions import (Fraction)
 
 
Line 1,988 ⟶ 2,690:
'''List of rows of Faulhaber fractions.'''
def go(rs, n):
xsdef =f(x, list(starmap(y):
lambda x, y:return Fraction(n, x) * y,
xs = ziplist(map(f, islice(count(2), m), rs))
))
return [Fraction(1 - sum(xs), 1)] + xs
 
return list(accumulate(
[[]] + list(islice(count(0), 1 + m)),
Line 2,004 ⟶ 2,706:
positive integers.
'''
def go(x, y):
return y * (n ** x)
 
return sum(
map(go, count(1), faulhaberTriangle(p)[-1])
list(starmap(
lambda x, y: y * (n ** x),
zip(count(1), faulhaberTriangle(p)[-1])
))
)
 
 
# TEST --------------------------- TEST -------------------------
def main():
'''Tests'''
Line 2,032 ⟶ 2,734:
 
 
# DISPLAY ------------------------- DISPLAY ------------------------
 
# fTable :: String -> (a -> String) ->
# (b -> String) -> (a -> b) -> [a] -> String
def fTable(s):
'''Heading -> x display function -> fx display function ->
fx display function -> f -> xs -> tabular string.
'''
def gogox(xShow, fxShow, f, xs):
ysdef = [xShowgofx(xfxShow) for x in xs]:
w = max(map(len, ys) def gof(f):
return s + '\n' + '\n'.join(map def goxs(xs):
lambda x, y: y.rjust(w, ' ') + ' ->ys '= + fxShow(f[xShow(x)), for x in xs]
xs w = max(map(len, ys))
 
))
def arrowed(x, y):
return lambda xShow: lambda fxShow: lambda f: lambda xs: go(
return y.rjust(w, ' ') + ' -> ' + (
xShow, fxShow, f, xs
fxShow(f(x))
)
)
return s + '\n' + '\n'.join(
map(arrowed, xs, ys)
)
return goxs
return gof
return gofx
return gox
 
 
# GENERIC ------------------------- GENERIC ------------------------
 
# compose (<<<) :: (b -> c) -> (a -> b) -> a -> c
Line 2,064 ⟶ 2,774:
def concat(xs):
'''The concatenation of all the elements
in a list or iterable.'''
'''
def f(ys):
zs = list(chain(*ys))
Line 2,081 ⟶ 2,792:
f lifted to a function over a list.
'''
returndef lambda xs: listgo(map(f, xs)):
return list(map(f, xs))
 
return go
 
 
Line 2,099 ⟶ 2,813:
representation of the ratio r.
'''
def go(n, r):
d =def f(r.denominator):
return str(r.numerator).rjust(m, ' ') +d (= r.denominator
('/' +return str(dr.numerator).ljustrjust(nm, ' ')) if 1 != d else+ (
(' /' *+ str(1d).ljust(n, +' n')) if 1 != d else (
' ' * (1 + n)
)
)
)return f
return lambda n: lambda r: go(n, r)
 
 
# MAIN ---
if __name__ == '__main__':
main()</langsyntaxhighlight>
{{Out}}
<pre>Faulhaber's triangle:
Line 2,130 ⟶ 2,846:
=={{header|Racket}}==
 
<langsyntaxhighlight lang="racket">#lang racket
(require math/number-theory)
 
Line 2,153 ⟶ 2,869:
(require rackunit)
(check-equal? (sum-k^p:formulaic 17 1000)
(for/sum ((k (in-range 1 (add1 1000)))) (expt k 17))))</langsyntaxhighlight>
{{out}}
<pre>'((1) (1/2 1/2) (1/6 1/2 1/3) (0 1/4 1/2 1/4) (-1/30 0 1/3 1/2 1/5) (0 -1/12 0 5/12 1/2 1/6) (1/42 0 -1/6 0 1/2 1/2 1/7) (0 1/12 0 -7/24 0 7/12 1/2 1/8) (-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9) (0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10))
56056972216555580111030077961944183400198333273050000</pre>
 
=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2017.05}}
{{trans|Sidef}}
 
<syntaxhighlight lang="raku" line># Helper subs
 
sub infix:<reduce> (\prev, \this) { this.key => this.key * (this.value - prev.value) }
 
sub next-bernoulli ( (:key($pm), :value(@pa)) ) {
$pm + 1 => [ map *.value, [\reduce] ($pm + 2 ... 1) Z=> 1 / ($pm + 2), |@pa ]
}
 
constant bernoulli = (0 => [1.FatRat], &next-bernoulli ... *).map: { .value[*-1] };
 
sub binomial (Int $n, Int $p) { combinations($n, $p).elems }
 
sub asRat (FatRat $r) { $r ?? $r.denominator == 1 ?? $r.numerator !! $r.nude.join('/') !! 0 }
 
 
# The task
sub faulhaber_triangle ($p) { map { binomial($p + 1, $_) * bernoulli[$_] / ($p + 1) }, ($p ... 0) }
 
# First 10 rows of Faulhaber's triangle:
say faulhaber_triangle($_)».&asRat.fmt('%5s') for ^10;
say '';
 
# Extra credit:
my $p = 17;
my $n = 1000;
say sum faulhaber_triangle($p).kv.map: { $^value * $n**($^key + 1) }</syntaxhighlight>
{{out}}
<pre> 1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
 
56056972216555580111030077961944183400198333273050000</pre>
 
=={{header|REXX}}==
<langsyntaxhighlight lang="rexx">Numeric Digits 100
Do r=0 To 20
ra=r-1
Line 2,257 ⟶ 3,018:
Parse Arg a,b
if b = 0 then return abs(a)
return gcd(b,a//b)</langsyntaxhighlight>
{{out}}
<pre> 1
Line 2,272 ⟶ 3,033:
56056972216555580111030077961944183400198333273050000
56056972216555580111030077961944183400198333273050000</pre>
 
=={{header|Ruby}}==
{{trans|D}}
<syntaxhighlight lang="ruby">class Frac
attr_accessor:num
attr_accessor:denom
 
def initialize(n,d)
if d == 0 then
raise ArgumentError.new('d cannot be zero')
end
 
nn = n
dd = d
if nn == 0 then
dd = 1
elsif dd < 0 then
nn = -nn
dd = -dd
end
 
g = nn.abs.gcd(dd.abs)
if g > 1 then
nn = nn / g
dd = dd / g
end
 
@num = nn
@denom = dd
end
 
def to_s
if self.denom == 1 then
return self.num.to_s
else
return "%d/%d" % [self.num, self.denom]
end
end
 
def -@
return Frac.new(-self.num, self.denom)
end
 
def +(rhs)
return Frac.new(self.num * rhs.denom + self.denom * rhs.num, rhs.denom * self.denom)
end
def -(rhs)
return Frac.new(self.num * rhs.denom - self.denom * rhs.num, rhs.denom * self.denom)
end
 
def *(rhs)
return Frac.new(self.num * rhs.num, rhs.denom * self.denom)
end
end
 
FRAC_ZERO = Frac.new(0, 1)
FRAC_ONE = Frac.new(1, 1)
 
def bernoulli(n)
if n < 0 then
raise ArgumentError.new('n cannot be negative')
end
 
a = Array.new(n + 1)
a[0] = FRAC_ZERO
 
for m in 0 .. n do
a[m] = Frac.new(1, m + 1)
m.downto(1) do |j|
a[j - 1] = (a[j - 1] - a[j]) * Frac.new(j, 1)
end
end
 
if n != 1 then
return a[0]
end
return -a[0]
end
 
def binomial(n, k)
if n < 0 then
raise ArgumentError.new('n cannot be negative')
end
if k < 0 then
raise ArgumentError.new('k cannot be negative')
end
if n < k then
raise ArgumentError.new('n cannot be less than k')
end
 
if n == 0 or k == 0 then
return 1
end
 
num = 1
for i in k + 1 .. n do
num = num * i
end
 
den = 1
for i in 2 .. n - k do
den = den * i
end
 
return num / den
end
 
def faulhaberTriangle(p)
coeffs = Array.new(p + 1)
coeffs[0] = FRAC_ZERO
q = Frac.new(1, p + 1)
sign = -1
for j in 0 .. p do
sign = -sign
coeffs[p - j] = q * Frac.new(sign, 1) * Frac.new(binomial(p + 1, j), 1) * bernoulli(j)
end
return coeffs
end
 
def main
for i in 0 .. 9 do
coeffs = faulhaberTriangle(i)
coeffs.each do |coeff|
print "%5s " % [coeff]
end
puts
end
end
 
main()</syntaxhighlight>
{{out}}
<pre> 1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10</pre>
 
=={{header|Scala}}==
{{trans|Java}}
<syntaxhighlight lang="scala">import java.math.MathContext
 
import scala.collection.mutable
 
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 {
val ZERO: Frac = Frac(0)
val ONE: Frac = Frac(1)
 
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 faulhaberTriangle(p: Int): List[Frac] = {
val coeffs = mutable.MutableList.fill(p + 1)(Frac.ZERO)
 
val q = Frac(1, p + 1)
var sign = -1
for (j <- 0 to p) {
sign *= -1
coeffs(p - j) = q * Frac(sign) * Frac(binomial(p + 1, j)) * bernoulli(j)
}
coeffs.toList
}
 
def main(args: Array[String]): Unit = {
for (i <- 0 to 9) {
val coeffs = faulhaberTriangle(i)
for (coeff <- coeffs) {
print("%5s ".format(coeff))
}
println()
}
println()
}
}</syntaxhighlight>
{{out}}
<pre> 1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 </pre>
=={{header|Scheme}}==
{{works with|Chez Scheme}}
<syntaxhighlight lang="scheme">; Return the first row-count rows of Faulhaber's Triangle as a vector of vectors.
(define faulhabers-triangle
(lambda (row-count)
; Calculate and store the value of the first column of a row.
; The value is one minus the sum of all the rest of the columns.
(define calc-store-first!
(lambda (row)
(vector-set! row 0
(do ((col-inx 1 (1+ col-inx))
(col-sum 0 (+ col-sum (vector-ref row col-inx))))
((>= col-inx (vector-length row)) (- 1 col-sum))))))
; Generate the Faulhaber's Triangle one row at a time.
; The element at row i >= 0, column j >= 1 (both 0-based) is the product
; of the element at i - 1, j - 1 and the fraction ( i / ( j + 1 ) ).
; The element at column 0 is one minus the sum of all the rest of the columns.
(let ((tri (make-vector row-count)))
(do ((row-inx 0 (1+ row-inx)))
((>= row-inx row-count) tri)
(let ((row (make-vector (1+ row-inx))))
(vector-set! tri row-inx row)
(do ((col-inx 1 (1+ col-inx)))
((>= col-inx (vector-length row)))
(vector-set! row col-inx
(* (vector-ref (vector-ref tri (1- row-inx)) (1- col-inx))
(/ row-inx (1+ col-inx)))))
(calc-store-first! row))))))
 
; Convert elements of a vector to a string for display.
(define vector->string
(lambda (vec)
(do ((inx 0 (1+ inx))
(str "" (string-append str (format "~7@a" (vector-ref vec inx)))))
((>= inx (vector-length vec)) str))))
 
; Display a Faulhaber's Triangle.
(define faulhabers-triangle-display
(lambda (tri)
(do ((inx 0 (1+ inx)))
((>= inx (vector-length tri)))
(printf "~a~%" (vector->string (vector-ref tri inx))))))
 
; Task..
(let ((row-count 10))
(printf "The first ~a rows of Faulhaber's Triangle..~%" row-count)
(faulhabers-triangle-display (faulhabers-triangle row-count)))
(newline)
(let ((power 17)
(sum-to 1000))
(printf "Sum over k=1..~a of k^~a using Faulhaber's Triangle..~%" sum-to power)
(let* ((tri (faulhabers-triangle (1+ power)))
(coefs (vector-ref tri power)))
(printf "~a~%" (do ((inx 0 (1+ inx))
(term-expt sum-to (* term-expt sum-to))
(sum 0 (+ sum (* (vector-ref coefs inx) term-expt))))
((>= inx (vector-length coefs)) sum)))))</syntaxhighlight>
{{out}}
<pre>
The first 10 rows of Faulhaber's Triangle..
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
 
Sum over k=1..1000 of k^17 using Faulhaber's Triangle..
56056972216555580111030077961944183400198333273050000
</pre>
 
=={{header|Sidef}}==
<langsyntaxhighlight lang="ruby">func faulhaber_triangle(p) {
{ binomial(p, _) * bernoulli(_) / p }.map(p ^.. 0)
}
Line 2,286 ⟶ 3,428:
 
say ''
say faulhaber_triangle(p+1).map_kv {|k,v| v * n**(k+1) }.sum</langsyntaxhighlight>
{{out}}
<pre>
Line 2,304 ⟶ 3,446:
 
Alternative solution:
<langsyntaxhighlight lang="ruby">func find_poly_degree(a) {
var c = 0
loop {
Line 2,325 ⟶ 3,467:
}
 
10.times { say faulhaber_triangle(_) }</langsyntaxhighlight>
{{out}}
<pre>
Line 2,338 ⟶ 3,480:
[-1/30, 0, 2/9, 0, -7/15, 0, 2/3, 1/2, 1/9]
[0, -3/20, 0, 1/2, 0, -7/10, 0, 3/4, 1/2, 1/10]
</pre>
 
=={{header|Visual Basic .NET}}==
{{trans|C#}}
<syntaxhighlight lang="vbnet">Module Module1
 
Class Frac
Private ReadOnly num As Long
Private ReadOnly denom As Long
 
Public Shared ReadOnly ZERO = New Frac(0, 1)
Public Shared ReadOnly ONE = 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")
End If
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
 
Private Shared Function Gcd(a As Long, b As Long) As Long
If b = 0 Then
Return a
Else
Return Gcd(b, a Mod b)
End If
End Function
 
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 Overrides Function ToString() As String
If denom = 1 Then
Return num.ToString
Else
Return String.Format("{0}/{1}", num, denom)
End If
End Function
 
Public Overrides 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
End Class
 
Function Bernoulli(n As Integer) As Frac
If n < 0 Then
Throw New ArgumentException("n may not be negative or zero")
End If
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 'first' Bernoulli number
If n <> 1 Then
Return a(0)
Else
Return -a(0)
End If
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
 
Function FaulhaberTriangle(p As Integer) As Frac()
Dim coeffs(p + 1) As Frac
For i = 1 To p + 1
coeffs(i - 1) = Frac.ZERO
Next
Dim q As New Frac(1, p + 1)
Dim sign = -1
For j = 0 To p
sign *= -1
coeffs(p - j) = q * New Frac(sign, 1) * New Frac(Binomial(p + 1, j), 1) * Bernoulli(j)
Next
Return coeffs
End Function
 
Sub Main()
For i = 1 To 10
Dim coeffs = FaulhaberTriangle(i - 1)
For Each coeff In coeffs
Console.Write("{0,5} ", coeff)
Next
Console.WriteLine()
Next
End Sub
 
End Module</syntaxhighlight>
{{out}}
<pre> 1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10</pre>
 
=={{header|V (Vlang)}}==
{{trans|Go}}
<syntaxhighlight lang="v (vlang)">import math.fractions
import math.big
 
fn bernoulli(n int) fractions.Fraction {
mut a := []fractions.Fraction{len: n+1}
for m,_ in a {
a[m] = fractions.fraction(1, i64(m+1))
for j := m; j >= 1; j-- {
mut d := a[j-1]
d = fractions.fraction(i64(j),i64(1)) * (d-a[j])
a[j-1]=d
}
}
// return the 'first' Bernoulli number
if n != 1 {
return a[0]
}
a[0] = a[0].negate()
return a[0]
}
fn binomial(n int, k int) i64 {
if n <= 0 || k <= 0 || n < k {
return 1
}
mut num, mut den := i64(1), i64(1)
for i := k + 1; i <= n; i++ {
num *= i64(i)
}
for i := 2; i <= n-k; i++ {
den *= i64(i)
}
return num / den
}
fn faulhaber_triangle(p int) []fractions.Fraction {
mut coeffs := []fractions.Fraction{len: p+1}
q := fractions.fraction(1, i64(p)+1)
mut t := fractions.fraction(1,1)
mut u := fractions.fraction(1,1)
mut sign := -1
for j,_ in coeffs {
sign *= -1
mut d := coeffs[p-j]
t=fractions.fraction(i64(sign),1)
u = fractions.fraction(binomial(p+1, j),1)
d=q*t
d*=u
d*=bernoulli(j)
coeffs[p-j]=d
}
return coeffs
}
fn main() {
for i in 0..10 {
coeffs := faulhaber_triangle(i)
for coeff in coeffs {
print("${coeff:5} ")
}
println('')
}
println('')
}</syntaxhighlight>
 
{{out}}
<pre>
1/1
1/2 1/2
1/6 1/2 1/3
0/1 1/4 1/2 1/4
-1/30 0/1 1/3 1/2 1/5
0/1 -1/12 0/1 5/12 1/2 1/6
1/42 0/1 -1/6 0/1 1/2 1/2 1/7
0/1 1/12 0/1 -7/24 0/1 7/12 1/2 1/8
-1/30 0/1 2/9 0/1 -7/15 0/1 2/3 1/2 1/9
0/1 -3/20 0/1 1/2 0/1 -7/10 0/1 3/4 1/2 1/10
</pre>
 
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-fmt}}
{{libheader|Wren-math}}
{{libheader|Wren-big}}
<syntaxhighlight lang="wren">import "./fmt" for Fmt
import "./math" for Int
import "./big" for BigRat
 
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] = BigRat.new(1, m+1)
var j = m
while (j >= 1) {
a[j-1] = (a[j-1] - a[j]) * BigRat.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 faulhaberTriangle = Fn.new { |p|
var coeffs = List.filled(p+1, null)
var q = BigRat.new(1, p+1)
var sign = -1
for (j in 0..p) {
sign = sign * -1
var b = BigRat.new(binomial.call(p+1, j), 1)
coeffs[p-j] = q * BigRat.new(sign, 1) * b * bernoulli.call(j)
}
return coeffs
}
 
BigRat.showAsInt = true
for (i in 0..9) {
var coeffs = faulhaberTriangle.call(i)
for (coeff in coeffs) Fmt.write("$5s ", coeff)
System.print()
}
System.print()
// get coeffs for (k + 1)th row
var k = 17
var cc = faulhaberTriangle.call(k)
var n = BigRat.new(1000, 1)
var np = BigRat.one
var sum = BigRat.zero
for (c in cc) {
np = np * n
sum = sum + np*c
}
System.print(sum)</syntaxhighlight>
 
{{out}}
<pre>
1
1/2 1/2
1/6 1/2 1/3
0 1/4 1/2 1/4
-1/30 0 1/3 1/2 1/5
0 -1/12 0 5/12 1/2 1/6
1/42 0 -1/6 0 1/2 1/2 1/7
0 1/12 0 -7/24 0 7/12 1/2 1/8
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10
 
56056972216555580111030077961944183400198333273050000
</pre>
 
Line 2,343 ⟶ 3,813:
{{libheader|GMP}} GNU Multiple Precision Arithmetic Library
Uses the code from [[Faulhaber's formula#zkl]]
<langsyntaxhighlight lang="zkl">foreach p in (10){
faulhaberFormula(p).apply("%7s".fmt).concat().println();
}
Line 2,351 ⟶ 3,821:
.walk() // -->(0, -3617/60 * 1000^2, 0, 595/3 * 1000^4 ...)
.reduce('+) // rat + rat + ...
.println();</langsyntaxhighlight>
{{out}}
<pre>
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/276292...353300"