import numpy as np
import matplotlib.pyplot as plt
from math import cos, sin, pi

# Define the golden ratio
phi = (1 + 5 ** 0.5) / 2

# Define the 12 vertices of the icosahedron
vertices = [
    (0, 1, phi), (0, 1, -phi),
    (0, -1, phi), (0, -1, -phi),
    (1, phi, 0), (1, -phi, 0),
    (-1, phi, 0), (-1, -phi, 0),
    (phi, 0, 1), (phi, 0, -1),
    (-phi, 0, 1), (-phi, 0, -1)
]
vertices = np.array(vertices) * (2 / (1 + 5 ** 0.5))

# Function to draw the icosahedron
def draw_icosahedron():
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    
    # Draw the 20 triangular faces
    face_indices = [
        (0, 1, 4), (0, 4, 3), (0, 3, 2), (0, 2, 1),
        (1, 2, 7), (1, 7, 6), (1, 6, 4),
        (2, 3, 11), (2, 11, 10), (2, 10, 7),
        (3, 4, 5), (3, 5, 11),
        (4, 6, 9), (4, 9, 5),
        (5, 9, 8), (5, 8, 11),
        (6, 7, 10), (6, 10, 9),
        (7, 10, 8), (7, 8, 6),
        (8, 10, 11), (8, 11, 9)
    ]
    
    for i, j, k in face_indices:
        ax.plot([vertices[i,0], vertices[j,0], vertices[k,0], vertices[i,0]],
               [vertices[i,1], vertices[j,1], vertices[k,1], vertices[i,1]],
               [vertices[i,2], vertices[j,2], vertices[k,2], vertices[i,2]], 'b-')
    
    # Set axis limits and labels
    ax.set_xlim3d(-1.2, 1.2)
    ax.set_ylim3d(-1.2, 1.2)
    ax.set_zlim3d(-1.2, 1.2)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Icosahedron')
    plt.show()

# Function to check prime number form and calculate modulo 6 residue
def prime_residue(n):
    if n % 6 in [1, 5]:
        return n % 6
    else:
        return None

# Visualize the prime number modular behavior
primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
residues = [prime_residue(p) for p in primes]

fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter([p for p, r in zip(primes, residues) if r == 1], [1] * sum(r == 1 for r in residues), color='blue', label='6k+1')
ax.scatter([p for p, r in zip(primes, residues) if r == 5], [5] * sum(r == 5 for r in residues), color='red', label='6k-1')
ax.set_xlabel('Prime Number')
ax.set_ylabel('Residue modulo 6')
ax.set_title('Prime Number Residues modulo 6')
ax.legend()
plt.show()

# Visualize the connection to icosahedral symmetry
draw_icosahedron()
print("The icosahedron has 36 rotational symmetries, which correspond to the 36 distinct residue combinations of primes modulo 6.")