linearAlgebra

Pure-C operations on double-precision vectors and matrices.

Matrices use contiguous row-major storage. Callers provide correctly sized, aligned buffers and representable dimensions. Inputs are preserved unless an output overlaps them; supported overlap and checked size limits are documented per operation. Functions without an explicit nonfinite policy use ordinary floating-point arithmetic. Norms sum squared components without range scaling, so extreme finite magnitudes can overflow or underflow. No function allocates heap memory or takes ownership of caller-provided storage. Operations needing scratch space use bounded local arrays.

Defines

DB0_EPS: Absolute small-value cutoff used by normalization and inverse checks.

LINEAR_ALGEBRA_MAX_ARRAY_SIZE: Maximum number of doubles in a local scratch array.

MXINDEX(dim2, row, col)

Compute a row-major matrix element index.

Parameters:

dim2 – Number of matrix columns.
row – Zero-based row index.
col – Zero-based column index.

Returns:

Linear index of the requested element.

Functions

void vElementwiseMult(double *v1, size_t dim, double *v2, double *result)

Multiply corresponding elements, allowing overlap with either input.

Note

Empty inputs cause no memory access; pointers may be null.

Note

If overlapping inputs require opposite traversal directions, the smaller input overlap with result must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Other layouts have no fixed size limit. Unrepresentable byte counts or excessive scratch requirements are logged and leave result unchanged. No heap memory is used.

Parameters:

v1 – First input vector.
dim – Number of elements.
v2 – Second input vector.
result – Output of the same shape; may overlap either or both inputs.

void vCopy(double *v, size_t dim, double *result)

Copy a vector, allowing overlapping input and output storage.

Note

If dim is zero, no memory is accessed and either pointer may be null.

Note

Unrepresentable byte counts are logged and leave result unchanged.

Parameters:

v – Input vector containing dim elements.
dim – Number of elements to copy.
result – Output vector containing dim elements; may overlap v.

void vSetZero(double *v, size_t dim)

Fill the output with zeros.

Note

An empty output causes no memory access; the output pointer may be null.

Parameters:

v – Output array.
dim – Number of elements.

void vSetOnes(double *v, size_t dim)

Fill the output with ones.

Note

An empty output causes no memory access; the output pointer may be null.

Parameters:

v – Output array.
dim – Number of elements.

void vAdd(double *v1, size_t dim, double *v2, double *result)

Add corresponding elements, allowing overlap with either input.

Note

Empty inputs cause no memory access; pointers may be null.

Note

If overlapping inputs require opposite traversal directions, the smaller input overlap with result must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Other layouts have no fixed size limit. Unrepresentable byte counts or excessive scratch requirements are logged and leave result unchanged. No heap memory is used.

Parameters:

v1 – First input vector.
dim – Number of elements.
v2 – Second input vector.
result – Output of the same shape; may overlap either or both inputs.

void vSubtract(double *v1, size_t dim, double *v2, double *result)

Subtract corresponding elements, allowing overlap with either input.

Note

Empty inputs cause no memory access; pointers may be null.

Note

If overlapping inputs require opposite traversal directions, the smaller input overlap with result must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Other layouts have no fixed size limit. Unrepresentable byte counts or excessive scratch requirements are logged and leave result unchanged. No heap memory is used.

Parameters:

v1 – First input vector.
dim – Number of elements.
v2 – Second input vector.
result – Output of the same shape; may overlap either or both inputs.

void vScale(double scaleFactor, double *v, size_t dim, double *result)

Scale every element, allowing overlapping input and output.

Note

Empty input causes no memory access; pointers may be null. Unrepresentable byte counts are logged and leave result unchanged. No heap storage or fixed scratch limit is used.

Parameters:

scaleFactor – Multiplier.
v – Input vector.
dim – Number of elements.
result – Output of the same shape; may overlap the input.

double vDot(double *v1, size_t dim, double *v2)

Compute the dot product of two vectors.

Note

An empty input returns zero without accessing either pointer.

Parameters:

v1 – First input vector.
dim – Number of elements.
v2 – Second input vector.

Returns:

Sum of corresponding component products.

void vOuterProduct(double *v1, size_t dim1, double *v2, size_t dim2, void *result)

Form a row-major outer product, allowing either input to overlap the output.

Note

A zero dimension causes no memory access; pointers may be null.

Note

Overlapping input vectors are saved in local scratch storage; their combined lengths must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Disjoint inputs and one-row or one-column outputs have no fixed size limit. Unrepresentable sizes or excessive scratch requirements are logged and leave result unchanged. No heap memory is used.

Parameters:

v1 – First vector with dim1 elements.
dim1 – Number of output rows.
v2 – Second vector with dim2 elements.
dim2 – Number of output columns.
result – Output matrix with dim1 * dim2 elements.

void vtMultM(double *v, void *mx, size_t dim1, size_t dim2, void *result)

Compute v transpose times mx.

Note

dim2 must not exceed LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid sizes are logged and leave result unchanged. An empty output causes no memory access. A zero inner dimension produces zeros. Pointers to empty arrays may be null.

Parameters:

v – Input vector with dim1 elements.
mx – Input matrix in contiguous row-major storage.
dim1 – Number of matrix rows.
dim2 – Number of matrix columns.
result – Output vector with dim2 elements; may overlap either input.

void vtMultMt(double *v, void *mx, size_t dim1, size_t dim2, void *result)

Compute v transpose times mx transpose.

Note

dim1 must not exceed LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid sizes are logged and leave result unchanged. An empty output causes no memory access. A zero inner dimension produces zeros. Pointers to empty arrays may be null.

Parameters:

v – Input vector with dim2 elements.
mx – Input matrix in contiguous row-major storage.
dim1 – Number of matrix rows.
dim2 – Number of matrix columns.
result – Output vector with dim1 elements; may overlap either input.

double vNorm(double *v, size_t dim)

Compute the Euclidean norm of a vector.

Note

The sum is not scaled; extreme finite magnitudes may overflow or underflow.

Note

An empty input returns zero without accessing the pointer.

Parameters:

v – Input vector.
dim – Number of elements.

Returns:

Square root of the sum of squared components.

double vMax(double *v, size_t dim)

Find the maximum element of a nonempty vector.

Parameters:

v – Input vector.
dim – Number of elements.

Returns:

Maximum element; any NaN element gives NaN, regardless of its position. Infinities retain their usual ordering. An empty input is logged and returns NaN without accessing v.

double vMaxAbs(double *v, size_t dim)

Find the maximum absolute element of a nonempty vector.

Parameters:

v – Input vector.
dim – Number of elements.

Returns:

Maximum absolute element; any NaN element gives NaN, regardless of its position. Either infinity gives positive infinity when no NaN is present. An empty input is logged and returns NaN without accessing v.

void vNormalize(double *v, size_t dim, double *result)

Normalize a vector, allowing overlapping input and output.

Note

A finite norm greater than DB0_EPS gives a unit vector; smaller norms give zeros. A nonfinite computed norm, including from NaN or infinite input, produces NaN in every output component. The norm computation does not scale extreme finite magnitudes.

Note

Empty input causes no memory access; pointers may be null. Unrepresentable byte counts are logged and leave result unchanged. No heap storage or fixed scratch limit is used.

Parameters:

v – Input vector.
dim – Number of elements.
result – Output vector; may overlap the input.

int vIsEqual(double *v1, size_t dim, double *v2, double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Note

An empty array compares equal for a valid tolerance, without accessing its storage.

Parameters:

v1 – Input array.
dim – Number of elements.
v2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

int vIsZero(double *v, size_t dim, double accuracy)

Compare every element with zero using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Note

An empty array compares equal for a valid tolerance, without accessing its storage.

Parameters:

v – Input array.
dim – Number of elements.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

void vPrint(FILE *pFile, const char *name, double *v, size_t dim)

Print the array name and values to a stream.

Parameters:

pFile – Destination stream.
name – Display label, passed as data rather than as a format string.
v – Input array.
dim – Number of elements.

void vSort(double *Input, double *Output, size_t dim)

Stably sort a vector in ascending order, with NaNs last.

Note

An empty input causes no memory access; pointers may be null. Unrepresentable byte counts are logged and leave Output unchanged. Equal values, including signed zeros and NaNs, retain their input order. Already sorted input takes linear time; worst-case work remains quadratic.

Parameters:

Input – Input vector containing dim elements.
Output – Output vector containing dim elements; may overlap Input.
dim – Number of elements to sort.

void v2Set(double v0, double v1, double result[2])

Set vector components in index order.

Parameters:

v0 – Value for component 0.
v1 – Value for component 1.
result – Output array.

void v2Copy(double v[2], double result[2])

Copy a two-component vector.

Parameters:

v – Input vector.
result – Output vector, which may overlap v.

void v2Scale(double scaleFactor, double v[2], double result[2])

Scale every element.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

scaleFactor – Multiplier.
v – Input array.
result – Output of the same shape; may overlap the input.

void v2SetZero(double v[2])

Fill the output with zeros.

Parameters:: v – Output array.

double v2Dot(double v1[2], double v2[2])

Compute the dot product of two vectors.

Parameters:

v1 – First input vector.
v2 – Second input vector.

Returns:

Sum of corresponding component products.

int v2IsEqual(double v1[2], double v2[2], double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

v1 – Input array.
v2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

int v2IsZero(double v[2], double accuracy)

Compare every element with zero using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

v – Input array.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

void v2Add(double v1[2], double v2[2], double result[2])

Add corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

v1 – First input.
v2 – Second input.
result – Output of the same shape; may overlap either input.

void v2Subtract(double v1[2], double v2[2], double result[2])

Subtract corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

v1 – First input.
v2 – Second input.
result – Output of the same shape; may overlap either input.

double v2Norm(double v1[2])

Compute the Euclidean norm of a vector.

Note

The sum is not scaled; extreme finite magnitudes may overflow or underflow.

Parameters:: v1 – Input vector.
Returns:: Square root of the sum of squared components.

void v2Normalize(double v1[2], double result[2])

Normalize a vector, allowing overlapping input and output.

Note

A finite norm greater than DB0_EPS gives a unit vector; smaller norms give zeros. A nonfinite computed norm, including from NaN or infinite input, produces NaN in every output component. The norm computation does not scale extreme finite magnitudes.

Parameters:

v1 – Input vector.
result – Output vector; may overlap the input.

void v3Set(double v0, double v1, double v2, double result[3])

Set vector components in index order.

Parameters:

v0 – Value for component 0.
v1 – Value for component 1.
v2 – Value for component 2.
result – Output array.

void v3Copy(double v[3], double result[3])

Copy a three-component vector.

Parameters:

v – Input vector.
result – Output vector, which may overlap v.

void v3SetZero(double v[3])

Fill the output with zeros.

Parameters:: v – Output array.

void v3Add(double v1[3], double v2[3], double result[3])

Add corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

v1 – First input.
v2 – Second input.
result – Output of the same shape; may overlap either input.

void v3Subtract(double v1[3], double v2[3], double result[3])

Subtract corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

v1 – First input.
v2 – Second input.
result – Output of the same shape; may overlap either input.

void v3Scale(double scaleFactor, double v[3], double result[3])

Scale every element.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

scaleFactor – Multiplier.
v – Input array.
result – Output of the same shape; may overlap the input.

double v3Dot(double v1[3], double v2[3])

Compute the dot product of two vectors.

Parameters:

v1 – First input vector.
v2 – Second input vector.

Returns:

Sum of corresponding component products.

void v3OuterProduct(double v1[3], double v2[3], double result[3][3])

Form the outer product of two three-component vectors.

Parameters:

v1 – First input vector, supplying the row factors.
v2 – Second input vector, supplying the column factors.
result – Output matrix with result[i][j] = v1[i] * v2[j]; may overlap either or both inputs.

void v3tMultM33(double v[3], double mx[3][3], double result[3])

Compute v transpose times mx.

Parameters:

v – Input vector.
mx – Input matrix.
result – Output array; may overlap either input.

void v3tMultM33t(double v[3], double mx[3][3], double result[3])

Compute v transpose times mx transpose.

Parameters:

v – Input vector.
mx – Input matrix.
result – Output array; may overlap either input.

double v3Norm(double v[3])

Compute the Euclidean norm of a vector.

Note

The sum is not scaled; extreme finite magnitudes may overflow or underflow.

Parameters:: v – Input vector.
Returns:: Square root of the sum of squared components.

void v3Normalize(double v[3], double result[3])

Normalize a vector, allowing overlapping input and output.

Note

A finite norm greater than DB0_EPS gives a unit vector; smaller norms give zeros. A nonfinite computed norm, including from NaN or infinite input, produces NaN in every output component. The norm computation does not scale extreme finite magnitudes.

Parameters:

v – Input vector.
result – Output vector; may overlap the input.

int v3IsEqual(double v1[3], double v2[3], double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

v1 – Input array.
v2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

int v3IsEqualRel(double v1[3], double v2[3], double accuracy)

Compare finite vectors relative to the Euclidean norm of the first vector.

Note

A zero reference vector or zero tolerance requires exact equality. NaN and infinity are rejected in either vector or the tolerance, including comparisons of equal infinities.

Parameters:

v1 – Reference vector with three components.
v2 – Vector to compare, in the same units as v1.
accuracy – Finite, nonnegative relative tolerance (dimensionless).

Returns:

1 if every component satisfies abs(v1[i] - v2[i]) / norm(v1) <= accuracy, otherwise 0.

int v3IsZero(double v[3], double accuracy)

Compare every element with zero using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

v – Input array.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

void v3Print(FILE *pFile, const char *name, double v[3])

Print the array name and values to a stream.

Parameters:

pFile – Destination stream.
name – Display label, passed as data rather than as a format string.
v – Input array.

void v3Cross(double v1[3], double v2[3], double result[3])

Compute the cross product of two vectors.

Parameters:

v1 – Input vector.
v2 – Input vector.
result – Output vector; may overlap any input.

void v3Perpendicular(double v[3], double result[3])

Compute a normalized vector perpendicular to the input.

Note

Zero input returns zero. Nonzero inputs must be finite and satisfy at least one of abs(v[0]) > DB0_EPS, abs(v[1]) > DB0_EPS, or v[2] != 0. The intermediate perpendicular components and their sum of squares must also remain finite. These restrictions apply even when the input norm is representable: for example, {1e-31, 0, 0} is unsupported and can raise a divide-by-zero exception and produce NaNs.

Parameters:

v – Input vector.
result – Output vector; may overlap any input.

void v3Tilde(double v[3], double result[3][3])

Form the skew-symmetric matrix whose product with a vector w is v cross w.

Parameters:

v – Input vector with three components.
result – Output matrix, which may overlap the input vector.

void v3Sort(double v[3], double result[3])

Stably sort three values in descending order, with NaNs last.

Note

Equal values, including signed zeros and NaNs, retain their input order.

Parameters:

v – Input vector.
result – Output vector; may overlap v.

void v3PrintScreen(const char *name, double v[3])

Print the array name and values to standard output.

Parameters:

name – Display label, passed as data rather than as a format string.
v – Input array.

void v4Set(double v0, double v1, double v2, double v3, double result[4])

Set vector components in index order.

Parameters:

v0 – Value for component 0.
v1 – Value for component 1.
v2 – Value for component 2.
v3 – Value for component 3.
result – Output array.

void v4Copy(double v[4], double result[4])

Copy a four-component vector.

Parameters:

v – Input vector.
result – Output vector, which may overlap v.

void v4SetZero(double v[4])

Fill the output with zeros.

Parameters:: v – Output array.

double v4Dot(double v1[4], double v2[4])

Compute the dot product of two vectors.

Parameters:

v1 – First input vector.
v2 – Second input vector.

Returns:

Sum of corresponding component products.

double v4Norm(double v[4])

Compute the Euclidean norm of a vector.

Note

The sum is not scaled; extreme finite magnitudes may overflow or underflow.

Parameters:: v – Input vector.
Returns:: Square root of the sum of squared components.

int v4IsEqual(double v1[4], double v2[4], double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

v1 – Input array.
v2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

int v4IsZero(double v[4], double accuracy)

Compare every element with zero using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

v – Input array.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

void v6Set(double v0, double v1, double v2, double v3, double v4, double v5, double result[6])

Set vector components in index order.

Parameters:

v0 – Value for component 0.
v1 – Value for component 1.
v2 – Value for component 2.
v3 – Value for component 3.
v4 – Value for component 4.
v5 – Value for component 5.
result – Output array.

void v6Copy(double v[6], double result[6])

Copy a six-component vector.

Parameters:

v – Input vector.
result – Output vector, which may overlap v.

double v6Dot(double v1[6], double v2[6])

Compute the dot product of two vectors.

Parameters:

v1 – First input vector.
v2 – Second input vector.

Returns:

Sum of corresponding component products.

void v6Scale(double scaleFactor, double v[6], double result[6])

Scale every element.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

scaleFactor – Multiplier.
v – Input array.
result – Output of the same shape; may overlap the input.

void v6OuterProduct(double v1[6], double v2[6], double result[6][6])

Form the outer product of two six-component vectors.

Parameters:

v1 – First input vector, supplying the row factors.
v2 – Second input vector, supplying the column factors.
result – Output matrix with result[i][j] = v1[i] * v2[j]; may overlap either or both inputs.

int v6IsEqual(double v1[6], double v2[6], double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

v1 – Input array.
v2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

void mLeastSquaresInverse(void *mx, size_t dim1, size_t dim2, void *result)

Compute the least-squares inverse using column-pivoted Householder QR.

Note

Requires dim1 >= dim2 and a finite, numerically full-column-rank matrix. Both dim1 * dim2 and dim2 * dim2 must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid dimensions are logged and leave result unchanged. An empty input causes no memory access; pointers may be null.

Note

Rank is rejected when a pivot column norm is at most DBL_EPSILON * dim1 times the largest original column norm. The squared product of the R diagonal must also exceed DB0_EPS, preserving the absolute Gram-determinant cutoff. Failure of these checks, nonfinite input or pivot norms, or nonfinite inverse entries is logged and zeros result. Overflow of the determinant-related product alone does not imply failure if the remaining checks pass. This routine does not compute a rank-deficient pseudoinverse.

Parameters:

mx – Input matrix in contiguous row-major storage.
dim1 – Number of input rows.
dim2 – Number of input columns.
result – Output matrix with dim2 rows and dim1 columns; may overlap mx.

void mMinimumNormInverse(void *mx, size_t dim1, size_t dim2, void *result)

Compute the minimum-norm inverse using column-pivoted Householder QR of the transpose.

Note

Requires dim1 <= dim2 and a finite, numerically full-row-rank matrix. Both dim1 * dim2 and dim1 * dim1 must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid dimensions are logged and leave result unchanged. An empty input causes no memory access; pointers may be null.

Note

The rank and determinant checks are those of mLeastSquaresInverse applied to the transpose, using dim2 for the rank tolerance’s dimension factor. Failure of these checks, nonfinite input or pivot norms, or nonfinite inverse entries is logged and zeros result. Overflow of the determinant-related product alone does not imply failure if the remaining checks pass. This routine does not compute a rank-deficient pseudoinverse.

Parameters:

mx – Input matrix in contiguous row-major storage.
dim1 – Number of input rows.
dim2 – Number of input columns.
result – Output matrix with dim2 rows and dim1 columns; may overlap mx.

void mCopy(void *mx, size_t dim1, size_t dim2, void *result)

Copy a matrix, allowing overlapping input and output storage.

Note

If either dimension is zero, no memory is accessed and either pointer may be null.

Note

Unrepresentable byte counts are logged and leave result unchanged.

Parameters:

mx – Input matrix stored contiguously in row-major order.
dim1 – Number of rows.
dim2 – Number of columns.
result – Output matrix of the same shape and layout; may overlap mx.

void mSetZero(void *result, size_t dim1, size_t dim2)

Fill the output with zeros.

Note

An empty output causes no memory access; the output pointer may be null.

Parameters:

result – Output array.
dim1 – Number of rows.
dim2 – Number of columns.

void mSetIdentity(void *result, size_t dim1, size_t dim2)

Set the diagonal to one and every other entry to zero.

Note

An empty output causes no memory access; the output pointer may be null.

Parameters:

result – Output array.
dim1 – Number of rows.
dim2 – Number of columns.

void mDiag(void *v, size_t dim, void *result)

Form a diagonal matrix from a vector, setting all off-diagonal entries to zero.

Note

If dim is zero, no memory is accessed and either pointer may be null.

Note

Unrepresentable byte counts are logged and leave result unchanged.

Parameters:

v – Input vector containing dim elements.
dim – Number of input elements and number of rows and columns in the output.
result – Output matrix with dim * dim elements in row-major order; may overlap v.

void mTranspose(void *mx, size_t dim1, size_t dim2, void *result)

Transpose a row-major matrix, allowing overlapping input and output.

Note

dim1 * dim2 must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid sizes are logged and leave result unchanged. An empty input causes no memory access; pointers may be null.

Parameters:

mx – Input matrix.
dim1 – Number of input rows.
dim2 – Number of input columns.
result – Output matrix with dim2 rows and dim1 columns.

void mAdd(void *mx1, size_t dim1, size_t dim2, void *mx2, void *result)

Add corresponding elements, allowing overlap with either input.

Note

Empty inputs cause no memory access; pointers may be null.

Note

If overlapping inputs require opposite traversal directions, the smaller input overlap with result must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Other layouts have no fixed size limit. Unrepresentable byte counts or excessive scratch requirements are logged and leave result unchanged. No heap memory is used.

Parameters:

mx1 – First input matrix.
dim1 – Number of rows.
dim2 – Number of columns.
mx2 – Second input matrix.
result – Output of the same shape; may overlap either or both inputs.

void mSubtract(void *mx1, size_t dim1, size_t dim2, void *mx2, void *result)

Subtract corresponding elements, allowing overlap with either input.

Note

Empty inputs cause no memory access; pointers may be null.

Note

If overlapping inputs require opposite traversal directions, the smaller input overlap with result must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Other layouts have no fixed size limit. Unrepresentable byte counts or excessive scratch requirements are logged and leave result unchanged. No heap memory is used.

Parameters:

mx1 – First input matrix.
dim1 – Number of rows.
dim2 – Number of columns.
mx2 – Second input matrix.
result – Output of the same shape; may overlap either or both inputs.

void mScale(double scaleFactor, void *mx, size_t dim1, size_t dim2, void *result)

Scale every element, allowing overlapping input and output.

Note

Empty input causes no memory access; pointers may be null. Unrepresentable byte counts are logged and leave result unchanged. No heap storage or fixed scratch limit is used.

Parameters:

scaleFactor – Multiplier.
mx – Input matrix.
dim1 – Number of rows.
dim2 – Number of columns.
result – Output of the same shape; may overlap the input.

void mMultM(void *mx1, size_t dim11, size_t dim12, void *mx2, size_t dim21, size_t dim22, void *result)

Compute mx1 times mx2.

Note

Requires dim12 == dim21, and dim11 * dim22 must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid sizes are logged and leave result unchanged. An empty output causes no memory access. A zero inner dimension produces zeros. Pointers to empty arrays may be null.

Parameters:

mx1 – First input matrix in contiguous row-major storage.
dim11 – Number of rows of mx1.
dim12 – Number of columns of mx1.
mx2 – Second input matrix in contiguous row-major storage.
dim21 – Number of rows of mx2.
dim22 – Number of columns of mx2.
result – Row-major output matrix, which may overlap either input.

void mtMultM(void *mx1, size_t dim11, size_t dim12, void *mx2, size_t dim21, size_t dim22, void *result)

Compute mx1 transpose times mx2.

Note

Requires dim11 == dim21, and dim12 * dim22 must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid sizes are logged and leave result unchanged. An empty output causes no memory access. A zero inner dimension produces zeros. Pointers to empty arrays may be null.

Parameters:

mx1 – First input matrix in contiguous row-major storage.
dim11 – Number of rows of mx1.
dim12 – Number of columns of mx1.
mx2 – Second input matrix in contiguous row-major storage.
dim21 – Number of rows of mx2.
dim22 – Number of columns of mx2.
result – Row-major output matrix, which may overlap either input.

void mMultMt(void *mx1, size_t dim11, size_t dim12, void *mx2, size_t dim21, size_t dim22, void *result)

Compute mx1 times mx2 transpose.

Note

Requires dim12 == dim22, and dim11 * dim21 must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid sizes are logged and leave result unchanged. An empty output causes no memory access. A zero inner dimension produces zeros. Pointers to empty arrays may be null.

Parameters:

mx1 – First input matrix in contiguous row-major storage.
dim11 – Number of rows of mx1.
dim12 – Number of columns of mx1.
mx2 – Second input matrix in contiguous row-major storage.
dim21 – Number of rows of mx2.
dim22 – Number of columns of mx2.
result – Row-major output matrix, which may overlap either input.

void mtMultMt(void *mx1, size_t dim11, size_t dim12, void *mx2, size_t dim21, size_t dim22, void *result)

Compute mx1 transpose times mx2 transpose.

Note

Requires dim11 == dim22, and dim12 * dim21 must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid sizes are logged and leave result unchanged. An empty output causes no memory access. A zero inner dimension produces zeros. Pointers to empty arrays may be null.

Parameters:

mx1 – First input matrix in contiguous row-major storage.
dim11 – Number of rows of mx1.
dim12 – Number of columns of mx1.
mx2 – Second input matrix in contiguous row-major storage.
dim21 – Number of rows of mx2.
dim22 – Number of columns of mx2.
result – Row-major output matrix, which may overlap either input.

void mMultV(void *mx, size_t dim1, size_t dim2, void *v, void *result)

Compute mx times v.

Note

dim1 must not exceed LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid sizes are logged and leave result unchanged. An empty output causes no memory access. A zero inner dimension produces zeros. Pointers to empty arrays may be null.

Parameters:

v – Input vector with dim2 elements.
mx – Input matrix in contiguous row-major storage.
dim1 – Number of matrix rows.
dim2 – Number of matrix columns.
result – Output vector with dim1 elements; may overlap either input.

void mtMultV(void *mx, size_t dim1, size_t dim2, void *v, void *result)

Compute mx transpose times v.

Note

dim2 must not exceed LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid sizes are logged and leave result unchanged. An empty output causes no memory access. A zero inner dimension produces zeros. Pointers to empty arrays may be null.

Parameters:

v – Input vector with dim1 elements.
mx – Input matrix in contiguous row-major storage.
dim1 – Number of matrix rows.
dim2 – Number of matrix columns.
result – Output vector with dim2 elements; may overlap either input.

double mTrace(void *mx, size_t dim)

Compute the trace of a square matrix.

Note

A zero dimension returns zero without accessing mx.

Parameters:

mx – Input matrix.
dim – Number of rows and columns.

Returns:

Sum of diagonal entries.

double mDeterminant(void *mx, size_t dim)

Compute the determinant of a square row-major matrix using pivoted LU.

Note

The (dim - 1) squared minor must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. For dim == 0, returns zero without accessing mx, preserving the existing convention. Scalar and 2-by-2 inputs use direct formulas. No singularity tolerance is applied. The LU pivot product uses a separate binary exponent to avoid intermediate range loss. Elimination arithmetic and a final determinant outside the double range can still lose range.

Parameters:

mx – Input matrix.
dim – Number of rows and columns.

Returns:

Determinant, or NaN after logging invalid dimensions.

void mCofactor(void *mx, size_t dim, void *result)

Compute the cofactor matrix, allowing overlapping input and output.

Note

dim * dim must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid dimensions are logged and leave result unchanged. For dim == 0, no memory is accessed and pointers may be null. Minors are evaluated independently, including for singular input matrices.

Parameters:

mx – Input square row-major matrix.
dim – Number of rows and columns.
result – Output square row-major matrix.

int mInverse(void *mx, size_t dim, void *result)

Invert a square row-major matrix using pivoted LU, allowing overlapping input and output.

Note

dim * dim must fit LINEAR_ALGEBRA_MAX_ARRAY_SIZE. Invalid dimensions, including dim == 0, are logged and leave result unchanged without accessing either pointer.

Note

The absolute determinant must exceed DB0_EPS; this existing cutoff is not scale invariant. A pivot is also rejected when its magnitude is at most DBL_EPSILON * dim times its row’s rounding scale. For P * mx = L * U, that scale is the sum of absolute entries in the original permuted row, plus sum(k < i, abs(L[i,k]) * sum(j >= k, abs(U[k,j]))) for pivot U[i,i]. This includes elimination growth while preserving graded diagonal matrices.

Note

A final condition check uses B[i,j] = mx[i,j] / r[i] / c[j], where r[i] is the largest absolute entry in row i and c[j] is the largest absolute entry in column j after row scaling. With X the computed inverse, the scaled inverse estimate is Y[i,j] = c[i] * X[i,j] * r[j]. The reciprocal condition estimate 1 / (norm_inf(B) * norm_inf(Y)) must exceed dim * DBL_EPSILON, where norm_inf is the maximum absolute row sum. Scaling is used only for this check; it does not change accepted inverse values. This policy can reject mathematically invertible matrices near numerical singularity and is not an exact rank test.

Note

Failed numerical checks, nonfinite input or inverse entries, and nonfinite condition-check norms are logged and zero result. Overflow of the determinant alone does not imply failure if the remaining checks pass.

Parameters:

mx – Input matrix.
dim – Number of rows and columns.
result – Output inverse; zeroed on numerical failure.

Returns:

Zero on success, one for numerical failure or invalid dimensions.

int mIsEqual(void *mx1, size_t dim1, size_t dim2, void *mx2, double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Note

An empty array compares equal for a valid tolerance, without accessing its storage.

Parameters:

mx1 – Input array.
dim1 – Number of rows.
dim2 – Number of columns.
mx2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

int mIsZero(void *mx, size_t dim1, size_t dim2, double accuracy)

Compare every element with zero using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Note

An empty array compares equal for a valid tolerance, without accessing its storage.

Parameters:

mx – Input array.
dim1 – Number of rows.
dim2 – Number of columns.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

void mPrintScreen(const char *name, void *mx, size_t dim1, size_t dim2)

Print the array name and values to standard output.

Parameters:

name – Display label, passed as data rather than as a format string.
mx – Input array.
dim1 – Number of rows.
dim2 – Number of columns.

void mPrint(FILE *pFile, const char *name, void *mx, size_t dim1, size_t dim2)

Print the array name and values to a stream.

Parameters:

pFile – Destination stream.
name – Display label, passed as data rather than as a format string.
mx – Input array.
dim1 – Number of rows.
dim2 – Number of columns.

void mGetSubMatrix(void *mx, size_t dim1, size_t dim2, size_t dim1Start, size_t dim2Start, size_t dim1Result, size_t dim2Result, void *result)

Extract a packed submatrix, allowing overlapping input and output.

Note

Invalid bounds or unrepresentable sizes are logged and leave result unchanged. An empty block causes no memory access; pointers may be null. Overlap is handled by directional row copies, without heap memory or a fixed scratch-buffer size limit.

Parameters:

mx – Input matrix in contiguous row-major storage.
dim1 – Number of input rows.
dim2 – Number of input columns.
dim1Start – First input row of the selected block.
dim2Start – First input column of the selected block.
dim1Result – Number of output rows.
dim2Result – Number of output columns.
result – Output matrix in contiguous row-major storage.

void mSetSubMatrix(void *mx, size_t dim1, size_t dim2, void *result, size_t dim1Result, size_t dim2Result, size_t dim1Start, size_t dim2Start)

Insert a packed submatrix, allowing overlapping input and output.

Note

Invalid bounds or unrepresentable sizes are logged and leave result unchanged. An empty block causes no memory access; pointers may be null. Overlap is handled by directional row copies, without heap memory or a fixed scratch-buffer size limit.

Parameters:

mx – Input matrix in contiguous row-major storage.
dim1 – Number of input rows.
dim2 – Number of input columns.
dim1Start – First output row of the selected block.
dim2Start – First output column of the selected block.
dim1Result – Number of output rows.
dim2Result – Number of output columns.
result – Output matrix in contiguous row-major storage; entries outside the block are preserved.

void m22Set(double m00, double m01, double m10, double m11, double m[2][2])

Set matrix entries in row-major order.

Parameters:

m00 – Value for row 0, column 0.
m01 – Value for row 0, column 1.
m10 – Value for row 1, column 0.
m11 – Value for row 1, column 1.
m – Output array.

void m22Copy(double mx[2][2], double result[2][2])

Copy a 2x2 matrix.

Parameters:

mx – Input matrix.
result – Output matrix, which may overlap mx.

void m22SetZero(double result[2][2])

Fill the output with zeros.

Parameters:: result – Output array.

void m22SetIdentity(double result[2][2])

Set the diagonal to one and every other entry to zero.

Parameters:: result – Output array.

void m22Transpose(double mx[2][2], double result[2][2])

Transpose a square matrix.

Parameters:

mx – Input matrix.
result – Output matrix; may overlap mx.

void m22Add(double mx1[2][2], double mx2[2][2], double result[2][2])

Add corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

mx1 – First input.
mx2 – Second input.
result – Output of the same shape; may overlap either input.

void m22Subtract(double mx1[2][2], double mx2[2][2], double result[2][2])

Subtract corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

mx1 – First input.
mx2 – Second input.
result – Output of the same shape; may overlap either input.

void m22Scale(double scaleFactor, double mx[2][2], double result[2][2])

Scale every element.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

scaleFactor – Multiplier.
mx – Input array.
result – Output of the same shape; may overlap the input.

void m22MultM22(double mx1[2][2], double mx2[2][2], double result[2][2])

Compute mx1 times mx2.

Parameters:

mx1 – Input matrix.
mx2 – Input matrix.
result – Output array; may overlap either input.

void m22tMultM22(double mx1[2][2], double mx2[2][2], double result[2][2])

Compute mx1 transpose times mx2.

Parameters:

mx1 – Input matrix.
mx2 – Input matrix.
result – Output array; may overlap either input.

void m22MultM22t(double mx1[2][2], double mx2[2][2], double result[2][2])

Compute mx1 times mx2 transpose.

Parameters:

mx1 – Input matrix.
mx2 – Input matrix.
result – Output array; may overlap either input.

void m22MultV2(double mx[2][2], double v[2], double result[2])

Compute mx times v.

Parameters:

mx – Input matrix.
v – Input vector.
result – Output array; may overlap either input.

void m22tMultV2(double mx[2][2], double v[2], double result[2])

Compute mx transpose times v.

Parameters:

mx – Input matrix.
v – Input vector.
result – Output array; may overlap either input.

double m22Trace(double mx[2][2])

Compute the trace of a square matrix.

Parameters:: mx – Input matrix.
Returns:: Sum of diagonal entries.

double m22Determinant(double mx[2][2])

Compute the determinant of a fixed-size square matrix.

Parameters:: mx – Input matrix.
Returns:: Determinant, without a singularity tolerance.

int m22IsEqual(double mx1[2][2], double mx2[2][2], double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

mx1 – Input array.
mx2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

int m22IsZero(double mx[2][2], double accuracy)

Compare every element with zero using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

mx – Input array.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

void m22Print(FILE *pFile, const char *name, double mx[2][2])

Print the array name and values to a stream.

Parameters:

pFile – Destination stream.
name – Display label, passed as data rather than as a format string.
mx – Input array.

int m22Inverse(double mx[2][2], double result[2][2])

Invert a fixed-size square matrix, allowing overlapping input and output.

Note

The absolute determinant must exceed both DB0_EPS and a roundoff bound: DBL_EPSILON times the matrix dimension squared times the sum of absolute determinant expansion terms. This rejects cancellation at roundoff while retaining graded diagonal matrices. Numerically singular matrices and nonfinite input, determinant, roundoff bound, or inverse entries are logged as failures. The direct inverse formula is retained.

Parameters:

mx – Input matrix.
result – Output inverse; zeroed on failure and allowed to overlap mx.

Returns:

Zero on success, one on failure.

void m22PrintScreen(const char *name, double mx[2][2])

Print the array name and values to standard output.

Parameters:

name – Display label, passed as data rather than as a format string.
mx – Input array.

void m33Set(double m00, double m01, double m02, double m10, double m11, double m12, double m20, double m21, double m22, double m[3][3])

Set matrix entries in row-major order.

Parameters:

m00 – Value for row 0, column 0.
m01 – Value for row 0, column 1.
m02 – Value for row 0, column 2.
m10 – Value for row 1, column 0.
m11 – Value for row 1, column 1.
m12 – Value for row 1, column 2.
m20 – Value for row 2, column 0.
m21 – Value for row 2, column 1.
m22 – Value for row 2, column 2.
m – Output array.

void m33Copy(double mx[3][3], double result[3][3])

Copy a 3x3 matrix.

Parameters:

mx – Input matrix.
result – Output matrix, which may overlap mx.

void m33SetZero(double result[3][3])

Fill the output with zeros.

Parameters:: result – Output array.

void m33SetIdentity(double result[3][3])

Set the diagonal to one and every other entry to zero.

Parameters:: result – Output array.

void m33Transpose(double mx[3][3], double result[3][3])

Transpose a square matrix.

Parameters:

mx – Input matrix.
result – Output matrix; may overlap mx.

void m33Add(double mx1[3][3], double mx2[3][3], double result[3][3])

Add corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

mx1 – First input.
mx2 – Second input.
result – Output of the same shape; may overlap either input.

void m33Subtract(double mx1[3][3], double mx2[3][3], double result[3][3])

Subtract corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

mx1 – First input.
mx2 – Second input.
result – Output of the same shape; may overlap either input.

void m33Scale(double scaleFactor, double mx[3][3], double result[3][3])

Scale every element.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

scaleFactor – Multiplier.
mx – Input array.
result – Output of the same shape; may overlap the input.

void m33MultM33(double mx1[3][3], double mx2[3][3], double result[3][3])

Compute mx1 times mx2.

Parameters:

mx1 – Input matrix.
mx2 – Input matrix.
result – Output array; may overlap either input.

void m33tMultM33(double mx1[3][3], double mx2[3][3], double result[3][3])

Compute mx1 transpose times mx2.

Parameters:

mx1 – Input matrix.
mx2 – Input matrix.
result – Output array; may overlap either input.

void m33MultM33t(double mx1[3][3], double mx2[3][3], double result[3][3])

Compute mx1 times mx2 transpose.

Parameters:

mx1 – Input matrix.
mx2 – Input matrix.
result – Output array; may overlap either input.

void m33MultV3(double mx[3][3], double v[3], double result[3])

Compute mx times v.

Parameters:

mx – Input matrix.
v – Input vector.
result – Output array; may overlap either input.

void m33tMultV3(double mx[3][3], double v[3], double result[3])

Compute mx transpose times v.

Parameters:

mx – Input matrix.
v – Input vector.
result – Output array; may overlap either input.

double m33Trace(double mx[3][3])

Compute the trace of a square matrix.

Parameters:: mx – Input matrix.
Returns:: Sum of diagonal entries.

double m33Determinant(double mx[3][3])

Compute the determinant of a fixed-size square matrix.

Parameters:: mx – Input matrix.
Returns:: Determinant, without a singularity tolerance.

int m33IsEqual(double mx1[3][3], double mx2[3][3], double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

mx1 – Input array.
mx2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

int m33IsZero(double mx[3][3], double accuracy)

Compare every element with zero using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

mx – Input array.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

void m33Print(FILE *pfile, const char *name, double mx[3][3])

Print the array name and values to a stream.

Parameters:

pfile – Destination stream.
name – Display label, passed as data rather than as a format string.
mx – Input array.

int m33Inverse(double mx[3][3], double result[3][3])

Invert a fixed-size square matrix, allowing overlapping input and output.

Note

The absolute determinant must exceed both DB0_EPS and a roundoff bound: DBL_EPSILON times the matrix dimension squared times the sum of absolute determinant expansion terms. This rejects cancellation at roundoff while retaining graded diagonal matrices. Numerically singular matrices and nonfinite input, determinant, roundoff bound, or inverse entries are logged as failures. The direct inverse formula is retained.

Parameters:

mx – Input matrix.
result – Output inverse; zeroed on failure and allowed to overlap mx.

Returns:

Zero on success, one on failure.

void m33SingularValues(double mx[3][3], double result[3])

Compute the three singular values in descending order.

Note

Non-finite input produces three NaNs. The output may alias a row of mx.

Parameters:

mx – Input matrix, preserved unless result overlaps its storage.
result – Three nonnegative singular values, in the same units as mx.

void m33EigenValues(double mx[3][3], double result[3])

Compute the real parts of the three eigenvalues in descending order.

Note

Exactly symmetric matrices use Jacobi rotations to retain accuracy near repeated eigenvalues. Other matrices use their characteristic polynomial. A complex-conjugate pair contributes its repeated real part; imaginary parts are not returned. Symmetry is checked exactly, without a tolerance.

Parameters:

mx – Finite input 3x3 matrix, preserved unless result overlaps its storage.
result – Three real parts, in the same units as mx; may overlap mx.

double m33ConditionNumber(double mx[3][3])

Compute the condition number from the largest and smallest singular values.

Note

A zero smallest singular value gives infinity when the largest is nonzero. The zero matrix and non-finite input give NaN.

Parameters:: mx – Input matrix, preserved by this function.
Returns:: The largest singular value divided by the smallest (dimensionless).

void m33PrintScreen(const char *name, double mx[3][3])

Print the array name and values to standard output.

Parameters:

name – Display label, passed as data rather than as a format string.
mx – Input array.

void m44Set(double m00, double m01, double m02, double m03, double m10, double m11, double m12, double m13, double m20, double m21, double m22, double m23, double m30, double m31, double m32, double m33, double m[4][4])

Set matrix entries in row-major order.

Parameters:

m00 – Value for row 0, column 0.
m01 – Value for row 0, column 1.
m02 – Value for row 0, column 2.
m03 – Value for row 0, column 3.
m10 – Value for row 1, column 0.
m11 – Value for row 1, column 1.
m12 – Value for row 1, column 2.
m13 – Value for row 1, column 3.
m20 – Value for row 2, column 0.
m21 – Value for row 2, column 1.
m22 – Value for row 2, column 2.
m23 – Value for row 2, column 3.
m30 – Value for row 3, column 0.
m31 – Value for row 3, column 1.
m32 – Value for row 3, column 2.
m33 – Value for row 3, column 3.
m – Output array.

void m44Copy(double mx[4][4], double result[4][4])

Copy a 4x4 matrix.

Parameters:

mx – Input matrix.
result – Output matrix, which may overlap mx.

void m44SetZero(double result[4][4])

Fill the output with zeros.

Parameters:: result – Output array.

void m44MultV4(double mx[4][4], double v[4], double result[4])

Compute mx times v.

Parameters:

mx – Input matrix.
v – Input vector.
result – Output array; may overlap either input.

double m44Determinant(double mx[4][4])

Compute the determinant of a fixed-size square matrix.

Parameters:: mx – Input matrix.
Returns:: Determinant, without a singularity tolerance.

int m44IsEqual(double mx1[4][4], double mx2[4][4], double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

mx1 – Input array.
mx2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

int m44Inverse(double mx[4][4], double result[4][4])

Invert a fixed-size square matrix, allowing overlapping input and output.

Note

The absolute determinant must exceed both DB0_EPS and a roundoff bound: DBL_EPSILON times the matrix dimension squared times the sum of absolute determinant expansion terms. This rejects cancellation at roundoff while retaining graded diagonal matrices. Numerically singular matrices and nonfinite input, determinant, roundoff bound, or inverse entries are logged as failures. The direct inverse formula is retained.

Parameters:

mx – Input matrix.
result – Output inverse; zeroed on failure and allowed to overlap mx.

Returns:

Zero on success, one on failure.

void m66Set(double m00, double m01, double m02, double m03, double m04, double m05, double m10, double m11, double m12, double m13, double m14, double m15, double m20, double m21, double m22, double m23, double m24, double m25, double m30, double m31, double m32, double m33, double m34, double m35, double m40, double m41, double m42, double m43, double m44, double m45, double m50, double m51, double m52, double m53, double m54, double m55, double m[6][6])

Set matrix entries in row-major order.

Parameters:

m00 – Value for row 0, column 0.
m01 – Value for row 0, column 1.
m02 – Value for row 0, column 2.
m03 – Value for row 0, column 3.
m04 – Value for row 0, column 4.
m05 – Value for row 0, column 5.
m10 – Value for row 1, column 0.
m11 – Value for row 1, column 1.
m12 – Value for row 1, column 2.
m13 – Value for row 1, column 3.
m14 – Value for row 1, column 4.
m15 – Value for row 1, column 5.
m20 – Value for row 2, column 0.
m21 – Value for row 2, column 1.
m22 – Value for row 2, column 2.
m23 – Value for row 2, column 3.
m24 – Value for row 2, column 4.
m25 – Value for row 2, column 5.
m30 – Value for row 3, column 0.
m31 – Value for row 3, column 1.
m32 – Value for row 3, column 2.
m33 – Value for row 3, column 3.
m34 – Value for row 3, column 4.
m35 – Value for row 3, column 5.
m40 – Value for row 4, column 0.
m41 – Value for row 4, column 1.
m42 – Value for row 4, column 2.
m43 – Value for row 4, column 3.
m44 – Value for row 4, column 4.
m45 – Value for row 4, column 5.
m50 – Value for row 5, column 0.
m51 – Value for row 5, column 1.
m52 – Value for row 5, column 2.
m53 – Value for row 5, column 3.
m54 – Value for row 5, column 4.
m55 – Value for row 5, column 5.
m – Output array.

void m66Copy(double mx[6][6], double result[6][6])

Copy a 6x6 matrix.

Parameters:

mx – Input matrix.
result – Output matrix, which may overlap mx.

void m66SetZero(double result[6][6])

Fill the output with zeros.

Parameters:: result – Output array.

void m66SetIdentity(double result[6][6])

Set the diagonal to one and every other entry to zero.

Parameters:: result – Output array.

void m66Transpose(double mx[6][6], double result[6][6])

Transpose a square matrix.

Parameters:

mx – Input matrix.
result – Output matrix; may overlap mx.

void m66Get33Matrix(size_t row, size_t col, double m[6][6], double mij[3][3])

Extract a 3x3 block from a 6x6 matrix.

Note

Invalid indices are logged; no input or output memory is accessed.

Parameters:

row – Block row index, either 0 or 1.
col – Block column index, either 0 or 1.
m – Input 6x6 matrix.
mij – Output packed 3x3 matrix; may overlap m.

void m66Set33Matrix(size_t row, size_t col, double mij[3][3], double m[6][6])

Insert a 3x3 block into a 6x6 matrix, preserving entries outside the selected block.

Note

Invalid indices are logged; no input or output memory is accessed.

Parameters:

row – Block row index, either 0 or 1.
col – Block column index, either 0 or 1.
mij – Input packed 3x3 matrix; may overlap m.
m – Output 6x6 matrix.

void m66Scale(double scaleFactor, double mx[6][6], double result[6][6])

Scale every element.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

scaleFactor – Multiplier.
mx – Input array.
result – Output of the same shape; may overlap the input.

void m66Add(double mx1[6][6], double mx2[6][6], double result[6][6])

Add corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

mx1 – First input.
mx2 – Second input.
result – Output of the same shape; may overlap either input.

void m66Subtract(double mx1[6][6], double mx2[6][6], double result[6][6])

Subtract corresponding elements.

Note

Uses fixed local storage; no heap allocation is required.

Parameters:

mx1 – First input.
mx2 – Second input.
result – Output of the same shape; may overlap either input.

void m66MultM66(double mx1[6][6], double mx2[6][6], double result[6][6])

Compute mx1 times mx2.

Parameters:

mx1 – Input matrix.
mx2 – Input matrix.
result – Output array; may overlap either input.

void m66tMultM66(double mx1[6][6], double mx2[6][6], double result[6][6])

Compute mx1 transpose times mx2.

Parameters:

mx1 – Input matrix.
mx2 – Input matrix.
result – Output array; may overlap either input.

void m66MultM66t(double mx1[6][6], double mx2[6][6], double result[6][6])

Compute mx1 times mx2 transpose.

Parameters:

mx1 – Input matrix.
mx2 – Input matrix.
result – Output array; may overlap either input.

void m66MultV6(double mx[6][6], double v[6], double result[6])

Compute mx times v.

Parameters:

mx – Input matrix.
v – Input vector.
result – Output array; may overlap either input.

void m66tMultV6(double mx[6][6], double v[6], double result[6])

Compute mx transpose times v.

Parameters:

mx – Input matrix.
v – Input vector.
result – Output array; may overlap either input.

int m66IsEqual(double mx1[6][6], double mx2[6][6], double accuracy)

Compare corresponding elements using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

mx1 – Input array.
mx2 – Second input.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

int m66IsZero(double mx[6][6], double accuracy)

Compare every element with zero using an absolute tolerance.

Note

NaN and infinity are rejected in the input or tolerance.

Parameters:

mx – Input array.
accuracy – Finite, nonnegative absolute tolerance, in the same units as the input.

Returns:

One if every absolute difference is at most accuracy, otherwise zero.

void m99SetZero(double result[9][9])

Fill the output with zeros.

Parameters:: result – Output array.

void cubicRoots(double a[3], double result[3])

Compute real parts of the roots of x^3 + a[2]*x^2 + a[1]*x + a[0] = 0.

Note

For one real root and a complex-conjugate pair, result[0] is the real root and result[1] and result[2] are the pair’s real part. Imaginary parts are not returned. All three outputs are roots when all roots are real, including repeated roots.

Parameters:

a – Three finite coefficients in ascending order, excluding the monic leading coefficient.
result – Three real parts; may overlap a. The output is not sorted.

double safeAcos(double x)

Compute an inverse cosine with its argument clamped to [-1, 1].

Parameters:: x – Input value.
Returns:: Angle in radians, or NaN for NaN input.

double safeAsin(double x)

Compute an inverse sine with its argument clamped to [-1, 1].

Parameters:: x – Input value.
Returns:: Angle in radians, or NaN for NaN input.

double safeSqrt(double x)

Take a square root after clamping negative inputs to zero.

Note

Finite negative inputs and quiet NaNs do not raise the invalid floating-point exception.

Parameters:: x – Input value.
Returns:: Zero for negative input, otherwise sqrt(x), preserving NaN and signed zero.