We study both analytically, using the renormalization group (RG) to two loop order, and numerically, using an exact polynomial algorithm, the disorder-induced glass phase of the two-dimensional XY model with quenched random symmetry-breaking fields and without vortices. In the super-rough glassy phase, i.e., below the critical temperature T-c, the disorder and thermally averaged correlation function B(r) of the phase field theta(x), B(r) = <(<[theta(x) - theta(x + r)](2)>)over bar> behaves, for r >> a, as B(r) similar or equal to A(tau)ln(2)(r/a) where r = vertical bar r vertical bar and a is a microscopic length scale. We derive the RG equations up to cubic order in tau = (T-c - T)/T-c and predict the universal amplitude A(tau) = 2 tau(2) - 2 tau(2) + O(tau(4)). The universality of A(tau) results from non-trivial cancellations between nonuniversal constants of RG equations. Using an exact polynomial algorithm on an equivalent dimer version of the model we compute A(tau) numerically and obtain a remarkable agreement with our analytical prediction, up to tau approximate to 0.5.