function index=decompose_index(k,n_grid,dim)
index=zeros(1,dim);
for i=dim:-1:1
    index(i)=floor((k-1)/n_grid^(i-1));
    k=k-index(i)*n_grid^(i-1);
end
index=index+1;
end