magic(n) := sum(i, i, 1, n^2)/n$ /* Defining it this way makes computing the inverse too slow */
magic(n) := n*(n^2+1)/2$ /* Therefore, we have to do it this way */

map(magic, makelist(i, i, 3, 22));
magic(1003);

inv_magic(lower) := block(n: 3, while(magic(n)<=lower) do(n: n+1), n)$
map(inv_magic, makelist(10^i, i, 1, 10));