윌런스 vs 신촌우왕: 소수수열 (feat. Claude AI)

[신촌우왕(paste0706)]

참조: 에라토스테네스-신촌우왕 체식

https://gall.dcinside.com/board/view/?id=mathematics&no=421526

수열수학(41): 에라토스테네스-신촌우왕 체식

"""
에라토스테네스-박찬우 소수 일반항 — 수식 충실 버전
Python 3.7 / 32비트 환경 호환 (NumPy, 외부 라이브러리 불필요)

기호 정의:
[b]^a = b // a (몫)
[b]_a = b % a (나머지)
δ(x, y) = 1 (x = y), 0 (x ≠ y)
τ(x, y) = 0 (x = y), 1 (x ≠ y)
ω(n, p) = δ(n, p) + τ([n]_p, 0)
W_x = τ(x, 1) · ∏_{p≤√x} ω(x, p) ← 소수 지시함수
π(k) = Σ_{x=1}^{k} W_x ← 소수 계수 함수

소수 일반항:
P_n = 1 + Σ_{k=1}^{π(k)<n} δ([π(k)]^n, 0)
"""

import math
import time


# ── 기본 기호 함수 ─────────────────────────────────────────

def delta(x, y):
"""δ(x, y) = 1 (x=y), 0 (x≠y)"""
return 1 if x == y else 0

def tau(x, y):
"""τ(x, y) = 0 (x=y), 1 (x≠y)"""
return 0 if x == y else 1

def omega(n, p):
"""ω(n, p) = δ(n, p) + τ([n]_p, 0)"""
return delta(n, p) + tau(n % p, 0)


# ── 소수 지시함수 ──────────────────────────────────────────

def W(x):
"""
소수 지시함수
W_x = τ(x, 1) · ∏_{p≤√x} ω(x, p)

x가 소수이면 1, 비소수이면 0
"""
if x <= 1:
return 0
val = tau(x, 1) # x=1 → 0, x>1 → 1
for p in range(2, int(math.sqrt(x)) + 1):
val *= omega(x, p)
if val == 0:
return 0
return val


# ── 소수 계수 함수 ─────────────────────────────────────────

def pi(k):
"""
소수 계수 함수
π(k) = Σ_{x=1}^{k} W_x

1부터 k까지 소수의 개수
"""
return sum(W(x) for x in range(1, k + 1))


# ── 소수 일반항 ────────────────────────────────────────────

def park_pn(n):
"""
에라토스테네스-박찬우 소수 일반항

P_n = 1 + Σ_{k=1}^{π(k)<n} δ([π(k)]^n, 0)

π(k) = Σ_{x=1}^{k} τ(x,1) · ∏_{p≤√k} ω(x,p)

Args:
n (int): 구하려는 소수의 순번 (1-based)

Returns:
int: n번째 소수
"""
if n < 1:
raise ValueError("n은 1 이상이어야 합니다.")

total = 0
k = 0
pi_k = 0 # π(k) 누적값

while pi_k < n: # 바운더리 컨디션: π(k) < n
k += 1
pi_k += W(k) # π(k) = Σ W_x 갱신
total += delta(pi_k // n, 0) # δ([π(k)]^n, 0)

return 1 + total


# ── 실행 예시 ──────────────────────────────────────────────
if __name__ == "__main__":

test_cases = [
(1, 2), (2, 3), (3, 5), (4, 7), (5, 11),
(6, 13), (7, 17), (8, 19), (9, 23), (10, 29),
(20, 71), (50, 229), (100, 541),
]

print("=" * 54)
print(" 에라토스테네스-박찬우 소수 일반항 (수식 충실 버전)")
print()
print(" P_n = 1 + Σ_{k=1}^{π(k)<n} δ([π(k)]^n, 0)")
print()
print(" δ(x,y)=1(x=y),0(x≠y) τ(x,y)=0(x=y),1(x≠y)")
print(" ω(n,p)=δ(n,p)+τ([n]_p,0)")
print(" W_x = τ(x,1)·∏_{p≤√x} ω(x,p)")
print(" π(k) = Σ_{x=1}^{k} W_x")
print("=" * 54)
print(f" {'n':>4} | {'P_n':>6} | {'시간':>10} | {'검증':>4}")
print("-" * 54)

for n, expected in test_cases:
t0 = time.perf_counter()
result = park_pn(n)
elapsed = time.perf_counter() - t0

def fmt(t):
if t < 0.001: return f"{t*1000:.3f}ms"
elif t < 1: return f"{t*1000:.2f}ms"
else: return f"{t:.3f}s"

mark = "✓" if result == expected else "✗"
print(f" {n:>4} | {result:>6} | {fmt(elapsed):>10} | {mark}")

print("=" * 54)

N = 1000
P_N = park_pn(N)
print(f"P_{N} = {P_N}")