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, (%)) import Data.List (intercalate)
-- FAULHABER ------------------------------------------------------------------- faulhaberTriangle :: Integer -> Ratio Integer faulhaberTriangle n = reverse $ faulHabers_ n
ft = faulhaberTriangle
faulHabers_ :: Integer -> Ratio Integer faulHabers_ 0 = 1 % 1 faulHabers_ n =
let xs = faulHabers_ (n - 1)
ys = zipWith (\nd j -> nd * (n % (j + 2))) (head xs) [0 ..]
in (((1 % 1) - sum ys) : ys) : xs
faulhaber :: Integer -> Ratio Integer -> Ratio Integer faulhaber k n =
sum $ zipWith (\nd i -> nd * (n ^ i)) (head $ faulHabers_ k) [1 ..]
-- DISPLAY --------------------------------------------------------------------- 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]
in if den == 1
then center n c (show num)
else let (q, r) = quotRem (n - 1) 2
in concat
[ justifyRight q c (show num)
, "/"
, justifyLeft (q + r) c (show den)
]
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)
, justifyRatio widths 0 ' ' (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=> FatRat.new(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
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