Faulhaber's triangle: Difference between revisions
m
→{{header|Wren}}: Minor tidy
m (→{{header|Wren}}: Minor tidy) |
|||
| (8 intermediate revisions by 3 users not shown) | |||
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}}==
Line 823 ⟶ 976:
{{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}}==
Line 2,127 ⟶ 2,322:
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}}==
Line 3,523 ⟶ 3,732:
{{libheader|Wren-math}}
{{libheader|Wren-big}}
<syntaxhighlight lang="
import "./math" for Int
import "./big" for BigRat
var bernoulli = Fn.new { |n|
| |||