# HW2 Companion — Machine Arithmetic, Matrix Multiplication, Linear Maps & Numerical Integration
> [!info]- How to use this companion
> This is a scaffold, not a solution key. For each problem you get three things:
> - **What to hand in** — a plain breakdown of what is a *written* calculation, what is a *code* calculation, and what to *report*.
> - **Hints** — a ladder of collapsed callouts. Open **Hint 1** first; only open the next if you are still stuck. The last hint gets you close, but the final answer stays yours.
> - **Code guidance** — a callout for each of the four languages (Python, R, MATLAB, Mathematica). Utility code you will reuse (random matrices, the trapezoid and Simpson rules) is given in full; the problem-specific code is a **skeleton with a few lines left for you** — the `# TODO` marks are where the real thinking goes.
>
> Everything you submit lives in one written PDF (calculations, tables, plots, short reflections), with your code files and generated figures turned in alongside it.
---
## Deliverables at a glance
| # | Topic | By hand | In code | Report |
|---|-------|---------|---------|--------|
| 1 | Machine arithmetic | the whole argument | optional spacing check | can it cross? (min / max), with justification |
| 2 | Matrix multiplication | complexity note | `MM`, benchmark vs. built-in | timing + speedup table, commentary |
| 3 | Action of a $2\times2$ matrix | the six line-plots, $\det$, eigenpairs | grid transform, eigen check | plots + a sentence per matrix |
| 4 | Normal quadrature | tie $F(1)$ to $\operatorname{erf}$ | trapezoid & Simpson | error table + points-needed |
| 5 | Normal vs. Cauchy | the two hand integrals | window experiments | what you notice + the money question |
Submit **one PDF** plus your **code files** and **figures**. Label every plot's axes; report numbers to the precision the problem asks.
---
## Building the test matrices (needed for Problem 2)
Problem 2 asks you to time your `MM` against the built-in product on **dense** matrices of three sizes (50, 200, 500) and on **sparse** versions of each. Below is the matrix-building code in all four languages — copy it into your driver script so the problem itself is just *multiply and time*, not *set up data*.
Three ideas to keep straight:
- **Seed the generator** (e.g. `default_rng(0)`, `set.seed(0)`, `rng(0)`, `SeedRandom[0]`) so your runs are reproducible and your two matrices are actually different draws.
- A [sparse matrix](https://en.wikipedia.org/wiki/Sparse_matrix) with $\sim5\%$ nonzeros stores only the nonzero entries and their positions (formats like [CSR/CSC](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_(CSR,_CRS_or_Yale_format))), so a $500\times500$ sparse matrix holds $\sim12{,}500$ numbers instead of $250{,}000$. The built-in sparse product **skips the zeros**; that is the whole reason it can be thousands of times faster.
- The density knob is *approximate* for the built-in random-sparse constructors (duplicate positions get merged), which is fine here — you want "about 5%," not exactly.
> [!example]- Python
>
> ```python
> import numpy as np
> from scipy import sparse
>
> rng = np.random.default_rng(0) # reproducible seed
> sizes = {"small": 50, "medium": 200, "large": 500}
>
> def dense_pair(n):
> """Two independent n x n dense matrices, entries in [0, 1)."""
> return rng.random((n, n)), rng.random((n, n))
>
> def sparse_pair(n, density=0.05):
> """Two n x n sparse matrices with ~5% nonzeros, CSR format."""
> A = sparse.random(n, n, density=density, format="csr", random_state=rng)
> B = sparse.random(n, n, density=density, format="csr", random_state=rng)
> return A, B
> ```
> *Look for:* `A.nnz` on a 500×500 sparse draw is about `0.05*500*500 = 12500`. The built-in dense product is `A @ B`; the sparse product is the same `A @ B` on CSR matrices.
> [!example]- R
>
> ```r
> set.seed(0) # reproducible seed
> library(Matrix) # sparse-matrix support
> sizes <- c(small = 50, medium = 200, large = 500)
>
> dense_pair <- function(n)
> list(A = matrix(runif(n * n), n, n),
> B = matrix(runif(n * n), n, n))
>
> sparse_pair <- function(n, density = 0.05) {
> nnz <- round(density * n * n)
> draw <- function()
> sparseMatrix(i = sample(n, nnz, replace = TRUE),
> j = sample(n, nnz, replace = TRUE),
> x = runif(nnz), dims = c(n, n))
> list(A = draw(), B = draw())
> }
> ```
> *Look for:* the built-in dense product is `A %*% B`; with `Matrix` sparse objects the same `A %*% B` dispatches to the sparse routine. `length(A@x)` reports the stored nonzeros.
> [!example]- MATLAB
>
> ```matlab
> rng(0); % reproducible seed
> sizes = [50 200 500];
>
> % dense: two independent n-by-n matrices, entries in (0,1)
> A = rand(n); B = rand(n);
>
> % sparse: ~5% nonzeros (sprand places them at random positions)
> As = sprand(n, n, 0.05);
> Bs = sprand(n, n, 0.05);
> ```
> *Look for:* `nnz(As)` is about `0.05*n*n`. The product is `A*B` in both cases — MATLAB picks the sparse algorithm automatically when the operands are sparse.
> [!example]- Mathematica *(numerical only)*
>
> ```wolfram
> SeedRandom[0]; (* reproducible seed *)
> sizes = {50, 200, 500};
>
> densePair[n_] := {RandomReal[{0, 1}, {n, n}], RandomReal[{0, 1}, {n, n}]};
>
> sparsePair[n_, density_: 0.05] := Module[{nnz = Round[density n n], draw},
> draw := SparseArray[
> Table[{RandomInteger[{1, n}], RandomInteger[{1, n}]} -> RandomReal[], nnz],
> {n, n}];
> {draw, draw}];
> ```
> *Look for:* `Length[SparseArray[...]["NonzeroValues"]]` is about `0.05 n^2`. Use `Dot[A, B]` for the product; on `SparseArray` operands it stays sparse. Keep everything numerical (no symbolic entries).
---
## Problem 1 — Machine arithmetic (the immortal bacterium)
**What to hand in.**
- **By hand:** the whole argument. This is a reasoning problem about the [floating-point](https://en.wikipedia.org/wiki/Floating-point_arithmetic) number line, not a coding one.
- **In code (optional):** one line to confirm the spacing at a chosen magnitude.
- **Report:** for the **minimum** and the **maximum** Earth–Mars separation, whether the simulated bacterium can traverse it — with the reasoning (the halt distance and the comparison) that gets you there.
> [!hint]- Hint 1 — what "gap" means
> Machine numbers are *not* evenly spaced. They cluster densely near zero and spread out as you move away, because $x=\pm(1+f)\times2^{e}$ keeps a fixed number of significand steps within each power-of-two band $[2^{e},2^{e+1})$. So the question "can one bacterial run change the stored number?" has a different answer near Earth than it does far out toward Mars.
> [!hint]- Hint 2 — the spacing formula
> Inside the band $[2^{e},2^{e+1})$ the distance between one machine number and the next is *constant*: with a 52-bit significand (double precision), that spacing is
> $\text{gap}(e)=2^{e}\cdot 2^{-52}.$
> This step size is the [unit in the last place](https://en.wikipedia.org/wiki/Unit_in_the_last_place), and $2^{-52}$ is [machine epsilon](https://en.wikipedia.org/wiki/Machine_epsilon). Notice the gap **doubles** every time the position doubles.
> [!hint]- Hint 3 — when does a run stop counting?
> A run of length $s$ moves the bacterium to a *new* machine number only if $s$ is at least the local gap. If $s<\text{gap}$, then (stored position) $+\,s$ rounds right back to the stored position — the bacterium runs in place forever. **Convert the run into the number line's units first:** one unit is a kilometer, and an *E. coli* run is $10$–$12$ micrometers, so a run is about $10^{-8}$ km. (Micrometers $\to$ meters $\to$ kilometers.)
> [!hint]- Hint 4 — find where it halts, then compare
> The gap grows with position while the run length stays fixed, so there is a first band where the gap first exceeds one run — that power of two is the farthest the bacterium can ever get, no matter how immortal it is. Solve $\text{gap}(e)\gtrsim s$ for the band, express that boundary in kilometers, and hold it against the two Earth–Mars separations (minimum $\approx 5.5\times10^{7}$ km, maximum $\approx 4.0\times10^{8}$ km). One of them sits below the wall; one sits above it.
> [!example]- Optional — check the spacing in code (all four languages)
>
> Each of these returns the distance to the next representable number above a value `x` — i.e. the gap at that magnitude. Compare it to your $\sim10^{-8}$ km run.
> ```python
> import numpy as np; np.spacing(5.0e7) # Python: gap just above 5e7
> ```
> ```r
> x <- 5.0e7; 2^(floor(log2(x)) - 52) # R: gap in the band containing x
> ```
> ```matlab
> eps(5.0e7) % MATLAB: gap just above 5e7
> ```
> ```wolfram
> x = 5.0*^7; 2^(Floor[Log2[x]] - 52) (* Mathematica: gap in x's band *)
> ```
> *Look for:* the gap at $5\times10^{7}$ is a few $\times10^{-9}$ km — just under a run — while at $4\times10^{8}$ it is several $\times10^{-8}$ km, comfortably larger than a run.
---
## Problem 2 — Matrix multiplication: your loops vs. the built-in
**What to hand in.**
- **By hand:** one or two sentences on why the naive product is $\mathcal{O}(n^{3})$ (count the multiplications in three nested loops).
- **In code:** `MM(A,B)` from explicit loops, with a dimension-mismatch guard; a check that `MM` matches the built-in on a small random pair; then timing of both routines across the three dense sizes and their sparse versions.
- **Report:** a table of times and **speedups** (built-in vs. `MM`) for each size, dense and sparse; a comment on how the gap grows with $n$ and why sparse is so dramatic.
> [!hint]- Hint 1 — the shape of MM
> Three nested loops: for each output row `i` and column `j`, accumulate the dot product $C_{ij}=\sum_k A_{ik}B_{kj}$. Before any of that, **guard the dimensions**: the number of columns of `A` must equal the number of rows of `B`, otherwise raise an error / call `error(...)` with a message that says the two shapes. Fail loudly, never silently.
> [!hint]- Hint 2 — verify before you time
> Multiply a small random pair (say $30\times30$) with both `MM` and the built-in and compare. Do **not** test for exact equality — floating-point reorderings mean the right check is that the largest entrywise difference is tiny, e.g. `max|MM - builtin| < 1e-10`. Only once this passes should you trust your timings.
> [!hint]- Hint 3 — timing honestly
> Wrap each call in the language's timer and record seconds. The **speedup** is $t_{\texttt{MM}}/t_{\text{built-in}}$. Expect it to climb steeply with $n$: the built-in uses cache-aware, [BLAS](https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms)-backed kernels while your triple loop does not. The naive $500\times500$ run may take a while; if it is painfully slow, time it once (or drop to $300$) and say so in your report — the *trend* is the deliverable.
> [!hint]- Hint 4 — the sparse contrast
> Build the sparse operands from the matrix section and multiply them with the **same** built-in operator; it skips the zeros, so it touches $\sim5\%$ of the work. Your naive `MM` does not know the zeros are special and grinds through all $n^3$ products anyway. Put the dense and sparse speedups side by side — that ratio is the story.
> [!example]- Python
>
> **`MM` skeleton** — fill the accumulation and the guard.
> ```python
> import numpy as np
>
> def MM(A, B):
> A = np.asarray(A, float); B = np.asarray(B, float)
> (m, n), (p, q) = A.shape, B.shape
> if n != p:
> raise ValueError(f"inner dims disagree: A is {m}x{n}, B is {p}x{q}")
> C = np.zeros((m, q))
> for i in range(m):
> for j in range(q):
> s = 0.0
> for k in range(n):
> s += 0.0 # TODO: accumulate A[i,k] * B[k,j]
> C[i, j] = s
> return C
> ```
> **Timing harness** (given in full):
> ```python
> import time
> def timed(fn, *args):
> t0 = time.perf_counter(); out = fn(*args)
> return out, time.perf_counter() - t0
>
> A, B = dense_pair(200)
> _, t_mine = timed(MM, A, B)
> _, t_builtin = timed(lambda X, Y: X @ Y, A, B)
> print(t_mine, t_builtin, t_mine / t_builtin) # speedup
> ```
> *Look for:* the verification `max(abs(MM(A,B) - A@B)) < 1e-10` before you report any speedup.
> [!example]- R
>
> **`MM` skeleton:**
> ```r
> MM <- function(A, B) {
> if (ncol(A) != nrow(B))
> stop(sprintf("inner dims disagree: A is %dx%d, B is %dx%d",
> nrow(A), ncol(A), nrow(B), ncol(B)))
> C <- matrix(0, nrow(A), ncol(B))
> for (i in seq_len(nrow(A)))
> for (j in seq_len(ncol(B))) {
> s <- 0
> for (k in seq_len(ncol(A)))
> s <- s + 0 # TODO: accumulate A[i,k] * B[k,j]
> C[i, j] <- s
> }
> C
> }
> ```
> **Timing** (given): `system.time( MM(A, B) )["elapsed"]` and `system.time( A %*% B )["elapsed"]`; speedup is their ratio.
> *Look for:* `max(abs(MM(A,B) - A %*% B)) < 1e-10`.
> [!example]- MATLAB
>
> **`MM` skeleton:**
> ```matlab
> function C = MM(A, B)
> [m, n] = size(A); [p, q] = size(B);
> if n ~= p
> error('inner dims disagree: A is %dx%d, B is %dx%d', m, n, p, q);
> end
> C = zeros(m, q);
> for i = 1:m
> for j = 1:q
> s = 0;
> for k = 1:n
> s = s + 0; % TODO: accumulate A(i,k) * B(k,j)
> end
> C(i, j) = s;
> end
> end
> end
> ```
> **Timing** (given): wrap calls in `tic; MM(A,B); t = toc;` and `tic; A*B; tb = toc;`; speedup `t/tb`.
> *Look for:* `max(abs(MM(A,B) - A*B), [], 'all') < 1e-10`.
> [!example]- Mathematica *(numerical only)*
>
> **`MM` skeleton:**
> ```wolfram
> MM[A_, B_] := Module[{m, n, p, q, C},
> {m, n} = Dimensions[A]; {p, q} = Dimensions[B];
> If[n =!= p, Return[$Failed]]; (* guard: print a dimension message too *)
> C = ConstantArray[0., {m, q}];
> Do[
> C[[i, j]] = Sum[0., {k, n}], (* TODO: accumulate A[[i,k]] B[[k,j]] *)
> {i, m}, {j, q}];
> C];
> ```
> **Timing** (given): `First @ AbsoluteTiming[ MM[A, B]; ]` and `First @ AbsoluteTiming[ A.B; ]`; speedup is their ratio.
> *Look for:* `Max @ Abs[MM[A, B] - A.B] < 10^-10`. Keep entries numerical so `.` uses the fast machine kernel.
---
## Problem 3 — How a $2\times2$ matrix acts on the plane
**What to hand in.**
- **By hand:** for each of the six matrices — (a) the two-line plot of $A\vec{x}=\vec{0}$ and its $\det$; (c) the eigenvalues and eigenvectors (a $2\times2$ eigenproblem is a quadratic you can solve on paper).
- **In code:** (b) transform the unit grid by each $A$ (use your `MM`); and a numerical eigen-decomposition to *check* your hand work.
- **Report:** the six line plots and the six grid-transform plots, plus one sentence per matrix connecting $\det A$ and the eigenstructure to the picture.
> [!hint]- Hint 1 — (a) the two lines and the determinant
> $A\vec{x}=\vec{0}$ is the pair of lines $ax+by=0$ and $cx+dy=0$, both through the origin. Their intersection is exactly the set of solutions. If $\det A\neq0$ the lines cross only at the origin (the unique solution $\vec{x}=\vec{0}$); if $\det A=0$ the two rows are proportional, so the lines coincide — a whole line of solutions. The determinant is your advance warning of which case you are in.
> [!hint]- Hint 2 — (b) the grid transform
> The columns of $A$ are the images of the standard basis vectors: $A\hat{\imath}$ is the first column, $A\hat{\jmath}$ the second. So $T$ sends the unit square to the parallelogram those two image-vectors span, and $|\det A|$ is that parallelogram's area (zero area means the plane collapsed onto a line). To draw it, stack your grid points as columns of a $2\times N$ array and hit them with `MM(A, P)`.
> [!hint]- Hint 3 — (c) eigenpairs by hand
> Solve $\det(A-\lambda I)=0$, a quadratic in $\lambda$; for each root solve $(A-\lambda I)\vec{v}=\vec{0}$ for the direction $\vec{v}$. Real, distinct $\lambda