Instantly share code, notes, and snippets.

@unaoya /satotate.py
Last active Dec 25, 2017

Embed
What would you like to do?
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')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment