Lagrange's Four Square theorem under modulus

Question

Given an integer n, print out all possible combinations of integer values of A,B,C and D which solve the equation A^2+B^2+C^2+D^2=N under modulus user given integer p

int sum(int n)
{
int a, b, c, d, na, nb, nc, nd;
    int count = 0;


    for (a = 0, na = n; a * a <= na; a++) 
    {

        for (b = a, nb = na - a * a; b * b <= nb; b++) 
        {


           for (c = b, nc = nb - b * b; c * c <= nc; c++) {
            nd = nc - c * c;
            d = sqrt (nd);
            if (d * d == nd)  
                {
                    cout<< a<< b<< c<< d;
                    count++;


                }
            }
        }
    }

   cout<<"Found  solutions :"<< count;
}

I want the values of a,b,c,d to be between 0 to (p-1).

I want to output the respective 4 numbers when squared and summed under modulo P make that respective product modulo P

For eg - 1 x 2 x 3 x 4 x 5 mod 100 = 120 mod 100 = ( 8^2 + 6^2 + 4^2 + 2^2 ) mod 100

Note :value of p can be either prime or non-prime.

How do I reflect it in the above code ?