## A.5 Mathematical Routines

This section describes a number of useful mathematical functions and classes that support basic operations in pbrt, such as solving small linear systems, manipulating matrices, and linear interpolation.

The Lerp() function linearly interpolates between the two provided values.

<<Global Inline Functions>>+=
inline Float Lerp(Float t, Float v1, Float v2) { return (1 - t) * v1 + t * v2; }

The Quadratic() function finds solutions of the quadratic equation ; the Boolean return value indicates whether solutions were found.

<<Global Inline Functions>>+=
inline bool Quadratic(Float a, Float b, Float c, Float *t0, Float *t1) { <<Find quadratic discriminant>>
double discrim = (double)b * (double)b - 4 * (double)a * (double)c; if (discrim < 0) return false; double rootDiscrim = std::sqrt(discrim);
double q; if (b < 0) q = -.5 * (b - rootDiscrim); else q = -.5 * (b + rootDiscrim); *t0 = q / a; *t1 = c / q; if (*t0 > *t1) std::swap(*t0, *t1); return true;
}

The implementation always uses double-precision floating-point values regardless of the type of Float in order to return a result with minimal floating-point error. If the discriminant () is negative, then there are no real roots and the function returns false.

double discrim = (double)b * (double)b - 4 * (double)a * (double)c; if (discrim < 0) return false; double rootDiscrim = std::sqrt(discrim);

The usual version of the quadratic equation can give poor numerical precision when due to cancellation error. It can be rewritten algebraically to a more stable form:

where

double q; if (b < 0) q = -.5 * (b - rootDiscrim); else q = -.5 * (b + rootDiscrim); *t0 = q / a; *t1 = c / q; if (*t0 > *t1) std::swap(*t0, *t1); return true;

### A.5.2 2 times 2 Linear Systems

There are a number of places throughout pbrt where we need to solve a linear system of the form

for values and . The SolveLinearSystem2x2() routine finds the closed-form solution to such a system. It returns true if it was successful and false if the determinant of is very small, indicating that the system is numerically ill-conditioned and either not solvable or likely to have unacceptable floating-point errors. In this case, no solution is returned.

<<Matrix4x4 Method Definitions>>=
bool SolveLinearSystem2x2(const Float A, const Float B, Float *x0, Float *x1) { Float det = A * A - A * A; if (std::abs(det) < 1e-10f) return false; *x0 = (A * B - A * B) / det; *x1 = (A * B - A * B) / det; if (std::isnan(*x0) || std::isnan(*x1)) return false; return true; }

### A.5.3 4 times 4 Matrices

The Matrix4x4 structure provides a low-level representation of matrices. It is an integral part of the Transform class.

<<Matrix4x4 Declarations>>=
struct Matrix4x4 { <<Matrix4x4 Public Methods>>
Matrix4x4() { m = m = m = m = 1.f; m = m = m = m = m = m = m = m = m = m = m = m = 0.f; } Matrix4x4(Float mat); Matrix4x4(Float t00, Float t01, Float t02, Float t03, Float t10, Float t11, Float t12, Float t13, Float t20, Float t21, Float t22, Float t23, Float t30, Float t31, Float t32, Float t33); bool operator==(const Matrix4x4 &m2) const { for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) if (m[i][j] != m2.m[i][j]) return false; return true; } bool operator!=(const Matrix4x4 &m2) const { for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) if (m[i][j] != m2.m[i][j]) return true; return false; } friend Matrix4x4 Transpose(const Matrix4x4 &); void Print(FILE *f) const { fprintf(f, "[ "); for (int i = 0; i < 4; ++i) { fprintf(f, " [ "); for (int j = 0; j < 4; ++j) { fprintf(f, "%f", m[i][j]); if (j != 3) fprintf(f, ", "); } fprintf(f, " ]\n"); } fprintf(f, " ] "); } static Matrix4x4 Mul(const Matrix4x4 &m1, const Matrix4x4 &m2) { Matrix4x4 r; for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) r.m[i][j] = m1.m[i] * m2.m[j] + m1.m[i] * m2.m[j] + m1.m[i] * m2.m[j] + m1.m[i] * m2.m[j]; return r; } friend Matrix4x4 Inverse(const Matrix4x4 &);
Float m; };

The default constructor, not shown here, sets the matrix to the identity matrix. The Matrix4x4 implementation also provides constructors that allow the user to pass an array of floats or 16 individual floats to initialize a Matrix4x4:

<<Matrix4x4 Public Methods>>=
Matrix4x4(Float mat); Matrix4x4(Float t00, Float t01, Float t02, Float t03, Float t10, Float t11, Float t12, Float t13, Float t20, Float t21, Float t22, Float t23, Float t30, Float t31, Float t32, Float t33);

The implementations of operators that test for equality and inequality are straightforward and not included in the text here.

The Matrix4x4 class supports a few low-level matrix operations. For example, Transpose() returns a new matrix that is the transpose of the original matrix.

<<Matrix4x4 Method Definitions>>+=
Matrix4x4 Transpose(const Matrix4x4 &m) { return Matrix4x4(m.m, m.m, m.m, m.m, m.m, m.m, m.m, m.m, m.m, m.m, m.m, m.m, m.m, m.m, m.m, m.m); }

The product of two matrices and is computed by setting the th element of the result to the inner product of the th row of with the th column of .

<<Matrix4x4 Public Methods>>+=
static Matrix4x4 Mul(const Matrix4x4 &m1, const Matrix4x4 &m2) { Matrix4x4 r; for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) r.m[i][j] = m1.m[i] * m2.m[j] + m1.m[i] * m2.m[j] + m1.m[i] * m2.m[j] + m1.m[i] * m2.m[j]; return r; }

Finally, Inverse() returns the inverse of the matrix. The implementation (not shown here) uses a numerically stable Gauss–Jordan elimination routine to compute the inverse.

<<Matrix4x4 Public Methods>>+=
friend Matrix4x4 Inverse(const Matrix4x4 &);