Faulhaber's triangle: Difference between revisions

Content added Content deleted

Inline

Revision as of 15:07, 7 June 2017

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

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" . fmap 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);

}

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

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

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)

}

1. First 10 rows of Faulhaber's triangle:

{ |p|

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

} << 1..10

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

@@ Line 35: / Line 35: @@
 * [https://en.wikipedia.org/wiki/Faulhaber%27s_formula Faulhaber's formula (Wikipedia)]
 * [http://www.ww.ingeniousmathstat.org/sites/default/files/Torabi-Dashti-CMJ-2011.pdf Faulhaber's triangle (PDF)]
+<br>
 =={{header|Haskell}}==