double xCoord = 0.;

double yCoord = 0.;

double zCoord = 0.;

{

return gUU * delta.uCoord * delta.uCoord + 2. * gUV * delta.uCoord * delta.vCoord +

gVV * delta.vCoord * delta.vCoord;

using Evaluate = void (*)(const void* context, const Vec2& uv, double& gUU, double& gUV, double& gVV);

Evaluate evaluate = nullptr;

const void* context = nullptr;

/// The 3D length squared spanned by the parametric displacement \a delta starting at \a uv.

double lengthSq(const Vec2& uv, const Vec2& delta) const

{

if (evaluate == nullptr) {

return delta.uCoord * delta.uCoord + delta.vCoord * delta.vCoord;

}

double gUU = 1.;

double gUV = 0.;

double gVV = 1.;

evaluate(context, uv, gUU, gUV, gVV);

return parametricLengthSq(gUU, gUV, gVV, delta);

}

/// The 3D distance squared between two nearby parametric points, with the form evaluated at \a from.

double distanceSq(const Vec2& from, const Vec2& to) const { return lengthSq(from, to - from); }

/// The largest 3D length a unit parametric displacement spans at \a uv: the square root of the larger eigenvalue.

double maxScale(const Vec2& uv) const

{

if (evaluate == nullptr) {

return 1.;

}

double gUU = 1.;

double gUV = 0.;

double gVV = 1.;

evaluate(context, uv, gUU, gUV, gVV);

const double trace = gUU + gVV;

const double determinant = gUU * gVV - gUV * gUV;

// the eigenvalues of a symmetric 2x2 form, guarded against a slightly negative discriminant

const double discriminant = std::max(0., trace * trace - 4. * determinant);

return std::sqrt(std::max(0., 0.5 * (trace + std::sqrt(discriminant))));

}

{

return {[](const void* context, const Vec2& uv, double& gUU, double& gUV, double& gVV) {

static_cast<const Surface*>(context)->parametricMetric(uv, gUU, gUV, gVV);

},

&surface};

{

return firstVector.xCoord * secondVector.xCoord + firstVector.yCoord * secondVector.yCoord +

firstVector.zCoord * secondVector.zCoord;

{

return {firstVector.yCoord * secondVector.zCoord - firstVector.zCoord * secondVector.yCoord,

firstVector.zCoord * secondVector.xCoord - firstVector.xCoord * secondVector.zCoord,

firstVector.xCoord * secondVector.yCoord - firstVector.yCoord * secondVector.xCoord};

{

const double vectorNorm = norm(vector);

if (vectorNorm <= kTolerance) {

return {};

}

return vector * (1. / vectorNorm);

{

if (dimension == 0) {

return vector.xCoord;

}

if (dimension == 1) {

return vector.yCoord;

}

return vector.zCoord;

{

if (dimension == 0) {

vector.xCoord = value;

} else if (dimension == 1) {

vector.yCoord = value;

} else {

vector.zCoord = value;

}

{

const double deltaU = firstPoint.uCoord - secondPoint.uCoord;

const double deltaV = firstPoint.vCoord - secondPoint.vCoord;

return deltaU * deltaU + deltaV * deltaV;

{

const Vec2 segmentVector = segmentEnd - segmentStart;

const double segmentLengthSq = segmentVector.uCoord * segmentVector.uCoord + segmentVector.vCoord * segmentVector.vCoord;

if (segmentLengthSq <= kToleranceSq) {

return distanceSq(point, segmentStart);

}

const double pointProjection = ((point.uCoord - segmentStart.uCoord) * segmentVector.uCoord +

(point.vCoord - segmentStart.vCoord) * segmentVector.vCoord) /

segmentLengthSq;

const double clampedProjection = std::max(0., std::min(1., pointProjection));

const Vec2 closestPoint{segmentStart.uCoord + clampedProjection * segmentVector.uCoord,

segmentStart.vCoord + clampedProjection * segmentVector.vCoord};

return distanceSq(point, closestPoint);

{

const Vec3 segmentVector = segmentEnd - segmentStart;

const double segmentLengthSq = normSq(segmentVector);

if (segmentLengthSq <= kToleranceSq) {

return distanceSq(point, segmentStart);

}

const double pointProjection = dot(point - segmentStart, segmentVector) / segmentLengthSq;

const double clampedProjection = std::max(0., std::min(1., pointProjection));

const Vec3 closestPoint = segmentStart + segmentVector * clampedProjection;

return distanceSq(point, closestPoint);

{

gUU = dot(axisU, axisU);

gUV = dot(axisU, axisV);

gVV = dot(axisV, axisV);

{

gUU = radius * radius;

gUV = 0.;

gVV = 1.;

{

gUU = radiusAtHeight * radiusAtHeight;

gUV = 0.;

gVV = 1. + slope * slope;

{

const double parallelRadius = radius * std::sin(theta);

gUU = parallelRadius * parallelRadius;

gUV = 0.;

gVV = radius * radius;

double& gVV)

{

const double ringRadius = majorRadius + minorRadius * std::cos(phiTube);

gUU = ringRadius * ringRadius;

gUV = 0.;

gVV = minorRadius * minorRadius;

{

return std::abs(firstDistance - secondDistance) <=

kIntersectionTolerance * std::max(1., std::max(std::abs(firstDistance), std::abs(secondDistance)));

double distance = 0.;

Vec3 normal;

/// The hit lies within the trim's on-boundary band, so its inside/outside side is a tie-break, not data.

bool onTrimBoundary = false;

Vec2 start;

Vec2 end;

Vec2 direction() const { return end - start; }

double lengthSq() const

{

const Vec2 delta = end - start;

return delta.uCoord * delta.uCoord + delta.vCoord * delta.vCoord;

}

bool degenerate() const { return lengthSq() <= kToleranceSq; }

/// Squared distance from a parametric point to this edge.

double distanceSq(const Vec2& point) const { return pointSegmentDistanceSq(point, start, end); }

/// Closest point on this edge to \a point. Returns the projected point and its clamped

/// parameter \a parameter in [0, 1] (0 at start, 1 at end). Degenerate edges return start.

Vec2 closestPoint(const Vec2& point, double& parameter) const

{

const Vec2 segmentVector = end - start;

const double segmentLengthSq = segmentVector.uCoord * segmentVector.uCoord +

segmentVector.vCoord * segmentVector.vCoord;

if (segmentLengthSq <= kToleranceSq) {

parameter = 0.;

return start;

}

const double projection = ((point.uCoord - start.uCoord) * segmentVector.uCoord +

(point.vCoord - start.vCoord) * segmentVector.vCoord) /

segmentLengthSq;

parameter = std::max(0., std::min(1., projection));

return {start.uCoord + parameter * segmentVector.uCoord, start.vCoord + parameter * segmentVector.vCoord};

}

/// Accumulate the edge endpoints into a parametric axis-aligned bounding box.

void extendBounds(Vec2& lower, Vec2& upper) const

{

lower.uCoord = std::min({lower.uCoord, start.uCoord, end.uCoord});

lower.vCoord = std::min({lower.vCoord, start.vCoord, end.vCoord});

upper.uCoord = std::max({upper.uCoord, start.uCoord, end.uCoord});

upper.vCoord = std::max({upper.vCoord, start.vCoord, end.vCoord});

}

Valid, ///< well-formed and already correctly oriented

Reversed, ///< well-formed but re-oriented to match its role (simple, logged repair)

NonFinite, ///< a vertex/edge contained a non-finite coordinate

Open, ///< an explicit edge list did not form a closed loop

TooFewVertices, ///< fewer than three distinct vertices after cleanup

DegenerateVertex, ///< a non-adjacent vertex coincided (self-touching / pinched loop)

ZeroArea ///< the loop encloses no area

{

const double scale = metric.maxScale(uv);

return scale > kTolerance ? kTolerance / scale : 0.;

std::vector<Vec2> vertices;

WireRole role = WireRole::Outer;

/// For each stored segment its input segment, or -1 once a vertex was dropped; sidecar v3 edge identities key on it.

std::vector<int> sourceEdge;

int edgeCount() const { return static_cast<int>(vertices.size()); }

/// The stored segment that came from input segment \a inputIndex, or -1 if there is none.

int storedIndexOfSource(int inputIndex) const

{

for (size_t index = 0; index < sourceEdge.size(); ++index) {

if (sourceEdge[index] == inputIndex) {

return static_cast<int>(index);

}

}

return -1;

}

SurfaceEdge edge(int index) const

{

const int count = edgeCount();

return {vertices[index % count], vertices[(index + 1) % count]};

}

/// Build and validate the wire from an implicitly closed vertex ring; \a metric turns separations into 3D lengths.

bool initialize(const std::vector<Vec2>& inputVertices, WireRole wireRole, WireStatus& status,

const ParametricMetric& metric = {})

{

role = wireRole;

vertices.clear();

vertices.reserve(inputVertices.size());

bool droppedAVertex = false;

for (const auto& vertex : inputVertices) {

if (!finite(vertex)) {

status = WireStatus::NonFinite;

return false;

}

if (vertices.empty() || metric.distanceSq(vertices.back(), vertex) > kToleranceSq) {

vertices.push_back(vertex);

} else {

droppedAVertex = true;

}

}

// drop an explicit closing duplicate (first == last)

if (vertices.size() > 1 && metric.distanceSq(vertices.front(), vertices.back()) <= kToleranceSq) {

vertices.pop_back();

droppedAVertex = true;

}

if (vertices.size() < 3) {

status = WireStatus::TooFewVertices;

return false;

}

// reject self-touching loops (non-adjacent coincident vertices)

for (size_t firstIndex = 0; firstIndex < vertices.size(); ++firstIndex) {

for (size_t secondIndex = firstIndex + 1; secondIndex < vertices.size(); ++secondIndex) {

if (metric.distanceSq(vertices[firstIndex], vertices[secondIndex]) <= kToleranceSq) {

status = WireStatus::DegenerateVertex;

return false;

}

}

}

const double area = signedArea();

if (std::abs(area) <= kAreaTolerance) {

status = WireStatus::ZeroArea;

return false;

}

// segment i is input segment i unless a vertex was dropped; then it is unknown

const int storedCount = static_cast<int>(vertices.size());

sourceEdge.assign(static_cast<size_t>(storedCount), -1);

if (!droppedAVertex) {

for (int index = 0; index < storedCount; ++index) {

sourceEdge[static_cast<size_t>(index)] = index;

}

}

// outer wires must wind CCW (positive area), inner wires CW (negative area)

const bool wantPositiveArea = (role == WireRole::Outer);

if ((area > 0.) != wantPositiveArea) {

std::reverse(vertices.begin(), vertices.end());

// reversing the ring maps old vertex k to new index n-1-k, so new segment j spans old

// vertices n-1-j and n-2-j, i.e. it is old segment n-2-j traversed backwards

std::vector<int> reversedSource(static_cast<size_t>(storedCount), -1);

for (int index = 0; index < storedCount; ++index) {

reversedSource[static_cast<size_t>(index)] =

sourceEdge[static_cast<size_t>((storedCount - 2 - index % storedCount + 2 * storedCount) % storedCount)];

}

sourceEdge.swap(reversedSource);

status = WireStatus::Reversed;

return true;

}

status = WireStatus::Valid;

return true;

}

/// Build and validate the wire from an ordered edge list, joining within \a joinTolerance through \a metric, as CurveWire does.

bool initializeFromEdges(const std::vector<SurfaceEdge>& edges, WireRole wireRole, WireStatus& status,

const ParametricMetric& metric = {}, double joinTolerance = kWireJoinTolerance)

{

if (edges.size() < 3) {

status = WireStatus::TooFewVertices;

return false;

}

for (size_t edgeIndex = 0; edgeIndex < edges.size(); ++edgeIndex) {

if (!finite(edges[edgeIndex].start) || !finite(edges[edgeIndex].end)) {

status = WireStatus::NonFinite;

return false;

}

const Vec2& nextStart = edges[(edgeIndex + 1) % edges.size()].start;

if (metric.distanceSq(edges[edgeIndex].end, nextStart) > joinTolerance * joinTolerance) {

status = WireStatus::Open;

return false;

}

}

std::vector<Vec2> ringVertices;

ringVertices.reserve(edges.size());

for (const auto& singleEdge : edges) {

ringVertices.push_back(singleEdge.start);

}

return initialize(ringVertices, wireRole, status, metric);

}

double signedArea() const

{

double area = 0.;

for (size_t vertexIndex = 0; vertexIndex < vertices.size(); ++vertexIndex) {

const auto& currentVertex = vertices[vertexIndex];

const auto& nextVertex = vertices[(vertexIndex + 1) % vertices.size()];

area += currentVertex.uCoord * nextVertex.vCoord - nextVertex.uCoord * currentVertex.vCoord;

}

return 0.5 * area;

}

/// Accumulate this wire's vertices into a parametric axis-aligned bounding box. This is

/// independent of any concrete surface so cylinders, spheres and cones can reuse it.

void parametricBounds(Vec2& lower, Vec2& upper) const

{

for (const auto& vertex : vertices) {

lower.uCoord = std::min(lower.uCoord, vertex.uCoord);

lower.vCoord = std::min(lower.vCoord, vertex.vCoord);

upper.uCoord = std::max(upper.uCoord, vertex.uCoord);

upper.vCoord = std::max(upper.vCoord, vertex.vCoord);

}

}

/// The de-duplicated vertex ring, closed back to its first vertex.

std::vector<Vec2> sampledBoundary() const

{

std::vector<Vec2> samples;

if (vertices.empty()) {

return samples;

}

samples.reserve(vertices.size() + 1);

samples.insert(samples.end(), vertices.begin(), vertices.end());

samples.push_back(vertices.front());

return samples;

}

/// Classify against the polygon with an on-boundary half-width of \a band, in parametric units.

WireClassification classify(const Vec2& point, double band) const

{

const double bandSq = band * band;

bool inside = false;

for (size_t vertexIndex = 0; vertexIndex < vertices.size(); ++vertexIndex) {

const auto& segmentStart = vertices[vertexIndex];

const auto& segmentEnd = vertices[(vertexIndex + 1) % vertices.size()];

if (pointSegmentDistanceSq(point, segmentStart, segmentEnd) <= bandSq) {

return WireClassification::Boundary;

}

const bool crossesScanline = (segmentStart.vCoord > point.vCoord) != (segmentEnd.vCoord > point.vCoord);

if (crossesScanline) {

const double intersectionU = segmentStart.uCoord + (point.vCoord - segmentStart.vCoord) *

(segmentEnd.uCoord - segmentStart.uCoord) /

(segmentEnd.vCoord - segmentStart.vCoord);

if (point.uCoord < intersectionU) {

inside = !inside;

}

}

}

return inside ? WireClassification::Inside : WireClassification::Outside;

}

/// \a metric sizes the band only: a polygon is exact, so its band is the length floor.

WireClassification classify(const Vec2& point, const ParametricMetric& metric = {}) const

{

return classify(point, trimLengthFloor(metric, point));

}

const Vec2& thirdVertex)

{

const double firstCross = cross2D(secondVertex - firstVertex, point - firstVertex);

const double secondCross = cross2D(thirdVertex - secondVertex, point - secondVertex);

const double thirdCross = cross2D(firstVertex - thirdVertex, point - thirdVertex);

return firstCross >= -kTolerance && secondCross >= -kTolerance && thirdCross >= -kTolerance;

{

std::vector<int> remainingIndices;

remainingIndices.reserve(wire.vertices.size());

if (wire.signedArea() >= 0.) {

for (size_t vertexIndex = 0; vertexIndex < wire.vertices.size(); ++vertexIndex) {

remainingIndices.push_back(static_cast<int>(vertexIndex));

}

} else {

for (size_t reverseIndex = wire.vertices.size(); reverseIndex > 0; --reverseIndex) {

remainingIndices.push_back(static_cast<int>(reverseIndex - 1));

}

}

std::vector<std::array<int, 3>> triangles;

size_t guardCounter = 0;

while (remainingIndices.size() > 3 && guardCounter++ < wire.vertices.size() * wire.vertices.size()) {

bool clippedEar = false;

for (size_t indexPosition = 0; indexPosition < remainingIndices.size(); ++indexPosition) {

const int previousIndex = remainingIndices[(indexPosition + remainingIndices.size() - 1) % remainingIndices.size()];

const int currentIndex = remainingIndices[indexPosition];

const int nextIndex = remainingIndices[(indexPosition + 1) % remainingIndices.size()];

const auto& previousVertex = wire.vertices[previousIndex];

const auto& currentVertex = wire.vertices[currentIndex];

const auto& nextVertex = wire.vertices[nextIndex];

if (cross2D(currentVertex - previousVertex, nextVertex - currentVertex) <= kTolerance) {

continue;

}

bool containsOtherVertex = false;

for (int candidateIndex : remainingIndices) {

if (candidateIndex == previousIndex || candidateIndex == currentIndex || candidateIndex == nextIndex) {

continue;

}

if (pointInTriangle(wire.vertices[candidateIndex], previousVertex, currentVertex, nextVertex)) {

containsOtherVertex = true;

break;

}

}

if (containsOtherVertex) {

continue;

}

triangles.push_back({previousIndex, currentIndex, nextIndex});

remainingIndices.erase(remainingIndices.begin() + indexPosition);

clippedEar = true;

break;

}

if (!clippedEar) {

break;

}

}

if (remainingIndices.size() == 3) {

triangles.push_back({remainingIndices[0], remainingIndices[1], remainingIndices[2]});

}

return triangles;

{

constexpr int fullTurnChunks = static_cast<int>(kTwoPi / kCoverChunkAngle); // eight

return std::max(1, std::min(fullTurnChunks, static_cast<int>(std::ceil(span / kCoverChunkAngle))));

{

const double atStart = a * std::cos(t0) + b * std::sin(t0);

const double atEnd = a * std::cos(t1) + b * std::sin(t1);

minimum = std::min(atStart, atEnd);

maximum = std::max(atStart, atEnd);

const double amplitude = std::hypot(a, b);

const double crest = std::atan2(b, a);

// shifted into [t0, t0 + 2pi), where a span of at most a full turn makes "<= t1" exactly the

// test for falling inside the interval

const double crestInRange = crest - kTwoPi * std::floor((crest - t0) / kTwoPi);

if (crestInRange <= t1) {

maximum = amplitude;

}

const double trough = crest + kPi;

const double troughInRange = trough - kTwoPi * std::floor((trough - t0) / kTwoPi);

if (troughInRange <= t1) {

minimum = -amplitude;

}

{

double minimum = 0.;

double maximum = 0.;

sinusoidRange(a, b, t0, t1, minimum, maximum);

return minimum;

{

double minimum = 0.;

double maximum = 0.;

sinusoidRange(a, b, t0, t1, minimum, maximum);

return maximum;

{

if (sweep >= kTwoPi - kTolerance) {

return true;

}

double delta = angle - start;

delta -= kTwoPi * std::floor(delta / kTwoPi); // wrap into [0, 2pi)

return delta <= sweep + tolerance || delta >= kTwoPi - tolerance;

{

nodes.assign(std::max(n, 1), 0.);

weights.assign(std::max(n, 1), 0.);

if (n < 1) {

return;

}

for (int i = 0; i < n; ++i) {

double root = std::cos(kPi * (i + 0.75) / (n + 0.5)); // asymptotic initial guess

double derivative = 1.;

for (int iteration = 0; iteration < 100; ++iteration) {

double previous = 1.;

double current = root;

for (int degreeIndex = 2; degreeIndex <= n; ++degreeIndex) {

const double next = ((2 * degreeIndex - 1) * root * current - (degreeIndex - 1) * previous) / degreeIndex;

previous = current;

current = next;

}

derivative = n * (root * current - previous) / (root * root - 1.);

const double delta = current / derivative;

root -= delta;

if (std::abs(delta) < 1.e-15) {

break;

}

}

nodes[i] = root;

weights[i] = 2. / ((1. - root * root) * derivative * derivative);

}

{

const double discriminant = coeffQ * coeffQ / 4. + coeffP * coeffP * coeffP / 27.;

if (!(coeffP < 0.) || discriminant > 0.) {

const double sqrtDiscriminant = std::sqrt(std::max(0., discriminant));

roots[0] = std::cbrt(-0.5 * coeffQ + sqrtDiscriminant) + std::cbrt(-0.5 * coeffQ - sqrtDiscriminant);

return 1;

}

// three real roots: coeffP < 0 here, so the trigonometric form is well defined

const double magnitude = 2. * std::sqrt(-coeffP / 3.);

const double cosineArgument = std::max(-1., std::min(1., 3. * coeffQ / (coeffP * magnitude)));

const double baseAngle = std::acos(cosineArgument);

for (int branch = 0; branch < 3; ++branch) {

roots[branch] = magnitude * std::cos((baseAngle - kTwoPi * branch) / 3.);

}

return 3;

NotAQuartic, ///< the leading coefficient vanishes; no roots are produced

Biquadratic, ///< the depressed quartic's odd term is zero, so y^4 + p y^2 + r = 0 is solved directly

Resolvent ///< Ferrari's general branch, through the resolvent cubic

std::array<double, 4> value{};

int count = 0;

void push_back(double root)

{

assert(count < 4 && "QuarticRoots holds at most four roots");

value[count++] = root;

}

double* begin() { return value.data(); }

double* end() { return value.data() + count; }

const double* begin() const { return value.data(); }

const double* end() const { return value.data() + count; }

size_t size() const { return static_cast<size_t>(count); }

bool empty() const { return count == 0; }

double operator[](size_t index) const { return value[index]; }

QuarticBranch* takenBranch = nullptr)

{

const auto note = [takenBranch](QuarticBranch branch) {

if (takenBranch) {

*takenBranch = branch;

}

};

note(QuarticBranch::NotAQuartic);

QuarticRoots roots;

// A genuine quartic needs only a non-zero leading coefficient. There is no scale to compare it

// against -- the normalisation below handles any coefficient ratio -- so the test is exact.

if (!(std::abs(a4) > 0.)) {

return roots; // the torus caller guarantees a4 = |dir|^4 > 0

}

// monic x^4 + b x^3 + c x^2 + d x + e

double coeffB = a3 / a4, coeffC = a2 / a4, coeffD = a1 / a4, coeffE = a0 / a4;

if (!std::isfinite(coeffB) || !std::isfinite(coeffC) || !std::isfinite(coeffD) || !std::isfinite(coeffE)) {

return roots; // a4 is denormal-small next to the rest, or an input was not finite

}

// Cauchy root bound rounded up to a power of two, so x = scale * y is exact; x^4 = 0 keeps scale = 1

const double rootBound = std::max({std::abs(coeffB), std::sqrt(std::abs(coeffC)),

std::cbrt(std::abs(coeffD)), std::sqrt(std::sqrt(std::abs(coeffE)))});

int boundExponent = 0;

std::frexp(rootBound, &boundExponent);

const double scale = std::ldexp(1., boundExponent);

coeffB /= scale;

coeffC /= scale * scale;

coeffD /= scale * scale * scale;

coeffE /= scale * scale * scale * scale;

// depress with y = z - b/4: z^4 + p z^2 + q z + r

const double termP = coeffC - 3. * coeffB * coeffB / 8.;

const double termQ = coeffD - coeffB * coeffC / 2. + coeffB * coeffB * coeffB / 8.;

const double termR =

coeffE - coeffB * coeffD / 4. + coeffB * coeffB * coeffC / 16. - 3. * coeffB * coeffB * coeffB * coeffB / 256.;

const double shift = -coeffB / 4.;

auto addQuadraticRoots = [&](double quadB, double quadC) {

const double discriminant = quadB * quadB - 4. * quadC;

if (discriminant < 0.) {

return; // complex pair

}

const double sqrtDiscriminant = std::sqrt(discriminant);

roots.push_back(shift + 0.5 * (-quadB - sqrtDiscriminant));

roots.push_back(shift + 0.5 * (-quadB + sqrtDiscriminant));

};

auto addBiquadraticRoots = [&]() {

// biquadratic z^4 + p z^2 + r = 0

const double discriminant = termP * termP - 4. * termR;

if (discriminant < 0.) {

return;

}

const double sqrtDiscriminant = std::sqrt(discriminant);

for (const double zSquared : {0.5 * (-termP + sqrtDiscriminant), 0.5 * (-termP - sqrtDiscriminant)}) {

if (zSquared >= 0.) {

const double z = std::sqrt(zSquared);

roots.push_back(shift + z);

roots.push_back(shift - z);

}

}

};

// q is zero to the precision of its terms, which normalisation bounds by 1: kQuarticEpsilon over the whole quartic, not over q's terms

bool biquadratic = std::abs(termQ) <= kQuarticEpsilon;

if (!biquadratic) {

note(QuarticBranch::Resolvent);

// resolvent cubic m^3 + p m^2 + (p^2/4 - r) m - q^2/8 = 0; its largest real root is > 0

const double cubicA2 = termP;

const double cubicA1 = termP * termP / 4. - termR;

const double cubicA0 = -termQ * termQ / 8.;

const double cubicP = cubicA1 - cubicA2 * cubicA2 / 3.;

const double cubicQ = 2. * cubicA2 * cubicA2 * cubicA2 / 27. - cubicA2 * cubicA1 / 3. + cubicA0;

std::array<double, 3> cubicRoots;

const int cubicCount = solveDepressedCubic(cubicP, cubicQ, cubicRoots);

double resolvent = 0.;

for (int index = 0; index < cubicCount; ++index) {

resolvent = std::max(resolvent, cubicRoots[index] - cubicA2 / 3.);

}

// a resolvent below the resolution of its cubic is noise; then the biquadratic branch is the better-conditioned answer

const double resolventScale = std::max({std::abs(cubicA2), std::sqrt(std::abs(cubicA1)),

std::cbrt(std::abs(cubicA0))});

if (resolvent > kQuarticEpsilon * resolventScale) {

const double sqrtTwoResolvent = std::sqrt(2. * resolvent);

const double linearTerm = sqrtTwoResolvent * termQ / (4. * resolvent);

addQuadraticRoots(-sqrtTwoResolvent, termP / 2. + resolvent + linearTerm);

addQuadraticRoots(sqrtTwoResolvent, termP / 2. + resolvent - linearTerm);

} else {

biquadratic = true;

}

}

if (biquadratic) {

note(QuarticBranch::Biquadratic);

addBiquadraticRoots();

}

// Newton polish against the monic quartic; a step longer than the Cauchy bound 2, or non-finite, is rejected

auto quartic = [&](double x) { return (((x + coeffB) * x + coeffC) * x + coeffD) * x + coeffE; };

auto quarticDerivative = [&](double x) { return ((4. * x + 3. * coeffB) * x + 2. * coeffC) * x + coeffD; };

for (double& root : roots) {

for (int iteration = 0; iteration < 2; ++iteration) {

const double step = quartic(root) / quarticDerivative(root);

if (std::isfinite(step) && std::abs(step) <= 2.) {

root -= step;

}

}

}

for (double& root : roots) {

root *= scale; // exact: scale is a power of two

}

return roots;

CurveKind kind = CurveKind::Line;

Vec2 lineStart; ///< line: start point (unused for arcs)

Vec2 lineEnd; ///< line: end point (unused for arcs)

Vec2 center; ///< arc: circle centre (unused for lines)

double radius = 0.; ///< arc: circle radius

double startAngle = 0.; ///< arc: start angle [rad]

double endAngle = 0.; ///< arc: end angle [rad] (sweep = endAngle - startAngle)

/// \name B-spline data (kind == BSpline): poles, optional weights and a clamped knot vector; the curve parameter runs on [0, 1]. @{

int degree = 0;

std::vector<Vec2> poles;

std::vector<double> weights;

std::vector<double> knots;

/// The flattened on-curve polyline, both ends included; CurveWire::initialize fills it and reversing clears it.

mutable std::vector<Vec2> bsplineCache;

/// @}

/// \name Loop-canonical endpoints: the seam vertices the curve's neighbours agree on, substituted at the polyline's ends

/// @{