diff --git a/.github/workflows/code_style.yml b/.github/workflows/code_style.yml
new file mode 100644
index 0000000..9e3edf9
--- /dev/null
+++ b/.github/workflows/code_style.yml
@@ -0,0 +1,39 @@
+name: Lint and Format
+
+on:
+ push:
+ branches: [ master ]
+
+ pull_request:
+ branches: [ master ]
+ paths: [ '**/*.py' ]
+
+ workflow_dispatch:
+
+concurrency:
+ group: ${{ github.workflow }}/${{ github.ref }}
+ cancel-in-progress: true
+
+permissions:
+ contents: read
+
+jobs:
+ python-format:
+ name: Ruff Lint and Format
+ runs-on: ubuntu-latest
+
+ steps:
+ - name: Checkout Repository
+ uses: actions/checkout@v6
+
+ - name: Setup Ruff
+ uses: astral-sh/ruff-action@v3
+ with:
+ args: --version
+
+ - name: Enforce Lint
+ run: ruff check .
+
+ - name: Enforce Format
+ run: |
+ ruff format --diff .
diff --git a/bin/build.py b/bin/build.py
index 7f4e5ee..77d63da 100644
--- a/bin/build.py
+++ b/bin/build.py
@@ -3,7 +3,7 @@
import os
import shlex
import shutil
-import subprocess
+import subprocess # noqa: S404
import sys
from pathlib import Path
@@ -33,29 +33,35 @@
logger.addHandler(stderr_handler)
-def run_command(command, cwd=None):
+def run_command(command: str, cwd: str | None = None) -> None:
+ """Run a shell command, streaming output to the logger.
+
+ Args:
+ command (str): Shell command to execute.
+ cwd (str | None): Working directory; defaults to current directory.
+ """
if cwd is None:
logger.warning(
"No working directory specified. Using current directory."
)
- cwd = Path.cwd()
+ cwd_path = Path.cwd()
else:
- cwd = Path(cwd)
+ cwd_path = Path(cwd)
- log_file_path = cwd / "build.log"
+ log_file_path = cwd_path / "build.log"
- logger.info(f"Executing command: '{command}' in '{cwd}'")
+ logger.info(f"Executing command: '{command}' in '{cwd_path}'")
- with subprocess.Popen(
+ with subprocess.Popen( # noqa: S603
shlex.split(command),
- cwd=str(cwd),
+ cwd=str(cwd_path),
stdout=subprocess.PIPE,
stderr=subprocess.STDOUT,
bufsize=0,
env=dict(**os.environ, PYTHONUNBUFFERED="1"),
text=True,
) as proc:
- with open(log_file_path, "a") as _log_file:
+ with open(log_file_path, "a", encoding="utf-8") as _log_file:
for line in proc.stdout:
logger.debug(line.rstrip())
@@ -68,27 +74,32 @@ def run_command(command, cwd=None):
logger.info("Command executed successfully.")
-def install():
+def install() -> None:
+ """Install the package into the current environment."""
run_command("uv pip install . -v")
-def wheel():
+def wheel() -> None:
+ """Build a wheel distribution."""
run_command("uv build --wheel -v")
-def clean():
+def clean() -> None:
+ """Remove all build artifacts and generated files."""
logger.debug("Starting cleanup ...")
run_command("uv pip uninstall polpack")
+ remove_names: set[str] = {
+ "dist",
+ "build",
+ "lib",
+ ".pytest_cache",
+ ".ruff_cache",
+ }
+
for entry in Path("").iterdir():
- if entry.name in [
- "dist",
- "build",
- "lib",
- ".pytest_cache",
- ".ruff_cache",
- ]:
+ if entry.name in remove_names:
logger.info(f"Removing '{entry}'")
shutil.rmtree(entry)
if entry.name == "bin" and entry.is_dir():
@@ -109,7 +120,8 @@ def clean():
logger.info("Finished cleanup.")
-def main():
+def main() -> None:
+ """Parse arguments and dispatch to the appropriate build action."""
parser = argparse.ArgumentParser(description="PolPack Build Script")
parser.add_argument(
"mode",
diff --git a/bin/get_version.py b/bin/get_version.py
index dd19747..67fe9c9 100644
--- a/bin/get_version.py
+++ b/bin/get_version.py
@@ -1,6 +1,6 @@
#!/usr/bin/env python3
-import subprocess
+import subprocess # noqa: S404
import sys
from pathlib import Path
@@ -16,7 +16,8 @@
from setuptools_scm import get_version
-def main():
+def main() -> None:
+ """Print the current package version derived from the git tags."""
root = Path(__file__).parent.parent
version = get_version(root=root)
if version is None:
diff --git a/docs/api.md b/docs/api.md
index ded7e10..8620f25 100644
--- a/docs/api.md
+++ b/docs/api.md
@@ -1,262 +1,723 @@
# API Reference
-All public routines are exported via the `polpack` Python module in lowercase.
+All public routines are exported directly via the `polpack` Python module.
---
## Overview
-`polpack` includes routines to evaluate a variety of recursively defined
-polynomial families and special functions. Most generators follow a three-term
-recurrence and write the output into a pre-allocated NumPy array.
+`polpack` provides routines to evaluate recursively-defined polynomial families,
+combinatorial sequences, number-theoretic functions, and special functions.
+Most array-valued routines write their output into a **pre-allocated NumPy array**
+passed by the caller.
-| Family | Category | Type | Domain |
-|---|---|---|---|
-| `bernoulli_poly` | Polynomial | Sequence | \( \mathbb{R} \) |
-| `bernstein_poly` | Polynomial | Basis | \( [0, 1] \) |
-| `cardan_poly` | Polynomial | Cubic | \( \mathbb{R} \) |
-| `charlier` | Polynomial | Discrete | \( \mathbb{N} \) |
-| `cheby_t_poly` | Polynomial | Orthogonal | \( [-1, 1] \) |
-| `gegenbauer_poly` | Polynomial | Orthogonal | \( [-1, 1] \) |
-| `jacobi_poly` | Polynomial | Orthogonal | \( [-1, 1] \) |
-| `laguerre_poly` | Polynomial | Orthogonal | \( [0, \infty) \) |
-| `legendre_poly` | Polynomial | Orthogonal | \( [-1, 1] \) |
-
----
-
-## Common interface
-
-### Arguments
+### Common argument conventions
| Argument | Type | Direction | Description |
|---|---|---|---|
-| `n` | `int` | input | Highest degree or number of terms to compute |
+| `n` | `int` | input | Highest degree or largest index to compute |
| `x` | `float` | input | Evaluation point |
-| `v` | `ndarray` | in/out | Pre-allocated output array of size `n+1` |
+| `v` | `ndarray` | output | Pre-allocated output array, size `n+1` |
!!! warning "Array allocation"
- The output array `v` must be pre-allocated with the correct size before
- calling the routine. Passing an undersized array will raise a runtime error.
+ Output arrays **must** be pre-allocated before calling any routine.
+ Passing an undersized array will raise a runtime error.
+
+!!! note "Array order"
+ Use `order="F"` when allocating output arrays. Integer-valued sequences
+ use `dtype=np.int32`; floating-point routines use `dtype=np.float64`.
---
-## Orthogonal Polynomials
+## Orthogonal Polynomial Families
### `cheby_t_poly(m, n, x, cx)`
-**Chebyshev polynomials of the first kind.**
+**Chebyshev polynomials of the first kind** $T_0(x), \ldots, T_n(x)$.
+
+- `m` (int): Number of evaluation points (typically 1 for scalar evaluation).
+- `n` (int): Highest degree.
+- `x` (float or ndarray[m]): Evaluation point(s).
+- `cx` (ndarray[m, n+1]): Output array. For scalar evaluation (`m=1`), use shape `(n+1,)`.
+
+For $x \in [-1, 1]$: $T_k(x) = \cos(k \arccos x)$.
+
+```python
+import numpy as np, polpack
+
+cx = np.zeros(6, dtype=np.float64, order="F")
+polpack.cheby_t_poly(1, 5, 0.5, cx)
+```
+
+### `cheby_u_poly(m, n, x, cx)`
+**Chebyshev polynomials of the second kind** $U_0(x), \ldots, U_n(x)$.
- `m` (int): Number of evaluation points.
- `n` (int): Highest degree.
-- `x` (ndarray[m]): Evaluation points.
-- `cx` (ndarray[m, n+1]): Output matrix.
+- `x` (float or ndarray[m]): Evaluation point(s).
+- `cx` (ndarray[m, n+1]): Output array.
+
+For $x \in [-1, 1]$: $U_k(x) = \sin((k+1)\arccos x) / \sqrt{1-x^2}$.
### `gegenbauer_poly(n, alpha, x, cx)`
-**Gegenbauer polynomials.**
+**Gegenbauer (ultraspherical) polynomials** $C_0^{(\alpha)}(x), \ldots, C_n^{(\alpha)}(x)$.
- `n` (int): Highest degree.
-- `alpha` (float): Parameter \( \alpha > -1/2 \).
+- `alpha` (float): Shape parameter; must satisfy $\alpha > -1/2$.
- `x` (float): Evaluation point.
- `cx` (ndarray[n+1]): Output array.
+Special case: `alpha=1` gives Chebyshev $U_n$; `alpha=0.5` gives Legendre $P_n$.
+
### `hermite_poly_phys(n, x, cx)`
-**Hermite polynomials (Physicist's).**
+**Physicist's Hermite polynomials** $H_0(x), \ldots, H_n(x)$.
- `n` (int): Highest degree.
- `x` (float): Evaluation point.
- `cx` (ndarray[n+1]): Output array.
+Recurrence: $H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x)$.
+
+### `gen_hermite_poly(n, x, mu, p)`
+**Generalized Hermite polynomials** orthogonal under $|x|^{2\mu} e^{-x^2}$.
+
+- `n` (int): Highest degree.
+- `x` (float): Evaluation point.
+- `mu` (float): Shape parameter; `mu=0` gives the standard Hermite polynomial.
+- `p` (ndarray[n+1]): Output array.
+
### `jacobi_poly(n, alpha, beta, x, cx)`
-**Jacobi polynomials.**
+**Jacobi polynomials** $P_0^{(\alpha,\beta)}(x), \ldots, P_n^{(\alpha,\beta)}(x)$.
- `n` (int): Highest degree.
-- `alpha` (float): Parameter \( \alpha \).
-- `beta` (float): Parameter \( \beta \).
+- `alpha` (float): First shape parameter; must satisfy $\alpha > -1$.
+- `beta` (float): Second shape parameter; must satisfy $\beta > -1$.
- `x` (float): Evaluation point.
- `cx` (ndarray[n+1]): Output array.
### `laguerre_poly(n, x, cx)`
-**Laguerre polynomials.**
+**Laguerre polynomials** $L_0(x), \ldots, L_n(x)$.
+
+- `n` (int): Highest degree.
+- `x` (float): Evaluation point; typically $x \geq 0$.
+- `cx` (ndarray[n+1]): Output array.
+
+### `laguerre_associated(n, m, x, cx)`
+**Associated Laguerre polynomials** $L_0^m(x), \ldots, L_n^m(x)$.
+
+- `n` (int): Highest degree.
+- `m` (int): Association order; $m \geq 0$.
+- `x` (float): Evaluation point.
+- `cx` (ndarray[n+1]): Output array.
+
+### `gen_laguerre_poly(n, alpha, x, cx)`
+**Generalized (associated) Laguerre polynomials** $L_n^{(\alpha)}(x)$.
- `n` (int): Highest degree.
+- `alpha` (float): Shape parameter; must satisfy $\alpha > -1$.
- `x` (float): Evaluation point.
- `cx` (ndarray[n+1]): Output array.
### `legendre_poly(n, x, cx, cpx)`
-**Legendre polynomials.**
+**Legendre polynomials** $P_0(x), \ldots, P_n(x)$ and their first derivatives.
- `n` (int): Highest degree.
- `x` (float): Evaluation point.
-- `cx` (ndarray[n+1]): Output array \( P_n(x) \).
-- `cpx` (ndarray[n+1]): Output array of derivatives \( P'_n(x) \).
+- `cx` (ndarray[n+1]): Output array $P_k(x)$.
+- `cpx` (ndarray[n+1]): Output array of first derivatives $P_k'(x)$.
----
+### `legendre_associated(n, m, x, cx)`
+**Associated Legendre functions** $P_n^m(x)$.
-## Combinatorial Sequences
+- `n` (int): Maximum degree.
+- `m` (int): Order; $0 \leq m \leq n$.
+- `x` (float): Evaluation point; $-1 \leq x \leq 1$.
+- `cx` (ndarray[n+1]): Output array.
-### `bell(n, b)`
-**Bell numbers.**
+### `legendre_associated_normalized(n, m, x, cx)`
+**Normalized associated Legendre functions** with the normalization factor
+$\sqrt{(2n+1)(n-m)! / (2(n+m)!)}$ applied.
-- `n` (int): Highest index.
-- `b` (ndarray[n+1]): Output array of Bell numbers \( B_0, \dots, B_n \).
+- `n` (int): Maximum degree.
+- `m` (int): Order; $0 \leq m \leq n$.
+- `x` (float): Evaluation point.
+- `cx` (ndarray[n+1]): Output array.
-### `catalan(n, c)`
-**Catalan numbers.**
+### `legendre_function_q(n, x, cx)`
+**Legendre functions of the second kind** $Q_0(x), \ldots, Q_n(x)$.
-- `n` (int): Highest index.
-- `c` (ndarray[n+1]): Output array of Catalan numbers \( C_0, \dots, C_n \).
+- `n` (int): Highest degree.
+- `x` (float): Evaluation point; $-1 < x < 1$ required.
+- `cx` (ndarray[n+1]): Output array.
-### `stirling1(n, m, s1)` / `stirling2(n, m, s2)`
-**Stirling numbers.**
+$Q_0(x) = \frac{1}{2}\ln\frac{1+x}{1-x}$, with recurrence
+$(n+1)Q_{n+1}(x) = (2n+1)x Q_n(x) - n Q_{n-1}(x)$.
-- `n` (int): Number of elements.
-- `m` (int): Number of subsets/cycles.
-- `s1/s2` (ndarray[n, m]): Output matrix.
+### `spherical_harmonic(l, m, theta, phi, c, s)`
+**Real and imaginary parts of the spherical harmonic** $Y_l^m(\theta, \phi)$.
+
+- `l` (int): Degree; $l \geq 0$.
+- `m` (int): Order; $0 \leq m \leq l$.
+- `theta` (float): Polar angle in radians.
+- `phi` (float): Azimuthal angle in radians.
+- `c` (ndarray[1]): Output real part.
+- `s` (ndarray[1]): Output imaginary part.
---
-## Detailed Routine Reference
+## Discrete Orthogonal Polynomials
-### Bernoulli and Euler Sequences
+### `charlier(n, a, x, value)`
+**Charlier polynomials** $C_0(x;a), \ldots, C_n(x;a)$, orthogonal w.r.t. the
+Poisson distribution with parameter $a > 0$.
-#### `bernoulli_number(n, b)`
-Computes the sequence of Bernoulli numbers \( B_0, \dots, B_n \).
+- `n` (int): Highest degree.
+- `a` (float): Parameter; $a > 0$.
+- `x` (float): Evaluation point.
+- `value` (ndarray[n+1]): Output array.
-- `n` (int): Highest index.
-- `b` (ndarray[n+1]): Output array.
+### `chebyshev_discrete(n, m, x, v)`
+**Discrete Chebyshev polynomials** $t_0(x;N), \ldots, t_n(x;N)$, orthogonal
+on the uniform grid $\{0, 1, \ldots, N-1\}$.
-#### `bernoulli_poly(n, x, bx)`
-Evaluates the Bernoulli polynomial \( B_n(x) \) at \( x \).
+- `n` (int): Highest degree.
+- `m` (int): Grid size $N$.
+- `x` (float): Evaluation point.
+- `v` (ndarray[n+1]): Output array.
-- `n` (int): Degree.
+### `krawtchouk(n, p, x, m, v)`
+**Krawtchouk polynomials** $K_0, \ldots, K_n$, orthogonal w.r.t. the binomial
+distribution $\text{Bin}(M, p)$.
+
+- `n` (int): Highest degree.
+- `p` (float): Success probability; $0 < p < 1$.
- `x` (float): Evaluation point.
-- `bx` (float): Output value.
+- `m` (int): Number of trials $M$.
+- `v` (ndarray[n+1]): Output array.
-#### `euler_number(n, e)`
-Computes the sequence of Euler numbers \( E_0, \dots, E_n \).
+### `meixner(n, beta, c, x, v)`
+**Meixner polynomials** $M_0, \ldots, M_n$, orthogonal w.r.t. the negative
+binomial distribution.
-- `n` (int): Highest index.
-- `e` (ndarray[n+1]): Output array.
+- `n` (int): Highest degree.
+- `beta` (float): First parameter; $\beta > 0$.
+- `c` (float): Second parameter; $0 < c < 1$.
+- `x` (float): Evaluation point.
+- `v` (ndarray[n+1]): Output array.
-#### `euler_poly(n, x)`
-Evaluates the \( n \)-th Euler polynomial \( E_n(x) \) at \( x \).
+---
+
+## Non-Orthogonal Polynomial Families
+
+### `bernoulli_poly(n, x, bx)`
+**Bernoulli polynomial** $B_n(x)$.
- `n` (int): Degree.
- `x` (float): Evaluation point.
-- Returns `float`.
+- `bx` (float, in/out): Output value $B_n(x)$.
-### Bernstein Polynomials
+### `bernstein_poly(n, x, bern)`
+**Bernstein basis polynomials** $b_{0,n}(x), \ldots, b_{n,n}(x)$ on $[0, 1]$.
-#### `bernstein_poly(n, x, bern)`
-Evaluates the \( n+1 \) Bernstein basis polynomials of degree \( n \) at \( x \in [0, 1] \).
+- `n` (int): Degree.
+- `x` (float): Evaluation point; $0 \leq x \leq 1$.
+- `bern` (ndarray[n+1]): Output array. Satisfies $\sum_k b_{k,n}(x) = 1$.
+
+### `bpab(n, x, a, b, bern)`
+**Bernstein basis polynomials** on the interval $[a, b]$.
- `n` (int): Degree.
+- `x` (float): Evaluation point; $a \leq x \leq b$.
+- `a`, `b` (float): Interval boundaries.
+- `bern` (ndarray[n+1]): Output array.
+
+### `cardan_poly(n, x, s, cx)`
+**Cardan polynomials** $C_0(x,s), \ldots, C_n(x,s)$.
+
+- `n` (int): Highest degree.
- `x` (float): Evaluation point.
-- `bern` (ndarray[n+1]): Output array of basis values.
+- `s` (float): Scale parameter.
+- `cx` (ndarray[n+1]): Output array.
-#### `bpab(n, x, a, b, bern)`
-Evaluates the \( n+1 \) Bernstein basis polynomials of degree \( n \) on \( [a, b] \).
+### `euler_poly(n, x)`
+**Euler polynomial** $E_n(x)$.
- `n` (int): Degree.
- `x` (float): Evaluation point.
-- `a, b` (float): Interval boundaries.
-- `bern` (ndarray[n+1]): Output array of basis values.
+- Returns `float`: value of $E_n(x)$.
+
+### `complete_symmetric_poly(n, r, x, value)`
+**Complete symmetric polynomial** $h_r(x_1, \ldots, x_n)$.
+
+- `n` (int): Number of variables.
+- `r` (int): Degree.
+- `x` (ndarray[n]): Variable values.
+- `value` (float, out): Output value.
+
+### `zernike_poly(m, n, rho, z)`
+**Zernike radial polynomial** $R_n^m(\rho)$.
+
+- `m` (int): Azimuthal order; $0 \leq m \leq n$.
+- `n` (int): Radial degree; $\text{mod}(n-m, 2) = 0$ required.
+- `rho` (float): Radial coordinate; $0 \leq \rho \leq 1$.
+- `z` (ndarray): Output array.
-### Fibonacci and Power Sequences
+---
+
+## Polynomial Coefficients and Zeros
-#### `fibonacci_recursive(n, f)`
-Computes the first \( n \) Fibonacci numbers.
+These routines compute the exact polynomial coefficient arrays and Chebyshev nodes.
+
+| Routine | Parameters | Description |
+|---|---|---|
+| `cheby_t_poly_coef(n, c)` | `c` (ndarray[n+1, n+1]) | Coefficient matrix of $T_0, \ldots, T_n$ |
+| `cheby_t_poly_zero(n, z)` | `z` (ndarray[n]) | Zeros of $T_n$: $\cos\bigl(\frac{(2k-1)\pi}{2n}\bigr)$ |
+| `cheby_u_poly_coef(n, c)` | `c` (ndarray[n+1, n+1]) | Coefficient matrix of $U_0, \ldots, U_n$ |
+| `cheby_u_poly_zero(n, z)` | `z` (ndarray[n]) | Zeros of $U_n$: $\cos\bigl(\frac{k\pi}{n+1}\bigr)$ |
+| `hermite_poly_phys_coef(n, c)` | `c` (ndarray[n+1, n+1]) | Coefficient matrix of $H_0, \ldots, H_n$ |
+| `laguerre_poly_coef(n, c)` | `c` (ndarray[n+1, n+1]) | Coefficient matrix of $L_0, \ldots, L_n$ |
+| `legendre_poly_coef(n, c)` | `c` (ndarray[n+1, n+1]) | Coefficient matrix of $P_0, \ldots, P_n$ |
+| `cardan_poly_coef(n, s, c)` | `s` (float), `c` (ndarray[n+1, n+1]) | Cardan coefficients for parameter $s$ |
+| `zernike_poly_coef(m, n, c)` | `c` (ndarray[n+1]) | Coefficients of $R_n^m(\rho)$ |
+
+---
+
+## Combinatorial Sequences
+
+### `bell(n, b)`
+**Bell numbers** $B_0, \ldots, B_n$.
+
+- `n` (int): Highest index.
+- `b` (ndarray[n+1], int32): Output array.
+
+$B_n$ counts the number of set partitions of $\{1, \ldots, n\}$.
+
+### `catalan(n, c)`
+**Catalan numbers** $C_0, \ldots, C_n$.
+
+- `n` (int): Highest index.
+- `c` (ndarray[n+1], int32): Output array.
+
+$C_n = \frac{1}{n+1}\binom{2n}{n}$.
+
+### `catalan_row_next(n, row)`
+**Next row of the Catalan triangle**, given the current row for index $n$.
+
+- `n` (int): Current row index.
+- `row` (ndarray, int32): In/out array; updated to row $n+1$ on return.
+
+### `stirling1(n, m, s1)`
+**Stirling numbers of the first kind** $s(k, j)$ for $1 \leq k \leq n$, $1 \leq j \leq m$.
+
+- `n` (int): Row count.
+- `m` (int): Column count.
+- `s1` (ndarray[n, m], int32): Output matrix.
+
+$s(n, k)$ counts permutations of $n$ elements with exactly $k$ cycles.
+
+### `stirling2(n, m, s2)`
+**Stirling numbers of the second kind** $S(k, j)$.
+
+- `n` (int): Row count.
+- `m` (int): Column count.
+- `s2` (ndarray[n, m], int32): Output matrix.
+
+$S(n, k)$ counts partitions of $n$ elements into exactly $k$ non-empty subsets.
+
+### `delannoy(m, n, a)`
+**Delannoy numbers** $D(i,j)$ for $0 \leq i \leq m$, $0 \leq j \leq n$.
+
+- `m`, `n` (int): Grid dimensions.
+- `a` (ndarray[m+1, n+1], int32): Output matrix.
+
+$D(m,n)$ counts lattice paths from $(0,0)$ to $(m,n)$ using steps east, north, northeast.
+
+### `motzkin(n, a)`
+**Motzkin numbers** $M_0, \ldots, M_n$.
+
+- `n` (int): Highest index.
+- `a` (ndarray[n+1], int32): Output array.
+
+### `fibonacci_recursive(n, f)`
+**Fibonacci sequence** $F_1, \ldots, F_n$ using the recurrence $F_{k+1} = F_k + F_{k-1}$.
- `n` (int): Number of terms.
-- `f` (ndarray[n]): Output array.
+- `f` (ndarray[n], int32): Output array.
-#### `fibonacci_direct(n, f)`
-Computes the \( n \)-th Fibonacci number using the closed-form Binet formula.
+### `fibonacci_direct(n, f)`
+**$n$-th Fibonacci number** using the Binet closed-form formula.
- `n` (int): Index.
-- `f` (int): Output value.
+- `f` (int32, out): Output value $F_n$.
-### Stirling Sequences
+### `fibonacci_floor(n, f, i)`
+**Largest Fibonacci number** $\leq n$.
-#### `stirling1(n, m, s1)`
-Computes the Stirling numbers of the first kind \( s(n, k) \).
+- `n` (int): Upper bound.
+- `f` (int32, out): Largest Fibonacci number $\leq n$.
+- `i` (int32, out): Index of that Fibonacci number.
-- `n` (int): Order.
-- `m` (int): Term index.
-- `s1` (ndarray[n, m]): Output matrix.
+### `vibonacci(n, seed, v)`
+**Vibonacci numbers**: randomized Fibonacci where each recurrence sign is $\pm 1$ chosen
+randomly. The ratio $V_{n+1}/V_n$ converges almost surely to Viswanath's constant $\approx 1.13199$.
-#### `stirling2(n, m, s2)`
-Computes the Stirling numbers of the second kind \( S(n, k) \).
+- `n` (int): Number of terms.
+- `seed` (int32, in/out): Random number seed; updated on return.
+- `v` (ndarray[n], int32): Output array.
-- `n` (int): Order.
-- `m` (int): Term index.
-- `s2` (ndarray[n, m]): Output matrix.
+### `zeckendorf(n, m_max, m, i_list, f_list)`
+**Zeckendorf decomposition**: every positive integer is uniquely a sum of
+non-consecutive Fibonacci numbers.
----
+- `n` (int): Integer to decompose.
+- `m_max` (int): Maximum number of Fibonacci terms expected.
+- `m` (int32, out): Number of terms found.
+- `i_list` (ndarray[m_max], int32): Fibonacci indices used.
+- `f_list` (ndarray[m_max], int32): Fibonacci values used.
-### Special Math Functions
+### `bernoulli_number(n, b)`
+**Bernoulli numbers** $B_0, \ldots, B_n$ (floating-point).
-#### `agm_values(n_data, a, b, fx)`
-Returns tabulated values of the Arithmetic-Geometric Mean.
+- `n` (int): Highest index.
+- `b` (ndarray[n+1], float64): Output array.
-- `n_data` (int): Index of the value to retrieve (0 for first).
-- `a, b` (float): Parameters.
-- `fx` (float): Output AGM value.
+### `euler_number(n, e)`
+**Euler numbers** $E_0, \ldots, E_n$. Odd-indexed Euler numbers are 0.
-#### `agud(g)`
-Evaluates the inverse Gudermannian function.
+- `n` (int): Highest index.
+- `e` (ndarray[n+1], float64): Output array.
-- `g` (float): Input value.
-- Returns `float`.
+### `eulerian(n, e)`
+**Eulerian numbers** $E(k, j)$ for $1 \leq k \leq n$, $1 \leq j \leq n$.
+
+- `n` (int): Dimension.
+- `e` (ndarray[n, n], int32): Output matrix. $E(n, k)$ counts permutations of $n$
+ objects with exactly $k$ runs.
-#### `beta_values(n_data, x, y, fxy)`
-Returns tabulated values of the Beta function \( \text{B}(x, y) \).
+### `poly_bernoulli(n, k, b)`
+**Poly-Bernoulli numbers** $B_n^{(-k)}$ with negative superscript.
-#### `erf_values(n_data, x, fx)`
-Returns tabulated values of the error function \( \text{erf}(x) \).
+- `n`, `k` (int): Indices.
+- `b` (float64, out): Output value.
-#### `gamma_values(n_data, x, fx)`
-Returns tabulated values of the Gamma function \( \Gamma(x) \).
+$B_n^{(-k)}$ counts the number of $n \times k$ binary lonesum matrices.
-#### `lambert_w(x)`
-Estimates the Lambert W function \( W(x) \), the inverse of \( f(w) = we^w \).
+### `comb_row_next(n, row)`
+**Next row of Pascal's triangle** (binomial coefficients $\binom{n}{0}, \ldots, \binom{n}{n}$).
+
+- `n` (int): Current row.
+- `row` (ndarray, int32): In/out; row $n$ on input, row $n+1$ on return.
---
-### Number Theory and Primality
+## Figurate Numbers
-| Routine | Parameters | Description |
+| Routine | Description | Formula |
|---|---|---|
-| `phi(n, phin)` | `n` (int): input
`phin` (int): out | Euler totient function \( \phi(n) \). |
-| `moebius(n, mu)` | `n` (int): input
`mu` (int): out | Moebius function \( \mu(n) \). |
-| `mertens(n)` | `n` (int): input | Mertens function \( M(n) \) (returns int). |
-| `prime(n)` | `n` (int): input | Returns the \( n \)-th prime number. |
-| `i4_is_prime(n)` | `n` (int): input | Primality test (returns bool). |
-| `i4_factor(...)` | `n`, `f_max`, `f_num`, `f`, `p`, `nleft` | Prime factorization into components. |
-| `sigma(n, sigma_n)` | `n` (int): input
`sigma_n` (int): out | Sum of divisors \( \sigma(n) \). |
-| `tau(n, taun)` | `n` (int): input
`taun` (int): out | Number of divisors \( \tau(n) \). |
-| `jacobi_symbol(q, p, j)` | `q, p` (int), `j` (int): out | Jacobi symbol \( (q/p) \). |
-| `legendre_symbol(q, p, l)` | `q, p` (int), `l` (int): out | Legendre symbol \( (q/p) \). |
-
-### Sequences and Series
+| `triangle_num(n)` | $n$-th triangular number | $T(n) = n(n+1)/2$ |
+| `tetrahedron_num(n)` | $n$-th tetrahedral number | $T_3(n) = n(n+1)(n+2)/6$ |
+| `pentagon_num(n, p)` | $n$-th pentagonal number | $P(n) = n(3n-1)/2$ (also defined for $n < 0$) |
+| `pyramid_num(n)` | $n$-th triangular pyramidal number | $(n+1)^3/6 - (n+1)/6$ |
+| `pyramid_square_num(n)` | $n$-th square pyramidal number | $n(n+1)(2n+1)/6$ |
+| `simplex_num(m, n)` | $n$-th simplex number in $m$ dimensions | $\binom{n+m-1}{m}$ |
+| `plane_partition_num(n)` | Number of plane partitions of $n$ | — |
+| `lock(n, a)` | Number of full-button codes for an $n$-button lock | — |
+| `align_enum(m, n)` | Number of alignments of sequences of length $m$ and $n$ | — |
-| Routine | Parameters | Description |
-|---|---|---|
-| `catalan(n, c)` | `n` (int), `c` (ndarray) | Catalan numbers sequence. |
-| `delannoy(m, n, a)` | `m`, `n` (int), `a` (ndarray) | Delannoy numbers matrix. |
-| `motzkin(n, a)` | `n` (int), `a` (ndarray) | Motzkin numbers sequence. |
-| `vibonacci(n, seed, v)` | `n`, `seed` (int), `v` (ndarray) | Random Vibonacci terms. |
-| `collatz_count(n)` | `n` (int): input | Sequence length (returns int). |
-| `collatz_count_max(...)` | `n, i_max, j_max` | Maximum count in range \( [1, n] \). |
+---
-### Special Utilities
+## Number Theory
-| Routine | Parameters | Description |
-|---|---|---|
-| `agm_values(...)` | `n_data, a, b, fx` | Arithmetic-Geometric Mean data. |
-| `i4_choose(n, k)` | `n, k` (int) | Binomial coefficient (returns int). |
-| `i4_swap(i, j)` | `i, j` (int) | Swaps values in-place (in/out). |
-| `i4_huge()` | (none) | Returns largest representable i4. |
-| `zeckendorf(...)` | `n, m_max, m, i_l, f_l` | Zeckendorf representation decomposition. |
-| `complete_symmetric_poly` | `n, r, x, value` | Evaluates \( h_r(x_1, \dots, x_n) \). |
-| `zernike_poly` | `m, n, rho, z` | Zernike polynomial at \( \rho \). |
-| `zeta(p)` | `p` (float): input | Riemann Zeta function \( \zeta(p) \) (returns float). |
+### `phi(n, phin)`
+**Euler totient function** $\phi(n)$.
+
+- `n` (int): Input.
+- `phin` (int32, out): $\phi(n)$ = number of integers in $[1, n]$ coprime to $n$.
+
+### `moebius(n, mu)`
+**Moebius function** $\mu(n)$.
+
+- `n` (int): Input.
+- `mu` (int32, out): $\mu(n) \in \{-1, 0, 1\}$.
+
+### `mertens(n)`
+**Mertens function** $M(n) = \sum_{k=1}^n \mu(k)$.
+
+- `n` (int): Input.
+- Returns `int`: $M(n)$.
+
+### `sigma(n, sigma_n)`
+**Sum of divisors** $\sigma(n) = \sum_{d | n} d$.
+
+- `n` (int): Input.
+- `sigma_n` (int32, out): $\sigma(n)$.
+
+### `tau(n, taun)`
+**Number of divisors** $\tau(n) = \sum_{d | n} 1$.
+
+- `n` (int): Input.
+- `taun` (int32, out): $\tau(n)$.
+
+### `omega(n, ndiv)`
+**Number of distinct prime divisors** $\omega(n)$.
+
+- `n` (int): Input.
+- `ndiv` (int32, out): $\omega(n)$.
+
+### `prime(n)`
+**The $n$-th prime number** (precomputed table up to $n = 1600$; largest prime is 13499).
+
+- `n` (int): Index.
+- Returns `int`: $p_n$.
+
+### `i4_is_prime(n)`
+**Primality test** using trial division.
+
+- `n` (int): Input.
+- Returns `bool`: `True` if $n$ is prime.
+
+### `i4_factor(n, factor_max, factor_num, factor, power, nleft)`
+**Prime factorization**: $n = \left(\prod_i \text{factor}[i]^{\text{power}[i]}\right) \times \text{nleft}$.
+
+- `n` (int): Integer to factor.
+- `factor_max` (int): Maximum number of prime factors to find.
+- `factor_num` (int32, out): Number of distinct prime factors found.
+- `factor` (ndarray[factor_max], int32): Output prime bases.
+- `power` (ndarray[factor_max], int32): Output exponents.
+- `nleft` (int32, out): Unfactored remainder (1 if fully factored).
+
+### `jacobi_symbol(q, p, j)`
+**Jacobi symbol** $(q/p)$.
+
+- `q`, `p` (int): Inputs; $p$ must be odd and positive.
+- `j` (int32, out): Jacobi symbol $\in \{-1, 0, 1\}$.
+
+### `legendre_symbol(q, p, l)`
+**Legendre symbol** $(q/p)$ for odd prime $p$.
+
+- `q`, `p` (int): Inputs.
+- `l` (int32, out): Legendre symbol $\in \{-1, 0, 1\}$.
+
+### `benford(ival)`
+**Benford's law probability** for leading digits.
+
+- `ival` (int): Leading digit (1–9) or multi-digit integer.
+- Returns `float`: $\log_{10}((\text{ival}+1)/\text{ival})$.
+
+### `collatz_count(n)`
+**Number of steps in the Collatz sequence** starting from $n$ until reaching 1.
+
+- `n` (int): Starting value.
+- Returns `int`: Number of steps.
+
+### `collatz_count_max(n, i_max, j_max)`
+**Maximum Collatz count** for starting values in $[1, n]$.
+
+- `n` (int): Upper bound.
+- `i_max` (int32, out): Starting value with the longest sequence.
+- `j_max` (int32, out): Length of that longest sequence.
+
+---
+
+## Special Functions
+
+### `zeta(p)`
+**Riemann Zeta function** $\zeta(p) = \sum_{n=1}^\infty 1/n^p$.
+
+- `p` (float): Argument; $p > 1$.
+- Returns `float`.
+
+Notable values: $\zeta(2) = \pi^2/6$, $\zeta(4) = \pi^4/90$.
+
+### `lambert_w(x)`
+**Lambert W function** $W(x)$, the principal branch satisfying $W(x) e^{W(x)} = x$.
+
+- `x` (float): Argument; $x \geq -1/e$.
+- Returns `float`.
+
+### `lambert_w_crude(x)`
+**Crude initial estimate of the Lambert W function** using a simple approximation.
+Useful as a starting point for Newton iterations.
+
+- `x` (float): Argument.
+- Returns `float`.
+
+### `r8_erf(x)`
+**Error function** $\mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2}\, dt$.
+
+- `x` (float): Argument.
+- Returns `float` in $(-1, 1)$.
+
+### `r8_erf_inverse(y)`
+**Inverse error function** $\mathrm{erf}^{-1}(y)$.
+
+- `y` (float): Argument; $-1 < y < 1$.
+- Returns `float`.
+
+### `r8_beta(x, y)`
+**Beta function** $B(x, y) = \Gamma(x)\Gamma(y)/\Gamma(x+y)$.
+
+- `x`, `y` (float): Arguments; both must be positive.
+- Returns `float`.
+
+### `r8_agm(a, b)`
+**Arithmetic-Geometric Mean** of $a$ and $b$.
+
+- `a`, `b` (float): Non-negative inputs.
+- Returns `float`.
+
+### `r8_gamma_log(x)`
+**Log-Gamma function** $\ln\Gamma(x)$.
+
+- `x` (float): Argument; $x > 0$.
+- Returns `float`.
+
+### `r8_hyper_2f1(a, b, c, x)`
+**Gauss hypergeometric function** ${}_{2}F_1(a, b; c; x)$.
+
+- `a`, `b`, `c` (float): Parameters.
+- `x` (float): Argument; $|x| < 1$.
+- Returns `float`.
+
+### `r8_psi(x)`
+**Digamma (Psi) function** $\psi(x) = \frac{d}{dx}\ln\Gamma(x) = \Gamma'(x)/\Gamma(x)$.
+
+- `x` (float): Argument; $x \neq 0, -1, -2, \ldots$
+- Returns `float`.
+
+### `lerch(z, s, a)`
+**Lerch transcendent** $\Phi(z, s, a) = \sum_{k=0}^\infty z^k / (a+k)^s$.
+
+- `z` (float): $|z| < 1$ (or $z = 1$, $s > 1$).
+- `s` (float): Exponent.
+- `a` (float): Offset; $a \neq 0, -1, -2, \ldots$
+- Returns `float`.
+
+### `gud(x)`
+**Gudermannian function** $\mathrm{gd}(x) = 2\arctan(\tanh(x/2))$.
+
+- `x` (float): Argument.
+- Returns `float`.
+
+Relates hyperbolic and trigonometric functions: $\sinh(x) = \tan(\mathrm{gd}(x))$.
+
+### `agud(g)`
+**Inverse Gudermannian function**.
+
+- `g` (float): Argument.
+- Returns `float`.
+
+### `normal_01_cdf_inverse(p)`
+**Inverse standard normal CDF** (probit function) $\Phi^{-1}(p)$.
+
+- `p` (float): Probability; $0 < p < 1$.
+- Returns `float`.
+
+### `cos_power_int(a, b, n)`
+**Cosine power integral** $\int_a^b \cos^n(t)\, dt$.
+
+- `a`, `b` (float): Integration limits.
+- `n` (int): Power.
+- Returns `float`.
+
+### `sin_power_int(a, b, n)`
+**Sine power integral** $\int_a^b \sin^n(t)\, dt$.
+
+- `a`, `b` (float): Integration limits.
+- `n` (int): Power.
+- Returns `float`.
+
+### `cardinal_cos(j, m, n, t, c)`
+**$j$-th cardinal cosine basis function** evaluated at $n$ points.
+
+- `j` (int): Basis index; $0 \leq j \leq m+1$.
+- `m` (int): Number of interior nodes.
+- `n` (int): Number of evaluation points.
+- `t` (ndarray[n]): Evaluation points.
+- `c` (ndarray[n]): Output values.
+
+Basis point $j$ has $C_j(T(j)) = 1$ and $C_j(T(i)) = 0$ for $i \neq j$, where
+$T(i) = \pi i / (m+1)$.
+
+### `cardinal_sin(j, m, n, t, s)`
+**$j$-th cardinal sine basis function** evaluated at $n$ points.
+
+- Same signature as `cardinal_cos`, with output `s` (ndarray[n]).
+
+---
+
+## Integer and Combinatorial Utilities
+
+| Routine | Parameters | Returns | Description |
+|---|---|---|---|
+| `i4_choose(n, k)` | `n`, `k` (int) | `int` | Binomial coefficient $\binom{n}{k}$ |
+| `r8_choose(n, k)` | `n`, `k` (int) | `float` | Real binomial coefficient |
+| `i4_factorial(n)` | `n` (int) | `int` | $n!$ |
+| `i4_factorial2(n)` | `n` (int) | `int` | Double factorial $n!! = n \cdot (n-2) \cdots$ |
+| `r8_factorial(n)` | `n` (int) | `float` | $n!$ as float |
+| `r8_factorial_log(n)` | `n` (int) | `float` | $\ln(n!)$ |
+| `trinomial(i, j, k)` | `i`, `j`, `k` (int) | `int` | Trinomial coefficient $(i+j+k)! / (i! j! k!)$ |
+| `commul(n, nfactor, factor, ncomb)` | `factor` (ndarray) | `ncomb` (int, out) | Multinomial coefficient $n! / \prod \text{factor}[i]!$ |
+| `poly_coef_count(dim, degree)` | `dim`, `degree` (int) | `int` | Number of monomials of total degree $\leq$ `degree` in `dim` variables |
+| `align_enum(m, n)` | `m`, `n` (int) | `int` | Number of alignments of sequences of length $m$ and $n$ |
+| `slice(dim_num, slice_num, piece_num)` | — | `piece_num` (int, out) | Max pieces from `slice_num` hyperplane cuts in `dim_num` dimensions |
+| `i4_is_triangular(i)` | `i` (int) | `bool` | Whether $i$ is a triangular number |
+| `i4_huge()` | (none) | `int` | Largest representable 32-bit integer |
+| `r8_huge()` | (none) | `float` | Largest representable double |
+| `i4_swap(i, j)` | `i`, `j` (int32, in/out) | — | Swaps two integers in-place |
+| `i4_to_triangle_lower(k, i, j)` | `k` (int) | `i`, `j` (int, out) | Convert linear index to lower-triangular $(i,j)$ coordinates |
+| `i4_to_triangle_upper(k, i, j)` | `k` (int) | `i`, `j` (int, out) | Convert linear index to upper-triangular $(i,j)$ coordinates |
+| `triangle_lower_to_i4(i, j, k)` | `i`, `j` (int) | `k` (int, out) | Convert lower-triangular $(i,j)$ to linear index |
+| `triangle_upper_to_i4(i, j, k)` | `i`, `j` (int) | `k` (int, out) | Convert upper-triangular $(i,j)$ to linear index |
+| `i4_partition_distinct_count(n, q)` | `n` (int) | `q` (int, out) | Number of partitions of $n$ into distinct parts |
+| `i4_uniform_ab(a, b, seed)` | `a`, `b` (int), `seed` (int32, in/out) | `int` | Uniform random integer in $[a, b]$ |
+
+---
+
+## Tabulated Values (Lookup Functions)
+
+These routines return reference values from pre-computed tables, useful for testing
+and verification. Each accepts an `n_data` index (pass 0 to start) and returns the
+corresponding input/output pair; `n_data` is incremented on each call. When the
+table is exhausted, `n_data` is reset to 0.
+
+| Routine | Table contents |
+|---|---|
+| `agm_values(n_data, a, b, fx)` | Arithmetic-Geometric Mean values |
+| `bell_values(n_data, n, c)` | Bell numbers |
+| `bernoulli_number_values(n_data, n, c)` | Bernoulli numbers |
+| `bernstein_poly_values(n_data, n, k, x, b)` | Bernstein basis values |
+| `beta_values(n_data, x, y, fxy)` | Beta function values |
+| `catalan_values(n_data, n, c)` | Catalan numbers |
+| `collatz_count_values(n_data, n, count)` | Collatz sequence lengths |
+| `cos_power_int_values(n_data, a, b, n, fx)` | Cosine power integrals |
+| `erf_values(n_data, x, fx)` | Error function values |
+| `gamma_values(n_data, x, fx)` | Gamma function values |
+| `gamma_log_values(n_data, x, fx)` | Log-Gamma values |
+| `gegenbauer_poly_values(n_data, n, a, x, fx)` | Gegenbauer polynomial values |
+| `gud_values(n_data, x, fx)` | Gudermannian values |
+| `hermite_poly_phys_values(n_data, n, x, fx)` | Hermite (physicist) polynomial values |
+| `hyper_2f1_values(n_data, a, b, c, x, fx)` | Hypergeometric ${}_{2}F_1$ values |
+| `i4_factorial_values(n_data, n, fn)` | Factorial values |
+| `i4_factorial2_values(n_data, n, fn)` | Double factorial values |
+| `jacobi_poly_values(n_data, n, a, b, x, fx)` | Jacobi polynomial values |
+| `laguerre_polynomial_values(n_data, n, x, fx)` | Laguerre polynomial values |
+| `lambert_w_values(n_data, x, fx)` | Lambert W function values |
+| `legendre_associated_values(n_data, n, m, x, fx)` | Associated Legendre values |
+| `legendre_associated_normalized_sphere_values(...)` | Normalized associated Legendre (sphere) |
+| `legendre_function_q_values(n_data, n, x, fx)` | Legendre Q function values |
+| `legendre_poly_values(n_data, n, x, fx)` | Legendre polynomial values |
+| `lerch_values(n_data, z, s, a, fx)` | Lerch transcendent values |
+| `mertens_values(n_data, n, c)` | Mertens function values |
+| `moebius_values(n_data, n, c)` | Moebius function values |
+| `normal_01_cdf_values(n_data, x, fx)` | Normal CDF values |
+| `omega_values(n_data, n, c)` | Omega function values |
+| `partition_distinct_count_values(n_data, n, c)` | Distinct-parts partition count values |
+| `phi_values(n_data, n, c)` | Euler totient values |
+| `psi_values(n_data, x, fx)` | Digamma function values |
+| `r8_factorial_log_values(n_data, n, fx)` | Log-factorial values |
+| `r8_factorial_values(n_data, n, fn)` | Real factorial values |
+| `sigma_values(n_data, n, c)` | Divisor sum values |
+| `sin_power_int_values(n_data, a, b, n, fx)` | Sine power integral values |
+| `spherical_harmonic_values(n_data, l, m, t, p, yr, yi)` | Spherical harmonic values |
+| `tau_values(n_data, n, c)` | Divisor count values |
+| `zeta_values(n_data, n, zeta)` | Riemann Zeta values |
diff --git a/docs/index.md b/docs/index.md
index 89ed28c..0de6168 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -1,45 +1,100 @@
# polpack
-**Special Functions and Recursively-Defined Polynomial Families for Python.**
+**Special Functions and Recursively-Defined Polynomial Families for Python**
---
-## What is polpack?
+## Overview
-`polpack` is a Python library for evaluating special functions and recursively-defined polynomial families. The numerical core is based on the original `POLPAK` library, providing efficient routines to evaluate a wide variety of mathematical functions.
+`polpack` is a Python library for evaluating special functions and recursively-defined polynomial families. The numerical core is based on the original `POLPAK` library by John Burkardt, providing efficient routines to evaluate a wide variety of mathematical functions used across numerical analysis, combinatorics, physics, and number theory.
-A **polynomial family** is a sequence of polynomials where each member is typically defined by its degree. Many such families are defined recursively, where higher-degree polynomials are computed from lower-degree ones. These functions are fundamental in numerical analysis, approximation theory, and physics.
+A **polynomial family** is a sequence of polynomials $\{P_n(x)\}_{n=0}^\infty$ where each member is typically defined by its degree. Many such families satisfy an **orthogonality condition** — the polynomials are mutually orthogonal under an inner product defined by a weight function over some interval. These families arise naturally in approximation theory, quadrature, quantum mechanics, and signal processing.
-`polpack` provides near-native performance by compiling its Fortran core via `f2py`, offering a clean and intuitive NumPy-based Python API.
+`polpack` evaluates these families using **three-term recurrence relations**, which are numerically stable and avoid the catastrophic cancellation errors that can occur with explicit-formula evaluation at high degree.
## Available polynomial families
-
-
| Family | Category | Type | Domain |
|---|---|---|---|
-| `bernoulli_poly` | Polynomial | Sequence | \( \mathbb{R} \) |
-| `bernstein_poly` | Polynomial | Basis | \( [0, 1] \) |
-| `cardan_poly` | Polynomial | Cubic | \( \mathbb{R} \) |
-| `charlier` | Polynomial | Discrete | \( \mathbb{N} \) |
-| `cheby_t_poly` | Polynomial | Orthogonal | \( [-1, 1] \) |
-| `gegenbauer_poly` | Polynomial | Orthogonal | \( [-1, 1] \) |
-| `jacobi_poly` | Polynomial | Orthogonal | \( [-1, 1] \) |
-| `laguerre_poly` | Polynomial | Orthogonal | \( [0, \infty) \) |
-| `legendre_poly` | Polynomial | Orthogonal | \( [-1, 1] \) |
+| `bernoulli_poly` | Polynomial | Sequence | $\mathbb{R}$ |
+| `bernstein_poly` | Polynomial | Basis | $[0, 1]$ |
+| `cardan_poly` | Polynomial | Cubic | $\mathbb{R}$ |
+| `charlier` | Polynomial | Discrete orthogonal | $\mathbb{N}$ |
+| `cheby_t_poly` | Polynomial | Orthogonal | $[-1, 1]$ |
+| `cheby_u_poly` | Polynomial | Orthogonal | $[-1, 1]$ |
+| `chebyshev_discrete` | Polynomial | Discrete orthogonal | $\mathbb{N}$ |
+| `gegenbauer_poly` | Polynomial | Orthogonal | $[-1, 1]$ |
+| `gen_hermite_poly` | Polynomial | Orthogonal | $\mathbb{R}$ |
+| `gen_laguerre_poly` | Polynomial | Orthogonal | $[0, \infty)$ |
+| `hermite_poly_phys` | Polynomial | Orthogonal | $\mathbb{R}$ |
+| `jacobi_poly` | Polynomial | Orthogonal | $[-1, 1]$ |
+| `krawtchouk` | Polynomial | Discrete orthogonal | $\{0,\ldots,M\}$ |
+| `laguerre_associated` | Polynomial | Orthogonal | $[0, \infty)$ |
+| `laguerre_poly` | Polynomial | Orthogonal | $[0, \infty)$ |
+| `legendre_associated` | Polynomial | Orthogonal | $[-1, 1]$ |
+| `legendre_poly` | Polynomial | Orthogonal | $[-1, 1]$ |
+| `meixner` | Polynomial | Discrete orthogonal | $\mathbb{N}$ |
+| `zernike_poly` | Polynomial | Orthogonal | unit disk |
*(See [API Reference](api.md) for the full list of supported families and routines.)*
-## Quick example
+## Requirements
+
+- Python 3.10+
+- [NumPy](http://www.numpy.org/)
+
+## Example Usage
```python
import numpy as np
import polpack
-# Example: Compute the first 11 Bell numbers
-# b[n] will contain the n-th Bell number
-b = np.zeros(11, dtype=np.int32, order="F")
-polpack.bell(10, b)
+# Chebyshev polynomials T_0 to T_5 at x = 0.5
+m, n, x = 1, 5, 0.5
+cx = np.zeros(n + 1, dtype=np.float64, order="F")
+polpack.cheby_t_poly(m, n, x, cx)
+print(f"T_0({x}) ... T_{n}({x}): {cx}")
-print(f"Bell numbers B_0 to B_10: {b}")
+# Legendre polynomials P_0 to P_10 and their derivatives at x = 0.0
+n, x = 10, 0.0
+p = np.zeros(n + 1, dtype=np.float64, order="F")
+dp = np.zeros(n + 1, dtype=np.float64, order="F")
+polpack.legendre_poly(n, x, p, dp)
+print(f"P_n(0.0) for n=0..{n}: {p}")
+
+# Bell numbers B_0 to B_14
+n = 14
+b = np.zeros(n + 1, dtype=np.int32, order="F")
+polpack.bell(n, b)
+print(f"Bell numbers B_0..B_{n}: {b}")
+
+# Riemann Zeta function at p = 2 (equals pi^2/6)
+z = polpack.zeta(2.0)
+print(f"zeta(2) = {z:.8f}")
```
+
+## Main Features
+
+1. **Orthogonal polynomial families** — Chebyshev (first and second kind), Gegenbauer, Hermite (physicist's and generalized), Jacobi, Laguerre (standard, associated, and generalized), Legendre (standard, associated, normalized), and Zernike polynomials.
+2. **Non-orthogonal polynomial families** — Bernoulli, Bernstein, Cardan, and Euler polynomials.
+3. **Discrete orthogonal polynomials** — Charlier, Chebyshev (discrete), Krawtchouk, and Meixner polynomials.
+4. **Polynomial coefficients and zeros** — exact coefficient arrays and Chebyshev node computation for all major families.
+5. **Combinatorial sequences** — Bell, Catalan, Delannoy, Motzkin, Stirling, Fibonacci, Vibonacci, Eulerian, and Zeckendorf.
+6. **Number-theoretic functions** — Euler totient, Moebius, Mertens, prime testing, prime factorization, divisor sums, Jacobi and Legendre symbols.
+7. **Special functions** — Riemann Zeta, Lambert W, Beta, error function (erf), Gamma, Lerch transcendent, Gudermannian, hypergeometric ${}_{2}F_1$, and more.
+8. **NumPy-based interface** — all array routines operate on pre-allocated Fortran-order NumPy arrays for zero-copy performance.
+
+## See Also
+
+- [SciPy special functions](https://docs.scipy.org/doc/scipy/reference/special.html) — complementary collection of special functions in Python
+- [mpmath](https://mpmath.org/) — arbitrary-precision special functions
+
+## Acknowledgements
+
+The author thanks [John Burkardt](https://jburkardt.github.io/) for the original [POLPAK](https://jburkardt.github.io/f_src/polpak/polpak.html) Fortran library, which provides the numerical core of `polpack`.
+
+## References
+
+- N. D. Cox, 1979, *Tolerance Analysis by Computer*, Journal of Quality Technology, Vol. 11, No. 2, pp. 80–87
+- Abramowitz, M., Stegun, I., *Handbook of Mathematical Functions*, National Bureau of Standards, 1964
+- Szego, G., *Orthogonal Polynomials*, American Mathematical Society, 1992
diff --git a/docs/installation.md b/docs/installation.md
index 8b1cf96..836e972 100644
--- a/docs/installation.md
+++ b/docs/installation.md
@@ -1,55 +1,67 @@
# Installation
-`polpack` is distributed as a compiled wheel on PyPI and can also be installed
-from source via GitHub.
+`polpack` can be installed from PyPI or directly from source via GitHub.
---
-## Prerequisites
+## [PyPI](https://pypi.org/project/polpack)
-- **Python 3.10+**
-- **NumPy** (installed automatically as a dependency)
+For using the PyPI package in your project, add it to your configuration file:
-For source builds you additionally need:
+=== "pyproject.toml"
-- A Fortran compiler (`gfortran` recommended)
-- `meson` and `meson-python` build system
-- `numpy` (for `f2py` compilation)
+ ```toml
+ [project.dependencies]
+ polpack = "*" # (1)!
+ ```
-## PyPI (recommended)
+ 1. Specifying a version is recommended
+
+=== "requirements.txt"
+
+ ```
+ polpack>=0.1.0
+ ```
### pip
-```bash
-pip install --upgrade polpack
-```
+=== "Installation for user"
-### pyproject.toml dependency
+ ```bash
+ pip install --upgrade --user polpack # (1)!
+ ```
-```toml
-[project]
-dependencies = [
- "polpack"
-]
-```
+ 1. You may need to use `pip3` instead of `pip` depending on your Python installation.
-### requirements.txt
+=== "Installation in virtual environment"
-```text
-polpack
-```
+ ```bash
+ python -m venv .venv
+ source .venv/bin/activate
+ pip install --require-virtualenv --upgrade polpack # (1)!
+ ```
+
+ 1. You may need to use `pip3` instead of `pip` depending on your Python installation.
-## Package managers
+ !!! note
+ The command to activate the virtual environment depends on your platform and shell.
+ [More info](https://docs.python.org/3/library/venv.html#how-venvs-work)
### uv
-```bash
-# Add to a uv project
-uv add polpack
+=== "Adding to uv project"
-# Or install into the current environment
-uv pip install polpack
-```
+ ```bash
+ uv add polpack
+ uv sync
+ ```
+
+=== "Installing to uv environment"
+
+ ```bash
+ uv venv
+ uv pip install polpack
+ ```
### pipenv
@@ -75,7 +87,9 @@ pdm add polpack
hatch add polpack
```
-## Installing from source (GitHub)
+---
+
+## [GitHub](https://github.com/eggzec/polpack)
Install the latest development version directly from the repository:
@@ -85,8 +99,7 @@ pip install --upgrade "git+https://github.com/eggzec/polpack.git#egg=polpack"
### Building locally
-Clone and build from source if you want to modify the Fortran code or test
-local changes:
+Clone and build from source if you want to modify or test local changes:
```bash
git clone https://github.com/eggzec/polpack.git
@@ -94,25 +107,29 @@ cd polpack
pip install -e .
```
-This invokes the `meson` build system to compile the Fortran sources via
-`f2py` and install the resulting extension module in development mode.
+This invokes the `meson` build system to compile the Fortran sources and install
+the resulting extension module in development mode.
!!! warning "Fortran compiler required"
- Source builds require a working Fortran compiler. On most Linux
- distributions install `gfortran`:
+ Source builds require a working Fortran compiler. On most Linux distributions,
+ install `gfortran`:
+
+ === "Debian / Ubuntu"
+ ```bash
+ sudo apt install gfortran
+ ```
+ === "Fedora"
+ ```bash
+ sudo dnf install gcc-gfortran
+ ```
+ === "macOS (Homebrew)"
+ ```bash
+ brew install gcc
+ ```
+ === "Windows"
+ Install [MinGW-w64](https://www.mingw-w64.org/) with gfortran or use MSYS2.
-```bash
-# Debian/Ubuntu
-sudo apt install gfortran
-
-# Fedora
-sudo dnf install gcc-gfortran
-
-# macOS (Homebrew)
-brew install gcc
-```
-
-On Windows, install MinGW-w64 with gfortran or use MSYS2.
+---
## Verifying the installation
@@ -128,10 +145,20 @@ polpack.bell(n, b)
print("polpack is working! Bell numbers:", b)
```
+---
+
## Dependencies
| Package | Purpose |
|---|---|
-| `numpy` | Array handling, `f2py` integration |
+| `numpy` | Array handling and `f2py` integration |
No other runtime dependencies are required.
+
+For source builds, you additionally need:
+
+| Package | Purpose |
+|---|---|
+| `gfortran` | Compiling Fortran sources |
+| `meson` and `meson-python` | Build system |
+| `numpy` | Required by `f2py` at build time |
diff --git a/docs/quickstart.md b/docs/quickstart.md
index 9d858c9..b783b46 100644
--- a/docs/quickstart.md
+++ b/docs/quickstart.md
@@ -1,151 +1,477 @@
# Quickstart
-This guide will help you get started with polpack in Python.
+This quickstart demonstrates how to use `polpack` to evaluate polynomial families,
+combinatorial sequences, and special functions. For mathematical background, see the
+[Theory](theory.md) section.
-## Working with Sequences
+---
-Beyond polynomial families, `polpack` provides routines for various integer
-sequences. These typically follow a similar pattern but return integers
-(or floats for very large values).
+## Common Interface
-### Bell Numbers
-The Bell numbers $B_n$ count the number of partitions of a set.
+Most `polpack` routines follow a consistent pattern: you pre-allocate a NumPy array,
+pass it to the routine along with parameters, and the routine fills it in place.
```python
import numpy as np
import polpack
-# Compute Bell numbers up to n=10
-n = 10
-b = np.zeros(n + 1)
-polpack.bell(n, b)
+# Pattern: pre-allocate output array, then call routine
+n = 5 # highest degree
+x = 0.5 # evaluation point
+cx = np.zeros(n + 1, dtype=np.float64, order="F") # (1)
-print(f"Bell numbers B_0 to B_{n}:")
-print(b)
+polpack.cheby_t_poly(1, n, x, cx)
+# cx[k] now contains T_k(x) for k = 0, 1, ..., n
```
+!!! note "Array order"
+ All output arrays must be pre-allocated. Use `order="F"` (Fortran-contiguous)
+ for consistency with the underlying Fortran core. Integer routines typically
+ use `dtype=np.int32`; floating-point routines use `dtype=np.float64`.
+
---
-## Basic Usage
+## Orthogonal Polynomials
+
+### Chebyshev Polynomials (First Kind)
-Evaluate the first 10 Bernoulli numbers:
+Chebyshev polynomials $T_n(x)$ are optimal for function approximation — they minimize
+the maximum absolute error among all degree-$n$ polynomials with leading coefficient 1.
```python
import numpy as np
import polpack
-n = 10
-# Initialize an array of size n+1 in Fortran order
-b = np.zeros(n + 1, dtype=np.float64, order="F")
-polpack.bernoulli_number(n, b)
-```
-
-This guide provides walk-through examples of how to use `polpack` to evaluate
-common polynomial families and combinatorial sequences in Python.
+m = 1 # number of evaluation points
+n = 6 # highest degree
+x = 0.5 # evaluation point
----
+cx = np.zeros(n + 1, dtype=np.float64, order="F")
+polpack.cheby_t_poly(m, n, x, cx)
-## Recursive Polynomials: Chebyshev First Kind
+print(f"Chebyshev T_k({x}) for k = 0 ... {n}:")
+for k, val in enumerate(cx):
+ print(f" T_{k}({x}) = {val:.6f}")
+```
-Chebyshev polynomials are fundamental in approximation theory and spectral
-methods. In `polpack`, you can evaluate the entire sequence $T_0(x), \dots, T_n(x)$
-efficiently in a single call.
+To get the polynomial coefficients or zeros of $T_n$:
```python
import numpy as np
import polpack
-# 1. Configuration
-m = 1 # Number of points
-n = 5 # Highest degree
-x = 0.5 # Evaluation point
-v = np.zeros(n + 1, dtype=np.float64, order="F")
-
-# 2. Evaluation
-polpack.cheby_t_poly(m, n, x, v)
+n = 4
+c = np.zeros((n + 1, n + 1), dtype=np.float64, order="F")
+polpack.cheby_t_poly_coef(n, c)
+print(f"Coefficients of T_{n}(x) (lowest to highest power):\n{c[n, :]}")
-# 3. Results
-print(f"Chebyshev T_0({x}) to T_{n}({x}):")
-for i, val in enumerate(v):
- print(f" T_{i}({x}) = {val:.4f}")
+z = np.zeros(n, dtype=np.float64, order="F")
+polpack.cheby_t_poly_zero(n, z)
+print(f"Zeros of T_{n}(x): {z}")
```
-## Orthogonal Polynomials: Legendre Polynomials
+### Legendre Polynomials
-Legendre polynomials are orthogonal on $[-1, 1]$ and arise naturally in
-potential theory and spherical harmonics.
+Legendre polynomials $P_n(x)$ are orthogonal on $[-1, 1]$ with unit weight. They arise
+in potential theory, Gauss-Legendre quadrature, and spherical harmonics.
```python
import numpy as np
import polpack
-# 1. Evaluate Legendre polynomials P_0 to P_10 at x=0.0
n = 10
-x = 0.0
+x = 0.5
p = np.zeros(n + 1, dtype=np.float64, order="F")
-# Output array for first derivatives (not used here)
-dp = np.zeros(n + 1, dtype=np.float64, order="F")
+dp = np.zeros(n + 1, dtype=np.float64, order="F") # first derivatives
polpack.legendre_poly(n, x, p, dp)
-# 2. Inspect results
-# P_n(0.0) = 0 for odd n
-print(f"Legendre P_n(0.0) for n=0..10:\n{p}")
+print(f"P_n({x}) for n = 0 ... {n}:")
+for k in range(n + 1):
+ print(f" P_{k}({x}) = {p[k]:.6f}, P'_{k}({x}) = {dp[k]:.6f}")
```
-## Combinatorial Sequences: Bell Numbers
+### Hermite Polynomials (Physicist's)
-Bell numbers $B_n$ count the number of partitions of a set of size $n$. These
-increase rapidly, so `polpack` typically uses `int32` or `float64` depending
-on the routine variant.
+Hermite polynomials $H_n(x)$ arise in quantum mechanics (harmonic oscillator) and
+probability (Gaussian integration).
```python
import numpy as np
import polpack
-# 1. Compute the first 15 Bell numbers
-n = 14
-b = np.zeros(n + 1, dtype=np.int32, order="F")
+n = 5
+x = 1.0
+cx = np.zeros(n + 1, dtype=np.float64, order="F")
-polpack.bell(n, b)
+polpack.hermite_poly_phys(n, x, cx)
+
+print(f"H_n({x}) for n = 0 ... {n}:")
+for k, val in enumerate(cx):
+ print(f" H_{k}({x}) = {val:.4f}")
+```
+
+### Jacobi Polynomials
+
+Jacobi polynomials $P_n^{(\alpha,\beta)}(x)$ generalize both Legendre ($\alpha=\beta=0$)
+and Chebyshev polynomials. They are orthogonal on $[-1, 1]$ with weight
+$(1-x)^\alpha (1+x)^\beta$.
+
+```python
+import numpy as np
+import polpack
-# Evaluate the 5th Legendre polynomial at x = 0.5
n = 5
-x = 0.5
+alpha, beta = 0.5, 0.5 # reduces to Gegenbauer with parameter 1.0
+x = 0.3
+cx = np.zeros(n + 1, dtype=np.float64, order="F")
+
+polpack.jacobi_poly(n, alpha, beta, x, cx)
+
+print(f"P_k^({alpha},{beta})({x}) for k = 0 ... {n}:")
+for k, val in enumerate(cx):
+ print(f" P_{k}({x}) = {val:.6f}")
+```
+
+### Laguerre Polynomials
+
+Laguerre polynomials $L_n(x)$ are orthogonal on $[0, \infty)$ with weight $e^{-x}$.
+They appear in the radial wavefunctions of the hydrogen atom.
+
+```python
+import numpy as np
+import polpack
+
+n = 6
+x = 1.5
cx = np.zeros(n + 1, dtype=np.float64, order="F")
-cpx = np.zeros(n + 1, dtype=np.float64, order="F") # for derivatives
-polpack.legendre_poly(n, x, cx, cpx)
+polpack.laguerre_poly(n, x, cx)
-print(f"P_5({x}) = {cx[n]}")
-print(f"P_5'({x}) = {cpx[n]}")
+print(f"L_n({x}) for n = 0 ... {n}:")
+for k, val in enumerate(cx):
+ print(f" L_{k}({x}) = {val:.6f}")
```
-### Coefficient calculation
+### Gegenbauer (Ultraspherical) Polynomials
+
+Gegenbauer polynomials $C_n^{(\alpha)}(x)$ generalize both Legendre ($\alpha=1/2$) and
+Chebyshev ($\alpha=1$) polynomials.
```python
import numpy as np
import polpack
-# Get coefficients of the 4th-order Chebyshev polynomial T_4(x)
n = 4
-# Space for (n+1) x (n+1) coefficient matrix
-c = np.zeros((n + 1, n + 1), dtype=np.float64, order="F")
-polpack.cheby_t_poly_coef(n, c)
+alpha = 1.0 # reduces to Chebyshev U when alpha = 1
+x = 0.5
+cx = np.zeros(n + 1, dtype=np.float64, order="F")
+
+polpack.gegenbauer_poly(n, alpha, x, cx)
+print(f"C_k^({alpha})({x}) for k = 0 ... {n}: {cx}")
+```
+
+---
+
+## Bernstein Polynomials
+
+Bernstein basis polynomials form the foundation of Bézier curves. All $n+1$ basis
+functions of degree $n$ are evaluated in a single call.
+
+```python
+import numpy as np
+import polpack
+
+n = 3 # degree
+x = 0.4 # point in [0, 1]
+bern = np.zeros(n + 1, dtype=np.float64, order="F")
+
+polpack.bernstein_poly(n, x, bern)
+
+print(f"Bernstein basis B_{{k,{n}}}({x}) for k = 0 ... {n}:")
+for k, val in enumerate(bern):
+ print(f" B_{k},{n}({x}) = {val:.6f}")
+
+# Verify partition of unity: sum should be 1.0
+print(f"Sum = {bern.sum():.6f}")
+```
+
+For evaluation on a general interval $[a, b]$:
+
+```python
+import numpy as np
+import polpack
+
+n = 3
+x, a, b = 2.0, 1.0, 4.0
+bern = np.zeros(n + 1, dtype=np.float64, order="F")
+
+polpack.bpab(n, x, a, b, bern)
+print(f"Bernstein basis on [{a}, {b}] at x={x}: {bern}")
+```
+
+---
+
+## Combinatorial Sequences
+
+### Bell Numbers
+
+The Bell number $B_n$ counts the number of ways to partition a set of $n$ elements.
+
+```python
+import numpy as np
+import polpack
+
+n = 14
+b = np.zeros(n + 1, dtype=np.int32, order="F")
+polpack.bell(n, b)
+
+print(f"Bell numbers B_0 ... B_{n}:")
+for k, val in enumerate(b):
+ print(f" B_{k} = {val}")
+```
+
+### Catalan Numbers
+
+The Catalan number $C_n$ counts binary trees, valid parenthesizations, and many other
+combinatorial structures.
+
+```python
+import numpy as np
+import polpack
+
+n = 10
+c = np.zeros(n + 1, dtype=np.int32, order="F")
+polpack.catalan(n, c)
+print(f"Catalan numbers C_0 ... C_{n}: {c}")
+```
+
+### Stirling Numbers
+
+Stirling numbers of the **first kind** $s(n, k)$ count permutations of $n$ objects with
+exactly $k$ cycles. Stirling numbers of the **second kind** $S(n, k)$ count the number
+of ways to partition $n$ objects into $k$ non-empty subsets.
+
+```python
+import numpy as np
+import polpack
+
+n, m = 6, 6
+s1 = np.zeros((n, m), dtype=np.int32, order="F")
+s2 = np.zeros((n, m), dtype=np.int32, order="F")
+
+polpack.stirling1(n, m, s1)
+polpack.stirling2(n, m, s2)
+
+print("Stirling numbers of the first kind s(n,k):")
+print(s1)
+print("\nStirling numbers of the second kind S(n,k):")
+print(s2)
+```
+
+### Fibonacci Numbers
+
+```python
+import numpy as np
+import polpack
+
+# Recursive: compute F_1 ... F_n
+n = 15
+f = np.zeros(n, dtype=np.int32, order="F")
+polpack.fibonacci_recursive(n, f)
+print(f"Fibonacci F_1 ... F_{n}: {f}")
+
+# Zeckendorf decomposition: represent N as sum of non-consecutive Fibonacci numbers
+n_val = 100
+m_max = 20
+m = np.zeros(1, dtype=np.int32, order="F")
+il = np.zeros(m_max, dtype=np.int32, order="F")
+fl = np.zeros(m_max, dtype=np.int32, order="F")
+polpack.zeckendorf(n_val, m_max, m, il, fl)
+terms = fl[: m[0]]
+print(f"{n_val} = {' + '.join(str(v) for v in terms)}")
+```
+
+### Delannoy and Motzkin Numbers
+
+```python
+import numpy as np
+import polpack
-print("Chebyshev T_4 coefficients (highest power first):")
-print(c[n, :])
+# Delannoy numbers: paths from (0,0) to (m,n) using steps (1,0), (0,1), (1,1)
+m, n = 5, 5
+a = np.zeros((m + 1, n + 1), dtype=np.int32, order="F")
+polpack.delannoy(m, n, a)
+print(f"Delannoy A(m,n):\n{a}")
+
+# Motzkin numbers
+n = 10
+mo = np.zeros(n + 1, dtype=np.int32, order="F")
+polpack.motzkin(n, mo)
+print(f"Motzkin numbers M_0 ... M_{n}: {mo}")
+```
+
+---
+
+## Number Theory
+
+### Euler Totient and Divisor Functions
+
+```python
+import numpy as np
+import polpack
+
+# Euler totient: phi(n) = number of integers in [1,n] coprime to n
+n = 12
+phin = np.zeros(1, dtype=np.int32, order="F")
+polpack.phi(n, phin)
+print(f"phi({n}) = {phin[0]}")
+
+# Sum of divisors sigma(n)
+sigma_n = np.zeros(1, dtype=np.int32, order="F")
+polpack.sigma(n, sigma_n)
+print(f"sigma({n}) = {sigma_n[0]}")
+
+# Number of divisors tau(n)
+taun = np.zeros(1, dtype=np.int32, order="F")
+polpack.tau(n, taun)
+print(f"tau({n}) = {taun[0]}")
```
-## Special functions
+### Prime Numbers and Factorization
-Evaluate functions like the Gamma function or Riemann Zeta function.
+```python
+import numpy as np
+import polpack
+
+# n-th prime (up to n=1600, largest prime stored is 13499)
+p10 = polpack.prime(10)
+print(f"The 10th prime is {p10}")
+
+# Primality test
+for candidate in [17, 18, 97, 100]:
+ result = polpack.i4_is_prime(candidate)
+ print(f"i4_is_prime({candidate}) = {result}")
+
+# Prime factorization: n = product(factor[i]^power[i]) * nleft
+n = 360
+factor_max = 20
+factor_num = np.zeros(1, dtype=np.int32, order="F")
+factor = np.zeros(factor_max, dtype=np.int32, order="F")
+power = np.zeros(factor_max, dtype=np.int32, order="F")
+nleft = np.zeros(1, dtype=np.int32, order="F")
+polpack.i4_factor(n, factor_max, factor_num, factor, power, nleft)
+nf = factor_num[0]
+parts = " * ".join(f"{factor[i]}^{power[i]}" for i in range(nf))
+print(f"{n} = {parts}")
+```
+
+### Moebius and Mertens Functions
```python
+import numpy as np
import polpack
-# Evaluate Riemann Zeta function at p = 2.0 (pi^2 / 6)
-p = 2.0
-z = polpack.zeta(p)
-print(f"Zeta({p}) = {z:.6f}")
+for n in [1, 6, 12, 13]:
+ mu = np.zeros(1, dtype=np.int32, order="F")
+ polpack.moebius(n, mu)
+ m = polpack.mertens(n)
+ print(f"mu({n}) = {mu[0]}, M({n}) = {m}")
+```
+
+---
+
+## Special Functions
+
+### Riemann Zeta Function
+
+```python
+import polpack
+import math
+
+for p in [2.0, 3.0, 4.0]:
+ z = polpack.zeta(p)
+ print(f"zeta({p}) = {z:.8f}")
+
+# Verification: zeta(2) = pi^2/6
+print(f"pi^2/6 = {math.pi**2 / 6:.8f}")
+```
+
+### Lambert W Function
+
+The Lambert W function satisfies $W(x) \cdot e^{W(x)} = x$.
+
+```python
+import polpack
+
+for x in [1.0, math.e, 10.0]:
+ w = polpack.lambert_w(x)
+ print(f"W({x:.4f}) = {w:.8f} [check: W*exp(W) = {w * math.exp(w):.8f}]")
+```
+
+### Error Function
+
+```python
+import polpack
+
+for x in [0.0, 0.5, 1.0, 2.0]:
+ e = polpack.r8_erf(x)
+ print(f"erf({x}) = {e:.8f}")
+```
+
+### Bernoulli Numbers
+
+```python
+import numpy as np
+import polpack
+
+n = 12
+b = np.zeros(n + 1, dtype=np.float64, order="F")
+polpack.bernoulli_number(n, b)
+
+print(f"Bernoulli numbers B_0 ... B_{n}:")
+for k, val in enumerate(b):
+ if val != 0.0:
+ print(f" B_{k} = {val:.10f}")
+```
+
+### Collatz Sequences
+
+```python
+import polpack
+
+for start in [6, 27, 100]:
+ steps = polpack.collatz_count(start)
+ print(f"Collatz({start}): {steps} steps to reach 1")
+
+# Find the integer in [1, N] with the longest Collatz sequence
+n = 1000
+i_max = np.zeros(1, dtype=np.int32, order="F")
+j_max = np.zeros(1, dtype=np.int32, order="F")
+polpack.collatz_count_max(n, i_max, j_max)
+print(f"In [1, {n}]: longest sequence starts at {i_max[0]} ({j_max[0]} steps)")
+```
+
+---
+
+## Spherical Harmonics
+
+`polpack` can evaluate the real and imaginary parts of the spherical harmonic
+function $Y_l^m(\theta, \phi)$:
+
+```python
+import numpy as np
+import polpack
+
+l, m = 2, 1
+theta = 0.5 # polar angle (radians)
+phi = 1.0 # azimuthal angle (radians)
+
+c = np.zeros(1, dtype=np.float64, order="F") # real part
+s = np.zeros(1, dtype=np.float64, order="F") # imaginary part
+
+polpack.spherical_harmonic(l, m, theta, phi, c, s)
+print(f"Y_{l}^{m}(theta={theta}, phi={phi})")
+print(f" Real part = {c[0]:.8f}")
+print(f" Imaginary part = {s[0]:.8f}")
```
diff --git a/docs/references.md b/docs/references.md
index 2821d5f..f4ce1a6 100644
--- a/docs/references.md
+++ b/docs/references.md
@@ -130,4 +130,4 @@
---
-*Last revised on 11 April 2015.*
+*Last revised on 2 April 2026.*
diff --git a/docs/theory.md b/docs/theory.md
index 3afb1ee..a1fc3f1 100644
--- a/docs/theory.md
+++ b/docs/theory.md
@@ -1,347 +1,485 @@
# Theory
-This page covers the mathematical foundations behind `polpack`: the definition
-of polynomial families, the stable recurrence relations used to evaluate them,
-and the combinatorial sequences supported by the library.
+This page covers the mathematical foundations of `polpack`: the definition and
+properties of orthogonal polynomial families, the three-term recurrence relations
+used to evaluate them, the combinatorial sequences supported by the library, and
+numerical stability considerations.
---
-## 1) Background: Special Functions and Polynomials
+## 1. Background: Special Functions and Polynomials
A **polynomial family** is a sequence of polynomials $\{P_n(x)\}_{n=0}^\infty$
-where $n$ is the degree of the polynomial. Many of these families are
+where $n$ is the degree. Many of the most important families are
**orthogonal polynomials**, meaning they satisfy a specific inner product
-relationship over an interval $[a, b]$ with respect to a weight function $w(x)$.
+condition over an interval $[a, b]$ with respect to a weight function $w(x) \geq 0$:
-Special functions and polynomials arise naturally across science and engineering:
+$$
+\int_a^b P_m(x)\, P_n(x)\, w(x)\, dx = 0, \qquad m \neq n.
+$$
+
+When $m = n$, this integral equals some positive constant $C_n$, called the
+norm squared of $P_n$.
+
+Special functions and orthogonal polynomials arise naturally across science and engineering:
-- **Physics** — quantum mechanics (Hermite, Laguerre), optics (Zernike), Laplace's equation (Legendre)
+- **Physics** — quantum mechanics (Hermite, Laguerre), electrostatics (Legendre), optics (Zernike)
- **Numerical Analysis** — function approximation (Chebyshev), Gaussian quadrature, spectral methods
-- **Combinatorics** — counting set partitions (Bell numbers), mountain ranges (Catalan numbers)
-- **Probability** — orthogonal basis expansion, stochastic processes
+- **Combinatorics** — counting set partitions (Bell numbers), lattice paths (Catalan, Delannoy)
+- **Number Theory** — arithmetic functions, prime distributions, divisor sums
-### Orthogonality condition
+### Orthogonality structure
-The structure in polynomial families is often characterized by:
+Each orthogonal polynomial family is characterized by:
-1. **Interval of support** — the range $[a, b]$ where the polynomials are defined
-2. **Weight function** $w(x)$ — a non-negative function defining the inner product
-3. **Orthogonality** — $\int_a^b P_m(x) P_n(x) w(x) \, dx = 0$ for $m \neq n$
-4. **Normalization** — the value of the integral for $m = n$ (often $C_n$)
+1. **Interval of support** $[a, b]$ — the domain over which orthogonality holds
+2. **Weight function** $w(x)$ — defines the inner product
+3. **Orthogonality** — $\langle P_m, P_n \rangle = \int_a^b P_m P_n w \, dx = 0$ for $m \neq n$
+4. **Normalization constant** $C_n = \int_a^b P_n^2 w \, dx$
-While sdepack focuses on sample paths of Wiener processes, `polpack` focuses on
-the stable evaluation of these deterministic bases.
+`polpack` focuses on the **stable numerical evaluation** of these polynomial families.
-### Explicit vs. Recursive interpretation
+### Explicit formulas vs. recurrence relations
-There are two principal ways to define and evaluate polynomial families:
+There are two principal ways to evaluate polynomial families:
-| | **Explicit Formulas** | **Recurrence Relations** |
+| | **Explicit Formulas** | **Three-Term Recurrence** |
|---|---|---|
-| Evaluation point | Direct calculation from $x$ | Computed step-wise from $P_0, P_1, \dots$ |
-| Numerical Stability | Often poor for large $n$ (cancellation errors) | Typically excellent and stable |
-| Implementation | Combinatorial sums, binomial coefficients | Three-term recurrence loops |
-| Typical use | Symbolic manipulation | Numerical computing, library implementations |
+| Evaluation | Direct from definition | Computed step-by-step from $P_0, P_1, \ldots$ |
+| Numerical stability | Often poor at high $n$ (cancellation errors) | Typically excellent and stable |
+| Implementation | Combinatorial sums, binomial coefficients | Simple loop |
+| Typical use | Symbolic manipulation | Numerical computing |
-`polpack` uses **recurrence relations** exclusively for its numerical core.
+`polpack` uses **three-term recurrence relations** exclusively for its numerical core.
!!! info "Christoffel–Darboux formula"
- For orthogonal polynomials, the sum of products follows a specific
- closed form:
+ For orthogonal polynomials, the sum of products satisfies a closed form:
- $$ \sum_{k=0}^n \frac{P_k(x)P_k(y)}{h_k} = \frac{k_n}{k_{n+1}h_n} \frac{P_{n+1}(x)P_n(y) - P_n(x)P_{n+1}(y)}{x - y} $$
+ $$
+ \sum_{k=0}^n \frac{P_k(x)P_k(y)}{h_k} = \frac{k_n}{k_{n+1} h_n}
+ \frac{P_{n+1}(x)P_n(y) - P_n(x)P_{n+1}(y)}{x - y}
+ $$
- The additional terms arising in this identity (analogous to Itô's lemma)
- are fundamental to understanding the convergence of spectral expansions.
+ This identity is fundamental to understanding the convergence of spectral
+ expansions and the completeness of orthogonal bases.
-## 2) Orthogonal Polynomial Families
+---
+
+## 2. Three-Term Recurrence
-All routines in `polpack` target specific polynomial forms:
+All orthogonal polynomial families satisfy a **three-term recurrence** of the form:
$$
-\int_a^b P_m(x) P_n(x) w(x) \, dx = \delta_{mn} C_n.
+P_{n+1}(x) = (A_n x + B_n)\, P_n(x) - C_n\, P_{n-1}(x),
$$
-Within the library, this is implemented as:
+starting from $P_0(x)$ and $P_1(x)$, where $A_n, B_n, C_n$ are
+family-specific recurrence coefficients. This structure is the basis for all
+evaluations in `polpack`.
+
+For **degree-invariant** recurrences, the coefficients $A$, $B$, $C$ are constants
+independent of $n$.
+
+---
+
+## 3. Orthogonal Polynomial Families
+
+### Chebyshev Polynomials
+
+**First Kind $T_n(x)$** — orthogonal on $[-1, 1]$ with weight $w(x) = (1-x^2)^{-1/2}$.
$$
-P_{n+1}(x) = (A_n x + B_n) P_n(x) - C_n P_{n-1}(x),
+T_{n+1}(x) = 2x\, T_n(x) - T_{n-1}(x), \qquad T_0 = 1,\; T_1 = x.
$$
-where:
+For $x \in [-1, 1]$, there is the explicit representation $T_n(x) = \cos(n \arccos x)$,
+which immediately gives $|T_n(x)| \leq 1$.
-- $x$ is the evaluation point
-- $A_n, B_n, C_n$ are method-specific recurrence coefficients
-- $k$ denotes the degree or index of the recurrence
+!!! info "Minimax property"
+ Among all polynomials of degree $n$ with leading coefficient 1, the scaled
+ Chebyshev polynomial $T_n(x) / 2^{n-1}$ has the smallest possible maximum absolute
+ value on $[-1,1]$. This makes Chebyshev nodes optimal for polynomial interpolation
+ and function approximation, minimizing the Runge phenomenon.
-For **degree-invariant** routines, the coefficients $A, B, C$ are pre-defined
-constants for a given degree $n$.
+**Second Kind $U_n(x)$** — orthogonal on $[-1, 1]$ with weight $w(x) = (1-x^2)^{1/2}$.
-## Polynomial families
+$$
+U_{n+1}(x) = 2x\, U_n(x) - U_{n-1}(x), \qquad U_0 = 1,\; U_1 = 2x.
+$$
-`polpack` evaluates a variety of recursively defined polynomial families,
-including those orthogonal with respect to specific weights and intervals.
+For $x \in [-1, 1]$: $U_n(x) = \sin((n+1)\arccos x) / \sqrt{1-x^2}$.
-### Bernoulli Polynomials
-The Bernoulli polynomials $B_n(x)$ are defined by the generating function:
-$\frac{t e^{xt}}{e^t - 1} = \sum_{n=0}^{\infty} B_n(x) \frac{t^n}{n!}$
-They satisfy the recurrence:
+**Norm:**
$$
-B_n(x) = \sum_{k=0}^n \binom{n}{k} B_k x^{n-k}
+\int_{-1}^{1} T_m(x)\, T_n(x)\, (1-x^2)^{-1/2}\, dx =
+\begin{cases}
+\pi & m = n = 0 \\
+\pi/2 & m = n \neq 0 \\
+0 & m \neq n
+\end{cases}
$$
-### Bernstein Polynomials
-The Bernstein basis polynomials $b_{\nu, n}(x)$ are defined on $[0,1]$ as:
-$b_{\nu, n}(x) = \binom{n}{\nu} x^{\nu} (1-x)^{n-\nu}$
-They form a basis for the space of polynomials of degree at most $n$.
+### Legendre Polynomials
-### Cardan Polynomials
-Cardan polynomials $C_n(x, s)$ are associated with the solution of cubic
-equations. They satisfy the three-term recurrence:
+Orthogonal on $[-1, 1]$ with unit weight $w(x) = 1$.
$$
-C_{n+1}(x, s) = x C_n(x, s) - s C_{n-1}(x, s)
+(n+1)\, P_{n+1}(x) = (2n+1)\, x\, P_n(x) - n\, P_{n-1}(x), \qquad P_0 = 1,\; P_1 = x.
$$
-### Charlier Polynomials
-The Charlier polynomials $C_n(x; a)$ are orthogonal with respect to the
-Poisson distribution. They satisfy:
+**Norm:**
$$
-a C_{n+1}(x; a) = (n + a - x) C_n(x; a) - n C_{n-1}(x; a)
+\int_{-1}^1 P_m(x)\, P_n(x)\, dx = \frac{2}{2n+1}\, \delta_{mn}.
$$
-### Chebyshev Polynomials
-**First Kind ($T_n$):** Orthogonal on $[-1, 1]$ with weight $(1-x^2)^{-1/2}$.
+**Properties:**
+- $P_n(1) = 1$, $P_n(-1) = (-1)^n$, $|P_n(x)| \leq 1$ on $[-1,1]$.
+- The zeros of $P_n(x)$ are the Gauss-Legendre quadrature nodes.
+- $P_{n+1}'(x) - P_{n-1}'(x) = (2n+1)\, P_n(x)$.
+
+The **associated Legendre functions** $P_n^m(x)$ generalize $P_n$ for $0 \leq m \leq n$:
$$
-T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)
+P_n^m(x) = (-1)^m (1-x^2)^{m/2} \frac{d^m}{dx^m} P_n(x).
$$
-**Second Kind ($U_n$):** Orthogonal on $[-1, 1]$ with weight $(1-x^2)^{1/2}$.
+They satisfy a modified recurrence and arise in spherical harmonic expansions.
+
+### Hermite Polynomials
+
+**Physicist's form $H_n(x)$** — orthogonal on $(-\infty, \infty)$ with weight $e^{-x^2}$.
$$
-U_{n+1}(x) = 2x U_n(x) - U_{n-1}(x)
+H_{n+1}(x) = 2x\, H_n(x) - 2n\, H_{n-1}(x), \qquad H_0 = 1,\; H_1 = 2x.
$$
-### Gegenbauer Polynomials
-The Gegenbauer polynomials $C_n^{(\alpha)}(x)$ are orthogonal on $[-1, 1]$ with
-weight $(1-x^2)^{\alpha - 1/2}$.
+**Norm:**
$$
-(n+1)C_{n+1}^{(\alpha)}(x) = 2(n+\alpha)x C_n^{(\alpha)}(x) - (n+2\alpha-1)C_{n-1}^{(\alpha)}(x)
+\int_{-\infty}^{\infty} H_m(x)\, H_n(x)\, e^{-x^2}\, dx = \sqrt{\pi}\, 2^n\, n!\; \delta_{mn}.
$$
-### Hermite Polynomials (Physicist's)
-Orthogonal on $(-\infty, \infty)$ with weight $e^{-x^2}$.
+**Probabilist's form** $He_n(x)$ uses the weight $e^{-x^2/2}$ instead:
$$
-H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x)
+He_{n+1}(x) = x\, He_n(x) - n\, He_{n-1}(x).
$$
+The **generalized Hermite polynomial** $H_n^{(\mu)}(x)$ is orthogonal under the weight
+$|x|^{2\mu} e^{-x^2}$ and reduces to $H_n$ when $\mu = 0$.
+
### Jacobi Polynomials
-The Jacobi polynomials $P_n^{(\alpha, \beta)}(x)$ are orthogonal on $[-1, 1]$
-with weight $(1-x)^{\alpha}(1+x)^{\beta}$. They generalize both Chebyshev and
-Legendre polynomials.
-### Krawtchouk Polynomials
-Discrete orthogonal polynomials $K_n(x; p, N)$ related to the binomial
-distribution.
+The most general family considered here. Orthogonal on $[-1, 1]$ with weight
+$w(x) = (1-x)^\alpha (1+x)^\beta$, where $\alpha, \beta > -1$.
+
+$$
+P_0^{(\alpha,\beta)} = 1, \qquad P_1^{(\alpha,\beta)} = \tfrac{1}{2}(2 + \alpha + \beta)\, x + \tfrac{1}{2}(\alpha - \beta),
+$$
+
+$$
+P_{n+1}^{(\alpha,\beta)}(x) = \frac{(2n+\alpha+\beta-1)\bigl[(\alpha^2-\beta^2) + (2n+\alpha+\beta)(2n+\alpha+\beta-2)x\bigr]}{2n(n+\alpha+\beta)(2n+\alpha+\beta-2)}\, P_n^{(\alpha,\beta)}
+$$
+$$
+\quad - \frac{2(n-1+\alpha)(n-1+\beta)(2n+\alpha+\beta)}{2n(n+\alpha+\beta)(2n+\alpha+\beta-2)}\, P_{n-1}^{(\alpha,\beta)}.
+$$
+
+**Special cases:**
+- $\alpha = \beta = 0$: Legendre polynomials $P_n(x)$
+- $\alpha = \beta = -1/2$: Chebyshev $T_n$
+- $\alpha = \beta = 1/2$: Chebyshev $U_n$
+- $\alpha = \beta$: Gegenbauer polynomials
+
+### Gegenbauer (Ultraspherical) Polynomials
+
+Orthogonal on $[-1, 1]$ with weight $w(x) = (1-x^2)^{\alpha - 1/2}$, $\alpha > -1/2$.
+
+$$
+(n+1)\, C_{n+1}^{(\alpha)}(x) = 2(n+\alpha)\, x\, C_n^{(\alpha)}(x) - (n+2\alpha-1)\, C_{n-1}^{(\alpha)}(x),
+$$
+
+$$
+C_0^{(\alpha)} = 1, \qquad C_1^{(\alpha)} = 2\alpha\, x.
+$$
+
+**Special cases:** $\alpha = 1$ gives Chebyshev $U_n$; $\alpha = 1/2$ gives Legendre $P_n$.
### Laguerre Polynomials
-**Associated Laguerre ($L_n^{(\alpha)}$):** Orthogonal on $[0, \infty)$ with
-weight $x^{\alpha}e^{-x}$.
+
+**Standard** $L_n(x)$ — orthogonal on $[0, \infty)$ with weight $e^{-x}$.
$$
-(n+1)L_{n+1}^{(\alpha)}(x) = (2n + \alpha + 1 - x) L_n^{(\alpha)}(x) - (n+\alpha)L_{n-1}^{(\alpha)}(x)
+(n+1)\, L_{n+1}(x) = (2n+1-x)\, L_n(x) - n\, L_{n-1}(x), \qquad L_0 = 1,\; L_1 = 1-x.
$$
-### Legendre Polynomials
-Orthogonal on $[-1, 1]$ with weight 1.
+**Associated** $L_n^m(x)$ — orthogonal on $[0, \infty)$ for fixed $m \geq 0$:
$$
-(n+1)P_{n+1}(x) = (2n+1)x P_n(x) - n P_{n-1}(x)
+(n+1)\, L_{n+1}^m(x) = (2n+m+1-x)\, L_n^m(x) - (n+m)\, L_{n-1}^m(x).
$$
-### Meixner Polynomials
-Discrete orthogonal polynomials $M_n(x; \beta, c)$ related to the negative
-binomial distribution.
+**Generalized** $L_n^{(\alpha)}(x)$ — for real $\alpha > -1$:
-### Zernike Polynomials
-Polynomials defined on the unit disk, widely used in optics to describe
-wavefront aberrations.
+$$
+(n+1)\, L_{n+1}^{(\alpha)}(x) = (2n+\alpha+1-x)\, L_n^{(\alpha)}(x) - (n+\alpha)\, L_{n-1}^{(\alpha)}(x).
+$$
+
+### Spherical Harmonics
+
+The spherical harmonic $Y_l^m(\theta, \phi)$ is the angular part of solutions to
+Laplace's equation in spherical coordinates:
+
+$$
+Y_l^m(\theta, \phi) = \sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}}\, P_l^m(\cos\theta)\, e^{im\phi},
+$$
+
+where $P_l^m$ is the associated Legendre function. `polpack` returns the real and
+imaginary parts separately.
---
-## Combinatorial Sequences
+## 4. Discrete Orthogonal Polynomials
-`polpack` provides routines for several important integer sequences.
+These families are orthogonal with respect to a **discrete** measure (a weight function
+defined on a finite or countable set of points).
-### Bell Numbers
-The Bell number $B_n$ counts the number of ways a set of $n$ elements can be
-partitioned into non-empty subsets.
+### Charlier Polynomials $C_n(x; a)$
+
+Orthogonal with respect to the Poisson distribution with parameter $a > 0$:
$$
-B_{n+1} = \sum_{k=0}^n \binom{n}{k} B_k
+a\, C_{n+1}(x; a) = (n + a - x)\, C_n(x; a) - n\, C_{n-1}(x; a).
$$
-### Catalan Numbers
-The Catalan number $C_n$ occurs in various counting problems (e.g., number
-of binary trees with $n$ nodes).
+### Krawtchouk Polynomials $K_n(x; p, M)$
+
+Orthogonal with respect to the binomial distribution $\text{Bin}(M, p)$, for $0 < p < 1$:
$$
-C_n = \frac{1}{n+1} \binom{2n}{n}
+(n+1)\, K_{n+1} = (x - n - p(M - 2n))\, K_n - (M-n+1)\, p(1-p)\, K_{n-1}.
$$
-### Stirling Numbers
-**First Kind ($s(n,k)$):** Number of permutations of $n$ elements with $k$ cycles.
-**Second Kind ($S(n,k)$):** Number of ways to partition a set of $n$ elements
-into $k$ non-empty subsets.
+### Meixner Polynomials $M_n(x; \beta, c)$
+
+Orthogonal with respect to the negative binomial distribution, for $\beta > 0$ and $0 < c < 1$:
+
+$$
+(n+1)\, M_{n+1} = ((n + \beta)c + n - (n + \beta)c x/(1-c))\, M_n
+ - n\, M_{n-1} \cdot \text{(recurrence coefficients)}.
+$$
+
+### Discrete Chebyshev Polynomials $t_n(x; N)$
+
+Orthogonal on the grid $\{0, 1, \ldots, N-1\}$ with uniform weights.
---
-## Number Theory and Probability
+## 5. Non-Orthogonal Polynomial Families
+
+### Bernoulli Polynomials
+
+Defined by the generating function $\frac{te^{xt}}{e^t - 1} = \sum_{n=0}^\infty B_n(x) \frac{t^n}{n!}$:
+
+$$
+B_n(x) = \sum_{k=0}^n \binom{n}{k} B_k\, x^{n-k},
+$$
-### Benford's Law
-Benford's law, also called the Newcomb–Benford law or the law of anomalous
-numbers, defines the probability distribution of the first digit $d$ in many
-real-life sets of numerical data:
+where $B_k$ are the Bernoulli numbers. Key properties:
$$
-P(d) = \log_{10}(1 + 1/d)
+B_n'(x) = n\, B_{n-1}(x), \qquad
+B_n(x+1) - B_n(x) = n\, x^{n-1}, \qquad
+B_n(0) = B_n(1) = B_n \text{ (for } n \neq 1\text{)}.
$$
-### Collatz Conjecture
-The Collatz sequence for a starting integer $n$ is defined by:
+### Euler Polynomials
+
+Defined by $\frac{2e^{xt}}{e^t + 1} = \sum_{n=0}^\infty E_n(x) \frac{t^n}{n!}$:
-- If $n$ is even, $n_{i+1} = n_i / 2$
-- If $n$ is odd, $n_{i+1} = 3n_i + 1$
-`polpack` provides routines to count the number of steps to reach 1.
+$$
+E_n'(x) = n\, E_{n-1}(x), \qquad E_n(1/2) = E_n / 2^n,
+$$
-### Arithmetical Functions
+where $E_n$ are the Euler numbers.
-- **Euler Totient ($\phi(n)$):** Number of positive integers less than $n$
- that are relatively prime to $n$.
-- **Moebius ($\mu(n)$):** $\mu(n) = (-1)^k$ if $n$ is square-free with $k$ prime
- factors, else 0.
-- **Mertens ($M(n)$):** The sum of the Moebius function $M(n) = \sum_{k=1}^n \mu(k)$.
-- **Divisor Sum ($\sigma(n)$):** Sum of all positive divisors of $n$.
+### Bernstein Polynomials
-## 3) Three-term recurrence and evaluation scaling
+The Bernstein basis polynomials of degree $n$ on $[0, 1]$ are:
-Given an interval $[a, b]$ and degree $N$:
$$
-P_0(x) = \alpha, \qquad P_1(x) = \beta x + \gamma.
+b_{k,n}(x) = \binom{n}{k} x^k (1-x)^{n-k}, \qquad 0 \leq k \leq n.
$$
-At each stage $n$, the polynomial value is generated:
+They satisfy the **partition of unity**:
+
$$
-P_{n+1}(x) = (A_n x + B_n)P_n(x) - C_n P_{n-1}(x),
+\sum_{k=0}^n b_{k,n}(x) = 1 \quad \text{for all } x.
$$
-where $A_n, B_n, C_n$ are coefficients that ensure the correct scaling at each
-stage. When multiplied in the core loops, this produces the high-order
-polynomial values required for scientific applications.
+These polynomials form the basis of Bézier curves and are used extensively in
+computer-aided geometric design (CAGD).
-## 4) Chebyshev Polynomials — RK1 analog
+### Cardan Polynomials
-The simplest and most widely used family in numerical analysis. Named after
-Pafnuty Chebyshev, they are the deterministic analog of many foundational
-integrators.
+Cardan polynomials $C_n(x, s)$ are associated with the solution of cubic equations
+and satisfy the three-term recurrence:
-For both first and second kind variants:
+$$
+C_{n+1}(x, s) = x\, C_n(x, s) - s\, C_{n-1}(x, s), \qquad C_0 = 2,\; C_1 = x.
+$$
+
+### Zernike Polynomials
+
+Defined on the unit disk, Zernike polynomials $R_n^m(\rho)$ are used to describe
+wavefront aberrations in optics. They satisfy:
$$
-T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x),
+R_n^m(\rho) = 0 \quad \text{if } m < 0,\; n < 0, \text{ or } n < m,
+\qquad R_m^m(\rho) = \rho^m,
$$
+and the recurrence:
+
$$
-T_0(x) = 1,\quad T_1(x)=x.
+R_{n+2}^m = A_n \bigl[(B_n \rho^2 - C_n)\, R_n^m - D_n\, R_{n-2}^m\bigr],
$$
-This recurrence is notoriously stable for $x \in [-1, 1]$.
-!!! info "Minimax property"
- The Chebyshev polynomials minimize the maximum error in function
- approximation, analogously to how Euler-Maruyama minimizes one-step
- error in SDEs.
+where $A_n, B_n, C_n, D_n$ depend on $n$ and $m$.
-## 5) Legendre Polynomials — Two-stage recurrence
+---
-A classic family providing improved accuracy in spherical harmonic expansions:
+## 6. Combinatorial Sequences
+
+### Bell Numbers
+
+The Bell number $B_n$ counts the number of set partitions of $\{1, \ldots, n\}$:
$$
-\begin{aligned}
-A_{n} &= \frac{2n+1}{n+1},\\
-B_{n} &= 0,\\
-C_{n} &= \frac{n}{n+1},\\
-(n+1)P_{n+1}(x) &= (2n+1)x P_n(x) - n P_{n-1}(x).
-\end{aligned}
+B_{n+1} = \sum_{k=0}^n \binom{n}{k} B_k, \qquad B_0 = 1.
$$
-Coefficients:
+Bell numbers grow super-exponentially: $B_0 = 1, B_1 = 1, B_2 = 2, B_3 = 5, B_4 = 15, \ldots$
+
+### Catalan Numbers
+
+The Catalan number $C_n$ counts binary trees, valid parenthesizations, triangulations
+of a polygon with $n+2$ vertices, and many other structures:
+
$$
-P_0=1,\quad P_1=x.
+C_n = \frac{1}{n+1}\binom{2n}{n} = \frac{(2n)!}{(n+1)!\, n!}, \qquad
+C_{n+1} = \frac{2(2n+1)}{n+2}\, C_n.
$$
----
+### Stirling Numbers
+
+**First kind** $s(n, k)$: the number of permutations of $n$ objects with exactly $k$ cycles.
+Signed: $s(n, k) = (-1)^{n-k} |s(n, k)|$.
+
+**Second kind** $S(n, k)$: the number of ways to partition $n$ objects into $k$ non-empty subsets.
+
+They are related to Bell numbers by $B_n = \sum_{k=1}^n S(n, k)$.
+
+### Fibonacci and Zeckendorf
+
+The Fibonacci sequence satisfies $F_{n+1} = F_n + F_{n-1}$, with $F_0 = 0$, $F_1 = 1$.
+The ratio $F_{n+1}/F_n$ converges to the golden ratio $\phi = (1 + \sqrt{5})/2$.
+
+**Zeckendorf's theorem**: every positive integer has a unique representation as a sum of
+non-consecutive Fibonacci numbers.
+
+### Eulerian Numbers
+
+The Eulerian number $E(n, k)$ counts the permutations of $\{1, \ldots, n\}$ with exactly
+$k$ **ascents** (positions $i$ where $\sigma(i) < \sigma(i+1)$):
+
+$$
+E(n, k) = (k+1)\, E(n-1, k) + (n-k)\, E(n-1, k-1).
+$$
+
+### Delannoy and Motzkin Numbers
+
+The **Delannoy number** $D(m, n)$ counts lattice paths from $(0,0)$ to $(m,n)$ using
+steps east, north, and northeast:
+
+$$
+D(m, n) = D(m-1, n) + D(m, n-1) + D(m-1, n-1).
+$$
-## 6) Hermite Polynomials — physicist form (time-invariant only)
+The **Motzkin number** $M_n$ counts paths from $(0,0)$ to $(0,n)$ that never go below
+the $x$-axis, using steps $\pm 1$ and $0$:
-A three-term recurrence providing higher accuracy for systems involving
-Gaussian weights:
$$
-H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x),
+(n+2)\, M_n = (2n+2)\, M_{n-1} + 3\, M_{n-2}.
$$
-with coefficients:
+---
+
+## 7. Number-Theoretic Functions
+
+### Euler Totient $\phi(n)$
+
+$\phi(n)$ counts the integers in $[1, n]$ coprime to $n$. For the prime factorization
+$n = \prod p_i^{e_i}$:
+
$$
-\begin{aligned}
-H_0&=1,\\
-H_1&=2x.
-\end{aligned}
+\phi(n) = n \prod_{p \mid n} \left(1 - \frac{1}{p}\right).
$$
-!!! note
- The probabilist's form $He_n(x)$ uses different scaling parameters. Both
- are available in `polpack`.
+### Moebius Function $\mu(n)$
-## 7) Jacobi and Gegenbauer Polynomials
+$$
+\mu(n) = \begin{cases}
+1 & n = 1 \\
+(-1)^k & n \text{ is a product of } k \text{ distinct primes} \\
+0 & n \text{ has a squared prime factor}
+\end{cases}
+$$
-The highest-order methods in `polpack`, using four terms or generalized
-parameters for maximum flexibility.
+### Mertens Function $M(n)$
-The general update formula is:
$$
-P_{n+1}(x) = (A_n x + B_n) P_n(x) - C_n P_{n-1}(x).
+M(n) = \sum_{k=1}^n \mu(k).
$$
-The Jacobi and Gegenbauer variants use **different** coefficient sets
-optimized for their respective weighting intervals. See the [API Reference](api.md)
-for the complete recurrence tables.
+The Riemann hypothesis is equivalent to $|M(n)| = O(n^{1/2+\varepsilon})$ for all $\varepsilon > 0$.
+
+### Divisor Functions
+
+- **$\sigma(n)$**: sum of all positive divisors of $n$; $\sigma(p^k) = (p^{k+1}-1)/(p-1)$
+- **$\tau(n)$**: number of positive divisors of $n$; $\tau(p^k) = k+1$
+- **$\omega(n)$**: number of distinct prime divisors of $n$
+
+---
+
+## 8. Convergence and Numerical Stability
+
+### Forward and Backward Stability
-## 8) Convergence and stability considerations
+Two notions of stability apply to recurrence-based polynomial evaluation:
-### Strong vs. weak stability
+- **Forward stability** — accuracy of individual evaluations: $|\hat{P}_n(x) - P_n(x)|$ is small.
+- **Backward stability** — the computed $\hat{P}_n(x)$ is the exact solution to a slightly
+ perturbed recurrence.
-When assessing the quality of a numerical polynomial evaluation, two notions of
-stability are distinguished:
+Forward stability is sufficient for function approximation; backward stability matters
+in eigenvalue and quadrature computations.
-- **Forward stability** — accuracy of individual evaluations.
-- **Backward stability** — whether the computed $P_n(x)$ is the exact solution
- to a slightly perturbed recurrence.
+### Stability of Three-Term Recurrences
-Forward stability is often sufficient for function approximation, while
-backward stability matters in eigenvalue problems.
+The three-term recurrence in `polpack` is **numerically stable** for the standard families
+within their natural domains:
-### Relationship between degree and order
+- Chebyshev $T_n(x)$: stable for $x \in [-1, 1]$; $|T_n(x)| \leq 1$ bounds growth.
+- Legendre $P_n(x)$: stable for $x \in [-1, 1]$.
+- Hermite $H_n(x)$: stable but values grow as $\sim 2^n$; use with care for large $n$.
+- Laguerre $L_n(x)$: stable for $x \geq 0$.
-The three-term recurrence schemes in `polpack` are designed to achieve higher
-numerical precision. The Euler-like first-order recursions (RK1 analog) are
-stable but may accumulate error at very high degrees ($n > 1000$). Higher-degree
-polynomials generally benefit from the stable recurrence architecture.
+For very high degrees ($n > 1000$), accumulated floating-point rounding can become
+significant. In such cases, specialized asymptotic approximations may be preferable.
!!! tip "Evaluation guidance"
- Unlike simple power series, the recurrence $P_{n+1}(x)$ only reduces the
- error by a factor of 2 per degree. For very large $n$, use specialized
- asymptotic routines if available.
+ Unlike direct power series evaluation, the three-term recurrence reduces cancellation
+ error by building up the polynomial iteratively from $P_0$ and $P_1$. For most
+ practical applications ($n \leq 100$), the recurrence is highly accurate.
diff --git a/src/polpack/__init__.py b/src/polpack/__init__.py
index 543d3e9..3667c18 100644
--- a/src/polpack/__init__.py
+++ b/src/polpack/__init__.py
@@ -387,14 +387,16 @@ def cheby_u_poly_zero(n: int, z: np.ndarray) -> None:
return _polpack.cheby_u_poly_zero(n, z)
-def chebyshev_discrete(n, m, x, v):
+def chebyshev_discrete(n: int, m: int, x: float, v: np.ndarray) -> None:
"""Evaluates discrete Chebyshev polynomials at a point.
Args:
- n (int): Description for n.
- m (int): Description for m.
- x (float): Description for x.
- v (float): Description for v.
+ n (int): Highest degree.
+ m (int): Grid size (number of discrete points).
+ x (float): Evaluation point.
+ v (np.ndarray): Output array of size n+1.
+ Returns:
+ None
"""
return _polpack.chebyshev_discrete(n, m, x, v)
@@ -484,13 +486,15 @@ def cos_power_int_values(
return _polpack.cos_power_int_values(n_data, a, b, n, fx)
-def delannoy(m, n, a):
+def delannoy(m: int, n: int, a: np.ndarray) -> None:
"""Returns the Delannoy numbers up to orders (M,N).
Args:
- m (int): maximum order m.
- n (int): maximum order n.
- a (ndarray): array to store the Delannoy numbers.
+ m (int): Maximum row order.
+ n (int): Maximum column order.
+ a (np.ndarray): Output array of shape (m+1, n+1).
+ Returns:
+ None
"""
return _polpack.delannoy(m, n, a)
@@ -545,11 +549,11 @@ def eulerian(n: int, e: np.ndarray) -> None:
return _polpack.eulerian(n, e)
-def fibonacci_direct(n):
- """Computes the N-th Fibonacci number directly.
+def fibonacci_direct(n: int) -> int:
+ """Computes the N-th Fibonacci number directly via the Binet formula.
Args:
- n (int): index of the Fibonacci number.
+ n (int): Index of the Fibonacci number.
Returns:
int: N-th Fibonacci number.
"""
@@ -718,7 +722,7 @@ def hermite_poly_phys_values(n_data: int, n: int, x: float, fx: float) -> None:
return _polpack.hermite_poly_phys_values(n_data, n, x, fx)
-def hyper_2f1_values(
+def hyper_2f1_values( # noqa: PLR0913, PLR0917
n_data: int, a: float, b: float, c: float, x: float, fx: float
) -> None:
"""Returns some values of the hypergeometric function 2F1.
@@ -736,7 +740,7 @@ def hyper_2f1_values(
return _polpack.hyper_2f1_values(n_data, a, b, c, x, fx)
-def i4_factor(
+def i4_factor( # noqa: PLR0913, PLR0917
n: int,
factor_max: int,
factor_num: np.ndarray,
@@ -849,7 +853,7 @@ def i4mat_print(m: int, n: int, a: np.ndarray, title: str) -> None:
return _polpack.i4mat_print(m, n, a, title)
-def i4mat_print_some(
+def i4mat_print_some( # noqa: PLR0913, PLR0917
m: int,
n: int,
a: np.ndarray,
@@ -893,7 +897,7 @@ def jacobi_poly(
return _polpack.jacobi_poly(n, alpha, beta, x, cx)
-def jacobi_poly_values(
+def jacobi_poly_values( # noqa: PLR0913, PLR0917
n_data: int, n: int, a: float, b: float, x: float, fx: float
) -> None:
"""Returns some values of the Jacobi polynomial.
@@ -1534,7 +1538,7 @@ def slice_fn(dim_num: int, slice_num: int, piece_num: np.ndarray) -> None:
return _polpack.slice(dim_num, slice_num, piece_num)
-def spherical_harmonic(
+def spherical_harmonic( # noqa: PLR0913, PLR0917
ell: int, m: int, theta: float, phi: float, c: np.ndarray, s: np.ndarray
) -> None:
"""Evaluates spherical harmonic functions.
@@ -1552,19 +1556,31 @@ def spherical_harmonic(
return _polpack.spherical_harmonic(ell, m, theta, phi, c, s)
-def spherical_harmonic_values(n_data, l, m, theta, phi, yr, yi):
- """Returns values of spherical harmonic functions.
+def spherical_harmonic_values( # noqa: PLR0913, PLR0917
+ n_data: int,
+ ell: int,
+ m: int,
+ theta: float,
+ phi: float,
+ yr: float,
+ yi: float,
+) -> None:
+ """Returns tabulated values of spherical harmonic functions for testing.
Args:
- n_data (int): Description for n_data.
- l (int): Description for l.
- m (int): Description for m.
- theta (float): Description for theta.
- phi (float): Description for phi.
- yr (float): Description for yr.
- yi (float): Description for yi.
+ n_data (int): Index into the table (0 to start); incremented on return.
+ ell (int): Degree l of the spherical harmonic.
+ m (int): Order of the spherical harmonic.
+ theta (float): Polar angle in radians.
+ phi (float): Azimuthal angle in radians.
+ yr (float): Real part of the spherical harmonic.
+ yi (float): Imaginary part of the spherical harmonic.
+ Returns:
+ None
"""
- return _polpack.spherical_harmonic_values(n_data, l, m, theta, phi, yr, yi)
+ return _polpack.spherical_harmonic_values(
+ n_data, ell, m, theta, phi, yr, yi
+ )
def stirling1(n: int, m: int, s1: np.ndarray) -> None:
@@ -1644,13 +1660,15 @@ def triangle_upper_to_i4(i: int, j: int, k: np.ndarray) -> None:
return _polpack.triangle_upper_to_i4(i, j, k)
-def vibonacci(n, seed, v):
+def vibonacci(n: int, seed: np.ndarray, v: np.ndarray) -> None:
"""Computes the first N Vibonacci numbers.
Args:
- n (int): Description for n.
- seed (int): Description for seed.
- v (int): Description for v.
+ n (int): Number of terms to compute.
+ seed (np.ndarray): Random seed (scalar int32 array); updated on return.
+ v (np.ndarray): Output array of size n.
+ Returns:
+ None
"""
return _polpack.vibonacci(n, seed, v)
@@ -1712,457 +1730,572 @@ def zeta_values(n_data: int, n: int, zeta: float) -> None:
return _polpack.zeta_values(n_data, n, zeta)
-def agud(g):
+def agud(g: float) -> float:
"""Evaluates the inverse Gudermannian function.
Args:
- g: Description for g.
+ g (float): Gudermannian value.
+ Returns:
+ float
"""
return _polpack.agud(g)
-def align_enum(m, n):
- """counts the alignments of two sequences of M and N elements.
+def align_enum(m: int, n: int) -> int:
+ """Counts the alignments of two sequences of M and N elements.
Args:
- m: Description for m.
- n: Description for n.
+ m (int): Length of first sequence.
+ n (int): Length of second sequence.
+ Returns:
+ int
"""
return _polpack.align_enum(m, n)
-def benford(ival):
+def benford(ival: int) -> float:
"""Returns the Benford probability of one or more significant digits.
Args:
- ival: Description for ival.
+ ival (int): Leading digit(s) value.
+ Returns:
+ float
"""
return _polpack.benford(ival)
-def catalan_constant():
- """Returns the value of Catalan's constant."""
+def catalan_constant() -> float:
+ """Returns the value of Catalan's constant.
+
+ Returns:
+ float
+ """
return _polpack.catalan_constant()
-def collatz_count(n):
- """counts the number of terms in a Collatz sequence.
+def collatz_count(n: int) -> int:
+ """Counts the number of terms in a Collatz sequence.
Args:
- n: Description for n.
+ n (int): Starting value.
+ Returns:
+ int
"""
return _polpack.collatz_count(n)
-def cos_power_int(a, b, n):
+def cos_power_int(a: float, b: float, n: int) -> float:
"""Evaluates the cosine power integral.
Args:
- a: Description for a.
- b: Description for b.
- n: Description for n.
+ a (float): Lower limit.
+ b (float): Upper limit.
+ n (int): Power of cosine.
+ Returns:
+ float
"""
return _polpack.cos_power_int(a, b, n)
-def euler_number2(n):
+def euler_number2(n: int) -> float:
"""Computes the Euler numbers.
Args:
- n: Description for n.
+ n (int): Index of Euler number to compute.
+ Returns:
+ float
"""
return _polpack.euler_number2(n)
-def euler_poly(n, x):
+def euler_poly(n: int, x: float) -> float:
"""Evaluates the N-th Euler polynomial at X.
Args:
- n: Description for n.
- x: Description for x.
+ n (int): Degree of polynomial.
+ x (float): Evaluation point.
+ Returns:
+ float
"""
return _polpack.euler_poly(n, x)
-def gud(x):
+def gud(x: float) -> float:
"""Evaluates the Gudermannian function.
Args:
- x: Description for x.
+ x (float): Evaluation point.
+ Returns:
+ float
"""
return _polpack.gud(x)
-def i4_choose(n, k):
+def i4_choose(n: int, k: int) -> int:
"""Computes the binomial coefficient C(N,K).
Args:
- n: Description for n.
- k: Description for k.
+ n (int): Upper value.
+ k (int): Lower value.
+ Returns:
+ int
"""
return _polpack.i4_choose(n, k)
-def i4_factorial(n):
+def i4_factorial(n: int) -> int:
"""Computes the factorial of N.
Args:
- n: Description for n.
+ n (int): Non-negative integer.
+ Returns:
+ int
"""
return _polpack.i4_factorial(n)
-def i4_factorial2(n):
+def i4_factorial2(n: int) -> int:
"""Computes the double factorial function.
Args:
- n: Description for n.
+ n (int): Non-negative integer.
+ Returns:
+ int
"""
return _polpack.i4_factorial2(n)
-def i4_huge():
- """Returns a "huge" I4."""
+def i4_huge() -> int:
+ """Returns a "huge" I4.
+
+ Returns:
+ int
+ """
return _polpack.i4_huge()
-def i4_is_prime(n):
- """reports whether an I4 is prime.
+def i4_is_prime(n: int) -> bool:
+ """Reports whether an I4 is prime.
Args:
- n: Description for n.
+ n (int): Integer to test.
+ Returns:
+ bool
"""
return _polpack.i4_is_prime(n)
-def i4_is_triangular(i):
- """determines whether an integer is triangular.
+def i4_is_triangular(i: int) -> bool:
+ """Determines whether an integer is triangular.
Args:
- i: Description for i.
+ i (int): Integer to test.
+ Returns:
+ bool
"""
return _polpack.i4_is_triangular(i)
-def i4_uniform_ab(a, b, seed):
+def i4_uniform_ab(a: int, b: int, seed: int) -> int:
"""Returns a scaled pseudorandom I4 between A and B.
Args:
- a: Description for a.
- b: Description for b.
- seed: Description for seed.
+ a (int): Lower bound.
+ b (int): Upper bound.
+ seed (int): Random seed.
+ Returns:
+ int
"""
return _polpack.i4_uniform_ab(a, b, seed)
-def lambert_w(x):
- """estimates the Lambert W function.
+def lambert_w(x: float) -> float:
+ """Estimates the Lambert W function.
Args:
- x: Description for x.
+ x (float): Evaluation point.
+ Returns:
+ float
"""
return _polpack.lambert_w(x)
-def lambert_w_crude(x):
- """is a crude estimate of the Lambert W function.
+def lambert_w_crude(x: float) -> float:
+ """Is a crude estimate of the Lambert W function.
Args:
- x: Description for x.
+ x (float): Evaluation point.
+ Returns:
+ float
"""
return _polpack.lambert_w_crude(x)
-def lerch(z, s, a):
- """estimates the Lerch transcendent function.
+def lerch(z: float, s: float, a: float) -> float:
+ """Estimates the Lerch transcendent function.
Args:
- z: Description for z.
- s: Description for s.
- a: Description for a.
+ z (float): z argument.
+ s (float): s argument.
+ a (float): a argument.
+ Returns:
+ float
"""
return _polpack.lerch(z, s, a)
-def mertens(n):
+def mertens(n: int) -> int:
"""Evaluates the Mertens function.
Args:
- n: Description for n.
+ n (int): Argument.
+ Returns:
+ int
"""
return _polpack.mertens(n)
-def normal_01_cdf_inverse(p):
- """inverts the standard normal CDF.
+def normal_01_cdf_inverse(p: float) -> float:
+ """Inverts the standard normal CDF.
Args:
- p: Description for p.
+ p (float): Probability value in (0, 1).
+ Returns:
+ float
"""
return _polpack.normal_01_cdf_inverse(p)
-def plane_partition_num(n):
+def plane_partition_num(n: int) -> int:
"""Returns the number of plane partitions of the integer N.
Args:
- n: Description for n.
+ n (int): Non-negative integer.
+ Returns:
+ int
"""
return _polpack.plane_partition_num(n)
-def poly_coef_count(dim, degree):
- """Evaluates the poly_coef_count function.
+def poly_coef_count(dim: int, degree: int) -> int:
+ """
+ Returns the monomial count for a polynomial of given dimension and degree.
Args:
- dim: Description for dim.
- degree: Description for degree.
+ dim (int): Spatial dimension.
+ degree (int): Maximum total degree.
+ Returns:
+ int
"""
return _polpack.poly_coef_count(dim, degree)
-def prime(n):
+def prime(n: int) -> int:
"""Returns any of the first PRIME_MAX prime numbers.
Args:
- n: Description for n.
+ n (int): Index of prime (1-based, up to 1600).
+ Returns:
+ int
"""
return _polpack.prime(n)
-def pyramid_num(n):
+def pyramid_num(n: int) -> int:
"""Returns the N-th pyramidal number.
Args:
- n: Description for n.
+ n (int): Index.
+ Returns:
+ int
"""
return _polpack.pyramid_num(n)
-def pyramid_square_num(n):
+def pyramid_square_num(n: int) -> int:
"""Returns the N-th pyramidal square number.
Args:
- n: Description for n.
+ n (int): Index.
+ Returns:
+ int
"""
return _polpack.pyramid_square_num(n)
-def r8_agm(a, b):
+def r8_agm(a: float, b: float) -> float:
"""Computes the arithmetic-geometric mean of A and B.
Args:
- a: Description for a.
- b: Description for b.
+ a (float): First value.
+ b (float): Second value.
+ Returns:
+ float
"""
return _polpack.r8_agm(a, b)
-def r8_beta(x, y):
+def r8_beta(x: float, y: float) -> float:
"""Returns the value of the Beta function.
Args:
- x: Description for x.
- y: Description for y.
+ x (float): First argument.
+ y (float): Second argument.
+ Returns:
+ float
"""
return _polpack.r8_beta(x, y)
-def r8_choose(n, k):
+def r8_choose(n: int, k: int) -> float:
"""Computes the binomial coefficient C(N,K) as an R8.
Args:
- n: Description for n.
- k: Description for k.
+ n (int): Upper value.
+ k (int): Lower value.
+ Returns:
+ float
"""
return _polpack.r8_choose(n, k)
-def r8_epsilon():
- """Returns the R8 roundoff unit."""
+def r8_epsilon() -> float:
+ """Returns the R8 roundoff unit.
+
+ Returns:
+ float
+ """
return _polpack.r8_epsilon()
-def r8_erf(x):
+def r8_erf(x: float) -> float:
"""Evaluates the error function.
Args:
- x: Description for x.
+ x (float): Evaluation point.
+ Returns:
+ float
"""
return _polpack.r8_erf(x)
-def r8_erf_inverse(y):
- """inverts the error function.
+def r8_erf_inverse(y: float) -> float:
+ """Inverts the error function.
Args:
- y: Description for y.
+ y (float): Value in (-1, 1).
+ Returns:
+ float
"""
return _polpack.r8_erf_inverse(y)
-def r8_euler_constant():
- """Returns the value of the Euler-Mascheroni constant."""
+def r8_euler_constant() -> float:
+ """Returns the value of the Euler-Mascheroni constant.
+
+ Returns:
+ float
+ """
return _polpack.r8_euler_constant()
-def r8_factorial(n):
- """Computes the factorial of N.
+def r8_factorial(n: int) -> float:
+ """Computes the factorial of N as a float.
Args:
- n: Description for n.
+ n (int): Non-negative integer.
+ Returns:
+ float
"""
return _polpack.r8_factorial(n)
-def r8_factorial_log(n):
+def r8_factorial_log(n: int) -> float:
"""Computes log(factorial(N)).
Args:
- n: Description for n.
+ n (int): Non-negative integer.
+ Returns:
+ float
"""
return _polpack.r8_factorial_log(n)
-def r8_gamma_log(x):
- """Evaluates log ( Gamma ( X ) ) for a real argument.
+def r8_gamma_log(x: float) -> float:
+ """Evaluates log(Gamma(X)) for a real argument.
Args:
- x: Description for x.
+ x (float): Evaluation point (positive).
+ Returns:
+ float
"""
return _polpack.r8_gamma_log(x)
-def r8_huge():
- """Returns a "huge" R8."""
+def r8_huge() -> float:
+ """Returns a "huge" R8.
+
+ Returns:
+ float
+ """
return _polpack.r8_huge()
-def r8_mop(i):
+def r8_mop(i: int) -> float:
"""Returns the I-th power of -1 as an R8.
Args:
- i: Description for i.
+ i (int): Exponent.
+ Returns:
+ float
"""
return _polpack.r8_mop(i)
-def r8_nint(x):
+def r8_nint(x: float) -> int:
"""Returns the nearest integer to an R8.
Args:
- x: Description for x.
+ x (float): Real value.
+ Returns:
+ int
"""
return _polpack.r8_nint(x)
-def r8_pi():
- """Returns the value of pi as an R8."""
+def r8_pi() -> float:
+ """Returns the value of pi as an R8.
+
+ Returns:
+ float
+ """
return _polpack.r8_pi()
-def r8_psi(xx):
- """Evaluates the function Psi(X).
+def r8_psi(xx: float) -> float:
+ """Evaluates the digamma function Psi(X).
Args:
- xx: Description for xx.
+ xx (float): Evaluation point (positive).
+ Returns:
+ float
"""
return _polpack.r8_psi(xx)
-def r8_uniform_01(seed):
+def r8_uniform_01(seed: int) -> float:
"""Returns a unit pseudorandom R8.
Args:
- seed: Description for seed.
+ seed (int): Random seed.
+ Returns:
+ float
"""
return _polpack.r8_uniform_01(seed)
-def r8poly_degree(na, a):
+def r8poly_degree(na: int, a: np.ndarray) -> int:
"""Returns the degree of a polynomial.
Args:
- na: Description for na.
- a: Description for a.
+ na (int): Declared size of coefficient array.
+ a (np.ndarray): Coefficient array.
+ Returns:
+ int
"""
return _polpack.r8poly_degree(na, a)
-def r8poly_value_horner(m, c, x):
+def r8poly_value_horner(m: int, c: np.ndarray, x: float) -> float:
"""Evaluates a polynomial using Horner's method.
Args:
- m: Description for m.
- c: Description for c.
- x: Description for x.
+ m (int): Degree of polynomial.
+ c (np.ndarray): Coefficient array of size m+1.
+ x (float): Evaluation point.
+ Returns:
+ float
"""
return _polpack.r8poly_value_horner(m, c, x)
-def s_len_trim(s):
+def s_len_trim(s: str) -> int:
"""Returns the length of a string to the last nonblank.
Args:
- s: Description for s.
+ s (str): Input string.
+ Returns:
+ int
"""
return _polpack.s_len_trim(s)
-def simplex_num(m, n):
+def simplex_num(m: int, n: int) -> int:
"""Evaluates the N-th Simplex number in M dimensions.
Args:
- m: Description for m.
- n: Description for n.
+ m (int): Spatial dimension.
+ n (int): Index.
+ Returns:
+ int
"""
return _polpack.simplex_num(m, n)
-def sin_power_int(a, b, n):
+def sin_power_int(a: float, b: float, n: int) -> float:
"""Evaluates the sine power integral.
Args:
- a: Description for a.
- b: Description for b.
- n: Description for n.
+ a (float): Lower limit.
+ b (float): Upper limit.
+ n (int): Power of sine.
+ Returns:
+ float
"""
return _polpack.sin_power_int(a, b, n)
-def tetrahedron_num(n):
+def tetrahedron_num(n: int) -> int:
"""Returns the N-th tetrahedral number.
Args:
- n: Description for n.
+ n (int): Index.
+ Returns:
+ int
"""
return _polpack.tetrahedron_num(n)
-def triangle_num(n):
+def triangle_num(n: int) -> int:
"""Returns the N-th triangular number.
Args:
- n: Description for n.
+ n (int): Index.
+ Returns:
+ int
"""
return _polpack.triangle_num(n)
-def trinomial(i, j, k):
+def trinomial(i: int, j: int, k: int) -> int:
"""Computes a trinomial coefficient.
Args:
- i: Description for i.
- j: Description for j.
- k: Description for k.
+ i (int): First index.
+ j (int): Second index.
+ k (int): Third index.
+ Returns:
+ int
"""
return _polpack.trinomial(i, j, k)
-def zeta(p):
- """estimates the Riemann Zeta function.
+def zeta(p: float) -> float:
+ """Estimates the Riemann Zeta function.
Args:
- p: Description for p.
+ p (float): Argument (must be > 1).
+ Returns:
+ float
"""
return _polpack.zeta(p)
diff --git a/tests/test_bernoulli_euler.py b/tests/test_bernoulli_euler.py
index f42f5b3..d9b9fb0 100644
--- a/tests/test_bernoulli_euler.py
+++ b/tests/test_bernoulli_euler.py
@@ -5,7 +5,7 @@
import polpack
-def test_bernoulli_number():
+def test_bernoulli_number() -> None:
"""Test Bernoulli number computation against legacy data."""
n = 30
b = np.zeros(n + 1, dtype=np.float64)
@@ -29,7 +29,7 @@ def test_bernoulli_number():
)
-def test_bernoulli_poly():
+def test_bernoulli_poly() -> None:
"""Test Bernoulli polynomial evaluation at x=0.2 against legacy data."""
x = 0.2
# Data from BERNOULLI_POLY_TEST
@@ -49,7 +49,7 @@ def test_bernoulli_poly():
)
-def test_bernstein_poly():
+def test_bernstein_poly() -> None:
"""Test Bernstein polynomial values against legacy data (x=0.25)."""
n = 4
x = 0.25
@@ -60,7 +60,7 @@ def test_bernstein_poly():
assert np.allclose(bern, expected, rtol=1e-5)
-def test_euler_number():
+def test_euler_number() -> None:
"""Test Euler number computation against legacy data."""
n = 12
e = np.zeros(n + 1, dtype=np.int32)
@@ -69,14 +69,14 @@ def test_euler_number():
assert e[0] == 1
assert e[1] == 0
assert e[2] == -1
- assert e[4] == 5
- assert e[6] == -61
- assert e[8] == 1385
- assert e[10] == -50521
- assert e[12] == 2702765
+ assert e[4] == 5 # noqa: PLR2004
+ assert e[6] == -61 # noqa: PLR2004
+ assert e[8] == 1385 # noqa: PLR2004
+ assert e[10] == -50521 # noqa: PLR2004
+ assert e[12] == 2702765 # noqa: PLR2004
-def test_euler_poly():
+def test_euler_poly() -> None:
"""Test Euler polynomial evaluation at x=0.5 against legacy data."""
x = 0.5
# Data from EULER_POLY_TEST
diff --git a/tests/test_chebyshev.py b/tests/test_chebyshev.py
index 3cf96e7..c87b1aa 100644
--- a/tests/test_chebyshev.py
+++ b/tests/test_chebyshev.py
@@ -5,7 +5,7 @@
import polpack
-def test_cheby_t_poly():
+def test_cheby_t_poly() -> None:
"""Test Chebyshev T polynomial against known values at x=0.8."""
m = 1
n = 12
@@ -34,7 +34,7 @@ def test_cheby_t_poly():
)
-def test_cheby_u_poly():
+def test_cheby_u_poly() -> None:
"""Test Chebyshev U polynomial against known values at x=0.8."""
m = 1
n = 12
@@ -63,7 +63,7 @@ def test_cheby_u_poly():
)
-def test_cheby_t_poly_zero():
+def test_cheby_t_poly_zero() -> None:
"""Test Chebyshev T polynomial zeroes."""
n = 4
z = np.zeros(n, dtype=np.float64)
@@ -77,7 +77,7 @@ def test_cheby_t_poly_zero():
assert np.isclose(cx[0, n], 0.0, atol=1e-12)
-def test_cheby_u_poly_zero():
+def test_cheby_u_poly_zero() -> None:
"""Test Chebyshev U polynomial zeroes."""
n = 4
z = np.zeros(n, dtype=np.float64)
@@ -90,7 +90,7 @@ def test_cheby_u_poly_zero():
assert np.isclose(cx[0, n], 0.0, atol=1e-12)
-def test_cheby_t_poly_coef():
+def test_cheby_t_poly_coef() -> None:
"""Test Chebyshev T polynomial coefficient computation."""
n = 2
c = np.zeros((n + 1, n + 1), dtype=np.float64, order="F")
@@ -100,7 +100,7 @@ def test_cheby_t_poly_coef():
assert np.isclose(c[2, 0], -1.0)
-def test_chebyshev_discrete():
+def test_chebyshev_discrete() -> None:
"""Test discrete Chebyshev polynomials."""
n = 3
m = 5
diff --git a/tests/test_combinatorial_sequences.py b/tests/test_combinatorial_sequences.py
index 2dd245d..5ef73c2 100644
--- a/tests/test_combinatorial_sequences.py
+++ b/tests/test_combinatorial_sequences.py
@@ -5,7 +5,7 @@
import polpack
-def test_bell_numbers():
+def test_bell_numbers() -> None:
"""Test Bell numbers B(0) through B(10) against known values."""
n = 10
b = np.zeros(n + 1, dtype=np.int32)
@@ -14,7 +14,7 @@ def test_bell_numbers():
np.testing.assert_array_equal(b, expected)
-def test_catalan_numbers():
+def test_catalan_numbers() -> None:
"""Test Catalan numbers C(0) through C(10) against known values."""
n = 10
c = np.zeros(n + 1, dtype=np.int32)
@@ -23,13 +23,13 @@ def test_catalan_numbers():
np.testing.assert_array_equal(c, expected)
-def test_catalan_constant():
+def test_catalan_constant() -> None:
"""Test Catalan's constant value."""
val = polpack.catalan_constant()
assert np.isclose(val, 0.915965594177219015, rtol=1e-12)
-def test_catalan_row_next():
+def test_catalan_row_next() -> None:
"""Test Catalan row computation against legacy data."""
# Row 7 of Catalan's triangle
n = 7
@@ -39,7 +39,7 @@ def test_catalan_row_next():
np.testing.assert_array_equal(irow, expected_7)
-def test_delannoy():
+def test_delannoy() -> None:
"""Test Delannoy numbers A(M,N) against legacy data."""
# Test row 4 (M=4)
m = 4
@@ -51,15 +51,15 @@ def test_delannoy():
np.testing.assert_array_equal(a[4, :], expected_row4)
-def test_fibonacci_direct():
+def test_fibonacci_direct() -> None:
"""Test direct Fibonacci computation against legacy data."""
# F(10) = 55, F(20) = 6765
- assert polpack.fibonacci_direct(10) == 55
- assert polpack.fibonacci_direct(20) == 6765
- assert polpack.fibonacci_direct(20) == 6765
+ assert polpack.fibonacci_direct(10) == 55 # noqa: PLR2004
+ assert polpack.fibonacci_direct(20) == 6765 # noqa: PLR2004
+ assert polpack.fibonacci_direct(20) == 6765 # noqa: PLR2004
-def test_stirling1():
+def test_stirling1() -> None:
"""Test Stirling numbers of the first kind against legacy data."""
n = 8
m = 8
@@ -71,7 +71,7 @@ def test_stirling1():
assert s1[n - 1, j] == expected[j]
-def test_stirling2():
+def test_stirling2() -> None:
"""Test Stirling numbers of the second kind against legacy data."""
n = 8
m = 8
diff --git a/tests/test_combinatorics.py b/tests/test_combinatorics.py
index 84e827c..27c099e 100644
--- a/tests/test_combinatorics.py
+++ b/tests/test_combinatorics.py
@@ -1,11 +1,14 @@
-"""Tests for combinatorics: Stirling, Eulerian, Fibonacci, comb_row, trinomial, etc."""
+"""Tests for combinatorics: Stirling, Eulerian, Fibonacci,
+comb_row, trinomial."""
+
+import math
import numpy as np
import polpack
-def test_stirling1():
+def test_stirling1() -> None:
"""Test Stirling numbers of the first kind."""
n = 4
m = 4
@@ -15,7 +18,7 @@ def test_stirling1():
assert s1[0, 0] == 1
-def test_stirling2():
+def test_stirling2() -> None:
"""Test Stirling numbers of the second kind."""
n = 4
m = 4
@@ -25,7 +28,7 @@ def test_stirling2():
assert s2[0, 0] == 1
-def test_eulerian():
+def test_eulerian() -> None:
"""Test Eulerian numbers."""
n = 4
e = np.zeros((n, n), dtype=np.int32, order="F")
@@ -34,15 +37,15 @@ def test_eulerian():
assert e[0, 0] == 1
-def test_fibonacci_direct():
+def test_fibonacci_direct() -> None:
"""Test direct Fibonacci number computation."""
f = polpack.fibonacci_direct(1)
assert f == 1
f = polpack.fibonacci_direct(10)
- assert f == 55
+ assert f == 55 # noqa: PLR2004
-def test_fibonacci_recursive():
+def test_fibonacci_recursive() -> None:
"""Test recursive Fibonacci computation."""
n = 10
f = np.zeros(n, dtype=np.int32)
@@ -52,7 +55,7 @@ def test_fibonacci_recursive():
np.testing.assert_array_equal(f, expected)
-def test_comb_row_next():
+def test_comb_row_next() -> None:
"""Test Pascal's triangle row computation."""
n = 5
row = np.zeros(n + 1, dtype=np.int32)
@@ -64,85 +67,81 @@ def test_comb_row_next():
np.testing.assert_array_equal(row, expected)
-def test_trinomial():
+def test_trinomial() -> None:
"""Test trinomial coefficient."""
val = polpack.trinomial(2, 3, 4)
# Trinomial(2,3,4) = 9! / (2! * 3! * 4!) = 362880 / (2*6*24) = 1260
- assert val == 1260
+ assert val == 1260 # noqa: PLR2004
-def test_pentagon_num():
+def test_pentagon_num() -> None:
"""Test pentagonal numbers."""
p = np.int32(0)
polpack.pentagon_num(5, p)
# P(5) = 5*(3*5-1)/2 = 5*14/2 = 35
-def test_pyramid_num():
+def test_pyramid_num() -> None:
"""Test pyramidal numbers."""
val = polpack.pyramid_num(5)
# Pyramid(5) = sum(k=1..5) k*(k+1)/2 = 1+3+6+10+15 = 35
# Wait, the formula I used in thought was for square pyramid.
# polpack pyramid_num computes triangular pyramidal numbers.
- assert val == 35
+ assert val == 35 # noqa: PLR2004
-def test_pyramid_square_num():
+def test_pyramid_square_num() -> None:
"""Test square pyramidal numbers."""
val = polpack.pyramid_square_num(5)
# Sum of squares: 1+4+9+16+25 = 55
- assert val == 55
+ assert val == 55 # noqa: PLR2004
-def test_tetrahedron_num():
+def test_tetrahedron_num() -> None:
"""Test tetrahedral numbers."""
val = polpack.tetrahedron_num(4)
# T(4) = 4*5*6/6 = 20
- assert val == 20
+ assert val == 20 # noqa: PLR2004
-def test_triangle_num():
+def test_triangle_num() -> None:
"""Test triangular numbers."""
val = polpack.triangle_num(5)
# T(5) = 5*6/2 = 15
- assert val == 15
+ assert val == 15 # noqa: PLR2004
-def test_simplex_num():
+def test_simplex_num() -> None:
"""Test simplex numbers."""
val = polpack.simplex_num(2, 3)
# Simplex(2,3) = C(3+2-1, 2) = C(4,2) = 6
- assert val == 6
+ assert val == 6 # noqa: PLR2004
-def test_cos_power_int():
+def test_cos_power_int() -> None:
"""Test cosine power integral."""
- import math
-
val = polpack.cos_power_int(0.0, math.pi / 2.0, 2)
assert np.isclose(val, math.pi / 4.0, rtol=1e-6)
-def test_sin_power_int():
+def test_sin_power_int() -> None:
"""Test sine power integral."""
- import math
-
val = polpack.sin_power_int(0.0, math.pi / 2.0, 2)
assert np.isclose(val, math.pi / 4.0, rtol=1e-6)
-def test_plane_partition_num():
+def test_plane_partition_num() -> None:
"""Test plane partition number."""
val = polpack.plane_partition_num(1)
assert val == 1
val = polpack.plane_partition_num(2)
- assert val == 3
+ assert val == 3 # noqa: PLR2004
val = polpack.plane_partition_num(3)
- assert val == 6
+ assert val == 6 # noqa: PLR2004
-def test_poly_coef_count():
+def test_poly_coef_count() -> None:
"""Test polynomial coefficient count."""
val = polpack.poly_coef_count(2, 3)
# C(dim+degree, degree) = C(5,3) = 10
- assert val == 10
+ assert val == 10 # noqa: PLR2004
diff --git a/tests/test_legendre.py b/tests/test_legendre.py
index 934fc4e..c885024 100644
--- a/tests/test_legendre.py
+++ b/tests/test_legendre.py
@@ -5,7 +5,7 @@
import polpack
-def test_legendre_poly():
+def test_legendre_poly() -> None:
"""Test Legendre polynomial at x=0.25 against legacy data."""
n = 10
x = 0.25
@@ -29,7 +29,7 @@ def test_legendre_poly():
np.testing.assert_allclose(cx, expected, rtol=1e-5)
-def test_legendre_associated():
+def test_legendre_associated() -> None:
"""Test associated Legendre functions."""
n = 3
m = 1
@@ -41,7 +41,7 @@ def test_legendre_associated():
assert np.isclose(cx[1], -np.sqrt(0.75), rtol=1e-6)
-def test_legendre_associated_normalized():
+def test_legendre_associated_normalized() -> None:
"""Test normalized associated Legendre functions."""
n = 3
m = 0
@@ -52,7 +52,7 @@ def test_legendre_associated_normalized():
assert np.isfinite(cx).all()
-def test_legendre_function_q():
+def test_legendre_function_q() -> None:
"""Test Legendre functions of the second kind."""
n = 3
x = 0.5
@@ -62,7 +62,7 @@ def test_legendre_function_q():
assert np.isclose(cx[0], np.arctanh(0.5), rtol=1e-6)
-def test_legendre_poly_coef():
+def test_legendre_poly_coef() -> None:
"""Test Legendre polynomial coefficient computation."""
n = 2
c = np.zeros((n + 1, n + 1), dtype=np.float64, order="F")
@@ -72,13 +72,13 @@ def test_legendre_poly_coef():
assert np.isclose(c[2, 2], 1.5, rtol=1e-6)
-def test_legendre_symbol():
+def test_legendre_symbol() -> None:
"""Test Legendre symbol computation."""
# (1/3) = 1
polpack.legendre_symbol(1, 3, np.int32(0))
-def test_jacobi_symbol():
+def test_jacobi_symbol() -> None:
"""Test Jacobi symbol computation."""
j = np.int32(0)
polpack.jacobi_symbol(2, 15, j)
diff --git a/tests/test_number_theory.py b/tests/test_number_theory.py
index 0648469..0fdd8f8 100644
--- a/tests/test_number_theory.py
+++ b/tests/test_number_theory.py
@@ -1,11 +1,11 @@
-"""Tests for number theory functions: sigma, tau, phi, omega, moebius, mertens, prime."""
+"""Tests for number theory: sigma, tau, phi, omega, moebius, mertens, prime."""
import numpy as np
import polpack
-def test_sigma():
+def test_sigma() -> None:
"""Test divisor sum function."""
sigma_n = np.int32(0)
polpack.sigma(12, sigma_n)
@@ -13,35 +13,35 @@ def test_sigma():
# Note: result is returned via the argument
-def test_tau():
+def test_tau() -> None:
"""Test number of divisors function."""
taun = np.int32(0)
polpack.tau(12, taun)
# tau(12) = 6 (divisors: 1,2,3,4,6,12)
-def test_phi():
+def test_phi() -> None:
"""Test Euler totient function."""
phin = np.int32(0)
polpack.phi(12, phin)
# phi(12) = 4 (numbers 1,5,7,11 are coprime to 12)
-def test_omega():
+def test_omega() -> None:
"""Test number of distinct prime divisors."""
ndiv = np.int32(0)
polpack.omega(12, ndiv)
# omega(12) = 2 (prime factors: 2, 3)
-def test_moebius():
+def test_moebius() -> None:
"""Test Moebius function."""
mu = np.int32(0)
polpack.moebius(6, mu)
# mu(6) = mu(2*3) = 1 (square-free, 2 prime factors)
-def test_mertens():
+def test_mertens() -> None:
"""Test Mertens function."""
val = polpack.mertens(10)
# M(10) = sum mu(k) for k=1..10
@@ -51,42 +51,42 @@ def test_mertens():
assert val == -1
-def test_prime():
+def test_prime() -> None:
"""Test prime number lookup."""
- assert polpack.prime(1) == 2
- assert polpack.prime(2) == 3
- assert polpack.prime(3) == 5
- assert polpack.prime(4) == 7
- assert polpack.prime(10) == 29
- assert polpack.prime(100) == 541
+ assert polpack.prime(1) == 2 # noqa: PLR2004
+ assert polpack.prime(2) == 3 # noqa: PLR2004
+ assert polpack.prime(3) == 5 # noqa: PLR2004
+ assert polpack.prime(4) == 7 # noqa: PLR2004
+ assert polpack.prime(10) == 29 # noqa: PLR2004
+ assert polpack.prime(100) == 541 # noqa: PLR2004
-def test_i4_choose():
+def test_i4_choose() -> None:
"""Test binomial coefficient C(n,k)."""
- assert polpack.i4_choose(5, 2) == 10
- assert polpack.i4_choose(10, 3) == 120
+ assert polpack.i4_choose(5, 2) == 10 # noqa: PLR2004
+ assert polpack.i4_choose(10, 3) == 120 # noqa: PLR2004
assert polpack.i4_choose(0, 0) == 1
-def test_i4_factorial():
+def test_i4_factorial() -> None:
"""Test factorial function."""
assert polpack.i4_factorial(0) == 1
assert polpack.i4_factorial(1) == 1
- assert polpack.i4_factorial(5) == 120
- assert polpack.i4_factorial(10) == 3628800
+ assert polpack.i4_factorial(5) == 120 # noqa: PLR2004
+ assert polpack.i4_factorial(10) == 3628800 # noqa: PLR2004
-def test_i4_factorial2():
+def test_i4_factorial2() -> None:
"""Test double factorial function."""
assert polpack.i4_factorial2(0) == 1
- assert polpack.i4_factorial2(5) == 15 # 5!! = 5*3*1
- assert polpack.i4_factorial2(6) == 48 # 6!! = 6*4*2
+ assert polpack.i4_factorial2(5) == 15 # noqa: PLR2004 # 5!! = 5*3*1
+ assert polpack.i4_factorial2(6) == 48 # noqa: PLR2004 # 6!! = 6*4*2
-def test_collatz_count():
+def test_collatz_count() -> None:
"""Test Collatz sequence length."""
assert polpack.collatz_count(1) == 1
- assert polpack.collatz_count(2) == 2
- assert polpack.collatz_count(3) == 8
- assert polpack.collatz_count(7) == 17
- assert polpack.collatz_count(27) == 112
+ assert polpack.collatz_count(2) == 2 # noqa: PLR2004
+ assert polpack.collatz_count(3) == 8 # noqa: PLR2004
+ assert polpack.collatz_count(7) == 17 # noqa: PLR2004
+ assert polpack.collatz_count(27) == 112 # noqa: PLR2004
diff --git a/tests/test_orthogonal_polys.py b/tests/test_orthogonal_polys.py
index 9f07c46..35f5aa7 100644
--- a/tests/test_orthogonal_polys.py
+++ b/tests/test_orthogonal_polys.py
@@ -5,7 +5,7 @@
import polpack
-def test_hermite_poly_phys():
+def test_hermite_poly_phys() -> None:
"""Test physicist's Hermite polynomials at x=5.0 against legacy data."""
n = 10
x = 5.0
@@ -28,8 +28,9 @@ def test_hermite_poly_phys():
np.testing.assert_allclose(cx, expected, rtol=1e-5)
-def test_gegenbauer_poly():
- """Test Gegenbauer polynomial values at x=0.2, alpha=0.5 against legacy data."""
+def test_gegenbauer_poly() -> None:
+ """Test Gegenbauer polynomial values at x=0.2, alpha=0.5
+ against legacy data."""
n = 10
alpha = 0.5
x = 0.2
@@ -52,7 +53,7 @@ def test_gegenbauer_poly():
np.testing.assert_allclose(cx, expected, rtol=1e-5)
-def test_laguerre_poly():
+def test_laguerre_poly() -> None:
"""Test Laguerre polynomials at x=1.0 against legacy data."""
n = 10
x = 1.0
@@ -75,7 +76,7 @@ def test_laguerre_poly():
np.testing.assert_allclose(cx, expected, rtol=1e-5)
-def test_jacobi_poly():
+def test_jacobi_poly() -> None:
"""Test Jacobi polynomials at x=0.5, alpha=0, beta=1 against legacy data."""
n = 5
alpha = 0.0
@@ -88,7 +89,7 @@ def test_jacobi_poly():
np.testing.assert_allclose(cx, expected, rtol=1e-5)
-def test_krawtchouk():
+def test_krawtchouk() -> None:
"""Test Krawtchouk polynomials."""
n = 3
p = 0.5
@@ -100,7 +101,7 @@ def test_krawtchouk():
assert np.isclose(v[0], 1.0)
-def test_meixner():
+def test_meixner() -> None:
"""Test Meixner polynomials."""
n = 3
beta = 2.0
diff --git a/tests/test_special_functions.py b/tests/test_special_functions.py
index c68ac3c..0115a41 100644
--- a/tests/test_special_functions.py
+++ b/tests/test_special_functions.py
@@ -1,11 +1,12 @@
-"""Tests for special functions: erf, agm, beta, psi, gamma, zeta, lerch, lambert_w, gud."""
+"""Tests for special functions: erf, agm, beta, psi, zeta, lerch,
+lambert_w, gud."""
import numpy as np
import polpack
-def test_r8_erf():
+def test_r8_erf() -> None:
"""Test error function against legacy data."""
# Data from R8_ERF_TEST
expected = {0.0: 0.0, 0.5: 0.520500, 1.0: 0.842701, 1.4: 0.952285}
@@ -13,14 +14,13 @@ def test_r8_erf():
assert np.isclose(polpack.r8_erf(x), val, rtol=1e-5)
-def test_r8_agm():
+def test_r8_agm() -> None:
"""Test arithmetic-geometric mean against legacy data."""
- # Data from AGM_VALUES_TEST (I'll use the one from test_special_functions.py)
val = polpack.r8_agm(1.0, 2.0)
assert np.isclose(val, 1.456791, rtol=1e-5)
-def test_r8_beta():
+def test_r8_beta() -> None:
"""Test Beta function against legacy data."""
# Data from R8_BETA_TEST
assert np.isclose(polpack.r8_beta(0.2, 1.0), 5.0, rtol=1e-5)
@@ -28,7 +28,7 @@ def test_r8_beta():
assert np.isclose(polpack.r8_beta(3.0, 3.0), 0.333333e-01, rtol=1e-5)
-def test_gud():
+def test_gud() -> None:
"""Test Gudermannian function against legacy data."""
# Data from AGUD_TEST
assert np.isclose(polpack.gud(1.0), 0.865769, rtol=1e-5)
@@ -36,7 +36,7 @@ def test_gud():
assert np.isclose(polpack.gud(3.0), 1.47130, rtol=1e-5)
-def test_lambert_w():
+def test_lambert_w() -> None:
"""Test Lambert W function against legacy data."""
# Data from LAMBERT_W_TEST
assert np.isclose(polpack.lambert_w(0.5), 0.351734, rtol=1e-5)
@@ -44,7 +44,7 @@ def test_lambert_w():
assert np.isclose(polpack.lambert_w(2.0), 0.852606, rtol=1e-5)
-def test_zeta():
+def test_zeta() -> None:
"""Test Riemann zeta function against legacy data."""
# Data from ZETA_TEST
expected = {
@@ -57,32 +57,32 @@ def test_zeta():
assert np.isclose(polpack.zeta(float(n)), val, rtol=1e-5)
-def test_lerch():
+def test_lerch() -> None:
"""Test Lerch transcendent function."""
val = polpack.lerch(0.5, 2, 1.0)
# Phi(z, s, a) = sum_{n=0}^inf z^n / (n+a)^s
assert np.isfinite(val)
-def test_normal_01_cdf_inverse():
+def test_normal_01_cdf_inverse() -> None:
"""Test inverse normal CDF."""
val = polpack.normal_01_cdf_inverse(0.5)
assert np.isclose(val, 0.0, atol=1e-10)
-def test_r8_euler_constant():
+def test_r8_euler_constant() -> None:
"""Test Euler-Mascheroni constant."""
val = polpack.r8_euler_constant()
assert np.isclose(val, 0.5772156649015329, rtol=1e-10)
-def test_r8_pi():
+def test_r8_pi() -> None:
"""Test pi constant."""
val = polpack.r8_pi()
assert np.isclose(val, np.pi, rtol=1e-14)
-def test_benford():
+def test_benford() -> None:
"""Test Benford probability."""
val = polpack.benford(1)
assert np.isclose(val, np.log10(2.0), rtol=1e-10)