Faulhaber's triangle

Faulhaber's triangle is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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:

$\sum _{k=1}^{n}k^{p}={1 \over p+1}\sum _{j=0}^{p}{p+1 \choose j}B_{j}n^{p+1-j}$

where $B_{n}$ 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:

$\sum _{k=1}^{n}k^{2}={1 \over 6}n+{1 \over 2}n^{2}+{1 \over 3}n^{3}$

Task

• show the first 10 rows of Faulhaber's triangle.
• using the 18th row of Faulhaber's triangle, compute the sum: $\sum _{k=1}^{1000}k^{17}$ (extra credit).

See also

• Bernoulli numbers
• Evaluate binomial coefficients
• Faulhaber's formula (Wikipedia)
• Faulhaber's triangle (PDF)

Haskell

Works with: GHC version 7.10.3

<lang haskell>import Data.Ratio (Ratio, numerator, denominator, (%)) import Data.List (intercalate)

faulhaberTriangle :: Integer -> Ratio Integer faulhaberTriangle n = reverse $ faulHabers_ n

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 ..]

-- TEST AND DISPLAY ------------------------------------------------------------ showRatio :: Ratio Integer -> String showRatio nd =

 if denominator nd /= 1
   then (intercalate "/" . fmap show . (<*>) [numerator, denominator] . pure) nd
   else show $ numerator nd

main :: IO () main =

 mapM_
   putStrLn
   [ unlines $ intercalate "\t" . (showRatio <$>) <$> faulhaberTriangle 9
   , showRatio $ 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);

}

1. First 10 rows of Faulhaber's triangle

foreach my $p (1 .. 10) {

   say map { sprintf("%6s", $_) } faulhaber_triangle($p);

}

1. 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

Works with: Rakudo version 2017.05

Translation of: 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=> 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 };

1. The task

sub faulhaber_triangle ($p) { map { binomial($p+1, $_) * bernoulli[$_] / ($p+1) }, ($p ... 0) }

1. First 10 rows of Faulhaber's triangle:

say faulhaber_triangle($_)».&asRat.fmt('%5s') for ^10; say ;

1. 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)

}

1. 1. First 10 rows of Faulhaber's triangle:

{ |p|

   say faulhaber_triangle(p).map{ '%6s' % .as_rat }.join

} << 1..10

1. 1. 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