# Polyhedral starting polygon for T(p,q): a "straight line" of slope (p,q) drawn on the
# boundary of the square-section, square-plan solid torus (all vertices rational).
# Author: David Niedbala Giraudin, ORCID 0009-0009-1526-1178
# Technical assistance: Claude (Anthropic)
from fractions import Fraction as Fr
from math import gcd

def grid_vertex(a, b, R, r):
    # a, b integers mod 4: plan angle a*90 deg, section angle b*90 deg
    c = [1, 0, -1, 0]; s = [0, 1, 0, -1]
    rho = R + r * c[b % 4]
    return (rho * c[a % 4], rho * s[a % 4], r * s[b % 4])

def p0_polygon(p, q, R, r, k):
    """Exact rational vertices.  The chart line alpha = p t, beta = q t (in units of a
    full turn, t in [0,1)) is cut at the grid lines alpha in Z/4, beta in Z/4; the
    cut points lie on edges of the trapezoidal faces; consecutive cut points lie on
    a common face, so the segment joining them lies in that face.  Each segment is
    then subdivided into k equal parts."""
    assert gcd(p, q) == 1
    ts = sorted(set([Fr(a, 4 * p) for a in range(4 * p)] + [Fr(b, 4 * q) for b in range(4 * q)]))
    pts = []
    for t in ts:
        al = 4 * p * t; be = 4 * q * t       # grid coordinates (in quarter turns)
        a0 = al.numerator // al.denominator; b0 = be.numerator // be.denominator
        fa = al - a0; fb = be - b0
        # point on the bilinear patch boundary: at least one of fa, fb is 0
        V00 = grid_vertex(a0, b0, R, r); V10 = grid_vertex(a0 + 1, b0, R, r)
        V01 = grid_vertex(a0, b0 + 1, R, r); V11 = grid_vertex(a0 + 1, b0 + 1, R, r)
        assert fa == 0 or fb == 0
        P = tuple((1 - fa) * (1 - fb) * V00[i] + fa * (1 - fb) * V10[i] + (1 - fa) * fb * V01[i] + fa * fb * V11[i] for i in range(3))
        pts.append(P)
    out = []
    n = len(pts)
    for i in range(n):
        A = pts[i]; B = pts[(i + 1) % n]
        for j in range(k):
            out.append(tuple(A[d] + (B[d] - A[d]) * Fr(j, k) for d in range(3)))
    return out, ts