gauss.py

def issol(x,y,z,p,k):
    return (pow(x,k,p) + pow(y,k,p)) % p == pow(z,k,p)

def count(p,k):
    c = 0
    for x in range(1,p):
        for y in range(1,p):
            for z in range(1,p):
                if issol(x,y,z,p,k):
                    c += 1
    return c

import math
def primes(n):
    result = set(range(2,n))
    for i in range(2,math.floor(math.sqrt(n))):
        for j in range(2,n//i):
            result.discard(i*j)
    return result

k=4
for p in primes(100):
    print('%d : %d' % (p, count(p,k)))

gauss2.py

import math
def primes(n):
    result = set(range(2,n))
    for i in range(2,math.floor(math.sqrt(n))+1):
        for j in range(2,n//i+1):
            result.discard(i*j)
    return result

import numpy as np
k = 6
n = 500
w = np.power(np.arange(1,n),k)
for p in primes(n):
    x = w[:p-1] % p
    print('%d : %d' % (p, (((((x + x.reshape((p-1,1))) - x.reshape(p-1,1,1)) % p==0).sum()))))