 Solution

Use schoolbook multiplication on the digit coefficients:

function Mul(a, b: number): number
    new(result, len(a)+len(b))
    for k := 0 to len(result)-1 do
        result[k] := 0
    for i := 0 to len(a)-1 do
        for j := 0 to len(b)-1 do
            result[i+j] := result[i+j] + a[i]*b[j]
    Normal(result)
    return result

This is discrete convolution of the digit arrays. Expanding the product showsbefore normalization. Normalization preserves the value and establishes the digit bounds. Allocating fresh storage makes input aliasing harmless. The convolution uses $O(LM)$ coefficient products for input lengths $L,M$. Intermediate products and sums need sufficiently wide or arbitrary-precision arithmetic; fitting each normalized digit into a small integer does not alone prevent accumulator overflow.