-
Notifications
You must be signed in to change notification settings - Fork 113
max_ik_solvers
title: Max IK Solvers description: published: true date: 2026-07-08T12:07:38.367Z tags: editor: markdown dateCreated: 2026-07-08T12:07:38.367Z
Gap tags. The following tag marks behavior not yet fully normative in this document:
- [CORE] — behavior implemented inside Max's binaries (geom/core DLLs), for which no source and no runnable Max installation exists. These are pinned by characterization against reference exports; candidate algorithms are given where strong leads exist.
Reference numerics are VS2008 x86, /fp:precise, x87 with control word 0x27F (53-bit precision, round-to-nearest), CRT transcendentals from msvcr90.dll. Any native Linux build intended to be bit-identical must be 32-bit -mfpmath=387, control word 0x27F, contraction and re-association disabled, with CRT-compatible sin/cos/acos/sqrt/fabs shims. Additionally, the SDK float-math primitives of §2.9 compile to x87 instruction sequences (not CRT calls) in the reference x86 build; a native build must reproduce those sequences as specified there rather than substituting libm.
Normative clarification. Under this environment, float subexpressions are evaluated on the x87 stack at 53-bit precision and are rounded to 32-bit float only at assignments, casts, and function-call boundaries taking float parameters. An inline expression like x*x+y*y+z*z therefore accumulates effectively in double. This is why Length() and FLength() differ (§2.6, §2.9) and why an SSE build cannot be bit-identical without shims. Spill points (where the compiler stores a float temporary to memory) are codegen-dependent; this is the class of divergence the §1.2 route eliminates.
Since the VS2008 toolchain runs under Wine, the cheapest path to bit-perfection is to compile the clean-room math kernel itself with cl.exe (x86, optimization and /fp switches taken from the shipped sample .vcproj) and run it under Wine inside the Linux pipeline, either in-process via a Wine bridge or as a subprocess oracle. This eliminates the entire class of "same C, different codegen" divergences (register spills, double-rounding sites, CRT transcendentals) in one move. A native Linux build can be introduced later and promoted only once it is bit-identical to the Wine build over the full verification corpus.
Matrices are 4×3 affine, row-vector convention: p′ = p·M; in A·B, A applies first. Rows 0–2 basis, row 3 translation. % is dot product, ^ is cross product. All quantities are 32-bit float except where noted. Degeneracy tests at call sites compare the LengthUnify() return against literal 0.0f with ==/!= exactly as written, except the spline solver's explicit 100.0f*FLT_EPSILON thresholds; note the compared value is produced by a [CORE] routine (§2.6).
inline double acos_safe(float x)
{
return x <= -1.0f ? PI :
x >= 1.0f ? 0.0 : acos((double)x);
}PI is the double constant π from the SDK trig header. The function returns double; every call site assigns the result to a float, so truncation to float happens at the assignment, not inside the function. The comparisons are float-literal comparisons against the float argument. A double overload (identical with -1.0/1.0 literals) exists but is not on the live paths.
inline Quat MakeRowQuat(const Point3& axis, float angle)
{
double halfAngle = 0.5 * angle; /* float promoted to double */
double c = cos(halfAngle); /* CRT double cos */
double s = -sin(halfAngle); /* CRT double sin, negated */
return Quat((float)s*axis.x, (float)s*axis.y, (float)s*axis.z, (float)c);
}Order of operations matters: s is truncated to float first, then multiplied by the float axis component (float × float), per component. The negated sine encodes the transposed (row-vector) rotation convention. The negation here also implies the library Quat(AngAxis) constructor does not negate — a useful constraint when characterizing that constructor (§2.8). Note these are the CRT double sin/cos, not the §2.9 float Sin/Cos — do not substitute one for the other.
The quaternion–pure-vector product is a solver-local operator with these exact component formulas (q = (x,y,z,w), v a 3-vector) — normative:
ret.w = - q.x*v.x - q.y*v.y - q.z*v.z;
ret.x = q.w*v.x + q.y*v.z - q.z*v.y;
ret.y = q.w*v.y + q.z*v.x - q.x*v.z;
ret.z = q.w*v.z + q.x*v.y - q.y*v.x;inline Point3 ApplyRowQuat(const Point3& p, const Quat& q)
{
Quat result = (q.Conjugate() * p) * q;
return Point3(result.x, result.y, result.z);
}The first product uses the local operator above. Conjugate() and the second product Quat::operator*(Quat) are declared in the public header but implemented in Max's core DLL, and are therefore [CORE]:
-
Conjugate()— documented contract: returns the conjugate. Near-certain candidate(-x,-y,-z,w); verify once with a single vector, then treat as pinned. -
operator*(Quat)— must be characterized before the swivel-angle paths can be trusted. Candidate: Max historically composes left-to-right (q1*q2 applies q1's rotation first), i.e. the reverse of the textbook Hamilton product. The SDK documentation's statement that the API rotation convention is the left-hand rule (opposite the UI) is consistent with this. Discriminating vectors: two 90° rotations about distinct axes; the two conventions produce opposite-handed compositions, so one vector pair settles it. Characterize jointly withQuat(AngAxis)andRotationValue::PreApplyTo, since only their composition is observable in exports.
RotAngleAxisMatrix(axis, angle), MatrixToEuler(M, eul, EULERTYPE_XYZ), Matrix3::Invert, Matrix3::operator*, Matrix3::PreRotateX/Y/Z, VectorTransform, and the Quat(AngAxis) conversion are [CORE]. Documented contracts and characterization leads:
-
PreRotateX/Y/Z(angle)is exactly multiplication on the left by the corresponding axis rotation matrix;RotateX/Y/Zmultiplies on the right. (Documented contract; the arithmetic remains [CORE].) -
RotAngleAxisMatrixexpects a normalized axis, and its angle direction is documented as opposite that ofAngAxisFromQ()— a free sign constraint for characterization. -
VectorTransformapplies the 3×3 part only (no translation), in row-vector convention; both argument orders are documented equivalent. -
MatrixToEuler: Max's Euler machinery derives from Shoemake's public Graphics Gems IVEulerAngles.c. Determine which Shoemake order code EULERTYPE_XYZ maps to, and whether the matrix is transposed on entry, using near-gimbal poses from reference exports as the discriminating vectors. -
RotationValue::PreApplyTo(Matrix3&)(spline twist path): functionally pre-applies the stored rotation to the matrix; characterize the exact matrix-from-quat formula and multiplication side together with the §2.3 product convention. NoteQuat::MakeMatrixhas a documented transpose flag whose default reproduces pre-Max-4 handedness — expect one of the two variants to appear here.
Chains evaluate in chain-root space (root at origin); rootLink.LinkMatrix(true) is the first child joint's transform in that space; the HI goal transform is relative to the chain root; DefaultZeroMap() maps an EE-axis direction to the zero-swivel reference vector and is [CORE] (goal-supplied). Both solvers are stateless per solve and re-derive the initial pose from initXYZ / initValue every evaluation.
Storage: three consecutive floats x, y, z; indexed access p[i] aliases (&x)[i] for read and write. Constructors: default performs no initialization; float constructor assigns directly; double and int constructors cast each component to float; array constructor takes elements 0–2.
Inline arithmetic, all component-wise in the written order, evaluated under §1.1 semantics:
a - b -> (a.x-b.x, a.y-b.y, a.z-b.z) /* + analogous */
a * b -> (a.x*b.x, a.y*b.y, a.z*b.z) /* element-wise; / analogous */
a * f, f * a -> (a.x*f, a.y*f, a.z*f) /* a/f divides; a+f adds f to each */
-=, +=, *=(float), /=(float), *=(Point3) /* in-place, component order x,y,z */
a % b -> (a.x*b.x + a.y*b.y + a.z*b.z) /* left-associative */
DotProd(a,b) -> (a.x*b.x + a.y*b.y + a.z*b.z) /* identical expression */
operator== / != -> exact per-component comparesLengths — the three differ, and the solvers' use of FLength specifically is normative:
Length(): (float)sqrt(x*x + y*y + z*z); /* CRT double sqrt; sum accumulated
at register precision per §1.1 */
FLength(): Sqrt(x*x + y*y + z*z); /* Sqrt per §2.9: float parameter,
x87 fsqrt */
LengthSquared(): (x*x + y*y + z*z);Because Sqrt takes a float parameter (§2.9), the register-precision sum is rounded to 32-bit float at the call boundary before the root is taken — a real numeric difference from Length(), which feeds the unrounded register-precision sum to the CRT double sqrt and truncates only on return. Both behaviors are normative.
[CORE] (declared with Point3 but implemented out-of-line): Normalize() (documented: each component divided by Length(); "more accurate than FNormalize"), FNormalize(), Unify() (in-place normalize), LengthUnify() (in-place normalize, returns the pre-normalization length), free Normalize(const Point3&), CrossProd(a,b) / operator^. Characterization requirements: which length function each normalizer uses; true division versus reciprocal-multiply; behavior at zero length (the solvers' == 0.0f tests imply the returned length is the raw computed value and no guard rewrites the vector — confirm). CrossProd candidate is the standard formula (a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x); since float products are exact in double, the result is association-insensitive up to the final double-rounding on store, so a small vector budget suffices.
Storage: float m[4][3] plus a flags DWORD with identity bits for position (row 3), rotation, and scale (POS_IDENT, ROT_IDENT, SCL_IDENT). Const row access returns row i as a Point3; GetRow(i) copies it; GetTrans() returns row 3. The default constructor leaves storage uninitialized and sets flags to 0 (i.e. not flagged identity); Matrix3(TRUE) builds the identity via a [CORE] routine (which presumably sets all identity flags). The four-row constructor and Set assign rows 0–3 via SetRow and then call ValidateFlags, which recomputes the identity bits. SetTrans clears POS_IDENT.
Flags are bit-exactness state. [CORE] routines (operator*, Invert, transforms) may branch on identity bits to skip work. Multiplication against a block that is exactly identity yields identical numbers on either path, and ValidateFlags presumably tests with exact compares — but this must be characterized, not assumed, because a fast path that skips arithmetic can round differently from one that performs it against nearly-identity values if the flag computation is tolerance-based. Implementers must reproduce flag propagation (constructors, SetRow, SetTrans, NoTrans, etc.) exactly as specified here; the characterization plan includes a ValidateFlags item (discriminator: matrices exactly identity vs identity ± 1 ulp in each block).
[CORE]: operator*, Invert (in-place and free Inverse), PreRotateX/Y/Z, RotateX/Y/Z, SetRow, NoTrans (zeros row 3), IdentityMatrix, ValidateFlags, RotAngleAxisMatrix, VectorTransform, PointTransform / operator*(Point3, Matrix3) (full affine).
Storage: floats x, y, z, w; vector part x/y/z, scalar w; q[i] maps 0→x … 3→w. Value constructors assign directly (double variant casts). The default constructor initializes to the identity quaternion (0,0,0,1) — this is normative.
[CORE]: Conjugate(), operator*(Quat) (both per §2.3), Quat(const Matrix3&), Quat(const AngAxis&), Quat(const Point3&, float) (documented to normalize), MakeMatrix(Matrix3&, bool transpose=false), Inverse. Live-path usage is limited to Conjugate, operator*, and Quat(AngAxis) (spline twist), plus whatever PreApplyTo uses internally. Quat(AngAxis) candidate: vector = axis·sin(a/2), w = cos(a/2), without the §2.2 negation (the solver-local negation exists precisely to compensate); sign convention is left-hand-rule per SDK documentation. Characterize jointly with the product and PreApplyTo as one observable unit.
The SDK float-math header supplies inline Sqrt, Sin, Cos, and SinCos. On the reference platform (x86, unmanaged) all four compile to x87 instruction sequences, not CRT calls; the CRT fallbacks that exist for other platforms are dead on the reference build and must not be used by a bit-exact replica.
Sqrt — the only §2.9 primitive on the live solver paths (via FLength, §2.6):
inline float Sqrt(float a) /* x86 reference build */
{
float r;
/* x87: load a (32-bit) onto the FP stack; execute fsqrt;
store ST(0) to r as 32-bit float. */
__asm { fld dword ptr a fsqrt fstp dword ptr [r] }
return r;
}Normative consequences:
- The parameter type is float: any register-precision sum passed in (e.g.
x*x+y*y+z*zfromFLength) is rounded to 32-bit float at the call boundary before the root is taken. -
fsqrtexecutes under the ambient control word 0x27F, i.e. rounds to 53-bit precision; the store then rounds to 24-bit. For a 24-bit input, a square root rounded first to 53 bits and then to 24 is provably identical to the directly rounded 24-bit result (the intermediate precision exceeds the 2p+2 threshold), so a native shim may implementSqrt(a)as(float)sqrt((double)a)provided the platform's doublesqrtis correctly rounded (IEEE requires this) — or, more simply, emitfsqrtvia inline asm on the 387 path. Either way, verify against the vector file; do not reason further from this note.
SinCos(angle, &s, &c) computes both via a single x87 fsincos on the 32-bit angle (cosine is at the top of the stack and is stored first, then sine); Sin/Cos call SinCos and discard one result. These are not on the live solver paths — the quaternion construction of §2.2 uses the CRT double sin/cos — but they are recorded here because Biped or node-TM code paths characterized later may turn out to use them, and fsincos differs numerically from the CRT polynomial implementations. If Stage 5/6 residual triage (§6) suspects a trig discrepancy, test the fsincos hypothesis before amending logic.
Accepted shapes: LinkCount() of 3 or 0; otherwise return bInvalidInitialValue. Missing HI goal interface: bInvalidArgument. Success: 0. Properties: swivel yes, EE rotation no, sliding joint no, analytic, single-chain, history-independent.
Set rootLink.rotXYZ = rootLink.initXYZ and each link's dofValue = initValue. Then, in chain space:
elbowTM = LinkOf(2).DofMatrix() * LinkOf(1).LinkMatrix(true)
* LinkOf(0).LinkMatrix(true) * rootLink.LinkMatrix(true);
elbowPos = elbowTM.GetTrans();
elbowUnit = Normalize(elbowPos);
M = rootLink.LinkMatrix(true);
LinkOf(0).ApplyLinkMatrix(M, true);
LinkOf(1).ApplyLinkMatrix(M, true);
LinkOf(2).ApplyLinkMatrix(M, true);
wristPos = M.GetTrans();
eeUnit = Normalize(wristPos); /* shoulder is at origin */
L1 = FLength(elbowPos);
L2 = FLength(wristPos - elbowPos);
Lsum = L1 + L2;
Ldiff = L1 - L2; /* signed; fabs applied at use */
d = FLength(wristPos);
Law-of-cosines angles, with exact float expression ordering (left-associative as written):
float elbowNum = L1*L1 + L2*L2 - d*d;
float cosElbow = elbowNum / (2 * L1 * L2); /* int 2 promoted; ((2*L1)*L2) */
float theta0 = acos_safe(cosElbow); /* double -> float at assignment */
float shoulderNum = d*d + L1*L1 - L2*L2;
float cosShoulder = shoulderNum / (2 * d * L1);
float phi0 = acos_safe(cosShoulder);Initial shoulder frame: row0 = eeUnit, row1 = DefaultZeroMap()(eeUnit) used as-is (deliberately not re-orthogonalized), row2 = Normalize(CrossProd(row0, row1)), row3 = (0,0,0).
Shoulder-local quantities: take M2 = rootLink.LinkMatrix(false); its translation is the local elbow position, normalized to elbowUnitLocal; apply the three link matrices to M2 as above to get the local wrist position, normalized to eeUnitLocal; the initial chain-plane normal is
initPlaneNormal = Normalize(CrossProd(eeUnitLocal, elbowUnitLocal));
g = Goal().GetTrans();
goalDist = FLength(g);
goalUnit = Normalize(g);
if (goalDist > Lsum) wrist = goalUnit * Lsum;
else if (goalDist < (float)fabs(Ldiff)) wrist = goalUnit * (float)fabs(Ldiff);
else wrist = g;(fabs is the double CRT function on the promoted float, cast back to float.) Recompute ee = Normalize(wrist), d = FLength(wrist), then elbowNum/cosElbow/theta and shoulderNum/cosShoulder/phi with the identical expressions of §3.1.
Reference vector p = DefaultZeroMap()(ee). Target normal t:
VH-target path (UseVHTarget()): t = Normalize(VHTarget()); then t -= (ee % t) * ee; if t.LengthUnify() == 0.0f then t = p, else t = ee ^ t.
Swivel in start joint (SwivelAngleParent() == kSAInStartJoint): t = ApplyRowQuat(DefaultZeroMap()(ee), MakeRowQuat(ee, SwivelAngle())).
Swivel in goal (kSAInGoal): set t = (0,0,0); eeInGoal = Inverse(Goal()).VectorTransform(ee); if eeInGoal.LengthUnify() != 0.0f: tInGoal = ApplyRowQuat(DefaultZeroMap()(eeInGoal), MakeRowQuat(eeInGoal, SwivelAngle())); t = Goal().VectorTransform(tInGoal); t -= (ee % t) * ee. Then, unconditionally on this path: if t.LengthUnify() == 0.0f, t = p.
Swivel scalar:
float swivel = p % t;
swivel = acos_safe(swivel); /* double -> float */
Point3 cross = p ^ t;
if (cross % ee < 0.0f)
swivel = 2.0 * PI - swivel; /* double expr, truncated */Frame update:
M1 = RotAngleAxisMatrix(p, phi - phi0);
M2 = RotAngleAxisMatrix(ee, swivel);
rx = VectorTransform(M2, VectorTransform(M1, ee));
ry = VectorTransform(M2, p);
rz = Normalize(CrossProd(rx, ry));Frame F has rows rx, ry, rz, zero translation, built over an identity matrix (note the flag consequences of §2.7 when rows are then overwritten). Shoulder rotation:
R = Inverse(initialShoulderFrame) * F;
R.PreRotateZ(rootLink.initXYZ.z);
R.PreRotateY(rootLink.initXYZ.y);
R.PreRotateX(rootLink.initXYZ.x);
MatrixToEuler(R, Eul, EULERTYPE_XYZ);
rootLink.rotXYZ = (Eul[0], Eul[1], Eul[2]);Behavioral note: this pre-rotation compensation is known-inexact for nonzero initial angles; the reference ships it anyway and the replica reproduces the same inexact result.
Elbow update:
E = RotAngleAxisMatrix(initPlaneNormal, theta - theta0);
E.PreRotateZ(LinkOf(0).initValue);
E.PreRotateY(LinkOf(1).initValue);
E.PreRotateX(LinkOf(2).initValue);
MatrixToEuler(E, elbowEul, EULERTYPE_XYZ);
LinkOf(0).dofValue = elbowEul[2];
LinkOf(1).dofValue = elbowEul[1];
LinkOf(2).dofValue = elbowEul[0];
return 0;The crosswise assignment (link 0 carries the Z component) is normative.
Initialization: initial EE unit initEEUnit = Normalize(rootLink.LinkMatrix(true).GetTrans()); initial frame built exactly as §3.1 (row1 from the zero map, not orthogonalized). Solve:
ee = Normalize(Goal().GetTrans());
ang = acos_safe(DotProd(initEEUnit, ee)); /* double -> float */
axis = Normalize(CrossProd(initEEUnit, ee)); /* degenerate when (anti)parallel —
inherited, not fixed */
M = RotAngleAxisMatrix(axis, ang);
rx = VectorTransform(M, initEEUnit);
ry = VectorTransform(M, initialFrame.row1);
rz = Normalize(CrossProd(rx, ry));Shoulder extraction, pre-rotations, and Euler write-back identical to §3.2. Return 0.
Solve queries the version-2 spline goal (IID_SPLINE_IKGOAL2); a null goal, goal node, or chain control returns bInvalidArgument. The initializer separately queries the version-1 interface and silently returns if absent. Properties: no swivel, no sliding joint, analytic, single chain, history-independent. Success returns 0, including the early-exit when no spline points were recorded — that early-exit is solver behavior (the chain is left untouched for that frame) and stays in the spec.
Clear the per-solve tables (bone lengths, axes, co-axes). With the joint iterator and InitJointAngles():
First pass: Begin(false), record the first pivot; iterate to the end; record the last pivot. chainAxis = last − first — never normalized; feed the raw vector to the zero map.
twistUp = DefaultZeroMap()(chainAxis);
chainUp = twistUp; /* copy of raw zero-map output */
twistUp = TwistParent().VectorTransform(twistUp);
twistUp = VectorTransform(twistUp, Inverse(linkChain.parentMatrix));
chainUp = VectorTransform(chainUp, Inverse(linkChain.parentMatrix));
twistUp.Unify();
chainUp.Unify();(twistUp is used in the §4.4 frame construction; chainUp feeds the per-bone pass below.)
Second pass, per joint:
mat = iter.ProximalFrame();
axisOrder = iter.GetJointAxes(); angles = iter.GetJointAngles();
for (i = 2; i >= 0; --i) {
if (axisOrder[i] == '_') break;
if (axisOrder[i] == 'x') mat.PreRotateX(angles[0]);
if (axisOrder[i] == 'y') mat.PreRotateY(angles[1]);
if (axisOrder[i] == 'z') mat.PreRotateZ(angles[2]);
}
mat.Invert();
boneAxis = iter.DistalEnd() * mat; boneAxis.Unify(); /* bone axis, local frame */Dominant index ei of boneAxis: start at 0 with maxAbs = fabs(boneAxis.x); if maxAbs < fabs(boneAxis.y) take 1; then if maxAbs < fabs(boneAxis.z) take 2 (strict <, so ties keep the earlier axis; fabs is the double CRT on promoted floats, compared as float — reproduce exactly).
boneUp = VectorTransform(chainUp, mat); /* mat is the INVERTED matrix */
coAxis = CrossProd(boneUp, boneAxis);
if (coAxis.LengthUnify() < 100.0f * FLT_EPSILON) {
coAxis = (0,0,0); coAxis[(ei + 2) % 3] = 1.0f;
coAxis = CrossProd(coAxis, boneAxis); coAxis.Unify();
}
coAxis = CrossProd(boneAxis, coAxis); /* final co-axis */Append boneAxis, coAxis; bone length FLength(DistalEnd() − Pivot()) appends and accumulates into the chain length.
Increments are exactly the float literals 0.001f, 0.0001f, 0.000001f (a fourth, 0.005f, exists off the live path). First-bone thresholds: eps1 = L0/7.0f, eps2 = L0/50.0f, eps3 = L0/100.0f. Cache prime: SplinePosAt(0.0f, 1, TRUE); thereafter every sample is SplinePosAt(u, 1) * toJointMat with toJointMat = goalNode->GetObjectTM(SolveTime()) * Inverse(parentMatrix). Spline evaluation itself is [CORE] — characterize the cached Bezier-spline position function against reference data before solver-level testing.
March from u = 0 while u <= 1.0f: sample, dist = FLength(sample − firstJointPos). Minimum tracking before advancing: if dist < distBest, record uBest = u, distBest, the sample point, uBeforeBest = uLast, and set a flag; then uLast = u. Advance: dist > eps1 → u += 0.001f; else dist > eps2 → u += 0.0001f; else dist > eps3 → u += 0.000001f; else append (u, point), set the running start point, and break. After advancing, if the flag is set, record uAfterBest = u (the u one step past the minimum) and clear the flag. u accumulates by repeated float addition in exactly this order — no closed-form reconstruction.
If the loop exits with u > 1.0f (no convergence): rescan for (u = uBeforeBest + 0.000001f; u < uBestInitial; u += 0.000001f) where uBestInitial is the pre-refinement uBest, updating uBest/distBest/best point only. If afterwards uBest >= uBestInitial, continue the same loop from the current u up to uAfterBest. Adopt u = uBest and the best point as the start; append. An offset from the true first-joint position is computed here and never used — dead code, do not apply it.
Open spline — while u <= 1.0f && j < boneCount, with dist = FLength(sample(u) − startPoint):
if (L[j] - dist > eps1) { uPrev = u; u += 0.001f; }
else if (L[j] - dist > eps2) { uPrev = u; u += 0.0001f; }
else if (L[j] - dist > eps3) { uPrev = u; u += 0.000001f; }
else if (dist - L[j] > eps3) { /* backtrack */
stepDelta = u - uPrev;
if (stepDelta == 0.001f) u = uPrev + 0.0001f; /* exact float equality */
else if (stepDelta == 0.0001f) u = uPrev + 0.000001f;
else { /* accept current point as the joint */
append (u, sample); startPoint = sample; uPrev = u; ++j;
if (j < boneCount) { eps1 = L[j]/5.0f; eps2 = L[j]/50.0f; eps3 = L[j]/100.0f; }
}
}
else { /* converged: accept */
append (u, sample); startPoint = sample; uPrev = u; ++j;
if (j < boneCount) { eps1 = L[j]/5.0f; eps2 = L[j]/50.0f; eps3 = L[j]/100.0f; }
}
Note the threshold divisor changes from 7 to 5 after the first bone, u is not advanced on acceptance (the next iteration measures zero distance and advances), and the exact-float-equality backtrack tests demand the same literals and the same repeated-addition accumulation of u. Under x87, force u to spill to a 32-bit variable each iteration (volatile float in the native build, or rely on the Wine/VS2008 kernel per §1.2) so the equality tests fire on the same branch as the reference.
Closed spline — identical structure with: termination du <= 1.0f && j < boneCount where du accumulates each increment and resets to 0.0f on acceptance; wrap after each advance if (u > 1.0) u = u - 1.00f (double comparison via promotion, float subtraction); backtrack delta stepDelta = (u > uPrev) ? u - uPrev : 1.0f + u - uPrev; backtrack retries adjust du += (0.0001f - 0.001f) and du += (0.000001f - 0.0001f) respectively and wrap with if (u > 1.0f) u -= 1.0f; the closed variant's acceptance and give-up branches do not update uPrev (the open variant's do) — reproduce this asymmetry.
If no points were recorded at all, return 0 without touching the chain.
For each bone index b in 0..boneCount−1:
twist = (b == 0) ? startTwist
: startTwist + (endTwist / (boneCount - 1)) * b;(For boneCount == 1 only b == 0 occurs, so the divisor (boneCount−1) is never evaluated; there is no division by zero on this path.)
Direction: only assigned when b <= uCount:
dir = (b == uCount - 1) ? Normalize(sample(1.00f) - point[b])
: (b == uCount) ? Normalize(sample(1.0f) - sample(0.99f))
: Normalize(point[b+1] - point[b]);For b > uCount the previous iteration's value persists unreassigned — replicate the carry-over.
Dominant index i of boneAxes[b] exactly as §4.1; j = (i+1)%3, k = (i+2)%3. Bone basis: an (uninitialized, flags = 0 per §2.7) matrix whose rows i, k, j are set to boneAxes[b], boneCoAxes[b], and CrossProd(boneCoAxes[b], boneAxes[b]) respectively — all three basis rows are covered — translation cleared, then inverted. Global frame over identity: row i = dir; norm = CrossProd(twistUp, dir); if norm.LengthUnify() < 100.0f*FLT_EPSILON, substitute norm = Unify(CrossProd(unitAxis(k), dir)); row j = norm; row k = CrossProd(dir, norm). Then global = boneBasis⁻¹ * global.
Twist application: q = Quat(AngAxis(boneAxes[b], twist)); the rotation value wraps q.Conjugate(). For b == 0 it pre-applies to global, and local = global. For b > 0: local = global * Inverse(frame[b-1]); pre-apply to local; global = local * frame[b-1]. Append global to the frame list; MatrixToEuler(local, Eul, EULERTYPE_XYZ); append the Euler triple. Finally iterate the joints from the chain start assigning the triples in order via SetJointAngles, and clear the per-solve tables.
The solvers are the small half of the problem. Everything they consume is built by Max's IK system, which is [CORE] end to end: construction of the LinkChain from nodes (chain-root space definition, LinkMatrix/DofMatrix/ApplyLinkMatrix semantics from bone offsets and rotation-controller axis ordering, parentMatrix), the HI goal transform and its animation sources (goal node TM, swivel angle track, VH target), the joint iterator (ProximalFrame, DistalEnd, GetJointAxes string, SetJointAngles write-back into node controllers), DefaultZeroMap, TwistParent, SolveTime, and the spline-goal position cache. Node-TM assembly must also track Matrix3 identity flags (§2.7), since they flow into [CORE] arithmetic.
The characterization corpus for all of these is the set of reference exports plus purpose-built test scenes evaluated by the original pipeline where any still exist in the asset history; each item gets its own vector file per Appendix A. The original NeL exporter sources in the Ryzom Core release remain the authoritative statement of what node-level data the pipeline sampled and are freely usable. That review must also settle the §7.6 question (whether the pipeline invoked the Biped export interface's scale removal or figure-mode calls before sampling), since the answer changes every sampled TM in the corpus.
Per-function, bottom-up, all comparisons on raw 32-bit patterns: CRT transcendental shims first, then §2 primitives (including a Sqrt vector file confirming the §2.9 shim equivalence on the target platform), then spline position sampling, then whole-solver joint Eulers per frame, then full-pipeline node transforms against reference exports. Every [CORE] item has a checked-in vector file with provenance (source export, frame, node).
Vector budget notes: MatrixToEuler near gimbal gets the largest budget. The Normalize/Unify/LengthUnify family needs vectors covering zero-length behavior and discriminating which length routine each uses; the §2.9 semantics sharpen this — FNormalize presumably rides FLength, whose call-boundary rounding is specified, so vectors whose squared sum rounds differently at float versus double discriminate cheaply. CrossProd needs only a small budget (association-insensitive per §2.6). Quat::Conjugate needs a single confirmation vector. The Quat product convention, Quat(AngAxis), and RotationValue::PreApplyTo are characterized as one jointly-observed unit with two-rotation discriminator vectors. Matrix3::ValidateFlags uses exact-identity versus identity ± 1 ulp discriminators. Sqrt needs a modest sweep (denormals, exact squares, values straddling representable midpoints) run identically on the Wine/VS2008 kernel and any native shim.
The Wine/VS2008 kernel build (§1.2) is the reference implementation; the native Linux build graduates only on corpus-wide bit-identity with it. Measure NeL track quantization early anyway — not to relax the target, but because knowing which ULP ranges are observable in exports tells you which characterization vectors can actually discriminate between candidate [CORE] implementations.