Remove recursion from Base.:^ (#618)

Co-authored-by: Martin Holters <martin.holters@hsu-hh.de>

Author: wheeheee

src/Filters/coefficients.jl

@@ -1,5 +1,7 @@
 # Filter types and conversions
 
+using Base: uabs
+
 abstract type FilterCoefficients{Domain} end
 
 Base.convert(::Type{T}, f::FilterCoefficients) where {T<:FilterCoefficients} = T(f)
@@ -45,11 +47,9 @@ Base.promote_rule(::Type{ZeroPoleGain{D,Z1,P1,K1}}, ::Type{ZeroPoleGain{D,Z2,P2,
 Base.inv(f::ZeroPoleGain{D}) where {D} = ZeroPoleGain{D}(f.p, f.z, inv(f.k))
 
 function Base.:^(f::ZeroPoleGain{D}, e::Integer) where {D}
-    if e < 0
-        return inv(f^-e)
-    else
-        return ZeroPoleGain{D}(repeat(f.z, e), repeat(f.p, e), f.k^e)
-    end
+    ae = uabs(e)
+    z, p = repeat(f.z, ae), repeat(f.p, ae)
+    return e < 0 ? ZeroPoleGain{D}(p, z, inv(f.k)^ae) : ZeroPoleGain{D}(z, p, f.k^ae)
 end
 
 #
@@ -184,11 +184,12 @@ end
 Base.inv(f::PolynomialRatio{D}) where {D} = PolynomialRatio{D}(f.a, f.b)
 
 function Base.:^(f::PolynomialRatio{D,T}, e::Integer) where {D,T}
+    ae = uabs(e)
+    b, a = f.b^ae, f.a^ae
     if e < 0
-        return PolynomialRatio{D}(f.a^-e, f.b^-e)
-    else
-        return PolynomialRatio{D}(f.b^e, f.a^e)
+        b, a = a, b
     end
+    return PolynomialRatio{D}(b, a)
 end
 
 coef_s(p::LaurentPolynomial) = p[end:-1:0]
@@ -304,7 +305,49 @@ Base.promote_rule(::Type{SecondOrderSections{D,T1,G1}}, ::Type{SecondOrderSectio
 
 SecondOrderSections{D,T,G}(f::SecondOrderSections) where {D,T,G} =
     SecondOrderSections{D,T,G}(f.biquads, f.g)
+SecondOrderSections{D,T,G}(f::Biquad{D}) where {D,T,G} = SecondOrderSections{D,T,G}([f], one(G))
+
 SecondOrderSections{D}(f::SecondOrderSections{D,T,G}) where {D,T,G} = SecondOrderSections{D,T,G}(f)
+SecondOrderSections{D}(f::Biquad{D,T}) where {D,T} = SecondOrderSections{D,T,Int}(f)
+
+*(f::SecondOrderSections{D}, g::Number) where {D} = SecondOrderSections{D}(f.biquads, f.g * g)
+*(g::Number, f::SecondOrderSections{D}) where {D} = SecondOrderSections{D}(f.biquads, f.g * g)
+*(f1::SecondOrderSections{D}, f2::SecondOrderSections{D}) where {D} =
+    SecondOrderSections{D}([f1.biquads; f2.biquads], f1.g * f2.g)
+*(f1::SecondOrderSections{D}, fs::SecondOrderSections{D}...) where {D} =
+    SecondOrderSections{D}(vcat(f1.biquads, map(f -> f.biquads, fs)...), f1.g * prod(f.g for f in fs))
+
+*(f1::Biquad{D}, f2::Biquad{D}) where {D} = SecondOrderSections{D}([f1, f2], 1)
+*(f1::Biquad{D}, fs::Biquad{D}...) where {D} = SecondOrderSections{D}([f1, fs...], 1)
+*(f1::SecondOrderSections{D}, f2::Biquad{D}) where {D} =
+    SecondOrderSections{D}([f1.biquads; f2], f1.g)
+*(f1::Biquad{D}, f2::SecondOrderSections{D}) where {D} =
+    SecondOrderSections{D}([f1; f2.biquads], f2.g)
+
+Base.inv(f::SecondOrderSections{D}) where {D} = SecondOrderSections{D}(inv.(f.biquads), inv(f.g))
+
+function Base.:^(f::SecondOrderSections{D}, e::Integer) where {D}
+    ae = uabs(e)
+    if e < 0
+        inv_f = inv(f)
+        return SecondOrderSections{D}(repeat(inv_f.biquads, ae), inv_f.g^ae)
+    else
+        return SecondOrderSections{D}(repeat(f.biquads, ae), f.g^ae)
+    end
+end
+
+function Base.:^(f::Biquad{D}, e::Integer) where {D}
+    ae = uabs(e)
+    return SecondOrderSections{D}(fill(e < 0 ? inv(f) : f, ae), 1)
+end
+
+function Biquad{D,T}(f::SecondOrderSections{D}) where {D,T}
+    if length(f.biquads) != 1
+        throw(ArgumentError("only a single second order section may be converted to a biquad"))
+    end
+    Biquad{D,T}(f.biquads[1] * f.g)
+end
+Biquad{D}(f::SecondOrderSections{D,T,G}) where {D,T,G} = Biquad{D,promote_type(T, G)}(f)
 
 function ZeroPoleGain{D,Z,P,K}(f::SecondOrderSections{D}) where {D,Z,P,K}
     z = Z[]
@@ -321,14 +364,6 @@ end
 ZeroPoleGain{D}(f::SecondOrderSections{D,T,G}) where {D,T,G} =
     ZeroPoleGain{D,complex(T),complex(T),G}(f)
 
-function Biquad{D,T}(f::SecondOrderSections{D}) where {D,T}
-    if length(f.biquads) != 1
-        throw(ArgumentError("only a single second order section may be converted to a biquad"))
-    end
-    Biquad{D,T}(f.biquads[1] * f.g)
-end
-Biquad{D}(f::SecondOrderSections{D,T,G}) where {D,T,G} = Biquad{D,promote_type(T,G)}(f)
-
 PolynomialRatio{D,T}(f::SecondOrderSections{D}) where {D,T} = PolynomialRatio{D,T}(ZeroPoleGain(f))
 PolynomialRatio{D}(f::SecondOrderSections{D}) where {D} = PolynomialRatio{D}(ZeroPoleGain(f))
 
@@ -447,39 +482,4 @@ end
 SecondOrderSections{D}(f::ZeroPoleGain{D,Z,P,K}) where {D,Z,P,K} =
     SecondOrderSections{D,promote_type(real(Z), real(P)), K}(f)
 
-
-SecondOrderSections{D,T,G}(f::Biquad{D}) where {D,T,G} = SecondOrderSections{D,T,G}([f], one(G))
-SecondOrderSections{D}(f::Biquad{D,T}) where {D,T} = SecondOrderSections{D,T,Int}(f)
 SecondOrderSections{D}(f::FilterCoefficients{D}) where {D} = SecondOrderSections{D}(ZeroPoleGain(f))
-
-*(f::SecondOrderSections{D}, g::Number) where {D} = SecondOrderSections{D}(f.biquads, f.g * g)
-*(g::Number, f::SecondOrderSections{D}) where {D} = SecondOrderSections{D}(f.biquads, f.g * g)
-*(f1::SecondOrderSections{D}, f2::SecondOrderSections{D}) where {D} =
-    SecondOrderSections{D}([f1.biquads; f2.biquads], f1.g * f2.g)
-*(f1::SecondOrderSections{D}, fs::SecondOrderSections{D}...) where {D} =
-    SecondOrderSections{D}(vcat(f1.biquads, map(f -> f.biquads, fs)...), f1.g * prod(f.g for f in fs))
-
-*(f1::Biquad{D}, f2::Biquad{D}) where {D} = SecondOrderSections{D}([f1, f2], 1)
-*(f1::Biquad{D}, fs::Biquad{D}...) where {D} = SecondOrderSections{D}([f1, fs...], 1)
-*(f1::SecondOrderSections{D}, f2::Biquad{D}) where {D} =
-    SecondOrderSections{D}([f1.biquads; f2], f1.g)
-*(f1::Biquad{D}, f2::SecondOrderSections{D}) where {D} =
-    SecondOrderSections{D}([f1; f2.biquads], f2.g)
-
-Base.inv(f::SecondOrderSections{D}) where {D} = SecondOrderSections{D}(inv.(f.biquads), inv(f.g))
-
-function Base.:^(f::SecondOrderSections{D}, e::Integer) where {D}
-    if e < 0
-        return inv(f)^-e
-    else
-        return SecondOrderSections{D}(repeat(f.biquads, e), f.g^e)
-    end
-end
-
-function Base.:^(f::Biquad{D}, e::Integer) where {D}
-    if e < 0
-        return inv(f)^-e
-    else
-        return SecondOrderSections{D}(fill(f, e), 1)
-    end
-end

test/filter_conversion.jl

@@ -219,7 +219,7 @@ end
         p = rand(ComplexF64, Npc) .- (0.5 + 0.5im);
         p = [p; conj(p); rand(Npr) .- 0.5; zeros(max(2Nzc+Nzr-2Npc-Npr, 0))]
         H′ = ZeroPoleGain(z, p,
-            (rand() + 0.5) * rand([-1, 1]), # non-zero gain with random sign
+            (rand() + 0.5) * rand((-1, 1)), # non-zero gain with random sign
         )
         maybe_biquad = length(z) ≤ 2 && length(p) ≤ 2 ? (Biquad,) : ()
         for T ∈ (PolynomialRatio, ZeroPoleGain, SecondOrderSections, maybe_biquad...)
@@ -236,6 +236,25 @@ end
             @test filt(H^0, x) ≈ x
         end
     end
+    # test that ^ doesn't cause stack overflow
+    let H = PolynomialRatio([1.0], [2.0])^typemin(Int8)
+        @test coefb(H) == [2.0^128]
+        @test coefa(H) == [1.0]
+    end
+    let zpg = ZeroPoleGain([1], [2], 3)^typemin(Int8)
+        @test length(zpg.z) == length(zpg.p) == 128
+        @test all(==(2), zpg.z)
+        @test all(==(1), zpg.p)
+        @test zpg.k ≈ inv(3)^128
+    end
+    let bq = Biquad(1:5...)
+        sos1 = bq^typemin(Int8)
+        sos2 = SecondOrderSections(bq)^typemin(Int8)
+        @test length(sos1.biquads) == length(sos2.biquads) == 128
+        @test all(==(inv(bq)), sos1.biquads)
+        @test all(==(inv(bq)), sos2.biquads)
+        @test sos1.g == sos2.g == 1
+    end
 end
 
 @testset "types" begin