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, (%))

-- 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);

}

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