#!/usr/bin/env python3

# Plot the density of primes

import numpy as np
import matplotlib.pyplot as plt
from sympy import primepi

# Define range for x
#x_values = np.linspace(1000, 10000000, 500)
x_values = np.logspace(1, 6, 250)  # 10^1 to 10^12, 250 data points
pi_x_values = np.array([primepi(int(x)) for x in x_values])  # Number of primes up to x
density_values = pi_x_values / x_values  # Density = π(x)/x
approx_density = 1 / np.log(x_values)  # Approximation by 1/ln(x)

# Plot the actual density and the approximation
plt.figure(figsize=(10, 6))
plt.semilogx(x_values, density_values, label=r'$\pi(n)/n$ (Actual Prime Density)', linewidth=2)
plt.semilogx(x_values, approx_density, linestyle='dashed', label=r'$1/\ln(n)$ (Approximation)', linewidth=2)

# Labels and legend
plt.xlabel(r'$n$', fontsize=12)
plt.ylabel('Density of Primes', fontsize=12)
plt.title('Density of Prime Numbers as n Increases', fontsize=14)
plt.legend()
plt.grid(True)

# Show the plot
plt.show()