We calculate the fermion propagator and the quark-antiquark Greens functions for a complete set of ultralocal fermion bilinears, ${{cal O}_Gamma}$ [$Gamma$: scalar (S), pseudoscalar (P), vector (V), axial (A) and tensor (T)], using perturbation theory up to one-loop and to lowest order in the lattice spacing. We employ the staggered action for fermions and the Symanzik Improved action for gluons. From our calculations we determine the renormalization functions for the quark field and for all ultralocal taste-singlet bilinear operators. The novel aspect of our calculations is that the gluon links which appear both in the fermion action and in the definition of the bilinears have been improved by applying a stout smearing procedure up to two times, iteratively. Compared to most other improved formulations of staggered fermions, the above action, as well as the HISQ action, lead to smaller taste violating effects. The renormalization functions are presented in the RI$$ scheme; the dependence on all stout parameters, as well as on the coupling constant, the number of colors, the lattice spacing, the gauge fixing parameter and the renormalization scale, is shown explicitly. We apply our results to a nonperturbative study of the magnetic susceptibility of QCD at zero and finite temperature. In particular, we evaluate the tensor coefficient, $tau$, which is relevant to the anomalous magnetic moment of the muon.