Faulhaber's triangle

From Rosetta Code
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:



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


FreeBASIC[edit]

Library: GMP
' version 12-08-2017
' compile with: fbc -s console
' uses GMP
 
#Include Once "gmp.bi"
 
#Define i_max 17
 
Dim As UInteger i, j, x
Dim As String s
Dim As ZString Ptr gmp_str : gmp_str = Allocate(100)
 
Dim As Mpq_ptr n, tmp1, tmp2, sum, one, zero
n = Allocate(Len(__mpq_struct)) : Mpq_init(n)
tmp1 = Allocate(Len(__mpq_struct)) : Mpq_init(tmp1)
tmp2 = Allocate(Len(__mpq_struct)) : Mpq_init(tmp2)
sum = Allocate(Len(__mpq_struct)) : Mpq_init(sum)
zero = Allocate(Len(__mpq_struct)) : Mpq_init(zero)
one = Allocate(Len(__mpq_struct)) : Mpq_init(one)
Mpq_set_ui(zero, 0, 0) ' 0/0 = 0
Mpq_set_ui(one , 1, 1) ' 1/1 = 1
 
Dim As Mpq_ptr Faulhaber_triangle(0 To i_max, 1 To i_max +1)
' only initialize the variables we need
For i = 0 To i_max
For j = 1 To i +1
Faulhaber_triangle(i, j) = Allocate(Len(__Mpq_struct))
Mpq_init(Faulhaber_triangle(i, j))
Next
Next
 
Mpq_set(Faulhaber_triangle(0, 1), one)
 
' we calculate the first 18 rows
For i = 1 To i_max
Mpq_set(sum, zero)
For j = i +1 To 2 Step -1
Mpq_set_ui(tmp1, i, j) ' i / j
Mpq_set(tmp2, Faulhaber_triangle(i -1, j -1))
Mpq_mul(Faulhaber_triangle(i, j), tmp2, tmp1)
Mpq_canonicalize(Faulhaber_triangle(i, j))
Mpq_add(sum, sum, Faulhaber_triangle(i, j))
Next
Mpq_sub(Faulhaber_triangle(i, 1), one, sum)
Next
 
Print "The first 10 rows"
For i = 0 To 9
For j = 1 To i +1
Mpq_get_str(gmp_str, 10, Faulhaber_triangle(i, j))
s = Space(6) + *gmp_str + Space(6)
x = InStr(s,"/")
If x = 0 Then x = 7 ' in case of 0 or 1
Print Mid(s, x -3, 7);
Next
Print
Next
print
 
' using the 17'the row
Mpq_set(sum, zero)
Mpq_set_ui(n, 1000, 1) ' 1000/1 = 1000
Mpq_set(tmp2, n)
For j = 1 To 18
Mpq_mul(tmp1, n, Faulhaber_triangle(17, j))
Mpq_add(sum, sum, tmp1)
Mpq_mul(n, n, tmp2)
Next
 
Mpq_get_str(gmp_str, 10, sum)
Print *gmp_str
 
' free memory
DeAllocate(gmp_str)
Mpq_clear(tmp1) : Mpq_clear(tmp2) : Mpq_clear(n)
Mpq_clear(zero) : Mpq_clear(one)  : Mpq_clear(sum)
 
For i = 0 To i_max
For j = 1 To i +1
Mpq_clear(Faulhaber_triangle(i, j))
Next
Next
 
' empty keyboard buffer
While Inkey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
Output:
The first 10 rows
   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

Haskell[edit]

Works with: GHC version 7.10.3
import Data.Ratio (Ratio, numerator, denominator, (%))
import Control.Arrow ((&&&))
 
-- FAULHABER -------------------------------------------------------------------
-- Infinite list of rows of Faulhaber's triangle
faulhaberTriangle :: [[Rational]]
faulhaberTriangle =
tail $
scanl
(\rs n ->
let xs = zipWith ((*) . (n %)) [2 ..] rs
in 1 - sum xs : xs)
[]
[0 ..]
 
-- 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 =
let widest f xs = maximum $ fmap (length . show . f) xs
in 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)
]
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[edit]

ES6[edit]

Translation of: Haskell

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, and good enough motivation to use a different language)

(() => {
 
// Order of Faulhaber's triangle -> rows of Faulhaber's triangle
// faulHaberTriangle :: Int -> [[Ratio Int]]
const faulhaberTriangle = n =>
map(x => tail(
scanl((a, x) => {
const ys = map((nd, i) =>
ratioMult(nd, Ratio(x, i + 2)), a);
return cons(ratioMinus(Ratio(1, 1), ratioSum(ys)), ys);
}, [], enumFromTo(0, x))
),
enumFromTo(0, 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(raise(n, i + 1), 1)),
last(faulhaberTriangle(p))
));
 
// 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));
};
 
// enumFromTo :: Int -> Int -> [Int]
const enumFromTo = (m, n) =>
Array.from({
length: Math.floor(n - m) + 1
}, (_, i) => m + i);
 
// 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;
 
// last :: [a] -> a
const last = xs => xs.length ? xs.slice(-1)[0] : undefined;
 
// 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];
 
// raise :: Num -> Int -> Num
const raise = (n, e) => Math.pow(n, e);
 
// replicate :: Int -> a -> [a]
const replicate = (n, x) =>
Array.from({
length: n
}, () => x);
 
// scanl :: (b -> a -> b) -> b -> [a] -> [b]
const scanl = (f, startValue, xs) =>
xs.reduce((a, x) => {
const v = f(a.acc, x);
return {
acc: v,
scan: cons(a.scan, v)
};
}, {
acc: startValue,
scan: [startValue]
})
.scan;
 
// show :: a -> String
const show = (...x) =>
JSON.stringify.apply(
null, x.length > 1 ? [x[0], null, x[1]] : x
);
 
// tail :: [a] -> [a]
const tail = xs => xs.length ? xs.slice(1) : undefined;
 
// 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)),
]
);
})();
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

Kotlin[edit]

Uses appropriately modified code from the Faulhaber's Formula task:

// version 1.1.2
 
import java.math.BigDecimal
import java.math.MathContext
 
val mc = MathContext(256)
 
fun gcd(a: Long, b: Long): Long = if (b == 0L) a else gcd(b, a % b)
 
class Frac : Comparable<Frac> {
val num: Long
val denom: Long
 
companion object {
val ZERO = Frac(0, 1)
val ONE = Frac(1, 1)
}
 
constructor(n: Long, d: Long) {
require(d != 0L)
var nn = n
var dd = d
if (nn == 0L) {
dd = 1
}
else if (dd < 0) {
nn = -nn
dd = -dd
}
val g = Math.abs(gcd(nn, dd))
if (g > 1) {
nn /= g
dd /= g
}
num = nn
denom = dd
}
 
constructor(n: Int, d: Int) : this(n.toLong(), d.toLong())
 
operator fun plus(other: Frac) =
Frac(num * other.denom + denom * other.num, other.denom * denom)
 
operator fun unaryMinus() = Frac(-num, denom)
 
operator fun minus(other: Frac) = this + (-other)
 
operator fun times(other: Frac) = Frac(this.num * other.num, this.denom * other.denom)
 
fun abs() = if (num >= 0) this else -this
 
override fun compareTo(other: Frac): Int {
val diff = this.toDouble() - other.toDouble()
return when {
diff < 0.0 -> -1
diff > 0.0 -> +1
else -> 0
}
}
 
override fun equals(other: Any?): Boolean {
if (other == null || other !is Frac) return false
return this.compareTo(other) == 0
}
 
override fun toString() = if (denom == 1L) "$num" else "$num/$denom"
 
fun toDouble() = num.toDouble() / denom
 
fun toBigDecimal() = BigDecimal(num).divide(BigDecimal(denom), mc)
}
 
fun bernoulli(n: Int): Frac {
require(n >= 0)
val a = Array(n + 1) { Frac.ZERO }
for (m in 0..n) {
a[m] = Frac(1, m + 1)
for (j in m downTo 1) a[j - 1] = (a[j - 1] - a[j]) * Frac(j, 1)
}
return if (n != 1) a[0] else -a[0] // returns 'first' Bernoulli number
}
 
fun binomial(n: Int, k: Int): Long {
require(n >= 0 && k >= 0 && n >= k)
if (n == 0 || k == 0) return 1
val num = (k + 1..n).fold(1L) { acc, i -> acc * i }
val den = (2..n - k).fold(1L) { acc, i -> acc * i }
return num / den
}
 
fun faulhaberTriangle(p: Int): Array<Frac> {
val coeffs = Array(p + 1) { Frac.ZERO }
val q = Frac(1, p + 1)
var sign = -1
for (j in 0..p) {
sign *= -1
coeffs[p - j] = q * Frac(sign, 1) * Frac(binomial(p + 1, j), 1) * bernoulli(j)
}
return coeffs
}
 
fun main(args: Array<String>) {
for (i in 0..9){
val coeffs = faulhaberTriangle(i)
for (coeff in coeffs) print("${coeff.toString().padStart(5)} ")
println()
}
println()
// get coeffs for (k + 1)th row
val k = 17
val cc = faulhaberTriangle(k)
val n = 1000
val nn = BigDecimal(n)
var np = BigDecimal.ONE
var sum = BigDecimal.ZERO
for (c in cc) {
np *= nn
sum += np * c.toBigDecimal()
}
println(sum.toBigInteger())
}
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[edit]

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);
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[edit]

Works with: Rakudo version 2017.05
Translation of: Sidef
# 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) }
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[edit]

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)
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[edit]

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

zkl[edit]

Uses the code from Faulhaber's formula#zkl and GMP (Gnu Multi Precision library).

foreach p in (10){ 
faulhaberFormula(p).apply("%7s".fmt).concat().println();
}
 
// each term of faulhaberFormula is BigInt/BigInt
[1..].zipWith(fcn(n,rat){ rat*BN(1000).pow(n) }, faulhaberFormula(17))
.walk() // -->(0, -3617/60 * 1000^2, 0, 595/3 * 1000^4 ...)
.reduce('+) // rat + rat + ...
.println();
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