The above formula doesn't seem to work. Let a=b=M=1. Give n=3, c=0, d=1. Then it is not true: 3^2=0-0+1^2.
In the case of a^2+b^2+ab, we can look at this in the form: \frac{a^3-b^3}{a-b}. So for a not equal to b, the case for 3, we have the form: a^3\equiv{b^3}\bmod{p} where the two integers are not the same. This shows that 3 is a divisor of p-1, and therefore p=3k+1. Since k is even, the form is 6k+1.
This then includes the primes 7, 13, 19, 31...As suggested by Ynaugh?
Whether there is a strictly positive integral form for the square of a prime is another matter. As ramsey2879 shows in his original entry: "Thus 3*3 +3*5 +5*5 = 49 and 9 - 3*8 + 8*8 and 25 - 5*8 + 84 also equal 49."
So here we have the number 7, of the form 3k+1, and 7=2^2+1^2+2x1, and so there is an all positve form for the square.
Also, I see that 13=3^2+1^1+3, and 13^2 = 8^2+7^2+56. Now, 19=3^2+2^2+6; 19^2=16^2+5^2+80.
31=5^2+5+1; 31^2=24^2+11^2+264.
There is, however, a key to this matter: We can decompose a^2+b^2+ab = (a-b\omega)(a-b\omega^2), where \omega =\frac{-1+\sqrt-3}{2}
For example: finding the case of 37 = 4^2+3^2+12. We decompose this into:(4-3\omega)(4-3\omega^2) then we square the first factor and simplify using the fact: 1+\omega+\omega^2=0. We then arrive at the first factor for the square is: (7-33\omega)
The square solution then is: 37^2=33^2+7^2+33x7
The use of this method then makes it clear, anytime we can find the above form for the prime, of the form 6k+1, we can build on this to any higher power as well as to any product of such primes. For example: 43^3=126^2+197^2+126x197.
However, if we are to have an all positive form then a and b must be of the different signs in the decomposition. Three is a problem in itself, 3=(1-\omega)(1-\omega^2), and there is no good form for 9, trivally, 3^2=3^2+3^3-3^2. Using the above ideas to get the form 3 times 37^2, we arrive at -26-63\omega, this produces 3x37^2=26^2+73^2-26x73