[linalg] NB comments

[mhoemmen]

Remove [linalg] from C++26?

FR-029-273 29.9 [linalg] Remove linalg from C++26 says:

Special math functions were added in the standard with the justification they were hard to implement and the expectation was that standard libraries maintainers would be the most qualified to do that work for the benefit of a small number of C++ users. But because standard library maintainers are not domain experts in all domains, and because their time is limited and their resources scarce....

(NVIDIA has an optimized implementation, btw)

Distinguish

• time / effort
• expertise

Everything requires time and effort. Special math functions require mathematical expertise.

Prof. Kahan (of IEEE 754 fame) circa 2004: the world creates one numerical analyst a year.

BLAS was designed as part of a separation of concerns between performance engineers and numerical analysts. I explained this in detail in my CppCon 2023 talk.

Special math function example: fast inverse square root

Depends on knowledge of

1. classical methods for approximating functions of real or complex numbers
2. floating-point representation on computers
3. computer hardware, e.g., relative performance of different instructions

Example of (1): Newton's method, a classical iteration $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$ for approximating a function using its derivative.

float Q_rsqrt(float number) {
  constexpr float threehalfs = 1.5f;

  float x2 = number * 0.5f;
  float y = number;
  std::int32_t i = std::bit_cast<std::int32_t>(y);
  i = 0x5f3759df - (i >> 1); // approx 1/sqrt
  y = std::bit_cast<float>(i);
  y = y * (threehalfs - (x2 * y * y)); // Newton's method, 1st iteration
  y = y * (threehalfs - (x2 * y * y)); // optional 2nd iteration

  return y;
}

Special math function example: Bessel functions of the first kind

$$ J_{\alpha}(x) = \sum_{m=0}^{\infty} \frac{(-1)^m}{m! \Gamma(m + \alpha + 1)} \left( \frac{x}{2} \right)^{2m + \alpha} $$

Solutions of the following differential equation:

$$ x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - \alpha^2) y = 0 $$

that "describe[s] the radial part of vibrations of a circular membrane":

BLAS example

Matrix-matrix multiply: $C = A * B$ where $C[r,c] = A[r,1] B[1,c] + \dots + A[r,K] B[K,c]$.

That's it! Normal plus and times! Elementary school math!