diff --git a/src/ArchimaxCopula.jl b/src/ArchimaxCopula.jl
index 1d6b6d6df..165aa6ebc 100644
--- a/src/ArchimaxCopula.jl
+++ b/src/ArchimaxCopula.jl
@@ -40,10 +40,45 @@ ArchimaxCopula(d, gen::Generator, ::NoTail) = ArchimedeanCopula(d, gen)
 ArchimaxCopula(d, ::IndependentGenerator, tail::Tail) = ExtremeValueCopula(d, tail) 
 Distributions.params(C::ArchimaxCopula) = (_as_tuple(Distributions.params(C.gen))..., _as_tuple(Distributions.params(C.tail))...)
 
-# Fast conditional distortion binding (bivariate)
-DistortionFromCop(C::ArchimaxCopula{2}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, ::Int) = BivArchimaxDistortion(C.gen, C.tail, Int8(js[1]), float(uⱼₛ[1]))
+function _cdf(C::ArchimaxCopula{d}, u) where {d}
+    # Basic support checks
+    T = eltype(u)
+    any(iszero, u) && return T(0)
+    all(isone, u) && return T(1)
+    # Compute x_i = ϕ⁻¹(u_i), S = ∑ x_i, ω = x/S, then ϕ(S·A(ω))
+    x = ϕ⁻¹.(C.gen, u)
+    S = sum(x)
+    S == 0 && return T(1)
+    ω = ntuple(i -> x[i] / S, d)
+    return ϕ(C.gen, S * A(C.tail, ω))
+end
+function Distributions._logpdf(C::ArchimaxCopula{d, TG, TT}, u) where {d, TG, TT}
+    @inbounds for ui in u
+        (0.0 < ui < 1.0) || return -Inf
+    end
+    val = _der(v -> Distributions.cdf(C, v), collect(u), ntuple(identity, d))
+    return log(max(val, 0))
+end
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{d, TG, TT}, X::AbstractMatrix{T}) where {T<:Real, d, TG, TT}
+    d == 2 && return @invoke Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{2, TG, TT}, X)
+    @assert size(X, 1) == d
+    U = rand(rng, Distributions.Uniform(), d, size(X, 2))
+    X .= inverse_rosenblatt(C, U)
+    return X
+end
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{d, TG, TT}, x::AbstractVector{T}) where {T<:Real, d, TG, TT}
+    d == 2 && return @invoke Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{2, TG, TT}, x)
+    u = rand(rng, Distributions.Uniform(), d)
+    x .= inverse_rosenblatt(C, u)
+    return x
+end
 
-# --- CDF ---
+
+
+
+
+
+###### Special methods for the bivariate cases
 function _cdf(C::ArchimaxCopula{2}, u)
     u1, u2 = u
     (0.0 ≤ u1 ≤ 1.0 && 0.0 ≤ u2 ≤ 1.0) || return 0.0
@@ -57,8 +92,6 @@ function _cdf(C::ArchimaxCopula{2}, u)
     t = _safett(y / S)                 # protect t≈0,1
     return ϕ(C.gen, S * A(C.tail, t))
 end
-
-# --- log-PDF stable ---
 function Distributions._logpdf(C::ArchimaxCopula{2, TG, TT}, u) where {TG, TT}
     T = promote_type(Float64, eltype(u))
     @assert length(u) == 2
@@ -90,17 +123,6 @@ function Distributions._logpdf(C::ArchimaxCopula{2, TG, TT}, u) where {TG, TT}
     base > 0 || return T(-Inf)
     return T(log(φpp) + log(base))
 end
-
-# --- Kendall τ: τ = τ_A + (1 - τ_A) τ_ψ ---
-τ(C::ArchimaxCopula) = begin
-    τA = τ(ExtremeValueCopula(2, C.tail))
-    τψ = τ(C.gen)
-    τA + (1 - τA) * τψ
-end
-
-
-# Use the matrix sampler for better efficiency
-# (if not working, maybe uncomment the vetor version ?)
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{2, TG, TT}, A::DenseMatrix{T}) where {T<:Real, TG, TT}
     evcop, frail = ExtremeValueCopula(2, C.tail), frailty(C.gen)
     Distributions._rand!(rng, evcop, A)
@@ -116,6 +138,12 @@ function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{
     x[2] = ϕ(C.gen, -log(v2)/M)
     return x
 end
+DistortionFromCop(C::ArchimaxCopula{2}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, ::Int) = BivArchimaxDistortion(C.gen, C.tail, Int8(js[1]), float(uⱼₛ[1]))
+τ(C::ArchimaxCopula{2, TG, TT})  where {TG, TT} = begin
+    τA = τ(ExtremeValueCopula(2, C.tail))
+    τψ = τ(C.gen)
+    τA + (1 - τA) * τψ
+end
 
 
 
diff --git a/src/ExtremeValueCopula.jl b/src/ExtremeValueCopula.jl
index 18e2d0ed4..e3b2d9869 100644
--- a/src/ExtremeValueCopula.jl
+++ b/src/ExtremeValueCopula.jl
@@ -45,7 +45,29 @@ end
 _cdf(C::ExtremeValueCopula{d, TT}, u) where {d, TT} = exp(-ℓ(C.tail, .- log.(u)))
 Distributions.params(C::ExtremeValueCopula) = Distributions.params(C.tail)
 
-#### Restriction to bivariate cases of the following methods: 
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ExtremeValueCopula{d, TT}, X::AbstractMatrix{T}) where {T<:Real, d, TT}
+    @assert size(X, 1) == d
+    U = rand(rng,  d, size(X, 2))
+    X .= inverse_rosenblatt(C, U)
+    return X
+end
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ExtremeValueCopula{d, TT}, x::AbstractVector{T}) where {T<:Real, d, TT}
+    u = rand(rng, d)
+    x .= inverse_rosenblatt(C, u)
+    return x
+end
+function Distributions._logpdf(C::ExtremeValueCopula{d, TT}, u) where {d, TT}
+    # domain checks
+    @inbounds for ui in u
+        (0.0 < ui < 1.0) || return -Inf
+    end
+    # Compute mixed partial ∂^d/∂u₁…∂u_d of the cdf via nested ForwardDiff
+    val = _der(v -> Distributions.cdf(C, v), collect(u), 1:d)
+    return log(max(val, 0))
+end
+
+
+###### Restriction to bivariate cases of the following methods: 
 function Distributions._logpdf(C::ExtremeValueCopula{2, TT}, u) where {TT}
     u1, u2 = u
     (0.0 < u1 ≤ 1.0 && 0.0 < u2 ≤ 1.0) || return -Inf 
@@ -86,4 +108,4 @@ function Distributions._rand!(rng::Distributions.AbstractRNG, C::ExtremeValueCop
     x[2] = exp(log(w)*(1-z)/a)
     return x
 end
-DistortionFromCop(C::ExtremeValueCopula{2, TT}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, ::Int) where TT = BivEVDistortion(C.tail, Int8(js[1]), float(uⱼₛ[1]))
+DistortionFromCop(C::ExtremeValueCopula{d, TT}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, ::Int) where {d, TT} = BivEVDistortion(C.tail, Int8(js[1]), float(uⱼₛ[1]))
\ No newline at end of file
diff --git a/src/Tail.jl b/src/Tail.jl
index f4cd8c5d7..06fde43d0 100644
--- a/src/Tail.jl
+++ b/src/Tail.jl
@@ -33,21 +33,18 @@ References:
 abstract type Tail end
 Base.broadcastable(tail::Tail) = Ref(tail)
 
-####### Functions you need to overload: 
+####### Main interface for any dimensions: 
 _is_valid_in_dim(tail::Tail, d::Int) = throw(ArgumentError("Validity of the tail type $(typeof(tail)) must be supplied by overwriting the function _is_valid_in_dim(tail::Tail, d::Int)"))
 A(::Tail, ω::NTuple{d,<:Real}) where {d} = throw(ArgumentError("Implement A(Tail{$d}, ω) en el simplex Δ_{d-1}"))
-
-####### Rest of the interface you can overload if more efficient:
-needs_binary_search(::Tail) = false
-# \ell function
 function ℓ(tail::Tail, x)
     s = sum(x)
     return s == 0 ? zero(eltype(x)) : s * A(tail, ntuple(i->x[i]/s, length(x)))
 end
 
 
-# A more friendly interface for models that are only bivariate: 
+####### A more friendly interface for models that are only bivariate: 
 abstract type Tail2 <: Tail end
+needs_binary_search(::Tail2) = false
 _is_valid_in_dim(::Tail2, d::Int) = (d==2)
 A(tail::Tail2, t::NTuple{2, <:Real}) = A(tail, t[1])
 dA(tail::Tail2, t::Real) = ForwardDiff.derivative(z -> A(tail, z), t)
diff --git a/src/Tail/GalambosTail.jl b/src/Tail/GalambosTail.jl
index a67b5daaa..97c0044ee 100644
--- a/src/Tail/GalambosTail.jl
+++ b/src/Tail/GalambosTail.jl
@@ -38,16 +38,20 @@ end
 const GalambosCopula{T} = ExtremeValueCopula{2, GalambosTail{T}}
 GalambosCopula(θ) =ExtremeValueCopula(2, GalambosTail(θ))
 Distributions.params(tail::GalambosTail) = (tail.θ,)
+_is_valid_in_dim(::GalambosTail, d::Int) = (d >= 2)
 
+function A(tail::GalambosTail, ω::NTuple{d,<:Real}) where {d}
+    θ = tail.θ
+    θ == 0 && return 1.0
+    isinf(θ) && return maximum(ω)
+    return -LogExpFunctions.expm1(-LogExpFunctions.logsumexp(-θ .* log.(ω))/θ)  # 1 - (∑ ω_i^{-θ})^{-1/θ}
+end
+
+#### Special bindings for dimension d == 2
 needs_binary_search(tail::GalambosTail) = (tail.θ > 19.5)
 function A(tail::GalambosTail, t::Real)
     tt = _safett(t)
-    θ  = tail.θ
-    if θ == 0
-        return 1.0
-    elseif isinf(θ)
-        return max(tt, 1-tt)
-    else
-        return -LogExpFunctions.expm1(-LogExpFunctions.logaddexp(-θ*log(tt), -θ*log(1-tt)) / θ)
-    end
-end
+    tail.θ == 0 && return 1.0
+    isinf(tail.θ) && return max(tt, 1-tt)
+    return -LogExpFunctions.expm1(-LogExpFunctions.logaddexp(-tail.θ*log(tt), -tail.θ*log(1-tt)) / tail.θ)
+end
\ No newline at end of file
diff --git a/src/Tail/LogTail.jl b/src/Tail/LogTail.jl
index cbc575a3f..2f612c49e 100644
--- a/src/Tail/LogTail.jl
+++ b/src/Tail/LogTail.jl
@@ -37,22 +37,230 @@ end
 const LogCopula{T} = ExtremeValueCopula{2, LogTail{T}}
 LogCopula(θ) = ExtremeValueCopula(2, LogTail(θ))
 Distributions.params(tail::LogTail) = (tail.θ,)
+_is_valid_in_dim(::LogTail, d::Int) = (d >= 2)
 
-function ℓ(tail::LogTail, t)
-    t₁, t₂ = t
+"""
+    A(tail::LogTail, ω::NTuple{d,Real}) where d
+
+Multi-dimensional Pickands dependence function on the simplex Δ_{d-1} for the
+logistic (a.k.a. Gumbel) extreme-value model with parameter θ ≥ 1.
+
+    A(ω) = (∑_{i=1}^d ω_i^θ)^{1/θ},   ω_i ≥ 0,  ∑ ω_i = 1.
+
+We implement this in a numerically stable manner using a log-sum-exp style
+aggregation on the log(ω_i) scale to mitigate underflow when some ω_i are very
+small.
+"""
+function A(tail::LogTail, ω::NTuple{d,<:Real}) where {d}
     θ = tail.θ
-    return (t₁^θ + t₂^θ)^(1/θ)
+    @inbounds begin
+        # Handle potential exact zeros (allowed at the boundary of the simplex)
+        # ω_i^θ = 0 if ω_i == 0 (θ ≥ 1). If all but one entry are zero the sum is 1.
+        # Fast path: if one coordinate is 1 (degenerate vertex) return 1.
+        for ωi in ω
+            if ωi == 1.0
+                return one(ωi)
+            end
+        end
+        # Collect scaled logs; skip zeros to avoid -Inf + later exp
+        logs = similar(ntuple(_->0.0, d))  # temporary tuple-like container
+        mx = -Inf
+        nz = 0
+        for i in 1:d
+            ωi = ω[i]
+            if ωi > 0
+                val = θ * log(ωi)
+                logs = Base.setindex(logs, val, i)
+                if val > mx; mx = val; end
+                nz += 1
+            else
+                logs = Base.setindex(logs, -Inf, i)
+            end
+        end
+        nz == 0 && return one(eltype(ω)) # degenerate (should not happen if ∑ ω_i =1)
+        s = 0.0
+        for i in 1:d
+            li = logs[i]
+            @inbounds if isfinite(li)
+                s += exp(li - mx)
+            end
+        end
+        return exp((mx + log(s)) / θ)
+    end
 end
 
-# A(t) for LogCopula (avec log-exp pour la stabilité)
-function A(tail::LogTail, t::Real)
-    θ = tail.θ
-    # log-sum-exp trick: log(t^θ + (1-t)^θ) = logsumexp(θ*log(t), θ*log1p(-t))
-    logB = LogExpFunctions.logaddexp(θ*log(t), θ*log1p(-t))
-    return exp(logB / θ)
+
+# Placeholder for future optimized logistic sampler (Dirichlet / positive stable approach)
+# Currently falls back to generic (inverse Rosenblatt) sampler above. Once validated
+# we can replace by an O(d) method.
+
+function Distributions._logpdf(C::ExtremeValueCopula{d, LogTail{Tθ}}, u) where {d, Tθ}
+    # Analytic logistic (multivariate Gumbel / LogTail) pdf for d ≥ 3
+    # ℓ(x) = (∑ x_i^θ)^{1/θ}, θ = C.tail.θ ≥ 1.
+    # Density: c(u) = exp(-ℓ(x)) * (1/∏ u_i) * (∏ x_i^{θ-1}) * (∑ x_i^θ)^{1/θ - d} * A_d(θ)
+    # where A_d(θ) obtained from Bell-type recursion:
+    #   F_k(θ) = falling factorial θ(θ-1)…(θ-k+1)
+    #   A[0]=1; A[n] = Σ_{k=1}^n (-1) * F_k(θ) * A[n-k] * binomial(n-1,k-1).
+    # For θ=1 (independence), return logpdf = 0.
+    # We guard against boundary (any u_i≈1 ⇒ x_i≈0). If θ>1 and some x_i == 0 ⇒ density → 0.
+    
+
+    d == 2 && return @invoke Distributions._logpdf(C::ExtremeValueCopula{2, LogTail{Tθ}}, u)
+    θ = C.tail.θ
+    # Domain check
+    @inbounds for ui in u
+        (0.0 < ui < 1.0) || return -Inf
+    end
+    if θ == 1
+        # Independence (limit case)
+        return 0.0
+    end
+    x = @inbounds (-log.(u))
+    # If any x_i == 0 (u_i==1) density -> 0
+    any(iszero, x) && return -Inf
+    # Core sums
+    xθ = map(y -> y^θ, x)
+    s = zero(eltype(xθ))
+    @inbounds for v in xθ; s += v; end
+    s == 0 && return -Inf
+    ℓx = s^(1/θ)
+    # product of x_i^{θ-1}
+    prod_r = zero(eltype(x)) + 1
+    @inbounds for xi in x
+        prod_r *= xi^(θ - 1)
+    end
+    # Recurrence for A_d(θ)
+    # falling factorial cache
+    F = Vector{eltype(x)}(undef, d)
+    F[1] = θ
+    @inbounds for k in 2:d
+        F[k] = F[k-1] * (θ - (k - 1))
+    end
+    A = Vector{eltype(x)}(undef, d+1)
+    A[1] = one(eltype(x)) # A[1] represents A_0
+    # Shifted indexing: A[n+1] stores A_n
+    for n in 1:d
+        acc = zero(eltype(x))
+        @inbounds for k in 1:n
+            w_k = F[k]
+            acc += (-1) * w_k * A[n - k + 1] * binomial(n - 1, k - 1)
+        end
+        A[n+1] = acc
+    end
+    A_dθ = A[d+1]
+    A_dθ <= 0 && return @invoke Distributions._logpdf(C::ExtremeValueCopula{d, Tail}, u) # fallback to AD if numerically unstable
+    logc = -ℓx - sum(log.(u)) + (θ - 1) * sum(log.(x)) + (1/θ - d) * log(s) + log(A_dθ)
+    return logc
 end
 
-# Première dérivée dA/dt (stable numériquement)
+## ---------------------------------------------------------------------------
+## Specialized sampler for logistic (LogTail) extreme value copula (any d)
+## Based on spectral measure representation: If E_i ~ Exp(1) independent and
+## G_i ~ Gamma(1/θ, 1) then set W_i = G_i / Σ G_i; define S = (Σ E_i / W_i)^{1/θ}.
+## Then U_i = exp(-(S * (E_i / W_i)^{1/θ})) has the logistic EV copula with θ.
+## (Derivation aligns with standard spectral construction for multivariate Gumbel.)
+## ---------------------------------------------------------------------------
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ExtremeValueCopula{d, LogTail{Tθ}}, X::AbstractMatrix{T}) where {d, Tθ, T<:Real}
+    @assert size(X,1) == d
+    θ = C.tail.θ
+    invθ = 1/θ
+    n = size(X,2)
+    @inbounds for col in 1:n
+        # Draw Dirichlet(1/θ,...,1/θ) via Gamma
+        gsum = zero(Float64)
+        for i in 1:d
+            gi = rand(rng, Distributions.Gamma(invθ))
+            X[i,col] = gi
+            gsum += gi
+        end
+        # Normalize to W_i
+        for i in 1:d
+            X[i,col] /= gsum
+        end
+        # Exponential draws
+        total = zero(Float64)
+        for i in 1:d
+            Ei = rand(rng)
+            while Ei <= 0.0; Ei = rand(rng); end # ensure >0
+            # store temporarily Ei in-place scaled later; reuse column
+            # Compute contribution: Ei / W_i
+            Wi = X[i,col]
+            Xi = Ei / Wi
+            X[i,col] = Xi  # hold Xi
+            total += Xi
+        end
+        S = total^invθ
+        for i in 1:d
+            Xi = X[i,col]
+            # Xi currently = Ei / W_i
+            X[i,col] = exp(-(S * Xi^invθ))
+        end
+    end
+    return X
+end
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ExtremeValueCopula{d, LogTail{Tθ}}, x::AbstractVector{T}) where {d, Tθ, T<:Real}
+    θ = C.tail.θ; invθ = 1/θ
+    # Dirichlet via gamma
+    gsum = 0.0
+    for i in 1:d
+        gi = rand(rng, Distributions.Gamma(invθ))
+        x[i] = gi
+        gsum += gi
+    end
+    for i in 1:d
+        x[i] /= gsum
+    end
+    total = 0.0
+    for i in 1:d
+        Ei = rand(rng)
+        while Ei <= 0.0; Ei = rand(rng); end
+        Wi = x[i]
+        Xi = Ei / Wi
+        x[i] = Xi
+        total += Xi
+    end
+    S = total^invθ
+    for i in 1:d
+        Xi = x[i]
+        x[i] = exp(-(S * Xi^invθ))
+    end
+    return x
+end
+
+
+## ---------------------------------------------------------------------------
+## Logistic (LogTail) multi-dimensional distortion specialization (p=1)
+## U_i | U_j = u_j has cdf F_{i|j}(u) = ∂_j C(u_i, u_j, 1, ...,1)/∂_j C(1, u_j,1,...,1)
+## For an extreme-value copula C(u) = exp(-ℓ(-log u)), with ℓ homogeneous of
+## order 1: ∂_j C = C * ( (∂_j ℓ)/u_j ). When fixing other arguments at 1,
+## x_k = -log 1 = 0 removes those coordinates. For the logistic model
+## ℓ(x) = (x_1^θ + x_2^θ)^{1/θ} if only coords i,j vary. We thus reduce to
+## the bivariate distortion with effective 2D tail (same θ), so we can reuse
+## the existing fast BivEVDistortion logic by mapping indices.
+## ---------------------------------------------------------------------------
+
+@inline function DistortionFromCop(C::ExtremeValueCopula{d, LogTail{Tθ}}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, i::Int) where {d, Tθ}
+    # For p=1 we only need a bivariate slice involving coordinates (i, j)
+    j = js[1]
+    if i == j
+        throw(ArgumentError("Cannot build distortion for identical conditioned and target index"))
+    end
+    # Reuse the existing BivEVDistortion with the same tail; orientation matters:
+    # BivEVDistortion expects j∈{1,2} as which coordinate is conditioned.
+    # We map (i,j) in original d-D space to a virtual 2D with ordering (free, conditioned).
+    # If original conditioned index is j, we set j_virtual = 2.
+    return BivEVDistortion(C.tail, Int8(2), float(uⱼₛ[1]))
+end
+
+
+##### Special binding for dim 2 
+## ---------------------------------------------------------------------------
+## Bivariate specializations (for performance / numerical stability)
+## ---------------------------------------------------------------------------
+A(tail::LogTail, t::Real) = begin
+    θ = tail.θ
+    return exp(LogExpFunctions.logaddexp(θ*log(t), θ*log1p(-t)) / θ)
+end
 function dA(tail::LogTail, t::Real)
     θ = tail.θ
 
@@ -72,8 +280,6 @@ function dA(tail::LogTail, t::Real)
 
     return Bpow * D
 end
-
-# Seconde dérivée d²A/dt² (stable numériquement)
 function d2A(tail::LogTail, t::Real)
     θ = tail.θ
 
@@ -107,4 +313,5 @@ function d2A(tail::LogTail, t::Real)
     term2 = exp((1 - θ) / θ * logB) * E
 
     return term1 + term2
-end
\ No newline at end of file
+end
+
diff --git a/test/HighDimEV_Archimax.jl b/test/HighDimEV_Archimax.jl
new file mode 100644
index 000000000..5d471576c
--- /dev/null
+++ b/test/HighDimEV_Archimax.jl
@@ -0,0 +1,56 @@
+@testitem "High-d EV and Archimax basics" tags=[:HighDim, :ExtremeValueCopula, :ArchimaxCopula] setup=[M] begin
+    using Copulas, Distributions, StableRNGs, Test
+    rng = StableRNG(42)
+
+    # 3D Logistic EV (aka Gumbel EV) via LogTail(θ)
+    θs = (1.2, 2.5)
+    for θ in θs
+        Cev = Copulas.ExtremeValueCopula(3, Copulas.LogTail(θ))
+        U = rand(rng, Cev, 200)
+        @test all(0 .<= U .<= 1)
+        @test 0.0 ≤ cdf(Cev, [0.5,0.5,0.5]) ≤ 1.0
+        @test isfinite(logpdf(Cev, [0.3,0.6,0.7]))
+        # small smoke on rosenblatt ∘ inverse
+        U2 = Copulas.inverse_rosenblatt(Cev, Copulas.rosenblatt(Cev, U))
+        @test isapprox(U, U2; atol=1e-2)
+    end
+
+    # 4D Galambos EV (negative logistic)
+    for θ in (0.7, 2.0)
+        Cev = Copulas.ExtremeValueCopula(4, Copulas.GalambosTail(θ))
+        U = rand(rng, Cev, 100)
+        @test all(0 .<= U .<= 1)
+        @test 0.0 ≤ cdf(Cev, fill(0.5, 4)) ≤ 1.0
+        @test isfinite(logpdf(Cev, fill(0.6, 4)))
+    end
+
+    # 3D Archimax: Clayton × Logistic tail
+    for (θg, θt) in ((1.5, 1.3), (3.0, 2.0))
+        Cax = Copulas.ArchimaxCopula(3, Copulas.ClaytonGenerator(θg), Copulas.LogTail(θt))
+        U = rand(rng, Cax, 150)
+        @test all(0 .<= U .<= 1)
+        @test 0.0 ≤ cdf(Cax, [0.2,0.7,0.5]) ≤ 1.0
+        @test isfinite(logpdf(Cax, [0.4,0.6,0.8]))
+    end
+end
+
+@testitem "Logistic EV specialized sampler & distortion" tags=[:ExtremeValueCopula, :LogisticSampler] setup=[M] begin
+    using Copulas, Distributions, StableRNGs, Test, Statistics
+    rng = StableRNG(84)
+    for d in (3,5)
+        θ = 2.5
+        C = Copulas.ExtremeValueCopula(d, Copulas.LogTail(θ))
+        n = 1500
+        U = rand(rng, C, n)
+        @test all(0 .<= U .<= 1)
+        # Basic marginal sanity: mean(U_i) should be close to 0.5 (uniform margins)
+        μ = mean(U[1,:])
+        @test 0.4 < μ < 0.6
+        # Distortion specialization check
+        j = 1; i = 2
+        u_j = U[j, 10]
+        D = Copulas.DistortionFromCop(C, (j,), (u_j,), i)
+        val = cdf(D, U[i,10])
+        @test 0.0 <= val <= 1.0
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 46fb7b62e..6e4f82176 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,4 +4,4 @@ using TestItemRunner
 # @run_package_tests filter=ti->(:GumbelBarnettCopula in ti.tags || :ArchimedeanCopula in ti.tags || :FrankCopula in ti.tags)
 # you can add verbose=true here 
 
-@run_package_tests
\ No newline at end of file
+@run_package_tests filter=ti->(:HighDim in ti.tags)
\ No newline at end of file