Faulhaber's triangle: Difference between revisions
m (→{{header|Haskell}}: Formatting lint, and a salute to Spoon!'s helpful Faulhaber edits) |
m (→{{header|Haskell}}: Top level tidy to avoid double `where`) |
||
| Line 47: | Line 47: | ||
where |
where |
||
fh_ :: [Rational] -> Integer -> [Rational] |
fh_ :: [Rational] -> Integer -> [Rational] |
||
fh_ x n = |
fh_ x n = |
||
let ys = zipWith (\nd j -> nd * (n % (j + 2))) x [0 ..] |
|||
where |
|||
in (1 - sum ys) : ys |
|||
-- p -> n -> Sum of the p-th powers of the first n positive integers |
-- p -> n -> Sum of the p-th powers of the first n positive integers |
||
Revision as of 05:25, 9 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:
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, (%))
-- FAULHABER ------------------------------------------------------------------- -- Infinite list of rows of Faulhaber's triangle faulhaberTriangle :: Rational faulhaberTriangle = tail $ scanl fh_ [] [0 ..]
where
fh_ :: [Rational] -> Integer -> [Rational]
fh_ x n =
let ys = zipWith (\nd j -> nd * (n % (j + 2))) x [0 ..]
in (1 - sum ys) : ys
-- p -> n -> Sum of the p-th powers of the first n positive integers faulhaber :: Int -> Rational -> Rational faulhaber p n =
sum $ zipWith (\nd i -> nd * (n ^ i)) (faulhaberTriangle !! p) [1 ..]
-- DISPLAY --------------------------------------------------------------------- -- (Max numerator+denominator widths) -> Column width -> Filler -> Ratio -> String justifyRatio :: (Int, Int) -> Int -> Char -> Rational -> 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 :: Rational -> (Int, Int) maxWidths xss = (n, d)
where widest f xs = maximum $ fmap (length . show . f) xs [n, d] = [widest numerator, widest denominator] <*> [concat xss]
-- 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)
]</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
JavaScript
ES6
JavaScript is probably not the right instrument to choose for this task, which requires both a ratio number type and arbitrary precision integers. JavaScript has neither – its only numeric datatype is the IEEE 754 double-precision floating-point format number, into which integers and all else must fit. (See the built-in JS name Number.MAX_SAFE_INTEGER)
This means that we can print Faulhaber's triangle (hand-coding some rudimentary ratio-arithmetic functions), but our only reward for evaluating faulhaber(17, 1000) is an integer overflow. With JS integers out of the box, we can get about as far as faulhaber(17, 8), or faulhaber(4, 1000).
(Further progress would entail implementing some hand-crafted representation of arbitrary precision integers – perhaps a bit beyond the intended scope of this task) <lang JavaScript>(() => {
// Order of Faulhaber's triangle -> rows of Faulhaber's triangle // faulHaberTriangle :: Int -> Ratio Int const faulhaberTriangle = n => reverse(fh_(n));
// p -> n -> Sum of the p-th powers of the first n positive integers
// faulhaber :: Int -> Ratio Int -> Ratio Int
const faulhaber = (p, n) =>
ratioSum(map((nd, i) =>
ratioMult(nd, Ratio(Math.pow(n, i + 1), 1)), head(fh_(p))));
// Order of Faulhaber's triangle -> reversed list of rows (last row first) // fh_ :: Int -> Ratio Int const fh_ = n => n === 0 ? ([ [Ratio(1, 1)] ]) : (() => { const xs = fh_(n - 1), ys = map((nd, j) => ratioMult(nd, Ratio(n, j + 2)), head(xs)); return cons( cons( ratioMinus(Ratio(1, 1), ratioSum(ys)), ys ), xs ); })();
// RATIOS -----------------------------------------------------------------
// (Max numr + denr widths) -> Column width -> Filler -> Ratio -> String
// justifyRatio :: (Int, Int) -> Int -> Char -> Ratio Integer -> String
const justifyRatio = (ws, n, c, nd) => {
const
w = max(n, ws.nMax + ws.dMax + 2),
[num, den] = [nd.num, nd.den];
return all(Number.isSafeInteger, [num, den]) ? (
den === 1 ? center(w, c, show(num)) : (() => {
const [q, r] = quotRem(w - 1, 2);
return concat([
justifyRight(q, c, show(num)),
'/',
justifyLeft(q + r, c, (show(den)))
]);
})()
) : "JS integer overflow ... ";
};
// Ratio :: Int -> Int -> Ratio
const Ratio = (n, d) => ({
num: n,
den: d
});
// ratioMinus :: Ratio -> Ratio -> Ratio
const ratioMinus = (nd, nd1) => {
const
d = lcm(nd.den, nd1.den);
return simpleRatio({
num: (nd.num * (d / nd.den)) - (nd1.num * (d / nd1.den)),
den: d
});
};
// ratioMult :: Ratio -> Ratio -> Ratio
const ratioMult = (nd, nd1) => simpleRatio({
num: nd.num * nd1.num,
den: nd.den * nd1.den
});
// ratioPlus :: Ratio -> Ratio -> Ratio
const ratioPlus = (nd, nd1) => {
const
d = lcm(nd.den, nd1.den);
return simpleRatio({
num: (nd.num * (d / nd.den)) + (nd1.num * (d / nd1.den)),
den: d
});
};
// ratioSum :: [Ratio] -> Ratio
const ratioSum = xs =>
simpleRatio(foldl((a, x) => ratioPlus(a, x), {
num: 0,
den: 1
}, xs));
// ratioWidths :: Ratio -> {nMax::Int, dMax::Int} const ratioWidths = xss => { return foldl((a, x) => { const [nw, dw] = ap( [compose(length, show)], [x.num, x.den] ), [an, ad] = ap( [curry(flip(lookup))(a)], ['nMax', 'dMax'] ); return { nMax: nw > an ? nw : an, dMax: dw > ad ? dw : ad }; }, { nMax: 0, dMax: 0 }, concat(xss)); };
// simpleRatio :: Ratio -> Ratio
const simpleRatio = nd => {
const g = gcd(nd.num, nd.den);
return {
num: nd.num / g,
den: nd.den / g
};
};
// GENERIC FUNCTIONS ------------------------------------------------------
// all :: (a -> Bool) -> [a] -> Bool const all = (f, xs) => xs.every(f);
// A list of functions applied to a list of arguments
// <*> :: [(a -> b)] -> [a] -> [b]
const ap = (fs, xs) => //
[].concat.apply([], fs.map(f => //
[].concat.apply([], xs.map(x => [f(x)]))));
// Size of space -> filler Char -> Text -> Centered Text
// center :: Int -> Char -> Text -> Text
const center = (n, c, s) => {
const [q, r] = quotRem(n - s.length, 2);
return concat(concat([replicate(q, c), s, replicate(q + r, c)]));
};
// compose :: (b -> c) -> (a -> b) -> (a -> c) const compose = (f, g) => x => f(g(x));
// concat :: a -> [a] | [String] -> String const concat = xs => xs.length > 0 ? (() => { const unit = typeof xs[0] === 'string' ? : []; return unit.concat.apply(unit, xs); })() : [];
// cons :: a -> [a] -> [a] const cons = (x, xs) => [x].concat(xs);
// 2 or more arguments
// curry :: Function -> Function
const curry = (f, ...args) => {
const go = xs => xs.length >= f.length ? (f.apply(null, xs)) :
function () {
return go(xs.concat(Array.from(arguments)));
};
return go([].slice.call(args, 1));
};
// flip :: (a -> b -> c) -> b -> a -> c const flip = f => (a, b) => f.apply(null, [b, a]);
// foldl :: (b -> a -> b) -> b -> [a] -> b const foldl = (f, a, xs) => xs.reduce(f, a);
// gcd :: Integral a => a -> a -> a
const gcd = (x, y) => {
const _gcd = (a, b) => (b === 0 ? a : _gcd(b, a % b)),
abs = Math.abs;
return _gcd(abs(x), abs(y));
};
// head :: [a] -> a const head = xs => xs.length ? xs[0] : undefined;
// intercalate :: String -> [a] -> String const intercalate = (s, xs) => xs.join(s);
// justifyLeft :: Int -> Char -> Text -> Text
const justifyLeft = (n, cFiller, strText) =>
n > strText.length ? (
(strText + cFiller.repeat(n))
.substr(0, n)
) : strText;
// justifyRight :: Int -> Char -> Text -> Text
const justifyRight = (n, cFiller, strText) =>
n > strText.length ? (
(cFiller.repeat(n) + strText)
.slice(-n)
) : strText;
// length :: [a] -> Int const length = xs => xs.length;
// lcm :: Integral a => a -> a -> a
const lcm = (x, y) =>
(x === 0 || y === 0) ? 0 : Math.abs(Math.floor(x / gcd(x, y)) * y);
// lookup :: Eq a => a -> [(a, b)] -> Maybe b
const lookup = (k, pairs) => {
if (Array.isArray(pairs)) {
let m = pairs.find(x => x[0] === k);
return m ? m[1] : undefined;
} else {
return typeof pairs === 'object' ? (
pairs[k]
) : undefined;
}
};
// map :: (a -> b) -> [a] -> [b] const map = (f, xs) => xs.map(f);
// max :: Ord a => a -> a -> a const max = (a, b) => b > a ? b : a;
// min :: Ord a => a -> a -> a const min = (a, b) => b < a ? b : a;
// quotRem :: Integral a => a -> a -> (a, a) const quotRem = (m, n) => [Math.floor(m / n), m % n];
// replicate :: Int -> a -> [a]
const replicate = (n, a) => {
let v = [a],
o = [];
if (n < 1) return o;
while (n > 1) {
if (n & 1) o = o.concat(v);
n >>= 1;
v = v.concat(v);
}
return o.concat(v);
};
// reverse :: [a] -> [a]
const reverse = xs =>
typeof xs === 'string' ? (
xs.split()
.reverse()
.join()
) : xs.slice(0)
.reverse();
// show :: a -> String
const show = (...x) =>
JSON.stringify.apply(
null, x.length > 1 ? [x[0], null, x[1]] : x
);
// unlines :: [String] -> String
const unlines = xs => xs.join('\n');
// TEST -------------------------------------------------------------------
const
triangle = faulhaberTriangle(9),
widths = ratioWidths(triangle);
return unlines(
map(row =>
concat(map(cell =>
justifyRatio(widths, 8, ' ', cell), row)), triangle)
) + '\n\n' + unlines(
[
'faulhaber(17, 1000)',
justifyRatio(widths, 0, ' ', faulhaber(17, 1000)),
'\nfaulhaber(17, 8)',
justifyRatio(widths, 0, ' ', faulhaber(17, 8)),
'\nfaulhaber(4, 1000)',
justifyRatio(widths, 0, ' ', faulhaber(4, 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 faulhaber(17, 1000) JS integer overflow ... faulhaber(17, 8) 2502137235710736 faulhaber(4, 1000) 200500333333300
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=> 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
REXX
<lang rexx>Numeric Digits 100 Do r=0 To 20
ra=r-1
If r=0 Then
f.r.1=1
Else Do
rsum=0
Do c=2 To r+1
ca=c-1
f.r.c=fdivide(fmultiply(f.ra.ca,r),c)
rsum=fsum(rsum,f.r.c)
End
f.r.1=fsubtract(1,rsum)
End
End
Do r=0 To 9
ol= Do c=1 To r+1 ol=ol right(f.r.c,5) End Say ol End
Say x=0 Do c=1 To 18
x=fsum(x,fmultiply(f.17.c,(1000**c))) End
Say k(x) s=0 Do k=1 To 1000
s=s+k**17 End
Say s Exit
fmultiply: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 res=(abs(ad)*abs(bd))'/'||(an*bn) Return s(ad,bd)k(res)
fdivide: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 res=s(ad,bd)(abs(ad)*bn)'/'||(an*abs(bd)) Return k(res)
fsum: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 n=an*bn d=ad*bn+bd*an res=d'/'n Return k(res)
fsubtract: Procedure Parse Arg a,b Parse Var a ad '/' an Parse Var b bd '/' bn If an= Then an=1 If bn= Then bn=1 n=an*bn d=ad*bn-bd*an res=d'/'n Return k(res)
s: Procedure Parse Arg ad,bd s=sign(ad)*sign(bd) If s<0 Then Return '-'
Else Return
k: Procedure Parse Arg a Parse Var a ad '/' an Select
When ad=0 Then Return 0 When an=1 Then Return ad Otherwise Do g=gcd(ad,an) ad=ad/g an=an/g Return ad'/'an End End
gcd: procedure Parse Arg a,b if b = 0 then return abs(a) return gcd(b,a//b)</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
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