Generate QAM coefficients?

Hi all, I am new here and I am not sure if this is the right topic to post under. Nonetheless, I was trying to write code that generates the coefficients of a quadrature amplitude modulated signal. The equation of a QAM signal is given below.

From the equation above, a and b are the coefficients which come from a set as shown in the image below. Also, a and b together comprise all possible combinations of the numbers in the set.

M is the constellation size, for instance, 4-QAM, 64-QAM. Therefore, an example of a and b for 4-QAM would be as seen below.

As you can see, the maximum value is +/- sqrt(M) - 1 and all possible combinations are shown. Restating my question, I would like help in writing a function that would generate these values for me. It is quite easy to manually do it for small values of M but obviously for values of M like 64 and 1024, it would be quite difficult. Thank you in advance for any assistance provided.