Moves lines, circles, arcs, points and polyline vertices so that the geometric and dimensional constraints placed on them are met.
How the precision is kept
The solver is an iteration in 64-bit numbers that stops when every equation is met within 1e-9 times the size of the geometry. That stopping tolerance, not the arithmetic, sets how exactly a constrained dimension is held. Derivatives are central differences, not exact derivatives.
Open item. Stop on a tolerance in each equation's own unit, scaled to the size of the constrained objects and not to their distance from the origin.
Calculations under this policy
9 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
The tolerance scales with the constrained objects' bounding box and with how far that box is from the origin. It is at least 1e-9. Circles and arcs count by center ± radius.
Each constraint is one or two functions of the current coordinates: coordinate differences, signed distance from a line (cross product ÷ length), differences of lengths or radii, or an angle from atan2 wrapped to (−π, π]. Parallel and perpendicular are the sine and cosine of the angle between the lines, so their residual has no length unit. Divisions by a length use 1e-300 in place of a zero length.
gsolve.js
lines 174 to 268
CadCraft 0.8.2
function equationsFor(c, G, D, ents, out) {
const r = c.refs;
if (!r.every((q) => ents.has(q.e) && refKind(ents.get(q.e), q.k))) return;
const deps = [...new Set(r.flatMap((q) => D.get(q.e)))];
const push = (...fs) => { for (const f of fs) out.push({ cid: c.id, f, deps }); };
const P = (i) => G.get(r[i].e).point(r[i].k);
const L = (i) => G.get(r[i].e).line(r[i].k);
const C = (i) => G.get(r[i].e).circle(r[i].k);
const K = (i) => refKind(ents.get(r[i].e), r[i].k);
const dir = (i) => { const [a, b] = L(i); return sub(b, a); };
const s = c.sign || 1;
switch (c.type) {
case 'coincident': {
const k2 = K(1);
if (k2 === 'P') push(() => P(0)[0] - P(1)[0], () => P(0)[1] - P(1)[1]);
else if (k2 === 'L') push(() => sd(P(0), L(1)));
else push(() => { const k = C(1); return len(sub(P(0), k.c)) - k.r; });
break;
}
case 'concentric': push(() => C(0).c[0] - C(1).c[0], () => C(0).c[1] - C(1).c[1]); break;
case 'collinear': push(() => sd(L(1)[0], L(0)), () => sd(L(1)[1], L(0))); break;
case 'fix': {
if (K(0) === 'P') { const v = c.value || [0, 0]; push(() => P(0)[0] - v[0], () => P(0)[1] - v[1]); break; }
const e = ents.get(r[0].e), fx = c.fixed || {};
const id = r[0].e;
for (const f of fieldsOf(e)) {
if (typeof fx[f] !== 'number') continue;
// whole entity: each field held at its value (an arc's angles too)
push(() => fieldValue(G.get(id), e, f) - fx[f]);
}
break;
}
case 'parallel': push(() => { const u = dir(0), v = dir(1); return cross(u, v) / ((len(u) * len(v)) || 1e-300); }); break;
case 'perpendicular': push(() => { const u = dir(0), v = dir(1); return dot(u, v) / ((len(u) * len(v)) || 1e-300); }); break;
case 'horizontal': case 'vertical': {
const k = c.type === 'horizontal' ? 1 : 0;
if (K(0) === 'L') push(() => { const [a, b] = L(0); return b[k] - a[k]; });
else push(() => P(1)[k] - P(0)[k]);
break;
}
case 'tangent': {
if (K(0) === 'L') push(() => s * sd(C(1).c, L(0)) - C(1).r);
else push(() => { const a = C(0), b = C(1), d = len(sub(b.c, a.c)); return c.kind === 'int' ? d - s * (a.r - b.r) : d - (a.r + b.r); });
break;
}
case 'equal':
if (K(0) === 'L') push(() => len(dir(0)) - len(dir(1)));
else push(() => C(0).r - C(1).r);
break;
case 'symmetric': {
// two points (or the matching ends of two lines, or two centres) mirror each other across line refs[2]
const pairs = [];
if (K(0) === 'P') pairs.push([() => P(0), () => P(1)]);
else if (K(0) === 'L') {
const sw = c.pair === 1;
pairs.push([() => L(0)[0], () => L(1)[sw ? 1 : 0]], [() => L(0)[1], () => L(1)[sw ? 0 : 1]]);
} else pairs.push([() => C(0).c, () => C(1).c]);
for (const [A, B] of pairs) {
push(() => { const a = A(), b = B(); return sd([(a[0] + b[0]) / 2, (a[1] + b[1]) / 2], L(2)); },
() => { const a = A(), b = B(), u = dir(2); return dot(sub(b, a), u) / (len(u) || 1e-300); });
}
if (K(0) === 'R') push(() => C(0).r - C(1).r);
break;
}
// ---- dimensional (value = the constraint's last resolved number)
case 'hdist': case 'vdist': {
const k = c.type === 'hdist' ? 0 : 1;
push(() => s * (P(1)[k] - P(0)[k]) - c.value);
break;
}
case 'dist': {
const k2 = K(1);
if (k2 === 'P') push(() => len(sub(P(1), P(0))) - c.value);
else if (K(0) === 'P') push(() => s * sd(P(0), L(1)) - c.value);
else push(() => s * sd(L(1)[0], L(0)) - c.value);
break;
}
case 'radius': push(() => C(0).r - c.value); break;
case 'diameter': push(() => 2 * C(0).r - c.value); break;
case 'angle': {
if (K(0) === 'R') {
// an arc's sweep (degrees)
push(() => { const gv = G.get(r[0].e), e = ents.get(r[0].e); return (wrap(fieldValue(gv, e, 'a1') - fieldValue(gv, e, 'a0') - c.value * Math.PI / 180)); });
break;
}
const [s1, s2] = c.dirs || [1, 1];
push(() => {
const u = dir(0).map((v) => v * s1), v = dir(1).map((w) => w * s2);
return wrap(Math.atan2(cross(u, v), dot(u, v)) - s * c.value * Math.PI / 180);
});
break;
}
default: break;
}
}
Circle of a polyline arc segment
64-bit
The center and radius of a polyline's arc segment as the solver sees it.
θ = 4·atan(b); r = |c ÷ (2·sin(θ ÷ 2))|; h = ± r·cos(θ ÷ 2); center = chord midpoint + normal × h
Recomputed from the two vertices at every evaluation; the bulge is not a variable, so the arc follows its ends. A zero chord is replaced by 1e-300 before dividing.
Central differences, not exact derivatives. f(x + h) − f(x − h) subtracts nearly equal numbers, so each entry carries rounding of the order of 2⁻⁵³ ÷ 10⁻⁶ ≈ 10⁻¹⁰ relative. This affects the steps taken and the rank analysis, not the residual the solve stops on. Entries that are exactly 0 or not finite are left out.
gsolve.js
lines 296 to 310
CadCraft 0.8.2
function jacobian(eqs, x) {
return eqs.map((q) => {
const row = new Map();
for (const k of q.deps) {
if (row.has(k)) continue;
const x0 = x[k], h = 1e-6 * Math.max(1, Math.abs(x0));
x[k] = x0 + h; const fp = q.f();
x[k] = x0 - h; const fm = q.f();
x[k] = x0;
const d = (fp - fm) / (2 * h);
if (d !== 0 && Number.isFinite(d)) row.set(k, d);
}
return row;
});
}
Cholesky solve
64-bit
Solves the damped normal equations for one step.
A = L·Lᵀ; L·y = b; Lᵀ·z = y
A dense factorization on a copy of the matrix, then forward and back substitution, all as running subtractions in 64-bit numbers. A pivot that is not greater than 0 returns null, and the caller raises the damping a hundredfold and tries again. There is no pivoting and no iterative refinement.
gsolve.js
lines 317 to 337
CadCraft 0.8.2
function cholSolve(A, b) {
const n = b.length;
const L = A.map((r) => r.slice());
for (let j = 0; j < n; j++) {
let s = L[j][j];
for (let k = 0; k < j; k++) s -= L[j][k] * L[j][k];
if (!(s > 0)) return null;
const d = Math.sqrt(s);
L[j][j] = d;
for (let i = j + 1; i < n; i++) {
let t = L[i][j];
for (let k = 0; k < j; k++) t -= L[i][k] * L[j][k];
L[i][j] = t / d;
}
}
const y = new Array(n);
for (let i = 0; i < n; i++) { let s = b[i]; for (let k = 0; k < i; k++) s -= L[i][k] * y[k]; y[i] = s / L[i][i]; }
const z = new Array(n);
for (let i = n - 1; i >= 0; i--) { let s = y[i]; for (let k = i + 1; k < n; k++) s -= L[k][i] * z[k]; z[i] = s / L[i][i]; }
return z;
}
Newton iteration
64-bit
Repeats damped least-norm steps until every equation is met.
(J·Jᵀ + μ·max(scale, 1e-12)·I)·y = −r; dx = Jᵀ·y; x ← x + α·dx, α = 1, ½, ¼ …
Stops when every |residual| ≤ tol, after 100 iterations, or when 12 tries with growing damping bring no decrease. μ starts at 1e-12, grows × 100 after a failed step (up to 1e6) and shrinks ÷ 10 after a good one; the line search halves the step up to 20 times and accepts any decrease of the sum of squares. Residuals of up to tol remain in the result.
gsolve.js
lines 339 to 378
CadCraft 0.8.2
function newton(eqs, x, tol, maxIter = 100) {
let r = residuals(eqs);
let f = norm2(r);
let mu = 1e-12;
let it = 0;
const met = () => r.every((v) => Math.abs(v) <= tol);
for (; it < maxIter; it++) {
if (met()) return { ok: true, it };
if (!x.length) break;
const J = jacobian(eqs, x);
const m = J.length;
const A = [];
let scale = 0;
for (let i = 0; i < m; i++) {
A.push(new Array(m));
for (let j = 0; j <= i; j++) { const v = sparseDot(J[i], J[j]); A[i][j] = v; if (j < i) A[j][i] = v; }
scale = Math.max(scale, A[i][i]);
}
let improved = false;
for (let tries = 0; tries < 12 && !improved; tries++) {
const B = A.map((row, i) => row.map((v, j) => (i === j ? v + mu * Math.max(scale, 1e-12) + 1e-300 : v)));
const y = cholSolve(B, r.map((v) => -v));
if (!y) { mu *= 100; continue; }
const dx = new Array(x.length).fill(0);
J.forEach((row, i) => { for (const [k, v] of row) dx[k] += v * y[i]; });
const x0 = x.slice();
let alpha = 1;
for (let ls = 0; ls < 20; ls++) {
for (let k = 0; k < x.length; k++) x[k] = x0[k] + alpha * dx[k];
const r2 = residuals(eqs), f2 = norm2(r2);
if (f2 < f) { r = r2; f = f2; improved = true; break; }
alpha /= 2;
}
if (!improved) { for (let k = 0; k < x.length; k++) x[k] = x0[k]; mu *= 100; if (mu > 1e6) break; }
else mu = Math.max(1e-12, mu / 10);
}
if (!improved) return { ok: met(), it };
}
return { ok: met(), it };
}
Writing the solution back
64-bit
The changed objects after a solve, and the remaining residual.
residual = √(Σ rᵢ²)
Solved values are written back only when the solve succeeded, and only for objects with a field that moved by more than 1e-12 × max(1, |value|). The values are the solver's own 64-bit results, not rounded.
gsolve.js
lines 394 to 415
CadCraft 0.8.2
/** Solve; on success the changed entities come back as new objects (nothing is mutated). */
export function solveSystem(constraints, get, { pinned = null, maxIter = 100 } = {}) {
const sys = buildSystem(constraints, get, pinned);
const tol = tolerance(sys.ents);
const res = newton(sys.eqs, sys.x, tol, maxIter);
const updates = new Map();
if (res.ok) {
for (const [id, e] of sys.ents) {
const vals = new Map();
let moved = false;
for (const f of fieldsOf(e)) {
const k = sys.vars.get(id + '|' + f);
if (k === undefined) continue;
vals.set(f, sys.x[k]);
if (Math.abs(sys.x[k] - getField(e, f)) > 1e-12 * Math.max(1, Math.abs(getField(e, f)))) moved = true;
}
if (moved) updates.set(id, withFields(e, vals));
}
}
const resid = residuals(sys.eqs);
return { ok: res.ok, updates, iterations: res.it, residual: Math.sqrt(norm2(resid)), worst: worstOf(sys.eqs, resid, tol) };
}
Radius sign and arc angles
64-bit
Puts solved values back into an object in its stored form.
A negative solved radius is made positive, and an arc then gets π added to both angles. Arc angles are reduced to [0, 2π) with the % operator.
gsolve.js
lines 45 to 62
CadCraft 0.8.2
/** A copy of `e` with new field values (Map field → number). Radii stay positive; arc angles in [0, 2π). */
export function withFields(e, vals) {
if (e.type === 'poly') {
const pts = e.pts.map((p, i) => {
const x = vals.get(i + 'x'), y = vals.get(i + 'y');
return x === undefined && y === undefined ? p : { ...p, x: x ?? p.x, y: y ?? p.y };
});
return { ...e, pts };
}
const out = { ...e };
for (const [f, v] of vals) out[f] = v;
if ((e.type === 'circle' || e.type === 'arc') && out.r < 0) {
out.r = -out.r;
if (e.type === 'arc') { out.a0 += Math.PI; out.a1 += Math.PI; }
}
if (e.type === 'arc') { out.a0 = norm(out.a0); out.a1 = norm(out.a1); }
return out;
}
Degrees of freedom and redundant constraints
64-bit
How many independent constraints there are, which ones add nothing, which conflict and which objects are fully fixed.
v = row ÷ |row|; v ← v − (v·b)·b for each kept row b, twice; redundant when |v| < 1e-7; degrees of freedom = variables − rank
Rows are normalized and orthogonalized twice against the rows kept so far, in Float64Array vectors. A row shorter than 1e-12, or with less than 1e-7 left, counts as redundant. A constraint is a conflict when its residual exceeds 1000 × tol at the current geometry; an object is fully determined when each of its variables has a projection above 1 − 1e-6.
gsolve.js
lines 429 to 466
CadCraft 0.8.2
export function analyze(constraints, get) {
const sys = buildSystem(constraints, get);
const { eqs, x, vars, ents } = sys;
const tol = tolerance(ents);
const J = jacobian(eqs, x);
const n = x.length;
const Q = [];
const redundant = new Set(), conflicts = new Set();
eqs.forEach((q, i) => {
if (Math.abs(q.f()) > tol * 1000) conflicts.add(q.cid);
const v = new Float64Array(n);
let vn = 0;
for (const [k, d] of J[i]) { v[k] = d; vn += d * d; }
vn = Math.sqrt(vn);
if (vn < 1e-12) { redundant.add(q.cid); return; }
for (let k = 0; k < n; k++) v[k] /= vn;
for (let pass = 0; pass < 2; pass++) {
for (const b of Q) {
let d = 0;
for (let k = 0; k < n; k++) d += b[k] * v[k];
if (d !== 0) for (let k = 0; k < n; k++) v[k] -= d * b[k];
}
}
let wn = 0;
for (let k = 0; k < n; k++) wn += v[k] * v[k];
wn = Math.sqrt(wn);
if (wn < 1e-7) { redundant.add(q.cid); return; }
for (let k = 0; k < n; k++) v[k] /= wn;
Q.push(v);
});
const fixedVar = new Array(n).fill(0);
for (const b of Q) for (let k = 0; k < n; k++) fixedVar[k] += b[k] * b[k];
const determined = new Set();
for (const [id, e] of ents) {
if (fieldsOf(e).every((f) => { const k = vars.get(id + '|' + f); return k === undefined || fixedVar[k] > 1 - 1e-6; })) determined.add(id);
}
return { dof: n - Q.length, rank: Q.length, vars: n, redundant: [...redundant], conflicts: [...conflicts], determined };
}
In detail
Method
Damped Newton with the least-norm step dx = Jᵀ(J·Jᵀ + μI)⁻¹(−r), a halving line search and Levenberg–Marquardt growth of μ, as the file's header states. J·Jᵀ is solved by a dense Cholesky factorization. Rank and redundant constraints by Gram–Schmidt over the Jacobian rows.
Number types
64-bit numbers; Float64Array for the Gram–Schmidt vectors.
Tolerances
tol = 1e-9 × extent, with extent = max(1, width, height, largest |coordinate|) of the constrained objects. The solve stops when every |residual| ≤ tol. Conflicts and the already-met test use 1000 × tol; the list of worst constraints uses 100 × tol. Finite-difference step 1e-6 × max(1, |x|). A row is redundant when less than 1e-7 of it is left after Gram–Schmidt. A field counts as moved above 1e-12 × max(1, |value|).
Limits
One tolerance is applied to residuals in drawing units (distances), in radians (angles) and to the unitless sine or cosine of the parallel and perpendicular equations. The extent includes the distance from the origin, so objects at coordinates near 10⁶ are accepted with a distance error of up to 0.001 units and an angle error of up to 0.001 rad. At most 100 iterations.
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.