Faulhaber's triangle
Named after 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:
where is the nth-Bernoulli number.
The first 5 rows of Faulhaber's triangle, are:
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
Using the third row of the triangle, we have:
- Task
- show the first 10 rows of Faulhaber's triangle.
- using the 18th row of Faulhaber's triangle, compute the sum: (extra credit).
- See also
- Bernoulli numbers
- Evaluate binomial coefficients
- Faulhaber's formula (Wikipedia)
- Faulhaber's triangle (PDF)
Haskell
<lang haskell>import Data.Ratio (Ratio, numerator, denominator, (%))
-- FAULHABER -------------------------------------------------------------------
-- Order of Faulhaber's triangle -> rows of Faulhaber's triangle faulhaberTriangle :: Integer -> Ratio Integer faulhaberTriangle = reverse . fh_
-- Order of Faulhaber's triangle -> reversed list of rows (last row first) fh_ :: Integer -> Ratio Integer fh_ 0 = 1 % 1 fh_ n =
let xs = fh_ (n - 1)
ys = zipWith (\nd j -> nd * (n % (j + 2))) (head xs) [0 ..]
in (((1 % 1) - sum ys) : ys) : xs
-- p -> n -> Sum of the p-th powers of the first n positive integers faulhaber :: Integer -> Ratio Integer -> Ratio Integer faulhaber p n = sum $ zipWith (\nd i -> nd * (n ^ i)) (head $ fh_ p) [1 ..]
-- DISPLAY ---------------------------------------------------------------------
-- (Max numerator+denominator widths) -> Column width -> Filler -> Ratio -> String justifyRatio :: (Int, Int) -> Int -> Char -> Ratio Integer -> String justifyRatio (wn, wd) n c nd =
let justifyLeft n c s = take n (s ++ replicate n c)
justifyRight n c s = drop (length s) (replicate n c ++ s)
center n c s =
let (q, r) = quotRem (n - length s) 2
in concat [replicate q c, s, replicate (q + r) c]
[num, den] = [numerator, denominator] <*> [nd]
w = max n (wn + wd + 2) -- Minimum column width, or more if specified.
in if den == 1
then center w c (show num)
else let (q, r) = quotRem (w - 1) 2
in concat
[ justifyRight q c (show num)
, "/"
, justifyLeft (q + r) c (show den)
]
-- List of Ratios -> (Max numerator width, Max denominator width) maxWidths :: Ratio Integer -> (Int, Int) maxWidths xss =
let widest f xs = maximum $ fmap (length . show . f) xs
[n, d] = [widest numerator, widest denominator] <*> [concat xss]
in (n, d)
-- TEST ------------------------------------------------------------------------
main :: IO ()
main = do
let triangle = faulhaberTriangle 9
widths = maxWidths triangle
mapM_
putStrLn
[ unlines ((justifyRatio widths 8 ' ' =<<) <$> triangle)
, (show . numerator) (faulhaber 17 1000)
]</lang>
- Output:
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
Perl
<lang perl>use 5.010; use List::Util qw(sum); use Math::BigRat try => 'GMP'; use ntheory qw(binomial bernfrac);
sub faulhaber_triangle {
my ($p) = @_;
map {
Math::BigRat->new(bernfrac($_))
* binomial($p, $_)
/ $p
} reverse(0 .. $p-1);
}
- First 10 rows of Faulhaber's triangle
foreach my $p (1 .. 10) {
say map { sprintf("%6s", $_) } faulhaber_triangle($p);
}
- Extra credit
my $p = 17; my $n = Math::BigInt->new(1000); my @r = faulhaber_triangle($p+1); say "\n", sum(map { $r[$_] * $n**($_ + 1) } 0 .. $#r);</lang>
- Output:
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
Perl 6
<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>
- Output:
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
Rexx
<lang rexx>Numeric Digits 100 Do r=0 To 20
ra=r-1
If r=0 Then
f.r.1=1
Else Do
rsum=0
Do c=2 To r+1
ca=c-1
f.r.c=fdivide(fmultiply(f.ra.ca,r),c)
rsum=fsum(rsum,f.r.c)
End
f.r.1=fsubtract(1,rsum)
End
End
Do r=0 To 9
ol= Do c=1 To r+1 ol=ol right(f.r.c,5) End Say ol End
Say x=0 Do c=1 To 18
x=fsum(x,fmultiply(f.17.c,(1000**c))) End
Say k(x) s=0 Do k=1 To 1000
s=s+k**17 End
Say s Exit
fmultiply: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 res=(abs(ad)*abs(bd))'/'||(an*bn) Return s(ad,bd)k(res)
fdivide: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 res=s(ad,bd)(abs(ad)*bn)'/'||(an*abs(bd)) Return k(res)
fsum: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 n=an*bn d=ad*bn+bd*an res=d'/'n Return k(res)
fsubtract: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 n=an*bn d=ad*bn-bd*an res=d'/'n Return k(res)
s: Procedure Parse Arg ad,bd s=sign(ad)*sign(bd) If s<0 Then Return '-'
Else Return
k: Procedure Parse Arg a Parse Var a ad '/' an Select
When ad=0 Then Return 0 When an=1 Then Return ad Otherwise Do g=gcd(ad,an) ad=ad/g an=an/g Return ad'/'an End End
gcd: procedure Parse Arg a,b if b = 0 then return abs(a) return gcd(b,a//b)</lang>
- Output:
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
56056972216555580111030077961944183400198333273050000
Sidef
<lang ruby>func faulhaber_triangle(p) {
{ binomial(p, _) * bernoulli(_) / p }.map(p ^.. 0)
}
-
- First 10 rows of Faulhaber's triangle:
{ |p|
say faulhaber_triangle(p).map{ '%6s' % .as_rat }.join
} << 1..10
-
- Extra credit:
const p = 17 const n = 1000
say say faulhaber_triangle(p+1).map_kv {|k,v| v * n**(k+1) }.sum</lang>
- Output:
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