function bernoulli(n)
    A = Vector{Rational{BigInt}}(undef, n + 1)
    for m = 0 : n
        A[m + 1] = 1 // (m + 1)
        for j = m : -1 : 1
            A[j] = j * (A[j] - A[j + 1])
        end
    end
    return A[1]
end

function display(n)
    B = map(bernoulli, 0 : n)
    pad = mapreduce(x -> ndigits(numerator(x)) + Int(x < 0), max, B)
    argdigits = ndigits(n)
    for i = 0 : n
        if numerator(B[i + 1]) & 1 == 1
            println(
                "B(", lpad(i, argdigits), ") = ",
                lpad(numerator(B[i + 1]), pad), " / ", denominator(B[i + 1])
            )
        end
    end
end

display(60)
