constant base = 1_000_000_000

function bcd9_mult(sequence a, sequence b)
sequence c
integer j
atom ci
    c = repeat(0,length(a)+length(b))
    for i=1 to length(a) do
        j = i+length(b)-1
        c[i..j] = sq_add(c[i..j],sq_mul(a[i],b))
    end for

    for i=1 to length(c) do
        ci = c[i]
        if ci>base then
            c[i+1] += floor(ci/base) -- carry
            c[i] = remainder(ci,base)
        end if
    end for

    if c[$]=0 then
        c = c[1..$-1]
    end if
    return c
end function

function atom_to_bcd9(atom a)
sequence s = {}
    while a>0 do
        s = append(s,remainder(a,base))
        a = floor(a/base)
    end while
    return s
end function

function bcd9_to_str(sequence a)
string s = sprintf("%d",a[$])
    for i=length(a)-1 to 1 by -1 do
        s &= sprintf("%09d",a[i])
    end for
    -- (might want to trim leading 0s here)
    return s
end function

sequence a, b, c

a = atom_to_bcd9(power(2,32))
printf(1,"a is %s\n",{bcd9_to_str(a)})

b = bcd9_mult(a,a)
printf(1,"a*a is %s\n",{bcd9_to_str(b)})

c = bcd9_mult(b,b)
printf(1,"a*a*a*a is %s\n",{bcd9_to_str(c)})