code wiki / (root) / nx_partloft_lib.nx
`nx_partloft_lib.nx` source

↩ module page · 1025 lines · 47789 B
1// nx_partloft_lib.nx -- THE LOFTED PART: Infinigen's NurbsPart, replicated in NishiLang from the first byte up
2// (procgen PG59 / anatomy AN19, 2026-09-19). Operator standing order the same day: the estate uses the best-in-class
3// published method as its BASELINE and improves it piece by piece; the ellipsoid-union skull is retired. This file is
4// the first piece. ONE part = a skeleton polyline with a normalised profile at every station, DECOMPOSED into the genes
5// the Infinigen creature system mutates, RECOMPOSED into a grid of NURBS control handles, EVALUATED as a cubic B-spline
6// surface, and WRITTEN as NXMSH2 through the estate's one mesh layout (nx_nxmesh_lib). Read from the primary source,
7// mirrored under knowledge/fetched/ and pinned on procgen.refs:
8//   infinigen/assets/creatures/parts/generic_nurbs.py   decompose_nurbs_handles / recompose_nurbs_handles /
9//                                                      NurbsPart.sample_params (the gene noise STRUCTURE)
10//   infinigen/assets/creatures/util/geometry/lofting.py  factorize_nurbs_handles, skeleton_to_tangents,
11//                                                      compute_profile_verts (frames by rotate_match_directions)
12//   infinigen/assets/creatures/util/geometry/nurbs.py    geomdl_nurbs: degree (3,3), clamped uniform u knots, cyclic v
13//                                                      wrapping three control points, sample_size from face_size
14//   infinigen/assets/creatures/util/part_util.py        nurbs_to_part: the attachment skeleton is the INTERIOR rows
15//   infinigen/assets/creatures/util/creature.py         apply_attach_transform: (u, v, rad) surface coordinates
16// CONVENTIONS, THEIRS: a part runs along +x (profiles face +x), y is lateral, z is up; skeleton_yoffs is the lateral
17// drift per segment and thetas the pitch in the xz plane (degrees there, VM_ONE radians here).
18// UNITS, OURS, and there are exactly three: POSITIONS in tenths of a millimetre (u10, NM_UNIT_PER_M, the NXMSH2
19// ruler's unit); FRACTIONS in PL_FQ = 65536 (Q16) -- gene proportions, lateral drifts, station parameters, normalised
20// profiles, surface parameters, the attachment coordinates, the gene variance, and the B-spline knots and basis; and
21// ANGLES in VM_ONE radians (Q12), the estate's one trigonometric unit, for the pitch genes and the yaw. A radius is a
22// position with sixteen fractional bits (Q16 u10), because a profile normalised by a radius rounded to 0.1 mm was
23// measured 0.25 percent off at 20 mm and would be 2.5 percent off on a finger: MEASURED on the first run of the gate,
24// where Q12 proportions truncated to whole u10 per segment lost 0.5 percent of a part's length (995 for 1000). Every
25// division that lands on a stored value rounds to nearest (pl_rdiv); truncation is how that half-percent was lost.
26// No float anywhere: the file's float32 is encoded once, at the door, by nx_nxmesh_lib / nx_vecmath.
27// DECLARED DEVIATIONS (each a rung on procgen, never a silent difference):
28//   1. np.random.normal is Irwin-Hall (twelve uniforms of a 31-bit LCG): same mean and variance, tails end at six
29//      sigma. The STRUCTURE of the noise is the source's exactly: one size draw, one length draw, one shared radius
30//      draw plus one per station, one per proportion, zero-meaned pitch noise, per-column times per-cell profile noise
31//      mirrored across the profile.
32//   2. The attachment ray-cast against a BVH is a parametric search: around the requested station the closed profile
33//      is searched for the point nearest the ray (least perpendicular distance, in front of the origin), coarse then
34//      bisected, over a one-segment window of u.
35//   3. The attachment skeleton is the clamped cubic B-spline of the interior rows, where the source runs two
36//      Catmull-Clark levels over the interior polyline (the limit of that subdivision IS the cubic B-spline).
37//   4. The emitted mesh carries two fan caps so the part is a closed solid for the lattice fusion; the source leaves
38//      the tube open and lets remesh close it. Winding is made outward per triangle here; the source runs
39//      normals_make_consistent.
40// LIB, no main. license_tier: ORIGINAL (a re-implementation from a read of BSD-3 source; no source bytes copied; the
41// original is an OUTSIDE oracle and never in the build path). No hw writes (Rule 26).
42import "nx_syscalls.nx"
43import "nx_vecmath.nx"
44import "nx_nxmesh_lib.nx"
45import "nx_itoa_lib.nx"
46
47const PL_ONE: i64 = 4096                // VM_ONE spelled once here, the ANGLE and unit-vector unit; the gate asserts the two agree
48const PL_FQ: i64 = 65536                // the FRACTION unit (Q16): genes, parameters, knots, basis
49const PL_FQ2: i64 = 4294967296          // PL_FQ squared: the divisor that takes a Q16 profile times a Q16 radius back to u10
50const PL_FQ_SHIFT: i64 = 16
51const PL_FQ_HALF_SHIFT: i64 = 8         // sqrt(x << 16) << 8 = sqrt(x) << 16: a length with sixteen fractional bits
52const PL_DEG: i64 = 3                   // degree_u, degree_v = (3, 3): nurbs.py geomdl_nurbs
53const PL_ORD: i64 = 4                   // PL_DEG + 1 nonzero basis functions on any span
54const PL_I64: i64 = 8
55const PL_V3: i64 = 3
56const PL_Q4: i64 = 4
57const PL_PART_BYTES: i64 = 256          // PlPart: 21 words, one page-rounded map
58const PL_MIN_RAD: i64 = 65536           // rads clipped to at least one u10 (Q16): np.clip(rads, 1e-3 m, 1e5 m)
59const PL_MIN_NU: i64 = 4                // a clamped cubic needs four control rows
60const PL_MIN_NV: i64 = 3                // and a periodic cubic three columns
61const PL_GAUSS_N: i64 = 12              // Irwin-Hall count: the variance of a sum of n uniforms is n/12, one at twelve
62const PL_LCG_MUL: i64 = 1103515245      // the ANSI C rand() recurrence, a published constant
63const PL_LCG_ADD: i64 = 12345
64const PL_LCG_MOD: i64 = 2147483648      // 2^31: state times multiplier stays under 2^62
65const PL_LCG_SHIFT: i64 = 16            // the low bits of an LCG are its weakest, discarded
66const PL_LCG_BITS: i64 = 32768          // 2^15 output range after the shift
67// THE GENE NOISE of generic_nurbs.py NurbsPart.sample_params, standard deviations in Q16 of the unit value; each is
68// the source's own figure, spelled in this file's unit (0.1 * 65536 = 6554, 0.15 * 65536 = 9830, 0.07 * 65536 = 4587,
69// 7 degrees = 7 * VM_PI / 180 = 500 in VM_ONE radians).
70const PL_SD_SIZE: i64 = 6554            // sz = N(1, 0.1)
71const PL_SD_LENGTH: i64 = 6554          // length *= sz * N(1, 0.1)
72const PL_SD_RADS: i64 = 6554            // rads *= sz * N(1, 0.1) ...
73const PL_SD_RAD_STATION: i64 = 9830     //              ... * N(1, 0.15) per station
74const PL_SD_PROP: i64 = 9830            // proportions *= N(1, 0.15) per segment
75const PL_SD_THETA: i64 = 500            // thetas += N(0, 7 deg) per segment, zero-meaned (VM_ONE radians)
76const PL_SD_PROF_COL: i64 = 4587        // profile noise N(1, 0.07) per column ...
77const PL_SD_PROF_CELL: i64 = 9830       //                ... times N(1, 0.15) per cell
78const PL_SEARCH_COARSE: i64 = 64        // attachment search: samples around the ring before the bisection
79const PL_SEARCH_REFINE: i64 = 12        // bisection steps after the coarse pass
80const PL_SEARCH_UHALF: i64 = 2          // u window: this many quarter-segments either side of the station
81const PL_SEARCH_UQUART: i64 = 4         // a segment is split in this many steps for the window
82const PL_FD_STEP: i64 = 256             // finite-difference step for tangents and normals, Q16 (1/256 of the domain)
83const PL_MAGIC_0: i64 = 78              // N
84const PL_MAGIC_1: i64 = 88              // X
85const PL_MAGIC_2: i64 = 77              // M
86const PL_MAGIC_3: i64 = 83              // S
87const PL_MAGIC_4: i64 = 72              // H
88const PL_MAGIC_5: i64 = 50              // 2
89const PL_LAYER_START_OFF: i64 = 16      // within a 24-byte layer record: 16 name bytes, then start u32, then count u32
90const PL_LAYER_COUNT_OFF: i64 = 20
91const PL_OFF_NORMALS: i64 = 36          // byte offset of the three per-vertex normals within a triangle record
92const PL_COL_BONE_R: i64 = 900          // the part's default colour, per-mille: a bone tint a viewer renders as bone
93const PL_COL_BONE_G: i64 = 880
94const PL_COL_BONE_B: i64 = 840
95const PL_MIN_RING: i64 = 3              // a cap fan needs three ring points
96const PL_SAVE_DIGITS: i64 = 16          // bytes reserved per integer on save: sign, 14 digits, separator
97const PL_SAVE_HEADER: i64 = 64
98const PL_MODE644: i64 = 420
99const PL_LOAD_HDR: i64 = 16             // the two-word length record sys_read_file fills
100const PL_FMT_VERSION: i64 = 2           // version 2: Q16 fractions and Q16 radii (version 1 was Q12, never shipped)
101const PL_HEADER_INTS: i64 = 7           // version nu nv root(3) length
102const PL_ASKEL_SPLINE: i64 = 1          // pl_askel took the cubic path (four or more interior rows)
103const PL_ASKEL_LINEAR: i64 = 2          // pl_askel took the linear path (fewer)
104const PL_POS_HUGE: i64 = 1000000000000000000
105
106struct PlPart { nu: i64, nv: i64, root: *i64, length: i64, prop: *i64, theta: *i64, yoff: *i64, rad: *i64, ts: *i64, prof: *i64, skel: *i64, tang: *i64, hnd: *i64, ku: *i64, kv: *i64, ka: *i64, bu: *i64, bv: *i64, bl: *i64, br: *i64 }
107
108func pl_abs(v: i64) -> i64 { if v < 0 { return 0 - v } return v }
109// division rounded to nearest, b > 0; the truncating form is how the first run lost half a percent of every length
110func pl_rdiv(a: i64, b: i64) -> i64 { if a >= 0 { return (a + b/2) / b } return (a - b/2) / b }
111func pl_len3(x: i64, y: i64, z: i64) -> i64 { return vm_isqrt(x*x + y*y + z*z) }
112// the same length with sixteen fractional bits (Q16 u10)
113func pl_len3_fx(x: i64, y: i64, z: i64) -> i64 { return vm_isqrt((x*x + y*y + z*z) << PL_FQ_SHIFT) << PL_FQ_HALF_SHIFT }
114// unit vector (Q12) of an integer vector in any unit; a null vector yields zeros and returns 0
115func pl_unit(o: *i64, a: *i64) -> i64 {
  let l: i64 = pl_len3(a[0], a[1], a[2])
  if l == 0 { o[0] = 0; o[1] = 0; o[2] = 0; return 0 }
  o[0] = pl_rdiv(a[0]*PL_ONE, l); o[1] = pl_rdiv(a[1]*PL_ONE, l); o[2] = pl_rdiv(a[2]*PL_ONE, l)
  return 1
120}
121// a station's radius in whole u10, rounded
122func pl_rad_u10(p: *PlPart, i: i64) -> i64 { return pl_rdiv(p.rad[i], PL_FQ) }
123
124// ---- the knot vectors, nurbs.py generate_knotvector ------------------------------------------------------------
125// u: clamped uniform over nu control rows, degree 3: four zeros, nu-4 interior knots k/(nu-3), four ones (Q16).
126// v: uniform unclamped over nv columns, k/(nv+3) for k = 0..nv+3, then the first three intervals wrapped past the end
127//    (kv[1..3] + kv[-1] - kv[0]), so a periodic cubic over nv+3 wrapped control columns has nv+7 knots.
128// a: the attachment skeleton's clamped knots over the nu-2 interior rows (filled only when four or more exist).
129func pl_knots(p: *PlPart) -> i64 {
  let nu: i64 = p.nu
  let nv: i64 = p.nv
  var i: i64 = 0
  while i < nu + PL_ORD {
      var k: i64 = 0
      if i >= PL_ORD { k = (i - PL_DEG) * PL_FQ / (nu - PL_DEG) }
      if i >= nu { k = PL_FQ }
      p.ku[i] = k
      i = i + 1
  }
  i = 0
  while i < nv + PL_ORD {
      p.kv[i] = i * PL_FQ / (nv + PL_DEG)
      i = i + 1
  }
  var w: i64 = 1
  while w <= PL_DEG {
      p.kv[nv + PL_DEG + w] = PL_FQ + w * PL_FQ / (nv + PL_DEG)
      w = w + 1
  }
  let na: i64 = nu - 2
  if na >= PL_ORD {
      i = 0
      while i < na + PL_ORD {
          var ka: i64 = 0
          if i >= PL_ORD { ka = (i - PL_DEG) * PL_FQ / (na - PL_DEG) }
          if i >= na { ka = PL_FQ }
          p.ka[i] = ka
          i = i + 1
      }
  }
  return 0
162}
163
164// a part with nu stations of nv profile points, every array owned by the part; arrays are page-rounded maps
165func pl_new(nu: i64, nv: i64) -> *PlPart {
  let p: *PlPart = (sys_mmap(PL_PART_BYTES)) as *PlPart
  p.nu = nu; p.nv = nv; p.length = 0
  p.root = sys_mmap(PL_V3*PL_I64) as *i64
  p.prop = sys_mmap(nu*PL_I64) as *i64
  p.theta = sys_mmap(nu*PL_I64) as *i64
  p.yoff = sys_mmap(nu*PL_I64) as *i64
  p.rad = sys_mmap(nu*PL_I64) as *i64
  p.ts = sys_mmap(nu*PL_I64) as *i64
  p.prof = sys_mmap(nu*nv*PL_V3*PL_I64) as *i64
  p.skel = sys_mmap(nu*PL_V3*PL_I64) as *i64
  p.tang = sys_mmap(nu*PL_V3*PL_I64) as *i64
  p.hnd = sys_mmap(nu*nv*PL_V3*PL_I64) as *i64
  p.ku = sys_mmap((nu+PL_ORD)*PL_I64) as *i64
  p.kv = sys_mmap((nv+PL_ORD+PL_DEG)*PL_I64) as *i64
  p.ka = sys_mmap((nu+PL_ORD)*PL_I64) as *i64
  p.bu = sys_mmap(PL_Q4*PL_I64) as *i64
  p.bv = sys_mmap(PL_Q4*PL_I64) as *i64
  p.bl = sys_mmap(PL_Q4*PL_I64) as *i64
  p.br = sys_mmap(PL_Q4*PL_I64) as *i64
  if nu >= PL_MIN_NU { if nv >= PL_MIN_NV { pl_knots(p) } }
  return p
187}
188
189// every gene and every derived array of src into dst (same nu, nv); returns 0, or -1 on a shape mismatch
190func pl_copy(dst: *PlPart, src: *PlPart) -> i64 {
  if dst.nu != src.nu { return 0 - 1 }
  if dst.nv != src.nv { return 0 - 1 }
  let nu: i64 = src.nu
  let nv: i64 = src.nv
  var i: i64 = 0
  while i < PL_V3 { dst.root[i] = src.root[i]; i = i + 1 }
  dst.length = src.length
  i = 0
  while i < nu { dst.prop[i] = src.prop[i]; dst.theta[i] = src.theta[i]; dst.yoff[i] = src.yoff[i]; dst.rad[i] = src.rad[i]; dst.ts[i] = src.ts[i]; i = i + 1 }
  i = 0
  while i < nu*PL_V3 { dst.skel[i] = src.skel[i]; dst.tang[i] = src.tang[i]; i = i + 1 }
  i = 0
  while i < nu*nv*PL_V3 { dst.prof[i] = src.prof[i]; dst.hnd[i] = src.hnd[i]; i = i + 1 }
  pl_knots(dst)
  return 0
206}
207
208// ---- rotation that carries unit vector a onto unit vector b (rotate_match_directions, one pair) ----------------
209// q = (a x b, 1 + a.b) normalised: the half-angle form needs no arccos and keeps a small angle to 1/4096 rad.
210// A null cross product (parallel OR antiparallel) yields the identity, as the source's |axis| <= 1e-4 branch does.
211func pl_rot_match(q: *i64, a: *i64, b: *i64) -> i64 {
  let cx: i64 = (a[1]*b[2] - a[2]*b[1]) / PL_ONE
  let cy: i64 = (a[2]*b[0] - a[0]*b[2]) / PL_ONE
  let cz: i64 = (a[0]*b[1] - a[1]*b[0]) / PL_ONE
  if cx*cx + cy*cy + cz*cz == 0 { vm_q_ident(q); return 0 }
  let d: i64 = (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) / PL_ONE
  q[0] = cx; q[1] = cy; q[2] = cz; q[3] = PL_ONE + d
  vm_q_norm(q, q)
  return 1
220}
221
222// skeleton_to_tangents: edge directions, averaged with the previous edge at interior stations, unit (Q12).
223// edge is caller scratch of n*3 words.
224func pl_tangents(skel: *i64, n: i64, tang: *i64, edge: *i64) -> i64 {
  var i: i64 = 0
  while i < n - 1 {
      edge[i*PL_V3+0] = skel[(i+1)*PL_V3+0] - skel[i*PL_V3+0]
      edge[i*PL_V3+1] = skel[(i+1)*PL_V3+1] - skel[i*PL_V3+1]
      edge[i*PL_V3+2] = skel[(i+1)*PL_V3+2] - skel[i*PL_V3+2]
      i = i + 1
  }
  i = 0
  while i < n {
      var ax: i64 = 0
      var ay: i64 = 0
      var az: i64 = 0
      var done: i64 = 0
      if i == 0 { ax = edge[0]; ay = edge[1]; az = edge[2]; done = 1 }
      if done == 0 {
          if i == n - 1 { ax = edge[(n-2)*PL_V3]; ay = edge[(n-2)*PL_V3+1]; az = edge[(n-2)*PL_V3+2]; done = 1 }
      }
      if done == 0 {
          ax = (edge[i*PL_V3] + edge[(i-1)*PL_V3]) / 2
          ay = (edge[i*PL_V3+1] + edge[(i-1)*PL_V3+1]) / 2
          az = (edge[i*PL_V3+2] + edge[(i-1)*PL_V3+2]) / 2
      }
      let l: i64 = pl_len3(ax, ay, az)
      if l == 0 { tang[i*PL_V3] = 0; tang[i*PL_V3+1] = 0; tang[i*PL_V3+2] = 0 } else {
          tang[i*PL_V3] = pl_rdiv(ax*PL_ONE, l); tang[i*PL_V3+1] = pl_rdiv(ay*PL_ONE, l); tang[i*PL_V3+2] = pl_rdiv(az*PL_ONE, l)
      }
      i = i + 1
  }
  return 0
254}
255
256// ---- DECOMPOSE: handles (nu x nv x 3, u10) -> genes. factorize_nurbs_handles + decompose_nurbs_handles. -------
257// Returns nu*nv, or a negative code for a shape the cubic surface cannot carry.
258func pl_decompose(hnd: *i64, nu: i64, nv: i64, p: *PlPart) -> i64 {
  if nu < PL_MIN_NU { return 0 - 1 }
  if nv < PL_MIN_NV { return 0 - 2 }
  if p.nu != nu { return 0 - 3 }
  if p.nv != nv { return 0 - 3 }
  var k: i64 = 0
  while k < nu*nv*PL_V3 { p.hnd[k] = hnd[k]; k = k + 1 }
  var i: i64 = 0
  while i < nu {
      var sx: i64 = 0
      var sy: i64 = 0
      var sz: i64 = 0
      var j: i64 = 0
      while j < nv { sx = sx + hnd[(i*nv+j)*PL_V3]; sy = sy + hnd[(i*nv+j)*PL_V3+1]; sz = sz + hnd[(i*nv+j)*PL_V3+2]; j = j + 1 }
      p.skel[i*PL_V3] = pl_rdiv(sx, nv); p.skel[i*PL_V3+1] = pl_rdiv(sy, nv); p.skel[i*PL_V3+2] = pl_rdiv(sz, nv)
      i = i + 1
  }
  let edge: *i64 = sys_mmap(nu*PL_V3*PL_I64) as *i64
  pl_tangents(p.skel, nu, p.tang, edge)
  let fwd: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  fwd[0] = PL_ONE; fwd[1] = 0; fwd[2] = 0
  let q: *i64 = sys_mmap(PL_Q4*PL_I64) as *i64
  let t3: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let d: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let r: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let lens: *i64 = sys_mmap(nu*PL_I64) as *i64
  i = 0
  while i < nu {
      t3[0] = p.tang[i*PL_V3]; t3[1] = p.tang[i*PL_V3+1]; t3[2] = p.tang[i*PL_V3+2]
      pl_rot_match(q, t3, fwd)
      var acc: i64 = 0
      var j: i64 = 0
      while j < nv {
          d[0] = hnd[(i*nv+j)*PL_V3] - p.skel[i*PL_V3]
          d[1] = hnd[(i*nv+j)*PL_V3+1] - p.skel[i*PL_V3+1]
          d[2] = hnd[(i*nv+j)*PL_V3+2] - p.skel[i*PL_V3+2]
          vm_q_rotate_v3(r, q, d)
          p.prof[(i*nv+j)*PL_V3] = r[0]; p.prof[(i*nv+j)*PL_V3+1] = r[1]; p.prof[(i*nv+j)*PL_V3+2] = r[2]
          acc = acc + pl_len3_fx(r[0], r[1], r[2])
          j = j + 1
      }
      var rad: i64 = pl_rdiv(acc, nv)
      if rad < PL_MIN_RAD { rad = PL_MIN_RAD }
      p.rad[i] = rad
      j = 0
      while j < nv {
          var c: i64 = 0
          while c < PL_V3 { p.prof[(i*nv+j)*PL_V3+c] = pl_rdiv(p.prof[(i*nv+j)*PL_V3+c] * PL_FQ2, rad); c = c + 1 }
          j = j + 1
      }
      i = i + 1
  }
  p.root[0] = p.skel[0]; p.root[1] = p.skel[1]; p.root[2] = p.skel[2]
  var total: i64 = 0
  i = 0
  while i < nu - 1 { lens[i] = pl_len3(edge[i*PL_V3], edge[i*PL_V3+1], edge[i*PL_V3+2]); total = total + lens[i]; i = i + 1 }
  p.length = total
  i = 0
  while i < nu - 1 {
      if total > 0 { p.prop[i] = pl_rdiv(lens[i] * PL_FQ, total) } else { p.prop[i] = 0 }
      p.theta[i] = vm_atan2(edge[i*PL_V3+2], edge[i*PL_V3])
      if lens[i] > 0 { p.yoff[i] = pl_rdiv(edge[i*PL_V3+1] * PL_FQ, lens[i]) } else { p.yoff[i] = 0 }
      i = i + 1
  }
  i = 0
  while i < nu { p.ts[i] = pl_rdiv(i * PL_FQ, nu - 1); i = i + 1 }
  pl_knots(p)
  sys_munmap(edge as *u8, nu*PL_V3*PL_I64); sys_munmap(fwd as *u8, PL_V3*PL_I64); sys_munmap(q as *u8, PL_Q4*PL_I64)
  sys_munmap(t3 as *u8, PL_V3*PL_I64); sys_munmap(d as *u8, PL_V3*PL_I64); sys_munmap(r as *u8, PL_V3*PL_I64)
  sys_munmap(lens as *u8, nu*PL_I64)
  return nu * nv
329}
330
331// ---- RECOMPOSE: genes -> skeleton -> handles. recompose_nurbs_handles + compute_profile_verts. -----------------
332func pl_recompose(p: *PlPart) -> i64 {
  let nu: i64 = p.nu
  let nv: i64 = p.nv
  if nu < PL_MIN_NU { return 0 - 1 }
  if nv < PL_MIN_NV { return 0 - 2 }
  p.skel[0] = p.root[0]; p.skel[1] = p.root[1]; p.skel[2] = p.root[2]
  var i: i64 = 0
  while i < nu - 1 {
      let len: i64 = pl_rdiv(p.length * p.prop[i], PL_FQ)
      let c: i64 = vm_cos(p.theta[i])
      let s: i64 = vm_sin(p.theta[i])
      p.skel[(i+1)*PL_V3] = p.skel[i*PL_V3] + pl_rdiv(len * c, PL_ONE)
      p.skel[(i+1)*PL_V3+1] = p.skel[i*PL_V3+1] + pl_rdiv(len * p.yoff[i], PL_FQ)
      p.skel[(i+1)*PL_V3+2] = p.skel[i*PL_V3+2] + pl_rdiv(len * s, PL_ONE)
      i = i + 1
  }
  let edge: *i64 = sys_mmap(nu*PL_V3*PL_I64) as *i64
  pl_tangents(p.skel, nu, p.tang, edge)
  let fwd: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  fwd[0] = PL_ONE; fwd[1] = 0; fwd[2] = 0
  let q: *i64 = sys_mmap(PL_Q4*PL_I64) as *i64
  let t3: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let tu: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let pos: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let d: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let r: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  i = 0
  while i < nu {
      // lerp_sample of the tangents and the skeleton at ts * (nu - 1)
      let f: i64 = p.ts[i] * (nu - 1)
      var idx: i64 = f / PL_FQ
      if idx < 0 { idx = 0 }
      if idx > nu - 1 { idx = nu - 1 }
      var rem: i64 = f - idx * PL_FQ
      if idx >= nu - 1 { rem = 0 }
      var c: i64 = 0
      while c < PL_V3 {
          let a0: i64 = p.tang[idx*PL_V3+c]
          let s0: i64 = p.skel[idx*PL_V3+c]
          var a1: i64 = a0
          var s1: i64 = s0
          if idx < nu - 1 { a1 = p.tang[(idx+1)*PL_V3+c]; s1 = p.skel[(idx+1)*PL_V3+c] }
          t3[c] = a0 + pl_rdiv((a1 - a0) * rem, PL_FQ)
          pos[c] = s0 + pl_rdiv((s1 - s0) * rem, PL_FQ)
          c = c + 1
      }
      pl_unit(tu, t3)
      pl_rot_match(q, fwd, tu)
      var j: i64 = 0
      while j < nv {
          d[0] = pl_rdiv(p.prof[(i*nv+j)*PL_V3] * p.rad[i], PL_FQ2)
          d[1] = pl_rdiv(p.prof[(i*nv+j)*PL_V3+1] * p.rad[i], PL_FQ2)
          d[2] = pl_rdiv(p.prof[(i*nv+j)*PL_V3+2] * p.rad[i], PL_FQ2)
          vm_q_rotate_v3(r, q, d)
          p.hnd[(i*nv+j)*PL_V3] = r[0] + pos[0]
          p.hnd[(i*nv+j)*PL_V3+1] = r[1] + pos[1]
          p.hnd[(i*nv+j)*PL_V3+2] = r[2] + pos[2]
          j = j + 1
      }
      i = i + 1
  }
  pl_knots(p)
  sys_munmap(edge as *u8, nu*PL_V3*PL_I64); sys_munmap(fwd as *u8, PL_V3*PL_I64); sys_munmap(q as *u8, PL_Q4*PL_I64)
  sys_munmap(t3 as *u8, PL_V3*PL_I64); sys_munmap(tu as *u8, PL_V3*PL_I64); sys_munmap(pos as *u8, PL_V3*PL_I64)
  sys_munmap(d as *u8, PL_V3*PL_I64); sys_munmap(r as *u8, PL_V3*PL_I64)
  return 0
398}
399
400// ---- THE GENES: NurbsPart.sample_params, deviation 1 in the header ----------------------------------------------
401func pl_lcg(st: *i64) -> i64 {
  st[0] = (st[0] * PL_LCG_MUL + PL_LCG_ADD) % PL_LCG_MOD
  return (st[0] >> PL_LCG_SHIFT) % PL_LCG_BITS
404}
405// uniform in [0, PL_FQ)
406func pl_urand(st: *i64) -> i64 { return pl_lcg(st) * PL_FQ / PL_LCG_BITS }
407// standard normal in Q16 by Irwin-Hall: twelve uniforms summed, centred
408func pl_nrand(st: *i64) -> i64 {
  var s: i64 = 0
  var k: i64 = 0
  while k < PL_GAUSS_N { s = s + pl_urand(st); k = k + 1 }
  return s - PL_GAUSS_N * PL_FQ / 2
413}
414// N(mean, sd * var) with sd and var_ in Q16 (a mean and sd in the caller's unit); a zero deviation returns the mean
415// exactly, so var_ = 0 is the identity by construction
416func pl_normal(st: *i64, mean: i64, sd: i64, var_: i64) -> i64 {
  let sdv: i64 = sd * var_ / PL_FQ
  if sdv == 0 { return mean }
  return mean + pl_rdiv(sdv * pl_nrand(st), PL_FQ)
420}
421// mutate the genes in place; var_ is the source's `var` in Q16 (PL_FQ = 1.0). Deterministic per (seed, var_).
422// The recomposed arrays are NOT refreshed: call pl_recompose after.
423func pl_mutate(p: *PlPart, seed: i64, var_: i64) -> i64 {
  let nu: i64 = p.nu
  let nv: i64 = p.nv
  let st: *i64 = sys_mmap(PL_I64) as *i64
  var s0: i64 = seed
  if s0 < 0 { s0 = 0 - s0 }
  st[0] = s0 % PL_LCG_MOD
  let sz: i64 = pl_normal(st, PL_FQ, PL_SD_SIZE, var_)
  p.length = pl_rdiv(pl_rdiv(p.length * sz, PL_FQ) * pl_normal(st, PL_FQ, PL_SD_LENGTH, var_), PL_FQ)
  let rshared: i64 = pl_normal(st, PL_FQ, PL_SD_RADS, var_)
  var i: i64 = 0
  while i < nu {
      var rr: i64 = pl_rdiv(p.rad[i] * sz, PL_FQ)
      rr = pl_rdiv(rr * rshared, PL_FQ)
      rr = pl_rdiv(rr * pl_normal(st, PL_FQ, PL_SD_RAD_STATION, var_), PL_FQ)
      if rr < PL_MIN_RAD { rr = PL_MIN_RAD }
      p.rad[i] = rr
      i = i + 1
  }
  i = 0
  while i < nu - 1 { p.prop[i] = pl_rdiv(p.prop[i] * pl_normal(st, PL_FQ, PL_SD_PROP, var_), PL_FQ); i = i + 1 }
  let dth: *i64 = sys_mmap(nu*PL_I64) as *i64
  var sum: i64 = 0
  i = 0
  while i < nu - 1 { dth[i] = pl_normal(st, 0, PL_SD_THETA, var_); sum = sum + dth[i]; i = i + 1 }
  let mean: i64 = sum / (nu - 1)
  i = 0
  while i < nu - 1 { p.theta[i] = p.theta[i] + dth[i] - mean; i = i + 1 }
  let col: *i64 = sys_mmap(nv*PL_I64) as *i64
  var j: i64 = 0
  while j < nv { col[j] = pl_normal(st, PL_FQ, PL_SD_PROF_COL, var_); j = j + 1 }
  let noise: *i64 = sys_mmap(nu*nv*PL_I64) as *i64
  i = 0
  while i < nu {
      j = 0
      while j < nv { noise[i*nv+j] = pl_rdiv(col[j] * pl_normal(st, PL_FQ, PL_SD_PROF_CELL, var_), PL_FQ); j = j + 1 }
      i = i + 1
  }
  // profile_noise[:, :m//2-1] = profile_noise[:, m//2:-1][:, ::-1]: column j takes column m-2-j for j < m/2-1
  i = 0
  while i < nu {
      j = 0
      while j < nv/2 - 1 { noise[i*nv+j] = noise[i*nv + (nv - 2 - j)]; j = j + 1 }
      i = i + 1
  }
  i = 0
  while i < nu {
      j = 0
      while j < nv {
          var c: i64 = 0
          while c < PL_V3 { p.prof[(i*nv+j)*PL_V3+c] = pl_rdiv(p.prof[(i*nv+j)*PL_V3+c] * noise[i*nv+j], PL_FQ); c = c + 1 }
          j = j + 1
      }
      i = i + 1
  }
  sys_munmap(st as *u8, PL_I64); sys_munmap(dth as *u8, nu*PL_I64); sys_munmap(col as *u8, nv*PL_I64)
  sys_munmap(noise as *u8, nu*nv*PL_I64)
  return 0
481}
482
483// ---- THE B-SPLINE SURFACE, degree (3,3): geomdl's evaluation in fixed point ----------------------------------
484// span k with U[k] <= x < U[k+1] over the valid domain U[3] .. U[ncp]; x at or past the end lands on the last span
485func pl_span(U: *i64, ncp: i64, x: i64) -> i64 {
  if x >= U[ncp] { return ncp - 1 }
  var k: i64 = PL_DEG
  while k < ncp - 1 {
      if x >= U[k+1] { k = k + 1 } else { return k }
  }
  return k
492}
493// the four nonzero cubic basis functions at x on span k (Piegl-Tiller A2.2), Q16; left/right are 4-word scratch
494func pl_basis(U: *i64, k: i64, x: i64, N: *i64, left: *i64, right: *i64) -> i64 {
  N[0] = PL_FQ
  var j: i64 = 1
  while j <= PL_DEG {
      left[j] = x - U[k+1-j]
      right[j] = U[k+j] - x
      var saved: i64 = 0
      var r: i64 = 0
      while r < j {
          let den: i64 = right[r+1] + left[j-r]
          var temp: i64 = 0
          if den != 0 { temp = N[r] * PL_FQ / den }
          N[r] = saved + right[r+1] * temp / PL_FQ
          saved = left[j-r] * temp / PL_FQ
          r = r + 1
      }
      N[j] = saved
      j = j + 1
  }
  return 0
514}
515// S(u, v) in u10 for u, v in [0, PL_FQ]; v maps onto the periodic domain [kv[3], kv[nv+3]], where the curve closes
516func pl_surface(p: *PlPart, u: i64, v: i64, out: *i64) -> i64 {
  let nu: i64 = p.nu
  let nv: i64 = p.nv
  var x: i64 = u
  if x < 0 { x = 0 }
  if x > PL_FQ { x = PL_FQ }
  var vv: i64 = v
  if vv < 0 { vv = 0 }
  if vv > PL_FQ { vv = PL_FQ }
  let vs: i64 = p.kv[PL_DEG]
  let ve: i64 = p.kv[nv + PL_DEG]
  let y: i64 = vs + (ve - vs) * vv / PL_FQ
  let ku: i64 = pl_span(p.ku, nu, x)
  let kv: i64 = pl_span(p.kv, nv + PL_DEG, y)
  pl_basis(p.ku, ku, x, p.bu, p.bl, p.br)
  pl_basis(p.kv, kv, y, p.bv, p.bl, p.br)
  var c: i64 = 0
  while c < PL_V3 {
      var acc: i64 = 0
      var a: i64 = 0
      while a < PL_ORD {
          let row: i64 = ku - PL_DEG + a
          var b: i64 = 0
          while b < PL_ORD {
              let colm: i64 = (kv - PL_DEG + b) % nv
              acc = acc + p.bu[a] * p.bv[b] / PL_FQ * p.hnd[(row*nv + colm)*PL_V3 + c]
              b = b + 1
          }
          a = a + 1
      }
      out[c] = pl_rdiv(acc, PL_FQ)
      c = c + 1
  }
  return 0
550}
551
552// ---- THE ATTACHMENT SKELETON: the interior rows, cubic where four exist (deviation 3), linear otherwise --------
553// t in [0, PL_FQ]; out = point (u10). Returns PL_ASKEL_SPLINE or PL_ASKEL_LINEAR.
554func pl_askel_point(p: *PlPart, t: i64, out: *i64) -> i64 {
  let na: i64 = p.nu - 2
  var tt: i64 = t
  if tt < 0 { tt = 0 }
  if tt > PL_FQ { tt = PL_FQ }
  if na >= PL_ORD {
      let k: i64 = pl_span(p.ka, na, tt)
      pl_basis(p.ka, k, tt, p.bu, p.bl, p.br)
      var c: i64 = 0
      while c < PL_V3 {
          var acc: i64 = 0
          var a: i64 = 0
          while a < PL_ORD {
              let row: i64 = 1 + k - PL_DEG + a
              acc = acc + p.bu[a] * p.skel[row*PL_V3 + c]
              a = a + 1
          }
          out[c] = pl_rdiv(acc, PL_FQ)
          c = c + 1
      }
      return PL_ASKEL_SPLINE
  }
  let f: i64 = tt * (na - 1)
  var idx: i64 = f / PL_FQ
  if idx > na - 2 { idx = na - 2 }
  if idx < 0 { idx = 0 }
  let rem: i64 = f - idx * PL_FQ
  var c2: i64 = 0
  while c2 < PL_V3 {
      let s0: i64 = p.skel[(1+idx)*PL_V3 + c2]
      let s1: i64 = p.skel[(2+idx)*PL_V3 + c2]
      out[c2] = s0 + pl_rdiv((s1 - s0) * rem, PL_FQ)
      c2 = c2 + 1
  }
  return PL_ASKEL_LINEAR
589}
590// the attachment skeleton's point and unit tangent (Q12) at t
591func pl_askel(p: *PlPart, t: i64, out: *i64, tan: *i64) -> i64 {
  let which: i64 = pl_askel_point(p, t, out)
  let a: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let b: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  var t0: i64 = t - PL_FD_STEP
  var t1: i64 = t + PL_FD_STEP
  if t0 < 0 { t0 = 0 }
  if t1 > PL_FQ { t1 = PL_FQ }
  pl_askel_point(p, t0, a)
  pl_askel_point(p, t1, b)
  a[0] = b[0] - a[0]; a[1] = b[1] - a[1]; a[2] = b[2] - a[2]
  pl_unit(tan, a)
  sys_munmap(a as *u8, PL_V3*PL_I64); sys_munmap(b as *u8, PL_V3*PL_I64)
  return which
605}
606
607// the squared perpendicular distance (u10^2) from the ray through o along the unit direction dw (Q12) to the surface
608// point at (u, v); a point behind the origin scores PL_POS_HUGE. Smaller is better: the search minimises it, and the
609// u10 quantisation of the surface bounds the search rather than the cosine's quantisation (a Q12 cosine reads 1.0 for
610// every direction within 22 milliradians -- MEASURED as a 4 u10 lateral miss on the first run of the gate).
611func pl_ray_perp2(p: *PlPart, u: i64, v: i64, o: *i64, dw: *i64, s: *i64) -> i64 {
  pl_surface(p, u, v, s)
  let dx: i64 = s[0] - o[0]
  let dy: i64 = s[1] - o[1]
  let dz: i64 = s[2] - o[2]
  let proj: i64 = (dx*dw[0] + dy*dw[1] + dz*dw[2]) / PL_ONE
  if proj <= 0 { return PL_POS_HUGE }
  return dx*dx + dy*dy + dz*dz - proj*proj
619}
620// ---- SURFACE COORDINATES -> WORLD: apply_attach_transform's ray, deviation 2 in the header --------------------
621// t in [0, PL_FQ] along the attachment skeleton; yaw in [-PL_FQ, PL_FQ], a half-turn per unit about the tangent from
622// straight down (-z in the part frame): 0 = down, +1/2 = +y, +1 = up, -1/2 = -y; rad in [0, PL_FQ] from the skeleton
623// (0) to the surface (1). out = location (u10), nrm = outward unit normal at the hit (Q12), tan = unit tangent (Q12),
624// info = [cosine of ray and hit (Q16), u of the hit, v of the hit, skeleton-to-hit distance u10, askel path]. Returns 0.
625func pl_surface_at(p: *PlPart, t: i64, yaw: i64, rad: i64, out: *i64, nrm: *i64, tan: *i64, info: *i64) -> i64 {
  let nu: i64 = p.nu
  let o: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let which: i64 = pl_askel(p, t, o, tan)
  let fwd: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  fwd[0] = PL_ONE; fwd[1] = 0; fwd[2] = 0
  let d0: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let ang: i64 = VM_PI * yaw / PL_FQ
  d0[0] = 0; d0[1] = vm_sin(ang); d0[2] = 0 - vm_cos(ang)
  let q: *i64 = sys_mmap(PL_Q4*PL_I64) as *i64
  pl_rot_match(q, fwd, tan)
  let dw: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  vm_q_rotate_v3(dw, q, d0)
  let s: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  // the station's u on the full surface: interior row 1 is t = 0, interior row nu-2 is t = 1
  var tt: i64 = t
  if tt < 0 { tt = 0 }
  if tt > PL_FQ { tt = PL_FQ }
  let uc: i64 = (PL_FQ + tt * (nu - 3)) / (nu - 1)
  let seg: i64 = PL_FQ / (nu - 1)
  let cell: i64 = PL_FQ / PL_SEARCH_COARSE
  var best: i64 = PL_POS_HUGE
  var bu: i64 = uc
  var bv: i64 = 0
  var kk: i64 = 0 - PL_SEARCH_UHALF
  while kk <= PL_SEARCH_UHALF {
      var uu: i64 = uc + kk * seg / PL_SEARCH_UQUART
      if uu < 0 { uu = 0 }
      if uu > PL_FQ { uu = PL_FQ }
      var m: i64 = 0
      while m < PL_SEARCH_COARSE {
          let vv: i64 = m * cell
          let sc: i64 = pl_ray_perp2(p, uu, vv, o, dw, s)
          if sc < best { best = sc; bu = uu; bv = vv }
          m = m + 1
      }
      kk = kk + 1
  }
  // bisection on v about the best coarse sample: the better of the two neighbouring midpoints, halving each step
  var lo: i64 = bv - cell
  var hi: i64 = bv + cell
  var step: i64 = 0
  while step < PL_SEARCH_REFINE {
      let ml: i64 = (lo + bv) / 2
      let mh: i64 = (bv + hi) / 2
      var vl: i64 = ml
      if vl < 0 { vl = vl + PL_FQ }
      var vh: i64 = mh
      if vh >= PL_FQ { vh = vh - PL_FQ }
      let sl: i64 = pl_ray_perp2(p, bu, vl, o, dw, s)
      let sh: i64 = pl_ray_perp2(p, bu, vh, o, dw, s)
      var moved: i64 = 0
      if sl < best { if sl <= sh { best = sl; hi = bv; bv = ml; moved = 1 } }
      if moved == 0 { if sh < best { best = sh; lo = bv; bv = mh; moved = 1 } }
      if moved == 0 { lo = ml; hi = mh }
      step = step + 1
  }
  if bv < 0 { bv = bv + PL_FQ }
  if bv >= PL_FQ { bv = bv - PL_FQ }
  pl_surface(p, bu, bv, s)
  let dist: i64 = pl_len3(s[0]-o[0], s[1]-o[1], s[2]-o[2])
  var cosq: i64 = 0
  if dist > 0 { cosq = pl_rdiv(((s[0]-o[0])*dw[0] + (s[1]-o[1])*dw[1] + (s[2]-o[2])*dw[2]) / PL_ONE * PL_FQ, dist) }
  out[0] = o[0] + pl_rdiv((s[0] - o[0]) * rad, PL_FQ)
  out[1] = o[1] + pl_rdiv((s[1] - o[1]) * rad, PL_FQ)
  out[2] = o[2] + pl_rdiv((s[2] - o[2]) * rad, PL_FQ)
  // the normal: central differences along u and v, crossed, oriented along the ray
  let su: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let sv: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  let tmp: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  var ua: i64 = bu - PL_FD_STEP
  var ub: i64 = bu + PL_FD_STEP
  if ua < 0 { ua = 0 }
  if ub > PL_FQ { ub = PL_FQ }
  pl_surface(p, ua, bv, tmp)
  pl_surface(p, ub, bv, su)
  su[0] = su[0] - tmp[0]; su[1] = su[1] - tmp[1]; su[2] = su[2] - tmp[2]
  var va: i64 = bv - PL_FD_STEP
  var vb: i64 = bv + PL_FD_STEP
  if va < 0 { va = va + PL_FQ }
  if vb >= PL_FQ { vb = vb - PL_FQ }
  pl_surface(p, bu, va, tmp)
  pl_surface(p, bu, vb, sv)
  sv[0] = sv[0] - tmp[0]; sv[1] = sv[1] - tmp[1]; sv[2] = sv[2] - tmp[2]
  tmp[0] = su[1]*sv[2] - su[2]*sv[1]
  tmp[1] = su[2]*sv[0] - su[0]*sv[2]
  tmp[2] = su[0]*sv[1] - su[1]*sv[0]
  pl_unit(nrm, tmp)
  if nrm[0]*dw[0] + nrm[1]*dw[1] + nrm[2]*dw[2] < 0 { nrm[0] = 0 - nrm[0]; nrm[1] = 0 - nrm[1]; nrm[2] = 0 - nrm[2] }
  info[0] = cosq; info[1] = bu; info[2] = bv; info[3] = dist; info[4] = which
  sys_munmap(o as *u8, PL_V3*PL_I64); sys_munmap(fwd as *u8, PL_V3*PL_I64); sys_munmap(d0 as *u8, PL_V3*PL_I64)
  sys_munmap(q as *u8, PL_Q4*PL_I64); sys_munmap(dw as *u8, PL_V3*PL_I64); sys_munmap(s as *u8, PL_V3*PL_I64)
  sys_munmap(su as *u8, PL_V3*PL_I64); sys_munmap(sv as *u8, PL_V3*PL_I64); sys_munmap(tmp as *u8, PL_V3*PL_I64)
  return 0
719}
720
721// ---- THE MESH: nurbs.py nurbs(..., method='geomdl', face_size) ---------------------------------------------
722// out = [ulength, vlength, sample_size, ring_points, ntris]. ulength is the longest u-polyline over the columns,
723// vlength the longest OPEN v-polyline over the rows (np.diff does not wrap), sample_size = L / face + 1 for both
724// parameters (one delta), ring_points = sample_size - 1 (the closing duplicate dropped), ntris counts the tube's two
725// triangles per quad plus the two cap fans. Returns ntris, or -1 when the face size leaves fewer than three ring
726// points (a cap fan cannot exist) -- a NAMED refusal, never a silent floor.
727func pl_mesh_counts(p: *PlPart, face: i64, out: *i64) -> i64 {
  let nu: i64 = p.nu
  let nv: i64 = p.nv
  var ul: i64 = 0
  var j: i64 = 0
  while j < nv {
      var acc: i64 = 0
      var i: i64 = 0
      while i < nu - 1 {
          let a: i64 = (i*nv+j)*PL_V3
          let b: i64 = ((i+1)*nv+j)*PL_V3
          acc = acc + pl_len3(p.hnd[b]-p.hnd[a], p.hnd[b+1]-p.hnd[a+1], p.hnd[b+2]-p.hnd[a+2])
          i = i + 1
      }
      if acc > ul { ul = acc }
      j = j + 1
  }
  var vl: i64 = 0
  var i2: i64 = 0
  while i2 < nu {
      var acc2: i64 = 0
      var j2: i64 = 0
      while j2 < nv - 1 {
          let a: i64 = (i2*nv+j2)*PL_V3
          let b: i64 = (i2*nv+j2+1)*PL_V3
          acc2 = acc2 + pl_len3(p.hnd[b]-p.hnd[a], p.hnd[b+1]-p.hnd[a+1], p.hnd[b+2]-p.hnd[a+2])
          j2 = j2 + 1
      }
      if acc2 > vl { vl = acc2 }
      i2 = i2 + 1
  }
  var L: i64 = ul
  if vl > L { L = vl }
  var f: i64 = face
  if f < 1 { f = 1 }
  let ss: i64 = L / f + 1
  let m: i64 = ss - 1
  out[0] = ul; out[1] = vl; out[2] = ss; out[3] = m
  if m < PL_MIN_RING { out[4] = 0 - 1; return 0 - 1 }
  let nt: i64 = 2 * (ss - 1) * m + 2 * m
  out[4] = nt
  return nt
769}
770// one triangle into the NXMSH2 buffer: winding flipped so the face normal points along `outward`
771func pl_put_tri(B: *u8, tb: i64, t: i64, V: *i64, ia: i64, ib: i64, ic: i64, ox: i64, oy: i64, oz: i64, xyz9: *i64) -> i64 {
  var b: i64 = ib
  var c: i64 = ic
  let e1x: i64 = V[ib*PL_V3] - V[ia*PL_V3]
  let e1y: i64 = V[ib*PL_V3+1] - V[ia*PL_V3+1]
  let e1z: i64 = V[ib*PL_V3+2] - V[ia*PL_V3+2]
  let e2x: i64 = V[ic*PL_V3] - V[ia*PL_V3]
  let e2y: i64 = V[ic*PL_V3+1] - V[ia*PL_V3+1]
  let e2z: i64 = V[ic*PL_V3+2] - V[ia*PL_V3+2]
  var nx: i64 = e1y*e2z - e1z*e2y
  var ny: i64 = e1z*e2x - e1x*e2z
  var nz: i64 = e1x*e2y - e1y*e2x
  if nx*ox + ny*oy + nz*oz < 0 { b = ic; c = ib; nx = 0 - nx; ny = 0 - ny; nz = 0 - nz }
  xyz9[0] = V[ia*PL_V3]; xyz9[1] = V[ia*PL_V3+1]; xyz9[2] = V[ia*PL_V3+2]
  xyz9[3] = V[b*PL_V3]; xyz9[4] = V[b*PL_V3+1]; xyz9[5] = V[b*PL_V3+2]
  xyz9[6] = V[c*PL_V3]; xyz9[7] = V[c*PL_V3+1]; xyz9[8] = V[c*PL_V3+2]
  nm_put_tri(B, tb, t, xyz9)
  let l: i64 = pl_len3(nx, ny, nz)
  var qx: i64 = 0
  var qy: i64 = 0
  var qz: i64 = 0
  if l > 0 { qx = pl_rdiv(nx*PL_ONE, l); qy = pl_rdiv(ny*PL_ONE, l); qz = pl_rdiv(nz*PL_ONE, l) }
  var v: i64 = 0
  while v < PL_V3 {
      let o: i64 = tb + t*NM_TRI_REC + PL_OFF_NORMALS + v*NM_VERT_STRIDE
      nm_put_u32(B, o, vm_int_to_f32(qx, PL_ONE))
      nm_put_u32(B, o + NM_F32_BYTES, vm_int_to_f32(qy, PL_ONE))
      nm_put_u32(B, o + NM_F32_BYTES*2, vm_int_to_f32(qz, PL_ONE))
      v = v + 1
  }
  nm_col_put(B, tb, t, PL_COL_BONE_R, PL_COL_BONE_G, PL_COL_BONE_B)
  return 0
803}
804// write the part as a one-layer NXMSH2 at `path`; rep = [ntris, bytes, sample_size, ring_points]. Returns ntris, or
805// -1 when the face size is refused by pl_mesh_counts, -2 when the file cannot be opened.
806func pl_mesh_write(p: *PlPart, face: i64, path: *u8, rep: *i64) -> i64 {
  let cnt: *i64 = sys_mmap(PL_I64*8) as *i64
  let nt: i64 = pl_mesh_counts(p, face, cnt)
  if nt < 0 { rep[0] = 0 - 1; rep[1] = 0; rep[2] = cnt[2]; rep[3] = cnt[3]; return 0 - 1 }
  let ss: i64 = cnt[2]
  let m: i64 = cnt[3]
  let V: *i64 = sys_mmap(ss*m*PL_V3*PL_I64) as *i64
  let cen: *i64 = sys_mmap(ss*PL_V3*PL_I64) as *i64
  let s: *i64 = sys_mmap(PL_V3*PL_I64) as *i64
  var i: i64 = 0
  while i < ss {
      let u: i64 = i * PL_FQ / (ss - 1)
      var cx: i64 = 0
      var cy: i64 = 0
      var cz: i64 = 0
      var j: i64 = 0
      while j < m {
          let v: i64 = j * PL_FQ / m
          pl_surface(p, u, v, s)
          V[(i*m+j)*PL_V3] = s[0]; V[(i*m+j)*PL_V3+1] = s[1]; V[(i*m+j)*PL_V3+2] = s[2]
          cx = cx + s[0]; cy = cy + s[1]; cz = cz + s[2]
          j = j + 1
      }
      cen[i*PL_V3] = pl_rdiv(cx, m); cen[i*PL_V3+1] = pl_rdiv(cy, m); cen[i*PL_V3+2] = pl_rdiv(cz, m)
      i = i + 1
  }
  let bytes: i64 = nm_file_bytes(1, nt)
  let B: *u8 = sys_mmap(bytes + PL_SAVE_HEADER)
  B[0] = PL_MAGIC_0 as u8; B[1] = PL_MAGIC_1 as u8; B[2] = PL_MAGIC_2 as u8; B[3] = PL_MAGIC_3 as u8
  B[4] = PL_MAGIC_4 as u8; B[5] = PL_MAGIC_5 as u8; B[6] = 0 as u8; B[7] = 0 as u8
  nm_put_u32(B, NM_OFF_NLAYERS, 1)
  nm_put_u32(B, NM_OFF_NTRIS, nt)
  B[NM_HDR] = 112 as u8; B[NM_HDR+1] = 97 as u8; B[NM_HDR+2] = 114 as u8; B[NM_HDR+3] = 116 as u8   // "part"
  nm_put_u32(B, NM_HDR + PL_LAYER_START_OFF, 0)
  nm_put_u32(B, NM_HDR + PL_LAYER_COUNT_OFF, nt)
  let tb: i64 = nm_tri_base(B)
  let xyz9: *i64 = sys_mmap(9*PL_I64) as *i64
  // the tube: quads (i,j) (i,j+1) (i+1,j+1) (i+1,j), each split at its first diagonal, cyclic in j
  var t: i64 = 0
  i = 0
  while i < ss - 1 {
      var j: i64 = 0
      while j < m {
          let jn: i64 = (j + 1) % m
          let a: i64 = i*m + j
          let b: i64 = i*m + jn
          let c: i64 = (i+1)*m + jn
          let d: i64 = (i+1)*m + j
          // outward for this quad: from the ring centre at row i to vertex a
          let ox: i64 = V[a*PL_V3] - cen[i*PL_V3]
          let oy: i64 = V[a*PL_V3+1] - cen[i*PL_V3+1]
          let oz: i64 = V[a*PL_V3+2] - cen[i*PL_V3+2]
          pl_put_tri(B, tb, t, V, a, b, c, ox, oy, oz, xyz9)
          t = t + 1
          pl_put_tri(B, tb, t, V, a, c, d, ox, oy, oz, xyz9)
          t = t + 1
          j = j + 1
      }
      i = i + 1
  }
  // the caps: fans from the ring centre, outward along the tube's own axis (deviation 4)
  let CV: *i64 = sys_mmap((m + 1)*PL_V3*PL_I64) as *i64
  var e: i64 = 0
  while e < 2 {
      var row: i64 = 0
      var ox: i64 = cen[0] - cen[PL_V3]
      var oy: i64 = cen[1] - cen[PL_V3+1]
      var oz: i64 = cen[2] - cen[PL_V3+2]
      if e == 1 {
          row = ss - 1
          ox = cen[(ss-1)*PL_V3] - cen[(ss-2)*PL_V3]
          oy = cen[(ss-1)*PL_V3+1] - cen[(ss-2)*PL_V3+1]
          oz = cen[(ss-1)*PL_V3+2] - cen[(ss-2)*PL_V3+2]
      }
      var j: i64 = 0
      while j < m {
          CV[j*PL_V3] = V[(row*m+j)*PL_V3]; CV[j*PL_V3+1] = V[(row*m+j)*PL_V3+1]; CV[j*PL_V3+2] = V[(row*m+j)*PL_V3+2]
          j = j + 1
      }
      CV[m*PL_V3] = cen[row*PL_V3]; CV[m*PL_V3+1] = cen[row*PL_V3+1]; CV[m*PL_V3+2] = cen[row*PL_V3+2]
      j = 0
      while j < m {
          pl_put_tri(B, tb, t, CV, m, j, (j + 1) % m, ox, oy, oz, xyz9)
          t = t + 1
          j = j + 1
      }
      e = e + 1
  }
  let fd: i64 = sys_openat_wr(path, PL_MODE644)
  if fd < 0 { rep[0] = 0 - 2; rep[1] = bytes; rep[2] = ss; rep[3] = m; return 0 - 2 }
  sys_write(fd, B, bytes)
  sys_close(fd)
  rep[0] = nt; rep[1] = bytes; rep[2] = ss; rep[3] = m
  sys_munmap(V as *u8, ss*m*PL_V3*PL_I64); sys_munmap(cen as *u8, ss*PL_V3*PL_I64); sys_munmap(s as *u8, PL_V3*PL_I64)
  sys_munmap(B, bytes + PL_SAVE_HEADER); sys_munmap(xyz9 as *u8, 9*PL_I64); sys_munmap(CV as *u8, (m + 1)*PL_V3*PL_I64)
  sys_munmap(cnt as *u8, PL_I64*8)
  return nt
903}
904
905// ---- THE PART FILE: the genes as whitespace-separated integers, so a part is DATA a board can carry -------------
906// NXPART <version>\n <nu> <nv>\n <root x y z>\n <length>\n <prop nu-1>\n <theta nu-1>\n <yoff nu-1>\n <rad nu>\n
907// <ts nu>\n <prof nu*nv*3>\n -- every value in this file's units (u10, Q16 fractions, Q16 u10 radii, VM_ONE radians).
908func pl_save_ints(fb: *u8, fo: i64, A: *i64, n: i64) -> i64 {
  var o: i64 = fo
  var i: i64 = 0
  while i < n {
      o = nxi_buf(fb, o, A[i])
      fb[o] = 32 as u8
      o = o + 1
      i = i + 1
  }
  fb[o] = 10 as u8
  return o + 1
919}
920func pl_save(path: *u8, p: *PlPart) -> i64 {
  let nu: i64 = p.nu
  let nv: i64 = p.nv
  let nints: i64 = PL_HEADER_INTS + 3*(nu-1) + 2*nu + nu*nv*PL_V3
  let cap: i64 = nints*PL_SAVE_DIGITS + PL_SAVE_HEADER
  let fb: *u8 = sys_mmap(cap)
  var fo: i64 = 0
  fb[0] = 78 as u8; fb[1] = 88 as u8; fb[2] = 80 as u8; fb[3] = 65 as u8; fb[4] = 82 as u8; fb[5] = 84 as u8; fb[6] = 32 as u8
  fo = 7
  fo = nxi_buf(fb, fo, PL_FMT_VERSION)
  fb[fo] = 10 as u8
  fo = fo + 1
  let hdr: *i64 = sys_mmap(PL_I64*2) as *i64
  hdr[0] = nu; hdr[1] = nv
  fo = pl_save_ints(fb, fo, hdr, 2)
  fo = pl_save_ints(fb, fo, p.root, PL_V3)
  hdr[0] = p.length
  fo = pl_save_ints(fb, fo, hdr, 1)
  fo = pl_save_ints(fb, fo, p.prop, nu - 1)
  fo = pl_save_ints(fb, fo, p.theta, nu - 1)
  fo = pl_save_ints(fb, fo, p.yoff, nu - 1)
  fo = pl_save_ints(fb, fo, p.rad, nu)
  fo = pl_save_ints(fb, fo, p.ts, nu)
  fo = pl_save_ints(fb, fo, p.prof, nu*nv*PL_V3)
  let fd: i64 = sys_openat_wr(path, PL_MODE644)
  if fd < 0 { return 0 - 1 }
  sys_write(fd, fb, fo)
  sys_close(fd)
  sys_munmap(fb, cap); sys_munmap(hdr as *u8, PL_I64*2)
  return fo
950}
951// the next integer at or after pos[0]; returns 1 and advances, or 0 at the end of the buffer
952func pl_tok(buf: *u8, n: i64, pos: *i64, out: *i64) -> i64 {
  var i: i64 = pos[0]
  while i < n {
      var c: i64 = buf[i] as i64
      var neg: i64 = 0
      if c == 45 {
          neg = 1
          i = i + 1
          if i < n { c = buf[i] as i64 } else { c = 0 }
      }
      var isd: i64 = 0
      if c >= 48 { if c <= 57 { isd = 1 } }
      if isd == 1 {
          var v: i64 = 0
          var run: i64 = 1
          while run == 1 {
              v = v*10 + (c - 48)
              i = i + 1
              if i < n { c = buf[i] as i64 } else { c = 0 }
              var still: i64 = 0
              if c >= 48 { if c <= 57 { still = 1 } }
              if still == 0 { run = 0 }
          }
          if neg == 1 { v = 0 - v }
          out[0] = v
          pos[0] = i
          return 1
      }
      i = i + 1
  }
  pos[0] = i
  return 0
984}
985func pl_load_ints(buf: *u8, n: i64, pos: *i64, A: *i64, cnt: i64, tok: *i64) -> i64 {
  var i: i64 = 0
  while i < cnt {
      if pl_tok(buf, n, pos, tok) == 0 { return 0 - 1 }
      A[i] = tok[0]
      i = i + 1
  }
  return cnt
993}
994// load a part file; returns the part, recomposed and ready to evaluate, or a null pointer when the file is absent,
995// carries another version, or ends before its declared count
996func pl_load(path: *u8) -> *PlPart {
  let ln: *i64 = sys_mmap(PL_LOAD_HDR) as *i64
  let buf: *u8 = sys_read_file(path, ln)
  if (buf as i64) == 0 { return 0 as *PlPart }
  let n: i64 = ln[0]
  let pos: *i64 = sys_mmap(PL_I64) as *i64
  let tok: *i64 = sys_mmap(PL_I64) as *i64
  pos[0] = 0
  if pl_tok(buf, n, pos, tok) == 0 { return 0 as *PlPart }
  if tok[0] != PL_FMT_VERSION { return 0 as *PlPart }
  if pl_tok(buf, n, pos, tok) == 0 { return 0 as *PlPart }
  let nu: i64 = tok[0]
  if pl_tok(buf, n, pos, tok) == 0 { return 0 as *PlPart }
  let nv: i64 = tok[0]
  if nu < PL_MIN_NU { return 0 as *PlPart }
  if nv < PL_MIN_NV { return 0 as *PlPart }
  let p: *PlPart = pl_new(nu, nv)
  if pl_load_ints(buf, n, pos, p.root, PL_V3, tok) < 0 { return 0 as *PlPart }
  if pl_tok(buf, n, pos, tok) == 0 { return 0 as *PlPart }
  p.length = tok[0]
  if pl_load_ints(buf, n, pos, p.prop, nu - 1, tok) < 0 { return 0 as *PlPart }
  if pl_load_ints(buf, n, pos, p.theta, nu - 1, tok) < 0 { return 0 as *PlPart }
  if pl_load_ints(buf, n, pos, p.yoff, nu - 1, tok) < 0 { return 0 as *PlPart }
  if pl_load_ints(buf, n, pos, p.rad, nu, tok) < 0 { return 0 as *PlPart }
  if pl_load_ints(buf, n, pos, p.ts, nu, tok) < 0 { return 0 as *PlPart }
  if pl_load_ints(buf, n, pos, p.prof, nu*nv*PL_V3, tok) < 0 { return 0 as *PlPart }
  pl_recompose(p)
  sys_munmap(ln as *u8, PL_LOAD_HDR); sys_munmap(pos as *u8, PL_I64); sys_munmap(tok as *u8, PL_I64)
  return p
1025}
nx_partloft_lib.nx source

`nx_partloft_lib.nx` source
