docs: standardize API contracts and numerical guides

.github/workflows/docs.yml

@@ -36,6 +36,14 @@ jobs:
         run: |
           python -m pip install -r requirements-docs.txt
           python tools/build_docs.py
+      - name: Verify Python documentation examples against the released package
+        run: |
+          python - <<'PY'
+          import pathlib, subprocess, sys, tomllib
+          version = tomllib.loads(pathlib.Path('pyproject.toml').read_text())['project']['version']
+          subprocess.run([sys.executable, '-m', 'pip', 'install', f'hmath-python[numpy,cffi]=={version}'], check=True)
+          PY
+          python tools/check_docs_examples.py
       - uses: actions/upload-pages-artifact@56afc609e74202658d3ffba0e8f6dda462b719fa # v3
         with:
           path: site

.github/workflows/python-wheels.yml

@@ -4,9 +4,9 @@ name: Python wheels
 on:
   push:
     branches: [master]
-    paths-ignore: ['docs/**', 'mkdocs.yml', 'requirements-docs.txt', 'tools/build_docs.py', 'tools/check_api_docs.py', 'tools/check_haskell_doc_examples.py', 'tools/check_docs_links.py', 'tools/haskell_reference.py', 'README.md']
+    paths-ignore: ['.github/workflows/docs.yml', 'docs/**', 'mkdocs.yml', 'requirements-docs.txt', 'tools/build_docs.py', 'tools/check_api_docs.py', 'tools/check_haskell_doc_examples.py', 'tools/check_docs_links.py', 'tools/haskell_reference.py', 'README.md']
   pull_request:
-    paths-ignore: ['docs/**', 'mkdocs.yml', 'requirements-docs.txt', 'tools/build_docs.py', 'tools/check_api_docs.py', 'tools/check_haskell_doc_examples.py', 'tools/check_docs_links.py', 'tools/haskell_reference.py', 'README.md']
+    paths-ignore: ['.github/workflows/docs.yml', 'docs/**', 'mkdocs.yml', 'requirements-docs.txt', 'tools/build_docs.py', 'tools/check_api_docs.py', 'tools/check_haskell_doc_examples.py', 'tools/check_docs_links.py', 'tools/haskell_reference.py', 'README.md']
   workflow_call:
     outputs:
       version:

.github/workflows/validate.yml

@@ -4,9 +4,9 @@ name: Validate
 on:
   push:
     branches: [master]
-    paths-ignore: ['docs/**', 'mkdocs.yml', 'requirements-docs.txt', 'tools/build_docs.py', 'tools/check_api_docs.py', 'tools/check_haskell_doc_examples.py', 'tools/check_docs_links.py', 'tools/haskell_reference.py', 'README.md']
+    paths-ignore: ['.github/workflows/docs.yml', 'docs/**', 'mkdocs.yml', 'requirements-docs.txt', 'tools/build_docs.py', 'tools/check_api_docs.py', 'tools/check_haskell_doc_examples.py', 'tools/check_docs_links.py', 'tools/haskell_reference.py', 'README.md']
   pull_request:
-    paths-ignore: ['docs/**', 'mkdocs.yml', 'requirements-docs.txt', 'tools/build_docs.py', 'tools/check_api_docs.py', 'tools/check_haskell_doc_examples.py', 'tools/check_docs_links.py', 'tools/haskell_reference.py', 'README.md']
+    paths-ignore: ['.github/workflows/docs.yml', 'docs/**', 'mkdocs.yml', 'requirements-docs.txt', 'tools/build_docs.py', 'tools/check_api_docs.py', 'tools/check_haskell_doc_examples.py', 'tools/check_docs_links.py', 'tools/haskell_reference.py', 'README.md']
   workflow_dispatch:
     inputs:
       full-matrix:

docs/contributing.md

@@ -22,7 +22,8 @@ put detailed instructions here and link to the authoritative explanation.
 
 ```sh
 python -m pip install -r requirements-docs.txt
-python -m mkdocs serve
+python tools/build_docs.py
+python -m mkdocs serve --config-file build-docs/mkdocs.generated.yml
 ```
 
 To build the site together with the complete Haskell API reference, install
@@ -36,3 +37,36 @@ CI builds with strict validation and deploys the generated `site/` directory
 to GitHub Pages after changes reach `master`. Do not commit generated HTML.
 Use runnable examples with expected results; distinguish available interfaces
 from future or Haskell-only features.
+
+## Reference writing standard
+
+Follow the [numpydoc section conventions](https://numpydoc.readthedocs.io/en/latest/format.html)
+for Python reference content and [Haddock conventions](https://haskell-haddock.readthedocs.io/latest/markup.html)
+for Haskell source documentation. The website uses Markdown; adopting these
+content conventions does not require changing its site generator.
+
+Each Python, C, or Wave entry needs a one-line purpose, exact signature,
+`Parameters` table with one row per argument in declaration order, a `Returns`
+section, and relevant `See also` links. Include types, defaults, units, shapes,
+bounds, ownership, and failure conditions where they apply. Use `None.` for
+an empty parameter list. Properties have no call parentheses. Explain errors
+with the actual exception/status/constructor, rather than “may fail”.
+
+Add examples for operations that need an explanation of usage. Show expected
+results, resource cleanup, and meaningful failure cases. For numerical
+algorithms, explain the equation, assumptions, stopping criteria, and accuracy
+limitations; keep implementation history out of user instructions. A simple
+getter does not need a paper citation or an artificial standalone example.
+
+Haskell module and area indexes are generated from declarations and authored
+Haddock prose, including constructors, record fields, and class methods.
+The API check includes symbolic function names as well as ordinary identifiers.
+Keep standalone task guides distinct from reference contracts and conceptual
+explanations, following [Diátaxis](https://diataxis.fr/start-here/).
+
+The documentation workflow checks source signatures and parameter order,
+reference coverage, local file/fragment links, Python examples against the
+released wheel, and complete Haskell example programs against this checkout.
+Run the C/Wave example programs with the native build when changing them.
+The home-page version is substituted from `pyproject.toml` during generation;
+documentation changes alone do not require a package release.

docs/explanation/numerical-conventions.md

@@ -0,0 +1,91 @@
+# Numerical conventions
+
+This page explains how to interpret numerical inputs and results. A function's
+reference entry takes precedence when it defines a more specific contract.
+
+## Numeric types and finite values
+
+Python, C, and Wave tensor kernels use IEEE 754 binary64 (`float`, `double`,
+and `f64`). Checked native kernels reject non-finite inputs/results. This
+does not make floating-point arithmetic exact: rounding and cancellation
+still occur for accepted finite values.
+
+Haskell functions have the numeric type shown in their signature. Generic
+`Num`, `Fractional`, or `Floating` functions inherit the selected type's
+arithmetic. For example, bounded integers can overflow. A generic type
+constraint alone does not promise finite-result validation.
+
+## Shapes and storage
+
+Row-major order means consecutive elements in a row are adjacent. For
+`[[1, 2], [3, 4]]`, flat values are `[1, 2, 3, 4]`.
+
+| Shape | Meaning | Logical element count |
+| --- | --- | --- |
+| `()` | Scalar | 1 |
+| `(n,)` | Vector | n |
+| `(m, n)` | Matrix | m × n |
+| `(0, n)` | Empty matrix with n columns | 0 |
+
+The native ABI bounds rank at 8 and logical element count at 1,048,576.
+These are native-interface limits, not universal limits for every Haskell
+module. Native strides count elements; NumPy strides count bytes. Storage
+span and logical element count can differ for a strided view.
+
+## Compare approximate results
+
+Use exact equality for dimensions, counts, and deliberately exact small
+examples. For general floating-point results, compare with explicit
+absolute and relative tolerances:
+
+$$
+|x-y| \leq \max(\mathrm{atol},\ \mathrm{rtol}\,\max(|x|,|y|)).
+$$
+
+This is the symmetric rule used by Python's `math.isclose`; a particular
+solver may use a different stopping rule. Choose tolerances for the scale
+of your data and the required accuracy rather than copying a universal
+constant. Absolute tolerance matters when the expected result is near zero.
+
+```python
+import math
+from hmath import Tensor
+
+with Tensor([0.1, 0.2]) as samples:
+    assert math.isclose(samples.sum(), 0.3, rel_tol=1e-12, abs_tol=1e-15)
+```
+
+## Residual, conditioning, and rank
+
+For a proposed solution of $A x = b$, the residual is
+$r = A x - b$. A small residual checks that the computed vector nearly
+satisfies the supplied equation. It does not by itself establish a small
+solution error: an ill-conditioned problem can amplify small perturbations
+in the data.
+
+QR rank is a numerical decision controlled by a relative pivot tolerance.
+`qrRank` is an estimate, not an exact algebraic proof. A rank-deficient
+factorization can be inspected, while the QR least-squares solver rejects
+it. See [QR contracts](../reference/qr.md) for the exact threshold and shapes.
+
+## Errors and stopping conditions
+
+| Interface | Failure reporting | What to inspect |
+| --- | --- | --- |
+| Python | `HMathError` or a documented Python exception | Exception type, `.status`, message |
+| C | Checked status and thread-local diagnostic | Status immediately, then error message on the same thread |
+| Wave | Wrapper/native status | Status; native message can be stale after wrapper-only rejection |
+| Haskell checked API | `Either MathError a` | `Left` constructor or accepted `Right` result |
+| Haskell legacy pure API | Function-specific precondition or pure error | Reference entry before calling |
+
+An accepted result is not a universal accuracy guarantee. For iterative
+methods, inspect the reported stopping reason and residual/error estimate.
+An evaluation or iteration budget limits work; it is not a requested error
+bound. Calculus `estimateError` values are numerical estimates, not rigorous
+bounds on every possible callback.
+
+## Further reading
+
+- [SciPy linear solve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html): shape contracts and singular/ill-conditioned systems.
+- [GNU Scientific Library linear algebra](https://www.gnu.org/software/gsl/doc/html/linalg.html): decompositions, rank, and numerical tolerances.
+- [Python `math.isclose`](https://docs.python.org/3/library/math.html#math.isclose): approximate comparison semantics.

docs/getting-started/choose-an-operation.md

@@ -0,0 +1,45 @@
+# Choose an operation
+
+Start with the problem you want to solve, then select an available interface.
+The Python wheel provides the native scalar/tensor interface. The Haskell
+library additionally supplies the broader numerical algorithms below.
+
+## Python, C, and Wave
+
+| Task | Python | C | Wave |
+| --- | --- | --- | --- |
+| Add scalars | `add(a, b)` | `hmath_scalar_add` | `scalar_add` |
+| Construct a matrix | `Tensor.matrix` | `hmath_matrix_copy` | `matrix_copy` |
+| Sum all elements | `tensor.sum()` | `hmath_tensor_sum` | `tensor_sum` |
+| Sum corresponding products | `a.dot(b)` | `hmath_tensor_dot` | `tensor_dot` |
+| Add corresponding elements | `a + b` | `hmath_tensor_add` | `tensor_add` |
+| Multiply two matrices | `a @ b` | `hmath_tensor_matmul` | `tensor_matmul` |
+| Exchange NumPy arrays | `Tensor.from_numpy`, `tensor.numpy` | Storage/DLPack ABI | Use the C ABI |
+| Reinterpret a contiguous shape | `tensor.reshape` | `hmath_tensor_reshape` | Use the C ABI |
+
+`dot` reduces every corresponding product to a scalar at any rank. Matrix
+multiplication requires rank two and matching inner dimensions. Elementwise
+addition requires identical shapes. Broadcasting and batched matrix products
+are not supported by these wrappers.
+
+Read the [Python](../reference/python.md), [C](../reference/c.md), or
+[Wave](../reference/wave.md) reference for complete contracts. The
+[matrix calculation guide](../guides/matrix-calculation.md) demonstrates the
+three operations on the same inputs.
+
+## Haskell numerical algorithms
+
+| Task | Start with | Guide / reference |
+| --- | --- | --- |
+| Solve a square linear system | `HMath.Matrix.LU.solve` | [Linear algebra API](../reference/haskell.md) |
+| Fit an overdetermined linear model | `HMath.Matrix.QR.leastSquares` | [Least squares](../guides/least-squares.md) |
+| Regularize a linear fit | `HMath.Matrix.QR.ridgeSolve` | [QR contracts](../reference/qr.md) |
+| Estimate a derivative | `HMath.Calculus.derivative` | [Calculus](../guides/calculus.md) |
+| Integrate a smooth finite-interval function | `HMath.Calculus.integrateComposite` | [Calculus](../guides/calculus.md) |
+| Request adaptive integration | `HMath.Calculus.adaptiveIntegrate` | [Calculus](../guides/calculus.md) |
+| Compute descriptive statistics | `HMath.Statistics` | [Summarize data](../guides/statistics.md) |
+| Explore optimization, probability, graphs, or transforms | Select a numerical area | [Haskell module index](../reference/haskell.md) |
+
+These functions require Haskell imports; they are not additional Python
+functions. Check [numerical conventions](../explanation/numerical-conventions.md)
+before selecting tolerances or interpreting convergence.

docs/guides/borrowed-storage.md

@@ -0,0 +1,84 @@
+# Borrow storage safely in C
+
+Use `hmath_tensor_borrow` when copying a read-only binary64 buffer is
+undesirable. The caller must guarantee valid aligned storage, its readable
+span, and immutability through **every alias** until the last storage owner
+is released. See [C integration](c.md) for compiler and linker setup.
+
+## Give the tensor a cleanup callback
+
+This complete example allocates three values, borrows them, reshapes the
+tensor without copying, and checks that cleanup happens only after both
+handles are released. Metadata and ownership operations need no runtime
+initialization.
+
+```c
+#include <hmath.h>
+#include <assert.h>
+#include <stdlib.h>
+
+struct storage_owner {
+    double *values;
+    unsigned releases;
+};
+
+static void release_storage(void *userdata, const double *data) {
+    struct storage_owner *owner = userdata;
+    assert(data == owner->values);
+    free(owner->values);
+    owner->values = NULL;
+    owner->releases++;
+}
+
+int main(void) {
+    struct storage_owner owner = {malloc(3 * sizeof(double)), 0};
+    hmath_tensor *vector = NULL, *matrix = NULL;
+    const size_t vector_shape[] = {3}, matrix_shape[] = {1, 3};
+    hmath_status status;
+    if (!owner.values) return 1;
+    owner.values[0] = 1; owner.values[1] = 2; owner.values[2] = 3;
+    status = hmath_tensor_borrow(owner.values, 3, vector_shape, NULL, 1,
+                                 release_storage, &owner, &vector);
+    if (status != HMATH_OK) {
+        free(owner.values); /* Failed construction did not consume storage. */
+        return 1;
+    }
+    status = hmath_tensor_reshape(vector, matrix_shape, 2, &matrix);
+    hmath_tensor_release(vector);
+    if (status != HMATH_OK) return 1; /* Final vector release already cleaned up. */
+    assert(owner.releases == 0);
+    hmath_tensor_release(matrix);
+    assert(owner.releases == 1);
+    return 0;
+}
+```
+
+The expected exit status is zero. The `userdata` object lives on the stack
+until both handles have been released. A real application must retain its
+owner object for at least as long as any handle or exported view can survive.
+
+## Check strided layouts
+
+For nonnegative element strides $s_i$ and dimensions $d_i$, a nonempty
+tensor reaches a storage span of
+
+$$
+1+\sum_i(d_i-1)s_i.
+$$
+
+`storage_count` must cover that span, including gaps; it is not simply the
+product of dimensions. C strides are measured in **doubles**, not bytes.
+Passing NULL strides selects contiguous row-major storage. An empty tensor
+does not read any values. Rank zero is a scalar and needs one value.
+
+## Avoid lifetime errors
+
+The deleter runs synchronously on the final releasing thread, exactly once
+for successfully borrowed storage. It must not throw across C or invoke
+runtime shutdown. Failed construction leaves ownership with the caller.
+`hmath_tensor_data` returns a borrowed pointer and does not create another
+owner; never free that pointer separately.
+
+See [the C reference](../reference/c.md),
+[architecture and ownership](../explanation/architecture.md), and
+[NumPy borrowing](numpy.md) for related contracts.

docs/guides/calculus.md

@@ -0,0 +1,62 @@
+# Differentiate and integrate in Haskell
+
+Use `HMath.Calculus` for checked finite differences and numerical quadrature
+of `Double` callbacks. Callbacks must terminate and be side-effect free.
+These operations are Haskell APIs; they are not Python tensor methods.
+
+## Use an explicit step and panel count
+
+Estimate the derivative of $x^2$ at 3 and its integral from 0 to 1:
+
+```haskell
+import qualified HMath.Calculus as Calculus
+
+main :: IO ()
+main = do
+    let derivative = Calculus.derivative Calculus.Central 1e-5 (\x -> x*x) 3
+        integral = Calculus.integrateComposite Calculus.Simpson 100 (\x -> x*x) 0 1
+        expect expected result = case result of
+            Left err -> fail (show err)
+            Right actual ->
+                if abs (actual - expected) <= 1e-8
+                    then pure ()
+                    else fail ("Unexpected estimate: " ++ show actual)
+    expect 6 derivative
+    expect (1/3) integral
+    putStrLn "derivative and integral checks passed"
+```
+
+The expected output is `derivative and integral checks passed`.
+The differentiation step must be positive, finite, and representable at
+the evaluation point. Making it arbitrarily small can worsen cancellation.
+Simpson's panel count must be positive and even; `ThreeEighths` requires a
+multiple of three. A fixed-panel calculation does not return an error bound.
+
+## Request adaptive refinement
+
+```haskell
+import qualified HMath.Calculus as Calculus
+
+main :: IO ()
+main = case Calculus.adaptiveIntegrate Calculus.defaultTolerance (\x -> x*x) 0 1 of
+    Left err -> fail (show err)
+    Right report -> do
+        if abs (Calculus.estimateValue report - 1/3) <= 1e-9
+            && Calculus.estimateError report >= 0
+            && Calculus.evaluations report <= 100000
+            then putStrLn "adaptive integration check passed"
+            else fail (show report)
+```
+
+The default adaptive tolerance uses `1e-10` absolute, `1e-8` relative, and
+100,000 evaluations. Its stopping target is
+$\max(\mathrm{absoluteTolerance},\ \mathrm{relativeTolerance}\,|\mathrm{estimateValue}|)$.
+Absolute tolerance must be positive; relative tolerance must be nonnegative;
+both must be finite. A nonempty interval needs at least 15 evaluations.
+
+If refinement cannot meet the target within the budget, the adaptive
+integrator returns `Left (NonConvergence ...)`. Non-finite sampled values
+also fail. `estimateError` is an estimate from the approximation rule,
+not a rigorous bound for arbitrary callbacks. Inspect the full
+[Calculus API](../reference/haskell.md) and
+[numerical conventions](../explanation/numerical-conventions.md).