Add dedicated epsilon page and reorganize internal utilities

- Created comprehensive "Step Size Selection" page covering:
  - Mathematical theory of optimal step sizes
  - Error analysis (truncation vs round-off)
  - Adaptive step sizing algorithms
  - Special cases for complex step and sparse Jacobians

- Renamed "Utilities" to "Internal Utilities" with complete coverage of:
  - Array utilities (_vec, _mat, setindex overloads)
  - Hessian utilities (_hessian_inplace, mutable_zeromatrix)
  - Sparse Jacobian internals (_make_Ji, _colorediteration!)
  - Performance optimizations and sparsity detection
  - Cache resize functions and error handling

- Updated pages.jl navigation structure
- Documentation builds cleanly with proper organization

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <[email protected]>

Author: ChrisRackauckas

docs/pages.jl

@@ -8,5 +8,6 @@ pages = [
     "Jacobians" => "jacobians.md",
     "Hessians" => "hessians.md",
     "Jacobian-Vector Products" => "jvp.md",
-    "Utilities" => "utilities.md"
+    "Step Size Selection" => "epsilons.md",
+    "Internal Utilities" => "utilities.md"
 ]

docs/src/epsilons.md

@@ -0,0 +1,71 @@
+# Step Size Selection (Epsilons)
+
+Functions and theory for computing optimal step sizes in finite difference approximations.
+
+## Theory
+
+The choice of step size (epsilon) in finite difference methods is critical for accuracy. Too large a step leads to truncation error, while too small a step leads to round-off error. The optimal step size balances these two sources of error.
+
+### Error Analysis
+
+For a function `f` with bounded derivatives, the total error in finite difference approximations consists of:
+
+1. **Truncation Error**: Comes from the finite difference approximation itself
+   - Forward differences: `O(h)` where `h` is the step size
+   - Central differences: `O(h²)`
+   - Hessian central differences: `O(h²)` for second derivatives
+
+2. **Round-off Error**: Comes from floating-point arithmetic
+   - Forward differences: `O(eps/h)` where `eps` is machine epsilon
+   - Central differences: `O(eps/h)`
+
+### Optimal Step Sizes
+
+Minimizing the total error `truncation + round-off` gives optimal step sizes:
+
+- **Forward differences**: `h* = sqrt(eps)` - balances `O(h)` truncation with `O(eps/h)` round-off
+- **Central differences**: `h* = eps^(1/3)` - balances `O(h²)` truncation with `O(eps/h)` round-off
+- **Hessian central**: `h* = eps^(1/4)` - balances `O(h²)` truncation for mixed derivatives
+- **Complex step**: `h* = eps` - no subtractive cancellation, only limited by machine precision
+
+## Adaptive Step Sizing
+
+The step size computation uses both relative and absolute components:
+
+```julia
+epsilon = max(relstep * abs(x), absstep) * dir
+```
+
+This ensures:
+- **Large values**: Use relative step `relstep * |x|` for scale-invariant accuracy
+- **Small values**: Use absolute step `absstep` to avoid underflow
+- **Direction**: Multiply by `dir` (±1) for forward differences
+
+## Implementation
+
+The step size computation is handled by internal functions:
+
+- **`compute_epsilon(fdtype, x, relstep, absstep, dir)`**: Computes the actual step size for a given finite difference method and input value
+- **`default_relstep(fdtype, T)`**: Returns the optimal relative step size for a given method and numeric type
+
+These functions are called automatically by all finite difference routines, but understanding their behavior can help with custom implementations or debugging numerical issues.
+
+## Special Cases
+
+### Complex Step Differentiation
+
+For complex step differentiation, the step size is simply machine epsilon since this method avoids subtractive cancellation entirely:
+
+⚠️ **Important**: The function `f` must be complex analytic when the input is complex!
+
+### Sparse Jacobians
+
+When computing sparse Jacobians with graph coloring, the step size is computed based on the norm of the perturbation vector to ensure balanced accuracy across all columns in the same color group.
+
+## Practical Considerations
+
+- **Default step sizes** are optimal for most smooth functions
+- **Custom step sizes** may be needed for functions with unusual scaling or near-discontinuities
+- **Relative steps** should scale with the magnitude of the input
+- **Absolute steps** provide a fallback for inputs near zero
+- **Direction parameter** allows for one-sided differences when needed (e.g., at domain boundaries)

docs/src/utilities.md

@@ -1,31 +1,73 @@
-# Utilities
+# Internal Utilities
 
-Internal utility functions and types that may be useful for advanced users.
+Internal utility functions and algorithms used throughout FiniteDiff.jl. These functions are primarily for advanced users who need to understand or extend the library's functionality.
 
-## Step Size Computation
+## Array Utilities
 
-While these functions are primarily internal, they can be useful for understanding and customizing finite difference computations:
+Internal utilities for handling different array types and ensuring compatibility across the Julia ecosystem:
 
-### Default Step Sizes
+### Type Conversion Functions
 
-Different finite difference methods use different optimal step sizes to balance truncation error and round-off error:
+These functions help FiniteDiff.jl work with various array types including `StaticArrays`, structured matrices, and GPU arrays:
 
-- **Forward differences**: `sqrt(eps(T))` - balances O(h) truncation error with round-off
-- **Central differences**: `cbrt(eps(T))` - balances O(h²) truncation error with round-off  
-- **Hessian central**: `eps(T)^(1/4)` - balances O(h²) truncation error for second derivatives
+- **`_vec(x)`**: Vectorizes arrays while preserving scalars unchanged
+- **`_mat(x)`**: Ensures matrix format, converting vectors to column matrices  
+- **`setindex(x, v, i...)`**: Non-mutating setindex operations for immutable arrays
 
-### Complex Step Differentiation
+### Hessian Utilities
 
-Complex step differentiation uses machine epsilon since it avoids subtractive cancellation:
+Helper functions specifically for Hessian computation:
 
-⚠️ **Important**: `f` must be a function of a real variable that is also complex analytic when the input is complex!
+- **`_hessian_inplace(x)`**: Determines whether to use in-place operations based on array mutability
+- **`__Symmetric(x)`**: Wraps matrices in `Symmetric` views for mathematical correctness
+- **`mutable_zeromatrix(x)`**: Creates mutable zero matrices compatible with input arrays
 
-## Array Utilities
+## Sparse Jacobian Internals
+
+Graph coloring and sparse matrix algorithms for efficient Jacobian computation:
+
+### Matrix Construction
+
+- **`_make_Ji(...)`**: Constructs Jacobian contribution matrices for both sparse and dense cases
+- **`_colorediteration!(...)`**: Core loop for sparse Jacobian assembly using graph coloring
+- **`_findstructralnz(A)`**: Finds structural non-zero patterns in dense matrices
+
+### Sparsity Detection
+
+- **`_use_findstructralnz(sparsity)`**: Determines when to use structural sparsity information
+- **`_use_sparseCSC_common_sparsity(J, sparsity)`**: Tests for common sparsity patterns between matrices
+
+### Performance Optimizations
+
+- **`fast_jacobian_setindex!(...)`**: Optimized index operations for sparse Jacobian assembly
+- **`void_setindex!(...)`**: Wrapper for setindex operations that discards return values
+
+## JVP Utilities
+
+- **`resize!(cache::JVPCache, i)`**: Resizes JVP cache arrays for dynamic problems
+- **`resize!(cache::JacobianCache, i)`**: Resizes Jacobian cache arrays for dynamic problems
+
+## Error Handling
+
+- **`fdtype_error(::Type{T})`**: Provides informative error messages for unsupported finite difference type combinations
+
+## Design Philosophy
+
+These internal utilities follow several design principles:
+
+1. **Type Stability**: All functions are designed for type-stable operations
+2. **Zero Allocation**: Internal utilities avoid allocations when possible
+3. **Genericity**: Support for multiple array types (dense, sparse, static, GPU)
+4. **Performance**: Optimized for the specific needs of finite difference computation
+5. **Safety**: Proper error handling and bounds checking where needed
+
+## Advanced Usage
 
-Internal utilities for handling different array types and ensuring compatibility:
+These functions are primarily internal, but advanced users may find them useful for:
 
-- `_vec(x)`: Vectorizes arrays while preserving scalars
-- `_mat(x)`: Ensures matrix format, converting vectors to column matrices
-- `setindex(x, v, i...)`: Non-mutating setindex for immutable arrays
+- **Custom finite difference implementations**
+- **Integration with other differentiation libraries**  
+- **Performance optimization in specialized applications**
+- **Understanding the implementation details of FiniteDiff.jl**
 
-These functions help FiniteDiff.jl work with various array types including `StaticArrays`, structured matrices, and GPU arrays.
+Most users should rely on the main API functions rather than calling these utilities directly.