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.
A.5.1 Solving Quadratic Equations
The Quadratic() function finds solutions of the quadratic equation ; the Boolean return value indicates whether solutions were found.
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.
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:
A.5.2 22 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.
A.5.3 44 Matrices
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:
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.
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 .
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.