factor with method selection

Hello,

Currently, `factor()` uses a hardcoded algorithm chain: trial division up to 2^32, then `pollardfactor()`. With ECM and eventually MPQS being added, users would benefit from:

1. Explicitly choosing a factorization algorithm
2. Tuning algorithm parameters (e.g., ECM smoothness bound)
3. A clean extension point for future algorithms

Proposed API

Following SciML conventions (cf. `solve(prob, alg`) in DifferentialEquations.jl ...):
```julia
abstract type FactorizationMethod end
struct TrialDivision <: FactorizationMethod end
struct PollardRho <: FactorizationMethod end
struct ECM <: FactorizationMethod
B1::Int
num_curves::Int
end
struct Auto <: FactorizationMethod end

# Keyword API (backwards compatible)
factor(n; method=Auto())

# Positional dispatch (Julian style)
factor(n, ::TrialDivision)
factor(n, ::PollardRho)
factor(n, ::ECM)
factor(n, ::Auto) # current polyalgorithm heuristic
```

`Auto()` would encode the current behavior (trial division → Pollard rho) and grow to include ECM/MPQS size-based heuristics.

Backward compatibility

- `factor(n)` continues to work unchanged (defaults to Auto())
- `factor(ContainerType, n)` unaffected
- `eachfactor(n)` would also accept a method keyword

Considerations

- Should strategy types be exported? Probably yes for discoverability.
- Should eachfactor support methods too, or only factor?
- Parameterless strategies (TrialDivision(), PollardRho()) could carry default tuning parameters later without breaking API.

Related:
- https://github.com/JuliaMath/Primes.jl/issues/159#issuecomment-2538065172