Bernoulli Power Series / Inverse Bernoulli Power Series

Forward: Function input c is the ordinary coefficient vector, real-valued . Function output d is the Bernoulli coefficient vector, real-valued. Use row vectors.
d = bernoulli_power_series(c)

Inverse: Function output c is the ordinary coefficient vector, real-valued . Function input d is the Bernoulli coefficient vector, real-valued. Use row vectors.
c = inverse_bernoulli_power_series(d)

Thorough theory can be found here. Methodology relies on definition of Bernoulli polynomials via inverted Dirichlet series.
https://qr.ae/pNvLNt

Quick explanation
Define an ordinary power series, where c_k | k=0,1,2... denotes the ordinary coefficient vector
y(x) = c_0 + xc_1 + x^2c_2...
Define a Bernoulli power series , where d_k | k=0,1,2... denotes the Bernoulli coefficient vector
y(x) = B_0(x)d_0 + B_1(x)d_1 + B_2(x)d_2...

Forward transform:
Function input is the finite ordinary coefficient vector c_k | k=0,1,2...K and function output is equal-length Bernoulli coefficient vector d_k | k=0,1,2...K

Inverse transform:
Function output is the finite ordinary coefficient vector c_k | k=0,1,2...K and function input is equal-length Bernoulli coefficient vector d_k | k=0,1,2...K

Example 1: Compute the Bernoulli power series of the 6th Bernoulli polynomial, B_6(x) = 1/42 - (1/2)x^2 + (5/2)x^4 - (3)x^5 + (1)x^6 then restore the ordinary coefficients.

c = [1/42 , 0 , -1/2 , 0 , 5/2 , -3 , 1];
d = bernoulli_power_series(c);
c = inverse_bernoulli_power_series(d);

answer: d = [0 , 0 , 0 , 0 , 0 , 0 , 1]

Example 2: Compute the Bernoulli power series of the 9th Partial sums polynomial, S_9(x) = -1/30x + (2/9)x^3 - (7/15)x^5 + (2/3)x^7 + (1/2)x^8 + (1/9)x^9 then restore the ordinary coefficients.

c = [0 , -1/30 , 0 , 2/9 , 0 , -7/15 , 0 , 2/3 , 1/2 , 1/9];
d = bernoulli_power_series(c);
c = inverse_bernoulli_power_series(d);

answer: d=[1/9 , 1 , 4 , 28/3 , 14 , 14 , 28/3 , 4 , 1 , 1/9]