# Exact-integer certification tools for C_22a (thickness of a smoothed polygon, # and embeddedness of a piecewise-linear path of polygons). # Author: David Niedbala Giraudin, ORCID 0009-0009-1526-1178 # Technical assistance: Claude (Anthropic) import numpy as np from fractions import Fraction as Fr from math import isqrt def to_int(X, S): return np.rint(np.asarray(X, dtype=np.float64) * S).astype(np.int64) # ---------- candidate pairs by an exact integer grid ---------- def grid_pairs(M2, R2): """M2: integer points. Returns all pairs (i= sqrt(R2)).""" c = isqrt(int(R2)) + 1 cell = M2 // c base = 1 << 20 key = ((cell[:, 0] + base) * (2 * base) + (cell[:, 1] + base)) * (2 * base) + (cell[:, 2] + base) order = np.argsort(key, kind='stable'); ks = key[order] out_i = []; out_j = [] n = len(M2) for dx in (-1, 0, 1): for dy in (-1, 0, 1): for dz in (-1, 0, 1): kk = ((cell[:, 0] + dx + base) * (2 * base) + (cell[:, 1] + dy + base)) * (2 * base) + (cell[:, 2] + dz + base) lo = np.searchsorted(ks, kk, 'left'); hi = np.searchsorted(ks, kk, 'right') cnt = hi - lo; tot = int(cnt.sum()) if tot == 0: continue ii = np.repeat(np.arange(n), cnt) offs = np.arange(tot) - np.repeat(np.cumsum(cnt) - cnt, cnt) jj = order[np.repeat(lo, cnt) + offs] m = ii < jj out_i.append(ii[m]); out_j.append(jj[m]) I = np.concatenate(out_i); J = np.concatenate(out_j) D = M2[I] - M2[J]; d2 = (D * D).sum(1) m = d2 <= R2 return np.stack([I[m], J[m]], 1) def sqrt_lo(n): # n: numpy int64 >= 0, < 2**53 ; float sqrt is correctly rounded return np.sqrt(n.astype(np.float64)) * (1 - 2.0 ** -50) def sqrt_hi(n): return np.sqrt(n.astype(np.float64)) * (1 + 2.0 ** -50) # ---------- embeddedness of the linear step C -> C2 ---------- def step_ok(C, C2): """C, C2: integer polygons (same N). True if every polygon on the segment [C, C2] is embedded, by Lemma 2 (linear steps): (i) non-adjacent edges at distance > 2 delta, delta = max vertex displacement; (ii) at each vertex 2 delta < min(|e0|,|e|) cos(theta/2). All tests are exact integer tests.""" assert int(np.abs(C).max()) < 2 ** 24 and int(np.abs(C2).max()) < 2 ** 24 # squares < 2**53 N = len(C) dv = C2 - C; d2 = (dv * dv).sum(1) Dmax = int(d2.max()) # delta**2 = Dmax exactly dU = isqrt(Dmax) + 1 # integer > delta E = np.roll(C, -1, 0) - C n2 = (E * E).sum(1) E0 = np.roll(E, 1, 0) dot = (E0 * E).sum(1) n0 = (E0 * E0).sum(1) # (ii) exactly: with s = sqrt(n0 n), cos(theta) = dot/s, m = min(n0, n), the condition # 4 delta^2 < m (1 + cos theta)/2 is implied by 8 Dmax (P + 1) < m (P + dot), P = isqrt(n0 n). for a0, a, d in zip(n0.tolist(), n2.tolist(), dot.tolist()): P = isqrt(a0 * a); m = min(a0, a) if not 8 * Dmax * (P + 1) < m * (P + d): return False, "adjacent" lhi = np.array([isqrt(int(x)) + 1 for x in n2.tolist()], dtype=np.int64) # integer > |e| M2 = C + np.roll(C, -1, 0) # doubled midpoints R = 4 * dU + 2 * int(lhi.max()) + 4 # in doubled units; pairs farther apart satisfy (i) pr = grid_pairs(M2, R * R) i, j = pr[:, 0], pr[:, 1] gap = np.abs(i - j); gap = np.minimum(gap, N - gap) msk = gap >= 2; i, j = i[msk], j[msk] D = M2[i] - M2[j]; dm2 = (D * D).sum(1) # = 4 |m_i - m_j|^2, exact # dist(e_i, e_j) >= |m_i - m_j| - (|e_i|+|e_j|)/2 > 2 delta if dm2 > (4 dU + lhi_i + lhi_j)^2 rhs = 4 * dU + lhi[i] + lhi[j] bad = np.nonzero(~(dm2 > rhs * rhs))[0] lim2 = 4 * Dmax # (2 delta)^2, exact Cl = C.tolist() for b in bad: # exact segment distance for the few close pairs a_, c_ = int(i[b]), int(j[b]) dd = seg_dist2(Cl[a_], Cl[(a_ + 1) % N], Cl[c_], Cl[(c_ + 1) % N]) if not dd > lim2: return False, ("pair", a_, c_, float(dd), Dmax) # pairs not returned by grid_pairs have 4|m_i-m_j|^2 > R^2 >= (4 dU + lhi_i + lhi_j)^2 return True, (len(bad), Dmax) # ---------- thickness of the arc-smoothed polygon ---------- def seg_dist2(a, b, c, d): """exact squared distance between segments ab and cd (Fractions / ints).""" u = [b[k] - a[k] for k in range(3)]; v = [d[k] - c[k] for k in range(3)]; w = [a[k] - c[k] for k in range(3)] dot = lambda x, y: x[0]*y[0] + x[1]*y[1] + x[2]*y[2] A = dot(u, u); B = dot(u, v); Cc = dot(v, v); Dd = dot(u, w); Ee = dot(v, w) cands = [] def pt(s, t): return [w[k] + s*u[k] - t*v[k] for k in range(3)] def clamp(x): return min(max(x, Fr(0)), Fr(1)) den = A*Cc - B*B if den != 0: s = clamp(Fr(B*Ee - Cc*Dd, den)); t = clamp(Fr(B*s + Ee, Cc) if Cc else Fr(0)) s = clamp(Fr(B*t - Dd, A) if A else Fr(0)) cands.append((s, t)) for s in (Fr(0), Fr(1)): cands.append((s, clamp(Fr(B*s + Ee, Cc)))) for t in (Fr(0), Fr(1)): cands.append((clamp(Fr(B*t - Dd, A)), t)) best = None for s, t in cands: p = pt(s, t); val = dot(p, p) if best is None or val < best: best = val return best def thickness_cert(F, S, slack=Fr(2, 5), verbose=True): """F: integer polygon (units 1/S). Certifies that the polygon smoothed by unit-radius circular arcs at the vertices has thickness >= 1. Returns (ok, L_up) with L_up an upper bound for its length (real units).""" N = len(F); P = [tuple(int(c) for c in row) for row in F] E = [tuple(P[(i+1) % N][k] - P[i][k] for k in range(3)) for i in range(N)] n = [e[0]*e[0] + e[1]*e[1] + e[2]*e[2] for e in E] llo = [isqrt(x) for x in n]; lhi = [isqrt(x) + 1 for x in n] tan_up = []; ok = True for v in range(N): e1 = E[v-1]; e2 = E[v] dot = e1[0]*e2[0] + e1[1]*e2[1] + e1[2]*e2[2] cr = (e1[1]*e2[2]-e1[2]*e2[1], e1[2]*e2[0]-e1[0]*e2[2], e1[0]*e2[1]-e1[1]*e2[0]) c2 = cr[0]**2 + cr[1]**2 + cr[2]**2 if dot <= 0: return False, "angle>=90 at vertex %d" % v tu = Fr(isqrt(c2) + 1, isqrt(n[v-1]*n[v]) + dot) tan_up.append(tu) if not min(llo[v-1], llo[v]) >= 2 * S * tu: # MinRad >= 1 return False, "minrad<1 at vertex %d" % v tmax = max(tan_up) if not max(lhi) < (Fr(333, 106) - slack) * S: # each edge shorter than pi - slack return False, "edge too long" if not 2 * tmax <= slack: return False, "tangent lengths too large" tau_v = [S * t * t / 2 for t in tan_up] # deviation bound (units) tau_e = [max(tau_v[i], tau_v[(i+1) % N]) for i in range(N)] # cumulative length bounds (units) cs_lo = [0]; cs_hi = [0] for i in range(N): cs_lo.append(cs_lo[-1] + llo[i]); cs_hi.append(cs_hi[-1] + lhi[i]) Lhi = cs_hi[-1]; Llo = cs_lo[-1] pi_lo = Fr(333, 106) thr = (pi_lo - slack) * S # candidate pairs M2 = np.asarray(F, dtype=np.int64) + np.roll(np.asarray(F, dtype=np.int64), -1, 0) taumax = max(tau_e); lmax = max(lhi) R = 2 * (2 * S + 2 * taumax + lmax) + 4 # exact (Fraction); pairs farther apart are fine Rc = -((-R.numerator) // R.denominator) if isinstance(R, Fr) else int(R) # ceiling pr = grid_pairs(M2, Rc * Rc) nexact = 0; nchk = 0; minmargin = None for i, j in pr: i = int(i); j = int(j) if i > j: i, j = j, i # upper bound of the shorter polygon arc between points of e_i and e_j fwd_hi = cs_hi[j+1] - cs_lo[i] back_hi = Lhi - (cs_lo[j] - cs_hi[i+1]) if min(fwd_hi, back_hi) < thr: continue nchk += 1 need = 2 * S + tau_e[i] + tau_e[j] D = M2[i] - M2[j] dm_lo = isqrt(int(D[0])**2 + int(D[1])**2 + int(D[2])**2) # doubled units if Fr(dm_lo, 2) - Fr(lhi[i] + lhi[j], 2) >= need: continue nexact += 1 d2 = seg_dist2(P[i], P[(i+1) % N], P[j], P[(j+1) % N]) if not d2 >= need * need: return False, "pair (%d,%d) too close" % (i, j) mg = d2 - need * need if minmargin is None or mg < minmargin: minmargin = mg if verbose: print(" thickness: %d vertex checks, %d pairs checked, %d exact; max tan %.4g" % (N, nchk, nexact, float(tmax)), flush=True) return True, Fr(Lhi, S)