diff --git a/src/sparse_ir/_util.py b/src/sparse_ir/_util.py
index 16ba422..def84dd 100644
--- a/src/sparse_ir/_util.py
+++ b/src/sparse_ir/_util.py
@@ -71,6 +71,72 @@ def check_range(x, xmin, xmax):
     return x
 
 
+def normalize_tau(statistics, tau, beta):
+    """Normalize τ to [0, β] with statistics-dependent periodicity.
+    
+    Handles boundary conditions based on statistics:
+    - Fermions ('F'): Anti-periodic G(τ + β) = -G(τ)
+    - Bosons ('B'): Periodic G(τ + β) = G(τ)
+    
+    This function maps τ values from the range [-β, β] to [0, β] with
+    appropriate sign factors, following the periodicity rules.
+    
+    Arguments:
+        statistics (str):
+            'F' for Fermionic or 'B' for Bosonic statistics.
+        tau (array_like):
+            Imaginary time value(s) in range [-β, β].
+        beta (float):
+            Inverse temperature.
+            
+    Returns:
+        tuple[ndarray, ndarray]:
+            (tau_normalized, sign) where:
+            - tau_normalized: τ values mapped to [0, β]
+            - sign: Sign factor (±1) for periodicity
+            
+    Raises:
+        ValueError: If tau is outside [-β, β] or statistics is invalid.
+        
+    Special cases:
+        - Negative zero (τ = -0.0) is treated as τ = β with appropriate sign
+        - For τ in [0, β]: returns (τ, +1)
+        - For τ in [-β, 0): returns (τ + β, sign) where sign depends on statistics
+        
+    .. versionadded:: 1.2
+    """
+    tau = np.asarray(tau, dtype=np.float64)
+    beta = float(beta)
+    
+    if statistics not in ('F', 'B'):
+        raise ValueError("statistics must be 'F' (Fermionic) or 'B' (Bosonic)")
+    
+    if np.any(tau < -beta) or np.any(tau > beta):
+        raise ValueError(f"τ must be in [-β, β] = [{-beta}, {beta}]")
+    
+    # Handle negative zero: τ = -0.0 → τ = β
+    is_neg_zero = (tau == 0.0) & np.signbit(tau)
+    
+    tau_normalized = np.where(is_neg_zero, beta, tau)
+    sign = np.ones_like(tau, dtype=np.float64)
+    
+    if statistics == 'F':
+        # Fermionic: anti-periodic
+        sign = np.where(is_neg_zero, -1.0, sign)
+    else:  # statistics == 'B'
+        # Bosonic: periodic
+        sign = np.where(is_neg_zero, 1.0, sign)
+    
+    # Normalize negative tau to [0, β]
+    mask_neg = tau_normalized < 0
+    tau_normalized = np.where(mask_neg, tau_normalized + beta, tau_normalized)
+    
+    if statistics == 'F':
+        sign = np.where(mask_neg, -sign, sign)
+    
+    return tau_normalized, sign
+
+
 def check_svd_result(svd_result, matrix_shape=None):
     """Checks that argument is a valid SVD triple (u, s, vH)"""
     u, s, vH = map(np.asarray, svd_result)
diff --git a/src/sparse_ir/augment.py b/src/sparse_ir/augment.py
index b26c50f..a8950af 100644
--- a/src/sparse_ir/augment.py
+++ b/src/sparse_ir/augment.py
@@ -111,14 +111,32 @@ def beta(self):
     def wmax(self):
         return self._basis.wmax
 
-    def default_tau_sampling_points(self, *, npoints=None):
+    def default_tau_sampling_points(self, *, npoints=None, use_positive_taus=True):
+        """Get default tau sampling points for augmented basis.
+        
+        Arguments:
+            npoints (int):
+                Minimum number of sampling points to return. If None, uses self.size.
+            use_positive_taus (bool):
+                If True, fold points to [0, β] range and sort them (default: True).
+                If False, points are in symmetric range.
+                
+                .. versionadded:: 1.2
+        """
         if npoints is None:
             npoints = self.size
 
         # Return the sampling points of the underlying basis, but since we
         # use the size of self, we add two further points.  One then has to
         # hope that these give good sampling points.
-        return basis_get_default_tau_sampling_points_ext(self._basis._ptr, npoints)
+
+        points = basis_get_default_tau_sampling_points_ext(self._basis._ptr, npoints)
+        
+        if use_positive_taus:
+            points = np.mod(points, self.beta)
+            points = np.sort(points)
+        
+        return points
 
     def default_matsubara_sampling_points(self, *, positive_only=False):
         """Get default Matsubara sampling points for augmented basis.
@@ -247,20 +265,38 @@ def hat(self, n):
 
 
 class TauConst(AbstractAugmentation):
-    """Constant in imaginary time/discrete delta in frequency"""
+    """Constant in imaginary time with statistics-dependent periodicity.
+    
+    Evaluates to a constant value in imaginary time with proper handling of
+    periodicity based on statistics:
+    - Fermions: Anti-periodic G(τ + β) = -G(τ)
+    - Bosons: Periodic G(τ + β) = G(τ)
+    
+    .. versionchanged:: 1.2
+        Added statistics parameter and support for [-β, β] range.
+    """
     @classmethod
     def create(cls, basis):
-        _check_bosonic_statistics(basis.statistics)
-        return cls(basis.beta)
+        return cls(basis.beta, basis.statistics)
 
-    def __init__(self, beta):
+    def __init__(self, beta, statistics='B'):
+        """
+        Arguments:
+            beta (float):
+                Inverse temperature.
+            statistics (str):
+                'F' for Fermionic or 'B' for Bosonic (default: 'B' for backward compatibility).
+        """
         if beta <= 0:
             raise ValueError("temperature must be positive")
+        if statistics not in ('F', 'B'):
+            raise ValueError("statistics must be 'F' or 'B'")
         self._beta = beta
+        self._statistics = statistics
 
     def __call__(self, tau):
-        tau = _util.check_range(tau, -self._beta, self._beta)
-        return np.broadcast_to(1 / np.sqrt(self._beta), tau.shape)
+        tau_normalized, sign = _util.normalize_tau(self._statistics, tau, self._beta)
+        return sign / np.sqrt(self._beta)
 
     def deriv(self, n=1):
         if n == 0:
@@ -269,27 +305,46 @@ def deriv(self, n=1):
             return lambda tau: np.zeros_like(tau)
 
     def hat(self, n):
-        n = _util.check_reduced_matsubara(n, zeta=0)
+        zeta = 1 if self._statistics == 'F' else 0
+        n = _util.check_reduced_matsubara(n, zeta=zeta)
         return np.sqrt(self._beta) * (n == 0).astype(complex)
 
 
 class TauLinear(AbstractAugmentation):
-    """Linear function in imaginary time, antisymmetric around beta/2"""
+    """Linear function in imaginary time with statistics-dependent periodicity.
+    
+    Evaluates to a linear function antisymmetric around β/2 with proper handling
+    of periodicity based on statistics:
+    - Fermions: Anti-periodic G(τ + β) = -G(τ)
+    - Bosons: Periodic G(τ + β) = G(τ)
+    
+    .. versionchanged:: 1.2
+        Added statistics parameter and support for [-β, β] range.
+    """
     @classmethod
     def create(cls, basis):
-        _check_bosonic_statistics(basis.statistics)
-        return cls(basis.beta)
+        return cls(basis.beta, basis.statistics)
 
-    def __init__(self, beta):
+    def __init__(self, beta, statistics='B'):
+        """
+        Arguments:
+            beta (float):
+                Inverse temperature.
+            statistics (str):
+                'F' for Fermionic or 'B' for Bosonic (default: 'B' for backward compatibility).
+        """
         if beta <= 0:
             raise ValueError("temperature must be positive")
+        if statistics not in ('F', 'B'):
+            raise ValueError("statistics must be 'F' or 'B'")
         self._beta = beta
+        self._statistics = statistics
         self._norm = np.sqrt(3/beta)
 
     def __call__(self, tau):
-        tau = _util.check_range(tau, -self._beta, self._beta)
-        x = 2/self._beta * tau - 1
-        return self._norm * x
+        tau_normalized, sign = _util.normalize_tau(self._statistics, tau, self._beta)
+        x = 2/self._beta * tau_normalized - 1
+        return sign * self._norm * x
 
     def deriv(self, n=1):
         if n == 0:
@@ -301,22 +356,43 @@ def deriv(self, n=1):
             return lambda tau: np.zeros_like(tau)
 
     def hat(self, n):
-        n = _util.check_reduced_matsubara(n, zeta=0)
+        zeta = 1 if self._statistics == 'F' else 0
+        n = _util.check_reduced_matsubara(n, zeta=zeta)
         inv_w = np.pi/self._beta * n
         inv_w = np.reciprocal(inv_w, out=inv_w, where=n.astype(bool))
         return self._norm * 2/1j * inv_w
 
 
 class MatsubaraConst(AbstractAugmentation):
-    """Constant in Matsubara, undefined in imaginary time"""
+    """Constant in Matsubara, undefined in imaginary time.
+    
+    This augmentation is constant in Matsubara frequency space and returns NaN
+    in imaginary time. The statistics parameter is accepted for type consistency
+    but does not affect the behavior.
+    
+    .. versionchanged:: 1.2
+        Accepts tau in [-β, β] range (previously [0, β]).
+        Added statistics parameter for consistency.
+    """
     @classmethod
     def create(cls, basis):
-        return cls(basis.beta)
+        return cls(basis.beta, basis.statistics)
 
-    def __init__(self, beta):
+    def __init__(self, beta, statistics=None):
+        """
+        Arguments:
+            beta (float):
+                Inverse temperature.
+            statistics (str, optional):
+                'F' for Fermionic or 'B' for Bosonic. Accepted for type consistency
+                but behavior is identical for both.
+        """
         if beta <= 0:
             raise ValueError("temperature must be positive")
+        if statistics is not None and statistics not in ('F', 'B'):
+            raise ValueError("statistics must be 'F' or 'B'")
         self._beta = beta
+        self._statistics = statistics
 
     def __call__(self, tau):
         tau = _util.check_range(tau, -self._beta, self._beta)
diff --git a/src/sparse_ir/basis.py b/src/sparse_ir/basis.py
index 4a7c006..4b6ce43 100644
--- a/src/sparse_ir/basis.py
+++ b/src/sparse_ir/basis.py
@@ -176,9 +176,23 @@ def shape(self):
         """Shape of the basis function set"""
         return self.s.shape
 
-    def default_tau_sampling_points(self, npoints=None):
-        """Get default tau sampling points."""
-        return basis_get_default_tau_sampling_points(self._ptr)
+    def default_tau_sampling_points(self, npoints=None, use_positive_taus=True):
+        """Get default tau sampling points.
+        
+        Arguments:
+            npoints (int):
+                Minimum number of sampling points to return (currently unused).
+            use_positive_taus (bool):
+                If True, fold points to [0, β] range and sort them (default: True).
+                If False, points are in symmetric range around β/2.
+                
+                .. versionadded:: 1.2
+        """
+        points = basis_get_default_tau_sampling_points(self._ptr)
+        if use_positive_taus:
+            points = np.mod(points, self.beta)
+            points = np.sort(points)
+        return points
 
     def default_omega_sampling_points(self, npoints=None):
         """Return default sampling points in imaginary time.
diff --git a/src/sparse_ir/dlr.py b/src/sparse_ir/dlr.py
index 65c6694..56c5ac2 100644
--- a/src/sparse_ir/dlr.py
+++ b/src/sparse_ir/dlr.py
@@ -209,8 +209,8 @@ def to_IR(self, g_dlr: np.ndarray, axis=0) -> np.ndarray:
             raise RuntimeError(f"Failed to convert DLR to IR: {ret}")
         return output
 
-    def default_tau_sampling_points(self):
-        return self._basis.default_tau_sampling_points()
+    def default_tau_sampling_points(self, **kwargs):
+        return self._basis.default_tau_sampling_points(**kwargs)
 
     def default_matsubara_sampling_points(self, **kwargs):
         return self._basis.default_matsubara_sampling_points(**kwargs)
diff --git a/src/sparse_ir/sampling.py b/src/sparse_ir/sampling.py
index f6ef891..b4e0a3f 100644
--- a/src/sparse_ir/sampling.py
+++ b/src/sparse_ir/sampling.py
@@ -43,9 +43,9 @@ def __init__(self, basis, sampling_points=None, use_positive_taus=True):
         self.basis = basis
 
         if sampling_points is None:
-            self.sampling_points = basis.default_tau_sampling_points()
-            if use_positive_taus:
-                self.sampling_points = np.mod(self.sampling_points, basis.beta)
+            self.sampling_points = basis.default_tau_sampling_points(
+                use_positive_taus=use_positive_taus
+            )
         else:
             self.sampling_points = np.asarray(sampling_points, dtype=np.float64)
 
diff --git a/tests/test_augment.py b/tests/test_augment.py
index 2019e21..6ae1802 100644
--- a/tests/test_augment.py
+++ b/tests/test_augment.py
@@ -4,6 +4,7 @@
 import sparse_ir
 from sparse_ir import poly
 from sparse_ir import augment
+from sparse_ir import _util
 
 import pytest
 
@@ -57,3 +58,196 @@ def test_vertex_basis(stat):
 
     np.testing.assert_allclose(giv, giv_reconst,
                                atol=np.abs(giv).max() * 1e-7, rtol=0)
+
+def test_normalize_tau_bosonic():
+    """Test normalize_tau for bosonic statistics"""
+    beta = 10.0
+    
+    # Test positive tau
+    tau_norm, sign = _util.normalize_tau('B', 5.0, beta)
+    assert tau_norm == 5.0
+    assert sign == 1.0
+    
+    # Test negative tau (periodic)
+    tau_norm, sign = _util.normalize_tau('B', -3.0, beta)
+    assert np.isclose(tau_norm, 7.0)
+    assert sign == 1.0
+    
+    # Test negative zero
+    tau_norm, sign = _util.normalize_tau('B', -0.0, beta)
+    assert tau_norm == beta
+    assert sign == 1.0
+    
+    # Test array input
+    taus = np.array([0.0, 5.0, -3.0, beta])
+    tau_norms, signs = _util.normalize_tau('B', taus, beta)
+    assert np.allclose(tau_norms, [0.0, 5.0, 7.0, beta])
+    assert np.allclose(signs, [1.0, 1.0, 1.0, 1.0])
+
+
+def test_normalize_tau_fermionic():
+    """Test normalize_tau for fermionic statistics"""
+    beta = 10.0
+    
+    # Test positive tau
+    tau_norm, sign = _util.normalize_tau('F', 5.0, beta)
+    assert tau_norm == 5.0
+    assert sign == 1.0
+    
+    # Test negative tau (anti-periodic)
+    tau_norm, sign = _util.normalize_tau('F', -3.0, beta)
+    assert np.isclose(tau_norm, 7.0)
+    assert sign == -1.0
+    
+    # Test negative zero
+    tau_norm, sign = _util.normalize_tau('F', -0.0, beta)
+    assert tau_norm == beta
+    assert sign == -1.0
+    
+    # Test array input
+    taus = np.array([0.0, 5.0, -3.0])
+    tau_norms, signs = _util.normalize_tau('F', taus, beta)
+    assert np.allclose(tau_norms, [0.0, 5.0, 7.0])
+    assert np.allclose(signs, [1.0, 1.0, -1.0])
+
+
+def test_normalize_tau_errors():
+    """Test normalize_tau error handling"""
+    beta = 10.0
+    
+    # Out of range
+    with pytest.raises(ValueError, match="τ must be in"):
+        _util.normalize_tau('B', -beta - 1, beta)
+    
+    with pytest.raises(ValueError, match="τ must be in"):
+        _util.normalize_tau('B', beta + 1, beta)
+    
+    # Invalid statistics
+    with pytest.raises(ValueError, match="statistics must be"):
+        _util.normalize_tau('X', 0.0, beta)
+
+
+@pytest.mark.parametrize("stat", ["F", "B"])
+def test_tau_const_periodicity(stat):
+    """Test TauConst with statistics-dependent periodicity"""
+    beta = 10.0
+    tc = augment.TauConst(beta, stat)
+    
+    # Test at tau=0
+    val0 = tc(0.0)
+    assert np.isclose(val0, 1.0 / np.sqrt(beta))
+    
+    # Test at tau=5
+    val_pos = tc(5.0)
+    assert np.isclose(val_pos, 1.0 / np.sqrt(beta))
+    
+    # Test at tau=-5 (should apply periodicity)
+    val_neg = tc(-5.0)
+    
+    if stat == 'F':
+        # Fermionic: anti-periodic
+        assert np.isclose(val_neg, -1.0 / np.sqrt(beta))
+    else:
+        # Bosonic: periodic
+        assert np.isclose(val_neg, 1.0 / np.sqrt(beta))
+
+
+@pytest.mark.parametrize("stat", ["F", "B"])
+def test_tau_linear_periodicity(stat):
+    """Test TauLinear with statistics-dependent periodicity"""
+    beta = 10.0
+    tl = augment.TauLinear(beta, stat)
+    
+    # Test at tau=0
+    val0 = tl(0.0)
+    norm = np.sqrt(3 / beta)
+    assert np.isclose(val0, -norm)  # x = 2*0/beta - 1 = -1
+    
+    # Test at tau=5 (middle)
+    val_mid = tl(5.0)
+    assert np.isclose(val_mid, 0.0)  # x = 2*5/10 - 1 = 0
+    
+    # Test at tau=-5
+    val_neg = tl(-5.0)
+    # tau_normalized = -5 + 10 = 5, x = 2*5/10 - 1 = 0
+    
+    if stat == 'F':
+        # Fermionic: anti-periodic, sign = -1
+        assert np.isclose(val_neg, 0.0)  # -1 * 0 = 0
+    else:
+        # Bosonic: periodic, sign = +1
+        assert np.isclose(val_neg, 0.0)  # +1 * 0 = 0
+
+
+def test_matsubara_const_range():
+    """Test MatsubaraConst accepts [-β, β] range"""
+    beta = 10.0
+    mc = augment.MatsubaraConst(beta)
+    
+    # Should accept negative tau and return NaN
+    assert np.isnan(mc(-5.0))
+    assert np.isnan(mc(5.0))
+    assert np.isnan(mc(0.0))
+    assert np.isnan(mc(-beta))
+    assert np.isnan(mc(beta))
+    
+    # Should reject out of range
+    with pytest.raises(ValueError):
+        mc(-beta - 1)
+    with pytest.raises(ValueError):
+        mc(beta + 1)
+
+
+@pytest.mark.parametrize("stat", ["F", "B"])
+def test_tau_const_with_statistics(stat):
+    """Test TauConst can be created with statistics parameter"""
+    beta = 10.0
+    basis = sparse_ir.FiniteTempBasis(stat, beta, wmax=2.0, eps=1e-6)
+    
+    # Test factory method
+    tc = augment.TauConst.create(basis)
+    assert tc._statistics == stat
+    
+    # Test direct creation
+    tc2 = augment.TauConst(beta, stat)
+    assert tc2._statistics == stat
+    
+    # Test evaluation works
+    val = tc(5.0)
+    assert np.isfinite(val)
+
+
+@pytest.mark.parametrize("stat", ["F", "B"])
+def test_tau_linear_with_statistics(stat):
+    """Test TauLinear can be created with statistics parameter"""
+    beta = 10.0
+    basis = sparse_ir.FiniteTempBasis(stat, beta, wmax=2.0, eps=1e-6)
+    
+    # Test factory method
+    tl = augment.TauLinear.create(basis)
+    assert tl._statistics == stat
+    
+    # Test direct creation
+    tl2 = augment.TauLinear(beta, stat)
+    assert tl2._statistics == stat
+    
+    # Test evaluation works
+    val = tl(5.0)
+    assert np.isfinite(val)
+
+
+def test_backward_compatibility():
+    """Test backward compatibility for augmentation constructors"""
+    beta = 10.0
+    
+    # Old API (without statistics) should default to Bosonic
+    tc = augment.TauConst(beta)
+    assert tc._statistics == 'B'
+    
+    tl = augment.TauLinear(beta)
+    assert tl._statistics == 'B'
+    
+    # MatsubaraConst can be created without statistics
+    mc = augment.MatsubaraConst(beta)
+    # Statistics is optional for MatsubaraConst
+    assert mc._statistics is None or mc._statistics in ('F', 'B')