sequence is1

function comb(atom res, sequence from, integer n, at=1, sequence chosen={})
    if length(chosen)=n then
        sequence digits = repeat(0,10)
        atom sumsq = 0
        for i=1 to length(chosen) do
            integer ci = chosen[i]
            sumsq += ci*ci
            digits[ci+1]+=1
        end for
        if sumsq=0 or is1[sumsq] then
            atom perms = factorial(length(chosen))
            for i=1 to 10 do
                if digits[i] then
                    perms /= factorial(digits[i])
                end if
            end for
            res += perms
        end if
    else
        for i=at to length(from) do
            res = comb(res,from,n,i,append(chosen,from[i]))
        end for
    end if
    return res
end function

procedure setis1(integer n)
    is1 = repeat(0,n*81)
    for i=1 to length(is1) do
        integer r, digit, w = i
        while 1 do
            r = 0
            while w!=0 do
                digit = mod(w,10)
                r += digit*digit
                w = floor(w/10)
            end while
            if r=89 then exit end if
            if r=1 then
                is1[i] = 1
                exit
            end if
            w = r
        end while
    end for
end procedure

constant nums = {0,1,2,3,4,5,6,7,8,9}
for i=1 to 16 do
    atom t0 = time()
    setis1(i)
    printf(1,"There are %d numbers from 1 to 10^%d ending with 89 (%3.2fs)\n",{power(10,i)-comb(0,nums,i),i,time()-t0})
end for