satotate.py

import numpy as np
import math
import matplotlib.pyplot as plt

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 list(result)
    
N = 20000
theta = []
ps = primes(N)
for p in ps:
    count = 1 #無限遠点
    for x in range(p):
        rhs = (pow(x,3,p) + pow(x,2,p) - x) % p
        if rhs == 0:
            count += 1
        elif pow(rhs,(p-1)//2,p) == 1:
            count += 2
    theta.append(math.acos((p+1-count)/(2*math.sqrt(p))))
    
plt.hist(theta, bins=18)
x = np.linspace(0,math.pi,100)
plt.plot(x, (np.sin(x)**2)*len(ps)/9)
plt.savefig('satotate.png')