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:








∑

k
=
1


n



k

p


=


1

p
+
1




∑

j
=
0


p





(



p
+
1

j


)




B

j



n

p
+
1
−
j




{\displaystyle \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




{\displaystyle 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:






∑

k
=
1


n



k

2


=


1
6


n
+


1
2



n

2


+


1
3



n

3




{\displaystyle \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: 




∑

k
=
1


1000



k

17




{\displaystyle \sum _{k=1}^{1000}k^{17}}

 (extra credit).
 See also

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


import Data.Ratio (Ratio, denominator, numerator, (%))
 
------------------------ FAULHABER -----------------------
 
faulhaber :: Int -> Rational -> Rational
faulhaber p n =
  sum $
    zipWith ((*) . (n ^)) [1 ..] (faulhaberTriangle !! p)
 
 
faulhaberTriangle :: [[Rational]]
faulhaberTriangle =
  tail $
    scanl
      ( \rs n ->
          let xs = zipWith ((*) . (n %)) [2 ..] rs
           in 1 - sum xs : xs
      )
      []
      [0 ..]
 
 
--------------------------- TEST -------------------------
main :: IO ()
main = do
  let triangle = take 10 faulhaberTriangle
      widths = maxWidths triangle
  mapM_
    putStrLn
    [ unlines
        ( (justifyRatio widths 8 ' ' =<<)
            <$> triangle
        ),
      (show . numerator) (faulhaber 17 1000)
    ]
 
------------------------- DISPLAY ------------------------
 
justifyRatio ::
  (Int, Int) -> Int -> Char -> Rational -> String
justifyRatio (wn, wd) n c nd =
  go $
    [numerator, denominator] <*> [nd]
  where
    -- Minimum column width, or more if specified.
    w = max n (wn + wd + 2)
    go [num, den]
      | 1 == den = center w c (show num)
      | otherwise =
        let (q, r) = quotRem (w - 1) 2
         in concat
              [ justifyRight q c (show num),
                "/",
                justifyLeft (q + r) c (show den)
              ]
 
justifyLeft :: Int -> a -> [a] -> [a]
justifyLeft n c s = take n (s <> replicate n c)
 
justifyRight :: Int -> a -> [a] -> [a]
justifyRight n c = (drop . length) <*> (replicate n c <>)
 
center :: Int -> a -> [a] -> [a]
center n c s =
  let (q, r) = quotRem (n - length s) 2
      pad = replicate q c
   in concat [pad, s, pad, replicate r c]
 
maxWidths :: [[Rational]] -> (Int, Int)
maxWidths xss =
  let widest f xs = maximum $ fmap (length . show . f) xs
   in ((,) . widest numerator <*> widest denominator) $
        concat xss