Thanks for all the responses to my question. Based on those, my own analyses, and lots of numerical work, this is what I’ve come up with. Let z be an ordinary complex number, z = x+iy = r exp(i t), where r and t are defined in the standard ways. I limit t to the branch [-pi, pi). Let n = a Gaussian integer; n = a+b i, where a and b are both integers and not both 0. For z to be an nth root of 1, then z^n = 1 = exp(i 2k pi), where k is an integer. And z^n = exp[n ln(z)], where ln(z) = ln(r)+i t, using only the principal value of the ln() function. So, assembling the parts, - ln(r) = 2k pi b / (a^2 + b^2), and - t = 2k pi a / (a^2 + b^2). Looking at the t equation, with the restriction that t is in [-pi,pi), - -pi <= 2k pi a / (a^2 + b^2) < pi, or - -(a^2 + b^2) <= 2k a < (a^2 + b^2). So, finding the number of roots for a given a+b i is just counting how many k’s will satisfy the inequality. Since everything in the last inequality is an integer, this can be performed exactly numerically. Here’s what I found: - For b = 0, a common factor of a drops out of the inequality and it reduces to the standard roots of 1 for integer n. - For a = 0, t = 0 and there are infinite real roots: z = exp(2k pi / b). - For |a| = |b|, the inequality reduces to: -|a| <= k < |a|, and there are 2|a| roots. - For a != 0, |a| can be factored out of each term, leaving -(a^2+b^2)/|a| <= 2k < (a^2+b^2)/|a|. Each side would contribute about (a^2+b^2)/|a| roots, but the 2k term cuts that in half, so the approximate number of roots is (a^2+b^2)/|a|. This result is exact for the previous three cases (assuming 1/0 is infinity). I've constructed a table of the numbers of roots for values of a & b. Is that of sufficient interest to be added to the OEIS? Kerry