fix(billing): bill the factorization numpy runs for method= and hermitian=True

[AndreasHad04]

Fixes #280.

What was wrong

Two options make numpy run a cheaper factorization than an SVD, and both were billed as the full SVD. Generator.multivariate_normal(..., method=) factors the covariance with an SVD, eigh or Cholesky, but multivariate_normal_flops always charged svd_cost(d, d, with_vectors=True): for the issue's 256x256 covariance every method cost 872,685,568, while the Cholesky factorization numpy actually runs bills 11,184,810 standalone. linalg.svd, linalg.pinv and linalg.matrix_rank with hermitian=True only forwarded the flag; numpy then runs eigh (or eigvalsh without vectors) instead of an SVD, yet the bill was identical to hermitian=False (872,415,232 for svd against 301,989,888 for the eigh it runs).

What it changes

multivariate_normal_flops takes the method, and its factorization term is what that factorization bills standalone: svd_cost(d, d, with_vectors=True) = 26d^3 for 'svd' (the default), eigh_cost(d) = 9d^3 for 'eigh' and cholesky_cost(d) = d^3/3 for 'cholesky'; the sampling term is unchanged. The Generator formula reads the keyword-only method=; RandomState and the module-level function, which always use SVD, bill as before. The issue's three methods now bill 872,685,568, 302,260,224 and 11,455,146.

For hermitian=True, the new hermitian_svd_cost(n, with_vectors=) in _decompositions.py bills the branch numpy runs: eigh (or eigvalsh) plus ordering the singular values by magnitude (sign, abs, an argsort and three take_along_axis gathers at weight 4, the sign folded into the vectors and the conjugate copy that forms Vh; values only, abs and a sort), each step at its op's own weight the way composite ops fold sub-tiers. svd bills it for a square matrix (302,797,824 in the issue, just above the 301,989,888 of eigh alone), keeping the truncated cost of a top-k request when that is lower, and pinv and matrix_rank use it for their SVD term (369,907,200 and 44,756,650 in the issue). hermitian=False, and non-square input, which numpy refuses under hermitian=True, bill as before. The registry notes for linalg.svd, linalg.pinv and Generator.multivariate_normal, ops.json (regenerated, three notes fields) and the cost-model rows describe the new terms.

Tests

New tests/test_mvn_method_hermitian_billing.py (22 tests): at production weights the issue's three methods bill exactly the standalone factorization plus the shared sampling term; under unit weights each method matches its closed form for two d/size combinations; draws equal numpy's for every method; RandomState and module-level sampling are unchanged. For hermitian=True, svd with and without vectors bills exactly the same as numpy's hermitian branch spelled out op by op with fnp, for n = 1, 2, 5, 64 and a batch, and returns numpy's values; the issue's rows land just above the eigh bill and keep pinv's and matrix_rank's other terms; closed forms match for svd, pinv and matrix_rank, with numpy's results; top-k requests bill the lower of the truncated cost and the hermitian path; hermitian=False is unchanged and non-square input still raises LinAlgError. 14 of the 22 fail without the source change.

Results

• New tests: 22 passed.
• Full default suite (uv run pytest): 10387 passed, 125 skipped, 2 xfailed.
• tests/parity fast gate: 136 passed.
• ruff check, ruff format --check, pyright on the changed files: clean.

Note: this branch and #291 both change matrix_rank in _properties.py and its cost-model row, so whichever merges second needs a small rebase.