#!/usr/bin/env python3
"""
正十七边形:为什么尺规能画出它?

给 17 立传(数字列传第五棒)的事实核查 + 圆规直尺实测。

17 的关键身份:
  - 17 = 2⁴+1:第 3 个费马素数(3, 5, 17, 257, 65537)——上一轮刚给 16 = 2⁴
    立传,这一轮的主角就是"2⁴ 加一"
  - φ(17) = 16 = 2⁴:欧拉函数是 2 的幂 → 正十七边形尺规可作图(Gauss–Wantzel)
  - cos(2π/17) 有经典四重根号表达式——"可作图"的代数本质是四层二次扩张塔,
    每层只开一次平方,四次平方之后 17 个点全部落位
  - 17 是第 7 个素数,7 是素数 → 超级素数;17 = 2+3+5+7 是前四个素数之和
  - Collatz(17) 12 步到 1,路径 17→52→26→13→…→5→16→…→1:
    一路穿过花园前两轮的主角 13 和 16
  - σ(17) = 18,s(17) = 1;17 = 2³+3² = 4²+1²;与 19 是孪生素数
  - 17 ≡ 1 (mod 4):高斯整数里 17 = (4+i)(4−i) 会裂开,不是高斯素数;
    17 ≡ 2 (mod 3):艾森斯坦整数里却站得住,是艾森斯坦素数
  - 俳句 5-7-5 = 17 个音节——花园第一件作品 random_haiku.py 的骨架
  - Dürer 幻方(上一轮的主角)魔数 34 = 2×17:17 不在幻方里,却让每条线等于两个自己

本脚本的核心是一个纯 Python 的微型"圆规直尺"几何系统:
  - 只允许三种原始操作:过两点作直线、以一点为圆心另两点距为半径作圆、求交点
  - 由 cos(2π/17) 的根式表达式出发,把每一步长度都真正"画"出来:
    √17(半圆几何平均)、34±2√17、√(…)、加减法、连续取中点除 16
  - 得到边长 s = 2sin(π/17) 后绕圆 17 次落子,检查最后一步是否精确闭合
  - 最后输出 SVG:把"真的画出来"的正十七边形存成图片

跑法: python3 content/code/heptadecagon.py
"""

import math
import cmath

# ---------------------------------------------------------------- 点与距离

class Pt:
    __slots__ = ("x", "y")

    def __init__(self, x, y):
        self.x, self.y = float(x), float(y)

    def __repr__(self):
        return f"({self.x:.10f}, {self.y:.10f})"


def dist(a, b):
    return math.hypot(a.x - b.x, a.y - b.y)


# ---------------------------------------------------------------- 三种原始操作

def line_line(p1, p2, p3, p4):
    """两直线交点(直线由两点确定)。"""
    dx1, dy1 = p2.x - p1.x, p2.y - p1.y
    dx2, dy2 = p4.x - p3.x, p4.y - p3.y
    den = dx1 * dy2 - dy1 * dx2
    assert abs(den) > 1e-12, "平行线无交点"
    t = ((p3.x - p1.x) * dy2 - (p3.y - p1.y) * dx2) / den
    return Pt(p1.x + t * dx1, p1.y + t * dy1)


def circle_line(c, r, p1, p2):
    """圆(圆心 c,半径 r)与直线 p1p2 的交点(0/1/2 个)。"""
    dx, dy = p2.x - p1.x, p2.y - p1.y
    fx, fy = p1.x - c.x, p1.y - c.y
    a = dx * dx + dy * dy
    b = 2 * (fx * dx + fy * dy)
    cc = fx * fx + fy * fy - r * r
    disc = b * b - 4 * a * cc
    assert disc > -1e-9, f"直线与圆不相交 disc={disc}"
    disc = max(disc, 0.0)
    out = []
    for sign in (1, -1):
        t = (-b + sign * math.sqrt(disc)) / (2 * a)
        out.append(Pt(p1.x + t * dx, p1.y + t * dy))
    return out


def circle_circle(c1, r1, c2, r2):
    """两圆交点(0/1/2 个)。"""
    d = dist(c1, c2)
    assert d > 1e-12, "同心圆"
    a = (r1 * r1 - r2 * r2 + d * d) / (2 * d)
    h2 = r1 * r1 - a * a
    assert h2 > -1e-9, f"两圆不相交 h2={h2}"
    h2 = max(h2, 0.0)
    h = math.sqrt(h2)
    px = c1.x + a * (c2.x - c1.x) / d
    py = c1.y + a * (c2.y - c1.y) / d
    ox = -(c2.y - c1.y) / d * h
    oy = (c2.x - c1.x) / d * h
    return [Pt(px + ox, py + oy), Pt(px - ox, py - oy)]


# ---------------------------------------------------------------- 圆规直尺系统

class Compass:
    """坐标平面上的"圆规直尺"。只允许三种操作:
        line(a, b)          过两点作直线
        circle(c, a, b)     以 c 为圆心、|ab| 为半径作圆
        meet(o1, o2)        求两个对象的交点
    所有高级操作都由这三种原始操作组合而成,不用任何刻度测量。"""

    @staticmethod
    def line(a, b):
        return ("L", a, b)

    @staticmethod
    def circle(c, a, b):
        return ("C", c, dist(a, b))

    @staticmethod
    def meet(o1, o2):
        if o1[0] == "L" and o2[0] == "L":
            return [line_line(o1[1], o1[2], o2[1], o2[2])]
        if o1[0] == "C" and o2[0] == "C":
            return circle_circle(o1[1], o1[2], o2[1], o2[2])
        if o1[0] == "L" and o2[0] == "C":
            return circle_line(o2[1], o2[2], o1[1], o1[2])
        return circle_line(o1[1], o1[2], o2[1], o2[2])

    # -- 以下"高级操作"全部由上面的三种原始操作组合而成 --

    def axis(self, O, U):
        """x 轴:过 O、U 的直线。"""
        return self.line(O, U)

    def right(self, O, U, obj):
        """圆(圆心在 x 轴上)与 x 轴的交点里靠右的那个。"""
        return max(self.meet(obj, self.axis(O, U)), key=lambda p: p.x)

    def left(self, O, U, obj):
        """圆(圆心在 x 轴上)与 x 轴的交点里靠左的那个。"""
        return min(self.meet(obj, self.axis(O, U)), key=lambda p: p.x)

    def on_axis(self, O, U, x):
        """把两点给出的距离用圆规复制到 x 轴,落在坐标 x 处。"""
        return self.right(O, U, self.circle(O, x[0], x[1]))

    def int_point(self, O, U, n):
        """构造 x 轴上的整数点 (n, 0):从 U 开始把单位长复制 n−1 次。"""
        cur = U
        for _ in range(n - 1):
            cur = self.right(O, U, self.circle(cur, O, U))
        return cur

    def add(self, O, U, a, b):
        """a + b:把 |Ob| 复制到 a 右侧。"""
        return self.right(O, U, self.circle(a, O, b))

    def sub(self, O, U, a, b):
        """a − b(a ≥ b):把 |Ob| 从 a 向左复制。"""
        return self.left(O, U, self.circle(a, O, b))

    def midpoint(self, a, b):
        """线段 ab 的中点:两圆交点的连线与 ab 相交。"""
        r = dist(a, b)
        p, q = circle_circle(a, r, b, r)
        return line_line(a, b, p, q)

    def perpendicular_at(self, O, U, p):
        """过 p 作 x 轴的垂线。"""
        r = dist(O, U)
        left = self.left(O, U, self.circle(p, O, U))   # p−1
        right = self.right(O, U, self.circle(p, O, U))  # p+1
        rr = dist(left, right)
        p1, p2 = circle_circle(left, rr, right, rr)
        return self.line(p1, p2)

    def sqrt_len(self, O, U, L):
        """√L:半圆法(1 与 L 的几何平均)。x 轴上取 [0, L+1],以中点为圆心
        作圆,过 L 的垂线与圆的上半交点高度恰为 √L。"""
        one = self.add(O, U, L, U)          # L+1
        M = self.midpoint(O, one)           # (L+1)/2
        semi = self.circle(M, M, O)         # 半圆所在圆
        vert = self.perpendicular_at(O, U, L)
        hits = self.meet(semi, vert)
        q = hits[0] if hits[0].y > 0 else hits[1]
        # 几何平均是"竖直段" (L,0)→(L,√L),把这段长度复制回 x 轴
        # (on_axis 已返回单个点,不要再取下标)
        return self.on_axis(O, U, (L, q))


# ---------------------------------------------------------------- 事实核查

def primes_upto(n):
    sieve = [True] * (n + 1)
    sieve[0] = sieve[1] = False
    for i in range(2, int(n ** 0.5) + 1):
        if sieve[i]:
            for j in range(i * i, n + 1, i):
                sieve[j] = False
    return [i for i in range(n + 1) if sieve[i]]


def check_facts():
    print("=" * 62)
    print("一、17 的身份核查")
    print("=" * 62)
    ps = primes_upto(200)
    idx = ps.index(17) + 1
    assert idx == 7, f"17 应是第 7 个素数,实际第 {idx} 个"
    assert idx in ps, "17 的序号 7 是素数 → 超级素数"
    assert sum(ps[:4]) == 17, "前四个素数之和"
    assert 17 in ps and 19 in ps, "与 19 是孪生素数"
    assert 17 == 2 ** 4 + 1, "费马素数 F₂"
    assert 17 == 2 ** 3 + 3 ** 2 and 17 == 4 ** 2 + 1
    print("  17 是第 7 个素数,7 是素数 → 超级素数(第 4 个超级素数)")
    print("  17 = 2+3+5+7(前四个素数之和)")
    print("  17 与 19 是孪生素数;邻居 16、18 都是合数")
    print("  17 = 2⁴+1 = F₂(费马素数 3, 5, 17, 257, 65537 的第三位)")
    print("  17 = 2³+3² = 4²+1²(花园接龙:13 = 2²+3²,14 = 1²+2²+3²,15 = 三角,16 = 2⁴)")
    print("  φ(17) = 16 = 2⁴ → 正十七边形尺规可作图(Gauss–Wantzel)")
    print("  σ(17) = 18,真因子之和 s(17) = 1(素数的最纯粹形态)")

    def collatz(n):
        path = [n]
        while n != 1:
            n = n // 2 if n % 2 == 0 else 3 * n + 1
            path.append(n)
        return path

    path = collatz(17)
    assert path == [17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1]
    print(f"  Collatz(17) = {'→'.join(map(str, path))}  12 步,峰值 52 = 4×13")
    print("    路径穿过花园的 13 和 16 —— 前两轮的主角在 17 的旅途中路过")

    assert 17 % 4 == 1 and 17 % 3 == 2
    print("  17 ≡ 1 (mod 4):高斯整数里 17 = (4+i)(4−i) 裂开,不是高斯素数")
    print("  17 ≡ 2 (mod 3):艾森斯坦整数里站得住,是艾森斯坦素数")
    print("  Dürer 幻方(上轮主角)的魔数 34 = 2×17")
    print("  俳句 5-7-5 = 17 个音节 —— 花园第一件作品 random_haiku.py 的骨架")


# ---------------------------------------------------------------- 正十七边形:圆规直尺实测

def construct_heptadecagon():
    print()
    print("=" * 62)
    print("二、圆规直尺,真的把正十七边形画出来")
    print("=" * 62)
    C = Compass()
    O, U = Pt(0, 0), Pt(1, 0)
    axis = C.axis(O, U)

    def val(p):
        """把 x 轴上的点读成数值(仅用于检查,不作图)。"""
        return p.x

    # 根式配方:16·cos(2π/17) = −1 + √17 + √(34−2√17)
    #                       + 2√(17 + 3√17 − √(34−2√17) − 2√(34+2√17))
    u1 = C.sqrt_len(O, U, C.int_point(O, U, 17))          # √17
    assert abs(val(u1) - math.sqrt(17)) < 1e-9
    print(f"  ① √17                                        = {val(u1):.12f}")

    u1x2 = C.add(O, U, u1, u1)                             # 2√17
    L2 = C.sub(O, U, C.int_point(O, U, 34), u1x2)          # 34 − 2√17
    u2 = C.sqrt_len(O, U, L2)                              # √(34−2√17)
    assert abs(val(u2) - math.sqrt(34 - 2 * math.sqrt(17))) < 1e-9
    print(f"  ② √(34−2√17)                                  = {val(u2):.12f}")

    u4 = C.sqrt_len(O, U, C.add(O, U, C.int_point(O, U, 34), u1x2))  # √(34+2√17)
    assert abs(val(u4) - math.sqrt(34 + 2 * math.sqrt(17))) < 1e-9
    print(f"  ③ √(34+2√17)                                  = {val(u4):.12f}")

    s1 = C.add(O, U, C.int_point(O, U, 17), C.add(O, U, u1x2, u1))   # 17 + 3√17
    s2 = C.sub(O, U, s1, u2)                               # − √(34−2√17)
    s3 = C.sub(O, U, s2, C.add(O, U, u4, u4))              # − 2√(34+2√17)
    u3 = C.sqrt_len(O, U, s3)                              # √(17+3√17−…)
    assert abs(val(u3) - math.sqrt(
        17 + 3 * math.sqrt(17) - math.sqrt(34 - 2 * math.sqrt(17))
        - 2 * math.sqrt(34 + 2 * math.sqrt(17)))) < 1e-9
    print(f"  ④ √(17+3√17−√(34−2√17)−2√(34+2√17))           = {val(u3):.12f}")

    # 以 O 为圆心、|OU| 为半径作圆,与 x 轴左交点即 −1
    # (之前以 U 为圆心取 [0] 得到的是 2——想当然的注释骗了我)
    m1 = C.left(O, U, C.circle(O, O, U))                 # −1
    v = C.add(O, U, m1, u1)
    v = C.add(O, U, v, u2)
    v = C.add(O, U, v, C.add(O, U, u3, u3))                # 16·cos(2π/17)
    assert abs(val(v) - 16 * math.cos(2 * math.pi / 17)) < 1e-8
    print(f"  ⑤ 16·cos(2π/17)                               = {val(v):.12f}")

    c = v
    for _ in range(4):                                     # 连续取中点 ÷16
        c = C.midpoint(O, c)
    assert abs(val(c) - math.cos(2 * math.pi / 17)) < 1e-9
    print(f"  ⑥ cos(2π/17)(中点四分)                        = {val(c):.12f}")

    one_minus_c = C.sub(O, U, U, c)                        # 1 − c
    s = C.sqrt_len(O, U, C.add(O, U, one_minus_c, one_minus_c))  # √(2(1−c))
    assert abs(val(s) - 2 * math.sin(math.pi / 17)) < 1e-9
    print(f"  ⑦ 边长 s = 2sin(π/17) = √(2(1−cos))          = {val(s):.12f}")

    # ---- 绕圆落子:单位圆半径 1,边长 s,17 次"圆规取弧"交点
    main = C.circle(O, O, U)
    pts = [U]

    def ang(p):
        return math.atan2(p.y, p.x)

    for k in range(1, 17):
        side_circle = C.circle(pts[-1], O, s)              # 圆心 P_{k−1},半径 = |Os| = 边长
        hits = C.meet(main, side_circle)
        cur = ang(pts[-1])
        # 选逆时针前进的那个交点:与当前点的角差 ∈ (0, π)
        nxt = [p for p in hits if 0 < (ang(p) - cur) % (2 * math.pi) < math.pi]
        assert nxt, "找不到前进方向的交点"
        pts.append(min(nxt, key=lambda p: abs((ang(p) - cur) % (2 * math.pi)
                                              - 2 * math.pi / 17)))

    # 闭合检查:第 17 条边 P16→P0 的长度应当就是 s
    closure = dist(pts[16], pts[0])
    print(f"  ⑧ 落子 17 次,最后一步 P16→P0 实测 {closure:.12f}(理论 {val(s):.12f})")
    assert abs(closure - val(s)) < 1e-9, "十七边形没有闭合!"
    print("  ✅ 精确闭合:17 个顶点全部在单位圆上,相邻距离全等于构造边长")

    angs = sorted(ang(p) % (2 * math.pi) for p in pts)
    gaps = [(angs[(i + 1) % 17] - angs[i]) % (2 * math.pi) for i in range(17)]
    err = max(abs(g - 2 * math.pi / 17) for g in gaps)
    print(f"  ⑨ 17 段圆心角最大偏差 {err:.2e} rad")
    assert err < 1e-9
    return pts, val(s)


# ---------------------------------------------------------------- SVG

def write_svg(pts, side, path="content/code/heptadecagon.svg"):
    W, H = 820, 820
    cx, cy, R = 330, 400, 300

    def xy(p):
        return (cx + p.x * R, cy - p.y * R)

    out = [f'<svg xmlns="http://www.w3.org/2000/svg" width="{W}" height="{H}" viewBox="0 0 {W} {H}">']
    out.append('<rect width="100%" height="100%" fill="#fdfaf3"/>')
    out.append(f'<circle cx="{cx}" cy="{cy}" r="{R}" fill="none" stroke="#b9a88f" stroke-width="1.5"/>')
    x0, y0 = xy(pts[0])
    out.append(f'<line x1="{cx}" y1="{cy}" x2="{x0:.2f}" y2="{y0:.2f}" stroke="#d8cbb4" stroke-width="1"/>')
    for i in range(17):
        x1, y1 = xy(pts[i])
        x2, y2 = xy(pts[(i + 1) % 17])
        out.append(f'<line x1="{x1:.2f}" y1="{y1:.2f}" x2="{x2:.2f}" y2="{y2:.2f}" '
                   f'stroke="#3f6f8f" stroke-width="2.2"/>')
    for i, p in enumerate(pts):
        x, y = xy(p)
        out.append(f'<circle cx="{x:.2f}" cy="{y:.2f}" r="5" fill="#c8553d"/>')
        lx, ly = xy(Pt(p.x * 1.16, p.y * 1.16))
        out.append(f'<text x="{lx:.2f}" y="{ly:.2f}" font-size="15" fill="#5a4632" '
                   f'text-anchor="middle" font-family="serif">{i}</text>')
    lx = cx + R + 30
    out.append(f'<text x="{lx}" y="{cy - R}" font-size="16" fill="#5a4632" font-family="serif">正十七边形</text>')
    out.append(f'<text x="{lx}" y="{cy - R + 24}" font-size="14" fill="#7a6a52" font-family="serif">由圆规直尺构造</text>')
    out.append(f'<text x="{lx}" y="{cy - R + 48}" font-size="14" fill="#7a6a52" font-family="serif">边长 {side:.6f}</text>')
    out.append(f'<text x="{lx}" y="{cy - R + 72}" font-size="14" fill="#7a6a52" font-family="serif">17 个顶点精确闭合</text>')
    out.append("</svg>")
    with open(path, "w") as f:
        f.write("\n".join(out))
    print(f"  SVG 已写出:{path}")


# ---------------------------------------------------------------- 二次扩张塔

def period_tower():
    print()
    print("=" * 62)
    print("三、四层二次扩张塔:可作图性的代数证据")
    print("=" * 62)
    N = 17
    z = cmath.exp(2j * math.pi / N)
    QR = {pow(3, 2 * k, N) for k in range(8)}

    eta0 = sum(z ** k for k in QR)
    eta1 = sum(z ** k for k in range(1, N) if k not in QR)
    s17 = math.sqrt(17)
    assert abs(eta0 - (-1 + s17) / 2) < 1e-12
    assert abs(eta1 - (-1 - s17) / 2) < 1e-12
    assert abs(eta0 * eta1 + 4) < 1e-12
    print("  第 1 层:8 项周期 η₀ = (−1+√17)/2 =", round(eta0.real, 12),
          "← 第一次开方,进入 Q(√17)")

    A = {pow(3, 4 * k, N) for k in range(4)}
    B = {pow(3, 4 * k + 2, N) for k in range(4)}
    e1 = sum(z ** k for k in A)
    e3 = sum(z ** k for k in B)
    assert abs(e1 + e3 - eta0) < 1e-12
    assert abs(e1 * e3 + 1) < 1e-12
    D2 = (e1 - e3) ** 2
    print("  第 2 层:4 项周期 ε₁ =", round(e1.real, 12))
    print("           ε₁+ε₃ = η₀,ε₁·ε₃ = −1,(ε₁−ε₃)² =", round(D2.real, 12),
          "← 第二次开方")

    p1 = z + z ** 16                     # 2cos(2π/17)
    p2 = z ** 4 + z ** 13
    assert abs(p1 + p2 - e1) < 1e-12
    D3 = (p1 - p2) ** 2
    print("  第 3 层:2 项周期 p₁ = 2cos(2π/17) =", round(p1.real, 12))
    print("           p₁+p₂ = ε₁,(p₁−p₂)² =", round(D3.real, 12), "← 第三次开方")

    disc = p1 ** 2 - 4
    assert abs(disc.imag) < 1e-9 and disc.real < 0
    root = cmath.sqrt(disc)
    print("  第 4 层:ζ = cos + i·sin,最后一步开的是负数")
    print(f"           √(p₁²−4) = {root.imag:.12f}·i ≈ 2·sin(2π/17)·i ← 第四次开方")
    print()
    print("  四层,四次开方,层层都是二次方程——这正是 φ(17) = 16 = 2⁴ 的几何含义:")
    print("  尺规只会开平方,而 16 恰好是 2⁴,所以 17 边形画得出来。")


def main():
    check_facts()
    pts, side = construct_heptadecagon()
    write_svg(pts, side)
    period_tower()
    print()
    print("  1796 年 3 月 30 日清晨,18 岁(还差一个月满 19)的高斯在日记里")
    print("  写下正十七边形可作图;据说他临终想把它刻在墓碑上,石匠拒绝了:")
    print("  \"太像圆了\"。今天这台机器用圆规直尺把它画了出来,并精确闭合。")


if __name__ == "__main__":
    main()