CGAC2022 Day 3: \$n\$-dimensional Chocolate Pyramid

Question

Part of Code Golf Advent Calendar 2022 event. See the linked meta post for details.

---

I've got an infinite supply of \$n\$-dimensional chocolate for some positive integer \$n\$. The shape of the chocolate is not important. You may assume that they are just \$n\$-dimensional hypercubes with side length \$1\$.

To celebrate the upcoming Christmas, I want to assemble them into a giant chocolate pyramid with \$x\$ layers. The base of the pyramid is an \$(n-1)\$-dimensional hypercube with side length \$x\$, which contains \$x^{n-1}\$ chocolates. The next layer is an \$(n-1)\$-dimensional hypercube with side length \$x-1\$, which contains \$(x-1)^{n-1}\$ chocolates. And so on. The top layer is a single chocolate.

For example, when \$n=3\$, the pyramid would contain \$1^2 + 2^2 + \dots + x^2 = \frac{1}{6}x(x+1)(2x+1)\$ chocolates.

Interestingly, for any dimension \$n\$, the number of chocolates in the pyramid is always a polynomial in \$x\$.

• When \$n=1\$, this is \$\sum_{k=1}^x k^0 = x\$.
• When \$n=2\$, this is \$\sum_{k=1}^x k^1 = \frac{1}{2}x^2+\frac{1}{2}x\$.
• When \$n=3\$, this is \$\sum_{k=1}^x k^2 = \frac{1}{3}x^3+\frac{1}{2}x^2+\frac{1}{6}x\$.
• When \$n=4\$, this is \$\sum_{k=1}^x k^3 = \frac{1}{4}x^4+\frac{1}{2}x^3+\frac{1}{4}x^2\$.
• When \$n=5\$, this is \$\sum_{k=1}^x k^4 = \frac{1}{5}x^5+\frac{1}{2}x^4+\frac{1}{3}x^3-\frac{1}{30}x\$.

The general formula for these polynomials is called the Faulhaber's Formula.

Your task is to find the coefficients of these polynomials.

Rules

The usual sequence rules apply. So you may:

• Output all the polynomials.
• Take an input \$n\$ and output the \$n\$-th polynomial.
• ake an input \$n\$ and output the first \$n\$ polynomial.

You may use \$0\$-indexing or \$1\$-indexing.

You may output the polynomials in any reasonable format. Here are some example formats:

• a list of coefficients, in descending order, e.g. \$\frac{1}{5}x^5+\frac{1}{2}x^4+\frac{1}{3}x^3-\frac{1}{30}x\$ is represented as [1/5,1/2,1/3,0,-1/30,0];
• a list of coefficients, in ascending order, e.g. \$\frac{1}{5}x^5+\frac{1}{2}x^4+\frac{1}{3}x^3-\frac{1}{30}x\$ is represented as [0,-1/30,0,1/3,1/2,1/5];
• a function that takes an input \$k\$ and gives the coefficient of \$x^k\$;
• a built-in polynomial object.

Since the coefficients are not integers, you may output them as rational numbers, floating-point numbers, or any other reasonable format.

You may also take two integers \$n, k\$, and output the coefficient of \$x^k\$ in \$n\$-th polynomial. You may assume that \$k\le n\$.

This is code-golf, so the shortest code in bytes wins.

Testcases

Here I output lists of coefficients in descending order.

1 -> [1, 0]
2 -> [1/2, 1/2, 0]
3 -> [1/3, 1/2, 1/6, 0]
4 -> [1/4, 1/2, 1/4, 0, 0]
5 -> [1/5, 1/2, 1/3, 0, -1/30, 0]
6 -> [1/6, 1/2, 5/12, 0, -1/12, 0, 0]
7 -> [1/7, 1/2, 1/2, 0, -1/6, 0, 1/42, 0]
8 -> [1/8, 1/2, 7/12, 0, -7/24, 0, 1/12, 0, 0]
9 -> [1/9, 1/2, 2/3, 0, -7/15, 0, 2/9, 0, -1/30, 0]
10 -> [1/10, 1/2, 3/4, 0, -7/10, 0, 1/2, 0, -3/20, 0, 0]
11 -> [1/11, 1/2, 5/6, 0, -1, 0, 1, 0, -1/2, 0, 5/66, 0]
12 -> [1/12, 1/2, 11/12, 0, -11/8, 0, 11/6, 0, -11/8, 0, 5/12, 0, 0]
13 -> [1/13, 1/2, 1, 0, -11/6, 0, 22/7, 0, -33/10, 0, 5/3, 0, -691/2730, 0]
14 -> [1/14, 1/2, 13/12, 0, -143/60, 0, 143/28, 0, -143/20, 0, 65/12, 0, -691/420, 0, 0]
15 -> [1/15, 1/2, 7/6, 0, -91/30, 0, 143/18, 0, -143/10, 0, 91/6, 0, -691/90, 0, 7/6, 0]
16 -> [1/16, 1/2, 5/4, 0, -91/24, 0, 143/12, 0, -429/16, 0, 455/12, 0, -691/24, 0, 35/4, 0, 0]

Answer 1

Python3, 207 bytes:

import math
R=range
C=math.comb
B=lambda n,c={0:1,1:0.5}:c[n]if n in c else 1-sum(B(k)/(n-k+1)*C(n,k)for k in R(n))
def r(N):
 d=[0]*-~N
 for i in R(N):d[N-i]+=(1/N)*(-1)**i*C(N,i)*[B(i),-.5][i==1]
 return d

Try it online!

Answer 2

JavaScript (ES6), 73 bytes

-7 bytes thanks to @Neil

Expects (n,k), with \$n\$ 0-indexed. Returns the coefficient of \$x^k\$ in the \$n\$-th polynomial as a floating-point number.

f=(n,k,p=k)=>p--?f(n,k,p)*(p-n)/~p:(b=j=>j--?b(j)-f(k,j):1)(k)/(n+1-k||1)

Try it online!

Less golfed

c = (n, k) => k-- ? c(n, k) * (k - n) / ~k : 1
b = (n, k) => k-- ? b(n, k)- f(n, k) : 1
f = (n, k) => k == n + 1 ? 0 : c(n, k) * b(k, k) / (n + 1 - k)

Try it online!

How?

As a starting point, I used:

• the first alternate expression given on Wikipedia for Faulhaber's formula:

  $$\color{grey}{\sum_{k=1}^n k^p=}\sum_{k=0}^p\dfrac{(-1)^{p-k}}{k+1}\binom{p}{k}B^-_{p-k}\color{grey}{n^{k+1}}$$

  (note that the convention used here is \$B^-_1=-1/2\$ rather than \$B^+_1=+1/2\$)

• the following recursive formula for Bernoulli numbers (also from Wikipedia):

  $$B^-_m=\delta_{m,0}-\sum_{k=0}^{m-1}\binom{m}{k}\dfrac{B^-_k}{m-k+1}$$

  (where \$\delta\$ is the Kronecker delta, which means that we have \$\delta_{m,0}=1\$ if \$m=0\$ and \$\delta_{m,0}=0\$ otherwise)

We can notice similarities between these formulae. So I've tried to factorize as much code as possible (especially the binomial and the division) and eventually came up with the following recursive functions:

$$b(n,k)=\begin{cases}1&\text{if }k=0\\ b(n,k-1)-f(n,k-1)&\text{if }k>0\end{cases}$$

$$f(n,k)=\begin{cases}0&\text{if }k=n+1\\ \dfrac{\binom{n}{k}\times b(k,k)}{n+1-k}&\text{if }k<n+1 \end{cases}$$

The coefficient of \$x^k\$ in the 0-indexed \$n\$-th polynomial is given by \$f(n,k)\$.

In the final JS code, the recursive computation of the binomial is merged into \$f\$ and the first argument of \$b\$ is actually taken from the scope of the caller.

Answer 3

Charcoal, 52 bytes

Ｎθ≔⟦¹⟧ηＦθ«⊞υ⁻¹↨Ｅυ∕×κ§ηλ⁻⊕ιλ¹⊞η⁰ＵＭη⁺κ§η⊖λ»Ｉ∕Ｅυ×ι§ηκθ0

Try it online! Link is to verbose version of code. Explanation:

Ｎθ

Input n.

≔⟦¹⟧η

Start with the first row of Pascal's Triangle.

Ｆθ«

Loop n times.

⊞υ⁻¹↨Ｅυ∕×κ§ηλ⁻⊕ιλ¹

Calculate the next Bernoulli number.

⊞η⁰ＵＭη⁺κ§η⊖λ

Calculate the next row of Pascal's triangle.

»Ｉ∕Ｅυ×ι§ηκθ

Output the first n coefficients.

0

Output the constant term.

55 bytes for a version with less floating-point inaccuracy using the newer version of Charcoal on ATO:

Ｎθ≔⟦¹⟧η≔Π…·¹θζＦθ«⊞υ⁻ζ↨Ｅ⮌×υη÷κ⁺²λ¹⊞η⁰ＵＭη⁺κ§η⊖λ»Ｉ∕×υη×θζ0

Attempt This Online! Link is to verbose version of code. Saves 8 bytes by being able to vectorised multiply two arrays. Explanation: Premultiplies the Bernoulli numbers by n! and then divides the final result by an extra n!.

Answer 4

Jelly,  24 23  15 bytes

c×>ɗþḶZæ*-ṪŻN-¦

A monadic Link that accepts a positive integer, \$n\$, and yields the coefficients in ascending order.

Try it online!

How?

c×>ɗþḶZæ*-ṪŻN-¦ - Link: positive integer, n
     Ḷ          - lowered range (n) -> [0 .. n-1]
    þ           - (i in [1 .. n]) table (j in [0 .. n-1]) with operator:
   ɗ            -   last three links as a dyad f(i,j):
c               -     (i) choose (j)
  >             -     (i) greater than (j)?
 ×              -     multiply -> iCj if i>j else 0
      Z         - transpose
         -      - -1
       æ*       - matrix power
          Ṫ     - tail
           Ż    - prepend a zero
              ¦ - sparse application...
             -  - ...to indices: -1 (penultimate)
            N   - ...action: negate

Answer 5

Desmos, 86 84 bytes

K=[0...m-1]
S(m)=join(nCr(m,K)/m∑_{k=0}^K∑_{v=0}^k(-1)^vnCr(k,v)(v+1)^K/(k+1),0)

My outputs quickly deviate from the actual solutions due to some weird floating point issues, though that is allowed by OP.

Try It On Desmos!

Try it On Desmos! - Prettified

Copying off formulas from Wikipedia go brrr :P

How?

I first started off from this formula from Wikipedia: $$S_m(n)=\frac1{m+1}\sum_{k=0}^m\binom{m+1}kB_k^+n^{m+1-k}$$ Where (formula also from Wikipedia): $$B_m^+=\sum_{k=0}^m\sum_{v=0}^k(-1)^v\binom kv\frac{(v+1)^m}{k+1}$$ So, $$S_m(n)=\frac{\displaystyle\sum_{k=0}^mn^{m+1-k}\binom{m+1}k\sum_{x=0}^k\sum_{v=0}^x(-1)^v\binom xv\frac{(v+1)^k}{x+1}}{m+1}$$ But this formula is shifted one to the left of what is expected; for example, \$S_3(n)\$ gives \$0.25n^4+0.5n^3+0.25n^2\$, which is actually the expected output for 4. To fix this, we will have to shift it back to the right; namely, we have to calculate \$S_{m-1}(n)\$ instead. Plugging this into our formula gives: $$S_{m-1}(n)=\frac{\displaystyle\sum_{k=0}^{m-1}n^{m-k}\binom mk\sum_{x=0}^k\sum_{v=0}^x(-1)^v\binom xv\frac{(v+1)^k}{x+1}}m$$ Notice that this is just a sum of monomials in terms of \$n\$, and that the sum goes from higher powers of \$n\$ to lower powers of \$n\$. So, the list of descending coefficients would be: $$\frac{\displaystyle\binom mk\sum_{x=0}^k\sum_{v=0}^x(-1)^v\binom xv\frac{(v+1)^k}{x+1}}m$$ for \$k=0,1,\ldots,m-1\$.

But looking back at the formula for \$S_{m-1}(n)\$, the outer sum only goes from monomials in terms of \$n^m\$ to \$n\$; specifically, we are missing the constant term. This means that the constant term is always \$0\$, which have the concatenate to the end of the list.

Answer 6

Python NumPy, 99 bytes

f=lambda n,i=1:q**~-n if i>n else(X:=f(n,i+1))+(q-1)**i/i*X[i-1]
from numpy import*
q=poly1d([1,0])

Attempt This Online!

Python, 103 bytes

def f(n):
 Q=8**n;D,N=1,Q**n
 while n:D*=n;s=(2*N//Q**n-~Q)//2%Q-Q//2;N=N*n+Q*(Q-1)**n*s;yield s/D;n-=1

Attempt This Online!

Pure Python version using this tip by @xnor.

Python NumPy, 104 bytes

lambda n:[Q:=q**~-n,*[Q:=Q+(q-1)**i/i*Q[i-1]for i in r_[n:0:-1]]][-1]
from numpy import*
q=poly1d([1,0])

Attempt This Online!

Now doing the Gauss elimination by hand.

Python NumPy, 115 bytes

lambda n:-q*([q**~-n]*mat([q**n+((q-1)**i-q**i)for i in r_[n:0:-1]])[:,1:].I).A1
from numpy import*
q=poly1d([1,0])

Attempt This Online!

Directly solves the linear system of equations that arises from checking terms between p(X) and p(X-1). Uses polynomial (X-1)^k to generate binomial coefficients as needed.

Returns actual polynomials (poly1d objects).

More explanation:

Consider the n+1 dimensional vector space of polynomials of degree <= n with monomials 1,X,X^2,... as a basis. Changing variables X |-> X+d for any constant d induces a linear map L_d from the vector space onto itself as evaluating a polynomial p(X) at X+d yields a new polynomial in X of the same degree. The coefficients of the matrix representing this map will be products of binomial coefficients and powers of d.

Now let us focus on the specific polynomial c(X) this challenge is asking for. This polynomial has the property c(j)-c(j-1)=j^(n-1) for any natural number j. As both sides of this equation are polynomials they must be identical c(X)-c(X-1) = X^(n-1). This can be written in terms of the substitution maps L_d, viz.:X^(n-1) =(L_0-L_-1) c.

What the code does is solving (or asking a library to solve) this equation.

Answer 7

05AB1E, 39 30 (or 29) bytes

LIÝδcθāε!©yL/yGü-ā*}®/}*I/`a)˜

-9 bytes and increased both speed and floating point accuracy (up to \$n\leq20\$) by porting @Neil's Charcoal answer instead.
Therefore also outputs as decimals.

Try it online or verify all test cases.

If we're allowed to output each inner value wrapped in a list, it's 1 byte less by removing the trailing ˜:
Try it online or verify all test cases.

Original 39 bytes answer:

ÎLD<δcÅl¾ζÐøU˜nO/тF©X2Føδ*O®}·s-}θ`s(s)

Port of @JonathanAllan's Jelly answer, except without the builtins.. :/
Therefore also outputs as decimals in reversed order.

Try it online or verify the first 6 test cases. Very slow. Also rounds the result to 5 decimal points in the footer (feel free to remove the 5.ò to see the actual unrounded output).

These kind of challenges aren't really 05AB1E's strong suit..

Explanation:

LIÝδcθ      # Get the input'th Pascal Triangle row:
L           #  Push a list in the range [1, (implicit) input]
 IÝ         #  Push another list in the range [0,input]
   δ        #  Apply double-vectorized over these two lists, creating a table of:
    c       #   Choose (binomial coefficient builtin)
     θ      #  Then only keep the last row of this matrix
āε!©yL/yGü-ā*}®/}
            # Get the first input amount of Bernoulli values:
ā           #  Push a list in the range [1,length] (without popping the list)
 ε          #  Map over each value:
  !         #   Push the factorial of this value
   ©        #   Store it in variable `®` (without popping)
    yL      #   Push a list in the range [1,value]
      /     #   Divide value! by each of these values
  yG        #   Loop the value-1 amount of times:
    ü       #    Map over each overlapping pair:
     -      #    Subtract them from one another
      ā     #    Push a list in the range [1,length] (basically [1,value])
       *    #    Multiply the values at the same positions in the lists
            #    (removing any access items of unequal length lists)
   }        #   Close the loop (resulting in a single wrapped value)
    ®/      #   Divide this wrapped result by the value! we've saved in variable `®`
 }          #  Close the map
*           # Multiply the values at the same positions in the lists together
 I/         # Divide each value by the input
   `a)      # Transform the trailing value to a 0:
   `        #  Pop the list and dump its values to the stack
    a       #  Transform the top value to 0 with an is_alphabetic builtin
     )      #  Wrap the values on the stack back into a list
      ˜     # (Optional: flatten the list of wrapped items to a list of items)
            # (after which the result is output implicitly)

!©yL/yGü-ā*}®/ calculates Bernoulli numbers, based on this CJam answer of @PeterTaylor.
(A 05AB1E port of his answer would start to give incorrect results from around \$n>20\$.)

Î           # Push 0 and the input-integer
 L          # Pop and push a list in the range [1,input]
  D<        # Duplicate, and decrease each value to range [0,input)
    δ       # Apply double-vectorized over these two lists, creating a table of:
     c      #  Choose (binomial coefficient builtin)
 Ål         # Get the lower triangle of this matrix
    ζ       # Zip/transpose this triangular matrix,
   ¾        # using 0 as filler for the unequal length rows
 ÐøU˜nO/тF©X2Føδ*O®}·s-}
            # Calculate the Inverse of this matrix:
 Ð          #  Triplicate this matrix
  ø         #  Zip/transpose the now square matrix back in the top copy
   U        #  Pop and put this transposed copy in variable `X`
  ˜         #  Flatten another copy
   n        #  Get the square of each value in this list
    O       #  Then sum it together
     /      #  Divide the values in the matrix by this sum
тF          #  Loop 100 times:
  ©         #   Store the current matrix in variable `®`
   X        #   Push the matrix from variable `X`
    2F      #   Inner loop 2 times:
      ø     #    Zip/transpose; swapping rows/columns
       δ    #    Apply double-vectorized on the two matrices:
        *   #     Multiply the values at the same positions in the matrices together
         O  #    Sum each inner row
          ® #    Push the matrix from variable `®`
     }·     #   After the inner loop: double each inner value in the top matrix
       s    #   Swap so the other matrix is at the top
        -   #   Subtract the values at the same positions of the matrices
 }          #  Close the outer loop
  θ         # Leave the last row of this Inverted matrix
   `        # Pop and dump its values to the stack
    s       # Swap so the second last value as at the top
     (      # Negate it
      s     # Swap it back
       )    # Wrap all values (including the initial 0) on the stack into a list
            # (which is output implicitly as result)

ÐøU˜nO/тF©X2Føδ*O®}·s- calculates the inverse of a matrix, which I've done before in this 05AB1E answer of mine.

Answer 8

APL(Dyalog Unicode), 25 24 22 20 bytes ^SBCS

-2 thanks to rak1507

-0,⊢⌷∘⌹(1∘+!⍥-⍳)⍤0∘⍳

Try it on APLgolf!

Outputs coefficients in ascending order.

Answer 9

Wolfram Mathematica, 22 bytes

Try it online!

HarmonicNumber[x,1-#]&

---

How?

Sum[k^(n - 1), {k, 1, x}] // FullSimplify reduces to HarmonicNumber[x, 1 - n].

Here HarmonicNumber[n, r] means the \$n\$-th harmonic number of order \$r\$: \$\ \ \ H_n^{(r)}=\sum_{i=1}^n 1 / i^r\$

Answer 10

PARI/GP, 21 bytes

n->sumformal(x^n--,x)

Attempt This Online!

A built-in.

---

PARI/GP, 36 bytes

n->(1/(1-1/matpascal(n))[,^n+1])[n,]

Attempt This Online!

Output a list of coefficients in ascending order.