diff --git a/lmo/_poly.py b/lmo/_poly.py
index 8d7e6a59..28af47cd 100644
--- a/lmo/_poly.py
+++ b/lmo/_poly.py
@@ -18,7 +18,7 @@
 
 
 if TYPE_CHECKING:
-    from .typing import np as lnpt
+    import lmo.typing.np as lnpt
 
 
 __all__ = (
diff --git a/lmo/contrib/scipy_stats.py b/lmo/contrib/scipy_stats.py
index bcf7911b..afbca912 100644
--- a/lmo/contrib/scipy_stats.py
+++ b/lmo/contrib/scipy_stats.py
@@ -28,10 +28,8 @@
 import lmo.typing as lmt
 import lmo.typing.np as lnpt
 import lmo.typing.scipy as lspt
-from lmo import (
-    inference,
-    l_moment as l_moment_est,
-)
+from lmo import inference
+from lmo._lm import l_moment as l_moment_est
 from lmo._utils import (
     clean_order,
     clean_orders,
@@ -1089,33 +1087,29 @@ def l_fit(
         self,
         data: lnpt.AnyVectorInt | lnpt.AnyVectorFloat,
         *args: float,
-        n_extra: int = 0,
-        trim: lmt.AnyTrim = 0,
+        n_extra: int = ...,
+        trim: lmt.AnyTrimInt = ...,
         full_output: Literal[True],
-        fit_kwargs: Mapping[str, Any] | None = None,
+        fit_kwargs: Mapping[str, Any] | None = ...,
         **kwds: Any,
-    ) -> tuple[float, ...]:
-        ...
-
+    ) -> tuple[float, ...]: ...
     @overload
     def l_fit(
         self,
         data: lnpt.AnyVectorInt | lnpt.AnyVectorFloat,
         *args: float,
-        n_extra: int = 0,
-        trim: lmt.AnyTrim = 0,
+        n_extra: int = ...,
+        trim: lmt.AnyTrimInt = ...,
         full_output: bool = ...,
-        fit_kwargs: Mapping[str, Any] | None = None,
+        fit_kwargs: Mapping[str, Any] | None = ...,
         **kwds: Any,
-    ) -> tuple[float, ...]:
-        ...
-
+    ) -> tuple[float, ...]: ...
     def l_fit(
         self,
         data: lnpt.AnyVectorInt | lnpt.AnyVectorFloat,
         *args: float,
         n_extra: int = 0,
-        trim: lmt.AnyTrim = 0,
+        trim: lmt.AnyTrimInt = 0,
         full_output: bool = False,
         fit_kwargs: Mapping[str, Any] | None = None,
         random_state: int | np.random.Generator | None = None,
diff --git a/lmo/ostats.py b/lmo/ostats.py
index fcd2804b..340d248d 100644
--- a/lmo/ostats.py
+++ b/lmo/ostats.py
@@ -22,7 +22,7 @@
 if TYPE_CHECKING:
     import optype.numpy as onpt
 
-    from .typing import np as lnpt
+    import lmo.typing.np as lnpt
 
 
 __all__ = 'weights', 'from_cdf'
diff --git a/lmo/pwm_beta.py b/lmo/pwm_beta.py
index 7df69d13..b1191727 100644
--- a/lmo/pwm_beta.py
+++ b/lmo/pwm_beta.py
@@ -24,13 +24,10 @@
 if TYPE_CHECKING:
     import optype.numpy as onpt
 
-    from .typing import np as lnpt
+    import lmo.typing.np as lnpt
 
 
-__all__ = (
-    'weights',
-    'cov',
-)
+__all__ = 'weights', 'cov'
 
 
 _F = TypeVar('_F', bound=np.floating[Any], default=np.float64)
diff --git a/lmo/special.py b/lmo/special.py
index 33b068cf..cff436df 100644
--- a/lmo/special.py
+++ b/lmo/special.py
@@ -8,21 +8,12 @@
 import optype.numpy as onpt
 import scipy.special as sc
 
+import lmo.typing as lmt
+import lmo.typing.np as lnpt
 from ._utils import clean_orders
-from .typing import (
-    AnyOrder,
-    AnyOrderND,
-    np as lnpt,
-)
 
 
-__all__ = (
-    'fpow',
-    'gamma2',
-    'harmonic',
-    'norm_sh_jacobi',
-    'fourier_jacobi',
-)
+__all__ = 'fpow', 'gamma2', 'harmonic', 'norm_sh_jacobi', 'fourier_jacobi'
 
 
 _DTYPE_CHARS: Final[str] = '?bBhHiIlLqQpP'
@@ -66,8 +57,6 @@ def fpow(
     /,
     out: onpt.Array[_T_shape, _T_float],
 ) -> onpt.Array[_T_shape, _T_float]: ...
-
-
 def fpow(
     x: lnpt.AnyScalarFloat | lnpt.AnyArrayFloat,
     n: lnpt.AnyScalarFloat | lnpt.AnyArrayFloat,
@@ -126,8 +115,6 @@ def gamma2(
     /,
     out: onpt.Array[_T_shape, _T_float],
 ) -> onpt.Array[_T_shape, _T_float]: ...
-
-
 def gamma2(
     a: lnpt.AnyScalarFloat,
     x: lnpt.AnyScalarFloat | lnpt.AnyArrayFloat,
@@ -224,17 +211,19 @@ def harmonic(
 
 
 @overload
-def norm_sh_jacobi(n: AnyOrder, alpha: float, beta: float) -> np.float64: ...
+def norm_sh_jacobi(
+    n: lmt.AnyOrder,
+    alpha: float,
+    beta: float,
+) -> np.float64: ...
 @overload
 def norm_sh_jacobi(
-    n: AnyOrderND,
+    n: lmt.AnyOrderND,
     alpha: float,
     beta: float,
 ) -> onpt.Array[onpt.AtLeast1D, np.float64]: ...
-
-
 def norm_sh_jacobi(
-    n: AnyOrder | AnyOrderND,
+    n: lmt.AnyOrder | lmt.AnyOrderND,
     alpha: float,
     beta: float,
 ) -> np.float64 | onpt.Array[onpt.AtLeast1D, np.float64]: