We present a probabilistic Las Vegas algorithm for computing the local zeta function of a genus-g hyperelliptic curve defined over F q with explicit real multiplication (RM) by an order Z[η] in a degree-g totally real number field. It is based on the approaches by Schoof and Pila in a more favorable case where we can split the-torsion into g kernels of endomorphisms, as introduced by Gaudry, Kohel, and Smith in genus 2. To deal with these kernels in any genus, we adapt a technique that the author, Gaudry, and Spaenlehauer introduced to model the-torsion by structured polynomial systems. Applying this technique to the kernels, the systems we obtain are much smaller and so is the complexity of solving them. Our main result is that there exists a constant c > 0 such that, for any fixed g, this algorithm has expected time and space complexity O((log q) c) as q grows and the characteristic is large enough. We prove that c ≤ 8 and we also conjecture that the result still holds for c = 6.