What is the general $n \times n$ form of the divisibility matrix? Background + Motivation
I had the following idea of using digits as vectors. Let us have say I have a three digit number $a$ in the basis $\lambda$:
$$ a = a_0 + a_1 \lambda + a_2 \lambda^2$$ 
Now, we want to convert $\lambda$ coefficients to $\lambda+1$ coefficients:  
$$ a = a_0 -a_1 +a_2 + (a_1 - 2a_2) (\lambda +1) + a_2(\lambda +1)^2 $$
We note if $a_0 -a_1 +a_2$ is divisible by $\lambda+1$ then so is $a$. And since we are only interested checking if the number is divisible, we introduce the divisibility matrix (for $3$ digits):
$$
(\lambda+1) |a   \implies (\lambda+ 1)|
\begin{pmatrix}
    1       & 0 & 0  \\
\end{pmatrix}
\begin{pmatrix}
    1       & -1 & 1  \\
    0       & 1 & -2 \\
    0       & 0 & 1
\end{pmatrix}
\begin{pmatrix}
    a_0 \\
    a_1 \\
    a_2
\end{pmatrix} 
$$ 
Similarly  if the $a$ is divisible by $\lambda+2$
$$
(\lambda+2) | a \implies (\lambda+ 2) |
\begin{pmatrix}
    1       & 0 & 0  \\
\end{pmatrix}
\begin{pmatrix}
    1       & -1 & 1  \\
    0       & 1 & -2 \\
    0       & 0 & 1
\end{pmatrix}^2
\begin{pmatrix}
    a_0 \\
    a_1 \\
    a_2
\end{pmatrix}
$$ 
Example
Let $\lambda = 10$ and $a=121$
Then we verify $1-2+1 = 0$ which is indeed divisible by $10 + 1 = 11$
Question
We have only done this for a $3$ digit number. What is general form of the divisibility matrix for a $n$ digit number?
 A: We make use of the following: 
$$a = a_0 + a_1 \lambda + a_2 \lambda^2 + \dots $$
Or with $\lambda + 1$ coefficients:
$$ a = b_0 + b_1 (\lambda +1) + b_2 (\lambda +1)^2 + \dots$$
We define $\lambda +1 = \beta $ and combine the above $2$ equations:
$$ a_0 + a_1(\beta -1) + a_2 (\beta -1)^2 + \dots = b_0 + b_1 \beta + b_2\beta^2 +\dots $$
Setting $\beta = 0$:
$$b_0 = a_0 -a_1 +a_2 -a_3 + \dots = \sum_{i=0}^\infty (-1)^i a_i $$
Differentiating and setting $\beta = 0$ again:
$$ b_1 = a_1 - 2 a_2 + 3 a_3 - 4 a_4 + \dots = \sum_{i=1}^\infty   a_i ( -1)^{i+1} i$$
Differentiating and setting $\beta = 0$ again:
$$  b_2 =  \frac{2!}{2! 0!} a_2 - \frac{3!}{1!2!} a_3 + \frac{4!}{2! 2!} a_4 - \frac{5!}{3! 2!} a_4 +  \dots = \sum_{i=2}^\infty a_i(-1)^{i} \text{ }{ }^i C_2  $$
Hence, in general:
$$ b_k = \sum_{i=k}^\infty a_i (-1)^{i-k} \text{ }{ }^i C_k  $$
with $k \neq 0$
Now, we can construct an $n \times n$ divisibility matrix:
$$D_{jk} = \begin{cases} 
      0 & j < k \\
      (-1)^k & j=0 \\
      (-1)^{j-k} \text{ }{ }^j C_k & \text{else}
   \end{cases}
$$
To write it some terms explicitly:
$$ D = 
\begin{pmatrix}
    1       & -1 & 1 & -1 & 1 &\dots  \\
    0       & 1 & -2  & 3 & -4 & \dots \\
    0       & 0 & 1 & -3 & 6 & \dots \\
    \vdots \\
    0       & 0 &  \dots &  & &1   &
\end{pmatrix}$$
Interestingly each column is related to the binomial tree.
A: I've made for me a set of matrix-functions in Pari/GP, that can easily be employed for this.
Assume a function on an argument $x$ which gives a row-vector as a result $$V(x)=[1,x,x^2,x^3,...,x^{n-1}]$$ 
Assume here the size/length $n$ as given as global variable. 
Now assume some rowvector with your coefficients $$A=[a_0,a_1,a_2,...,a_{n-1}]$$
Let's use the Pari/GP-notation "~" for transposes. Then the matrix-product
$$ f(x) = V(x) \cdot A\sim $$
gives a polynomial in $x$ of order $n-1$ .        
Now define the upper-triangular Pascal-/Biomialmatrix of size $n\times n$
$$P = \small \begin{bmatrix} 
 1 & 1 & 1 & 1 & 1 & 1 & \cdots \\ 
 . & 1 & 2 & 3 & 4 & 5 & \cdots \\ 
 . & . & 1 & 3 & 6 & 10 & \cdots \\ 
 . & . & . & 1 & 4 & 10 & \cdots \\ 
 . & . & . & . & 1 & 5 & \cdots \\ 
 . & . & . & . & . & 1& \cdots \\
 \vdots &\vdots &\vdots &\vdots &\vdots &\vdots & \ddots\end{bmatrix}$$
Then, by the binomial-theorem
$$ V(x) \cdot P = V(x+1)$$ 
To have now again $f(x)$ by the matrixproduct with $V(x+1)$ and $A$ you need the inverse of $P$ to formally write
$$ \begin{array} {} f(x) &=& V(x) \cdot A\sim   \\
&=&V(x) \cdot I \cdot A\sim & \text{writing $I$ for the identity matrix} \\
&=& V(x) \cdot (P \cdot P^{-1} ) \cdot A\sim  \\
&=&(V(x) \cdot P) \cdot (P^{-1}  \cdot A\sim ) \\
&=&V(x+1)  \cdot B_1\sim  \\
\end{array}$$
where I write $B_1$ for the vector of your $[b_0,b_1,...,b_{n-1}]$.       
For finite size $n$ this is pretty much generalizable - you can use integer powers $h$ of $P$ to note
$$ f(x) =  (V(x) \cdot P^h) \cdot (P^{-h} \cdot A\sim ) = V(x+h) \cdot B_h\sim $$ 
and even fractional powers of $P$ are definable (via the matrixlogarithm and -exponential).    
For the case of infinite size (thus using formal powerseries instead of polynomials of order $n-1$) there are concerns of convergence (or at least summability in the sense of divergent summation) of the right part of the product $B_h$; but in many cases this can be done without much additional ado.
For instance by this ansatz it is very easy to understand (and implement) Cesaro- and Euler-summation of adjustable order.

Note also, that $P$ can be understand as matrix-exponential of the differentiation-operator, when given as a matrix:
$$ P = \exp(L) \implies L=\log(P)$$ and numerically the matrix $L$ looks like
$$L=\small \begin{bmatrix} 
 . & 1 & . & . & . & . \\ 
 . & . & 2 & . & . & . \\ 
 . & . & . & 3 & . & . \\ 
 . & . & . & . & 4 & . \\ 
 . & . & . & . & . & 5 \\ 
 . & . & . & . & . & .
 \end{bmatrix}
$$ (here of size $6 \times 6$)
With this you can also reproduce your above equations which employ the derivatives by
$$ f'(x) = V(x) \cdot L \cdot A \sim  \\
 = a_1 + 2 a_2 x + 3 a_3 x^2 + 4 a_4 x^3 + 5 a_5 x^4 $$
(for size $n=6$).
