Use Half{Odd{Int}} in Jacobi-Ultraspherical conversions

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.28"
+version = "0.6.29"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
@@ -11,8 +11,10 @@ DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OddEvenIntegers = "8d37c425-f37a-4ca2-9b9d-a61bc06559d2"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
@@ -30,8 +32,10 @@ DomainSets = "0.5, 0.6"
 FastGaussQuadrature = "0.4, 0.5"
 FastTransforms = "0.12, 0.13, 0.14, 0.15.1"
 FillArrays = "0.11, 0.12, 0.13, 1"
+HalfIntegers = "1.5"
 IntervalSets = "0.5, 0.6, 0.7"
 LazyArrays = "0.22, 1"
+OddEvenIntegers = "0.1"
 Reexport = "0.2, 1"
 SpecialFunctions = "0.10, 1.0, 2"
 Static = "0.8"

src/ApproxFunOrthogonalPolynomials.jl

@@ -85,6 +85,9 @@ using StaticArrays: SVector
 
 import LinearAlgebra: isdiag, eigvals, eigen
 
+using OddEvenIntegers
+using HalfIntegers
+
 export bandmatrices_eigen, SymmetricEigensystem, SkewSymmetricEigensystem
 
 points(d::IntervalOrSegmentDomain{T},n::Integer) where {T} =
@@ -134,6 +137,44 @@ compare_op() = ≈
 compare_op(::Any, args...) = compare_op(args...)
 compare_orders(a::Number, b::Number) = compare_op(a, b)(a, b)
 
+# work around type promotions to preserve types for StepRanges involving HalfOddIntegers with a unit step
+const HalfOddInteger{T<:Integer} = Half{Odd{T}}
+decreasingunitsteprange(start, stop) = start:-1:stop
+decreasingunitsteprange(start::HalfOddInteger, stop::Integer) = start:-1:oftype(start, stop - half(1))
+decreasingunitsteprange(start::Integer, stop::HalfOddInteger) = start:-1:oftype(start, stop - half(1))
+
+
+# return 1/2, possibly preserving types but not being too fussy
+_onehalf(x) = onehalf(x)
+_onehalf(::Union{StaticInt, StaticFloat64}) = static(0.5)
+_onehalf(::Integer) = half(Odd(1))
+
+# return -1/2, possibly preserving types but not being too fussy
+_minonehalf(x) = -onehalf(x)
+_minonehalf(@nospecialize(_::Union{StaticInt, StaticFloat64})) = static(-0.5)
+_minonehalf(::Integer) = half(Odd(-1))
+
+isapproxminhalf(a) = a ≈ _minonehalf(a)
+isapproxminhalf(@nospecialize(a::StaticInt)) = isequalminhalf(a)
+isapproxminhalf(::Integer) = false
+
+isequalminhalf(x) = x == _minonehalf(x)
+isequalminhalf(@nospecialize(a::StaticInt)) = false
+isequalminhalf(::Integer) = false
+
+isequalhalf(x) = x == _onehalf(x)
+isequalhalf(@nospecialize(a::StaticInt)) = false
+isequalhalf(x::Integer) = false
+
+isapproxinteger_addhalf(a) = !isapproxinteger(a) && isapproxinteger(a+_onehalf(a))
+isapproxinteger_addhalf(a::HalfOddInteger) = true
+isapproxinteger_addhalf(@nospecialize(a::StaticInt)) = false
+function isapproxinteger_addhalf(a::StaticFloat64)
+    x = mod(a, static(1))
+    x == _onehalf(x) || dynamic(x) ≈ 0.5
+end
+isapproxinteger_addhalf(::Integer) = false
+
 include("bary.jl")
 
 

src/Spaces/Jacobi/Jacobi.jl

@@ -1,6 +1,5 @@
 export Jacobi, NormalizedJacobi, Legendre, NormalizedLegendre
 
-
 """
 `Jacobi(b,a)` represents the space spanned by Jacobi polynomials `P_k^{(a,b)}`,
 which are orthogonal with respect to the weight `(1+x)^β*(1-x)^α`
@@ -17,8 +16,15 @@ Legendre(domain = ChebyshevInterval()) = Jacobi(0,0,Domain(domain)::Domain)
 Legendre(s::PolynomialSpace) = Legendre(Jacobi(s))
 Legendre(s::Jacobi) = s.a == s.b == 0 ? s : throw(ArgumentError("can't convert $s to Legendre"))
 Jacobi(b::Number,a::Number,d=ChebyshevInterval()) = Jacobi(promote(b, a)..., Domain(d)::Domain)
-Jacobi(A::Ultraspherical) = Jacobi(order(A)-0.5,order(A)-0.5,domain(A))
-Jacobi(A::Chebyshev) = Jacobi(-0.5,-0.5,domain(A))
+function Jacobi(A::Ultraspherical)
+    m = order(A)
+    n = m + half(Odd(-1))
+    Jacobi(n,n,domain(A))
+end
+function Jacobi(A::Chebyshev)
+    n = half(Odd(-1))
+    Jacobi(n,n,domain(A))
+end
 Jacobi(A::Jacobi) = A
 
 const NormalizedJacobi{D<:Domain,R,T} = NormalizedPolynomialSpace{Jacobi{D,R,T},D,R}
@@ -44,7 +50,7 @@ end
 
 function Ultraspherical(J::Jacobi)
     if J.a == J.b
-        Ultraspherical(J.a+0.5,domain(J))
+        Ultraspherical(J.a+_onehalf(J.a),domain(J))
     else
         error("Cannot construct Ultraspherical with a=$(J.a) and b=$(J.b)")
     end
@@ -62,8 +68,6 @@ convert(::Type{Jacobi{D,R1,T1}},J::Jacobi{D,R2,T2}) where {D,R1,R2,T1,T2} =
 compare_orders((Aa, Ba)::NTuple{2,Number}, (Ab, Bb)::NTuple{2,Number}) = compare_orders(Aa, Ba) && compare_orders(Ab, Bb)
 spacescompatible(a::Jacobi, b::Jacobi) = compare_orders((a.a, b.a), (a.b, b.b)) && domainscompatible(a,b)
 
-isapproxinteger_addhalf(a) = !isapproxinteger(a) && isapproxinteger(a+0.5)
-isapproxinteger_addhalf(::Integer) = false
 function canonicalspace(S::Jacobi)
     if isapproxinteger_addhalf(S.a) && isapproxinteger_addhalf(S.b)
         Chebyshev(domain(S))

src/Spaces/Jacobi/JacobiOperators.jl

@@ -158,8 +158,8 @@ function Conversion(L::Jacobi,M::Jacobi)
             # Conversion(L, M) == Conversion(J, M) * Conversion(L, J)
             # Conversion(L, J) = Conversion(Jacobi(L.b, L.a, dm), Jacobi(M.b, L.a, dm))
             # Conversion(J, M) = Conversion(Jacobi(M.b, L.a, dm), Jacobi(M.b, M.a, dm))
-            CLJ = [ConcreteConversion(Jacobi(b-1,L.a,dm), Jacobi(b, L.a, dm)) for b in M.b:-1:L.b+1]
-            CJM = [ConcreteConversion(Jacobi(M.b,a-1,dm), Jacobi(M.b, a, dm)) for a in M.a:-1:L.a+1]
+            CLJ = [ConcreteConversion(Jacobi(b-1,L.a,dm), Jacobi(b, L.a, dm)) for b in decreasingunitsteprange(M.b, L.b+1)]
+            CJM = [ConcreteConversion(Jacobi(M.b,a-1,dm), Jacobi(M.b, a, dm)) for a in decreasingunitsteprange(M.a, L.a+1)]
             C = [CJM; CLJ]
             return ConversionWrapper(TimesOperator(C))
         end
@@ -221,102 +221,98 @@ end
 
 # Assume m is compatible
 
-function Conversion(A::PolynomialSpace,B::Jacobi)
+function Conversion(A::PolynomialSpace, B::Jacobi)
     @assert domain(A) == domain(B)
     J = Jacobi(A)
     J == B ? ConcreteConversion(A,B) :
-             ConversionWrapper(TimesOperator(Conversion(J,B),Conversion(A,J)))
+             ConversionWrapper(SpaceOperator(TimesOperator(Conversion(J,B),Conversion(A,J)), A, B))
 end
 
-function Conversion(A::Jacobi,B::PolynomialSpace)
+function Conversion(A::Jacobi, B::PolynomialSpace)
     @assert domain(A) == domain(B)
     J = Jacobi(B)
     J == A ? ConcreteConversion(A,B) :
-             ConversionWrapper(TimesOperator(Conversion(J,B),Conversion(A,J)))
+             ConversionWrapper(SpaceOperator(TimesOperator(Conversion(J,B),Conversion(A,J)), A, B))
 end
 
-isequalminhalf(x) = x == -0.5
-isequalminhalf(@nospecialize ::Integer) = false
-
-function Conversion(A::Jacobi,B::Chebyshev)
+function Conversion(A::Jacobi, B::Chebyshev)
     @assert domain(A) == domain(B)
     if isequalminhalf(A.a) && isequalminhalf(A.b)
         ConcreteConversion(A,B)
     elseif A.a == A.b == 0
-        ConversionWrapper(
-            SpaceOperator(
-                ConcreteConversion(Ultraspherical(1//2),B),
-                A,B))
+        ConversionWrapper(SpaceOperator(ConcreteConversion(Ultraspherical(A), B), A, B))
     elseif A.a == A.b
         US = Ultraspherical(A)
-        ConversionWrapper(Conversion(US,B)*ConcreteConversion(A,US))
+        ConversionWrapper(SpaceOperator(TimesOperator(Conversion(US,B), ConcreteConversion(A,US)), A, B))
     else
         J = Jacobi(B)
         ConcreteConversion(J,B)*Conversion(A,J)
     end
 end
 
-function Conversion(A::Chebyshev,B::Jacobi)
+function Conversion(A::Chebyshev, B::Jacobi)
     @assert domain(A) == domain(B)
     if isequalminhalf(B.a) && isequalminhalf(B.b)
         ConcreteConversion(A,B)
     elseif B.a == B.b == 0
-        ConversionWrapper(
-            SpaceOperator(
-                ConcreteConversion(A,Ultraspherical(1//2,domain(B))),
-                A,B))
+        ConversionWrapper(SpaceOperator(ConcreteConversion(A, Ultraspherical(B)), A, B))
     elseif B.a == B.b
         US = Ultraspherical(B)
-        ConcreteConversion(US,B) * Conversion(A,US)
+        ConversionWrapper(SpaceOperator(TimesOperator(ConcreteConversion(US,B), Conversion(A,US)), A, B))
     else
         J = Jacobi(A)
         Conversion(J,B)*ConcreteConversion(A,J)
     end
 end
 
 
-function Conversion(A::Jacobi,B::Ultraspherical)
+@inline function _Conversion(A::Jacobi, B::Ultraspherical)
     @assert domain(A) == domain(B)
     if isequalminhalf(A.a) && isequalminhalf(A.b)
-        ConversionWrapper(Conversion(Chebyshev(domain(A)),B)*
-            ConcreteConversion(A,Chebyshev(domain(A))))
+        C = Chebyshev(domain(A))
+        ConversionWrapper(SpaceOperator(
+            TimesOperator(Conversion(C,B), ConcreteConversion(A,C)), A, B))
     elseif isequalminhalf(A.a - order(B)) && isequalminhalf(A.b - order(B))
         ConcreteConversion(A,B)
     elseif A.a == A.b == 0
-        ConversionWrapper(
-            SpaceOperator(
-                Conversion(Ultraspherical(1//2),B),
-                A,B))
+        ConversionWrapper(SpaceOperator(Conversion(Ultraspherical(A), B), A, B))
     elseif A.a == A.b
         US = Ultraspherical(A)
-        ConversionWrapper(Conversion(US,B)*ConcreteConversion(A,US))
+        ConversionWrapper(SpaceOperator(
+            TimesOperator(Conversion(US,B), ConcreteConversion(A,US)), A, B))
     else
         J = Jacobi(B)
         ConcreteConversion(J,B)*Conversion(A,J)
     end
 end
 
-function Conversion(A::Ultraspherical,B::Jacobi)
+@inline function _Conversion(A::Ultraspherical, B::Jacobi)
     @assert domain(A) == domain(B)
     if isequalminhalf(B.a) && isequalminhalf(B.b)
-        ConversionWrapper(ConcreteConversion(Chebyshev(domain(A)),B)*
-            Conversion(A,Chebyshev(domain(A))))
+        C = Chebyshev(domain(B))
+        ConversionWrapper(SpaceOperator(
+            TimesOperator(ConcreteConversion(C, B), Conversion(A, C)), A, B))
     elseif isequalminhalf(B.a - order(A)) && isequalminhalf(B.b - order(A))
         ConcreteConversion(A,B)
     elseif B.a == B.b == 0
-        ConversionWrapper(
-            SpaceOperator(
-                Conversion(A,Ultraspherical(1//2,domain(B))),
-                A,B))
+        ConversionWrapper(SpaceOperator(Conversion(A, Ultraspherical(B)), A, B))
     elseif B.a == B.b
         US = Ultraspherical(B)
-        ConversionWrapper(ConcreteConversion(US,B)*Conversion(A,US))
+        ConversionWrapper(SpaceOperator(
+            TimesOperator(ConcreteConversion(US,B), Conversion(A,US)), A, B))
     else
         J = Jacobi(A)
         Conversion(J,B)*ConcreteConversion(A,J)
     end
 end
 
+@static if VERSION >= v"1.8"
+    Base.@constprop :aggressive Conversion(A::Jacobi, B::Ultraspherical) = _Conversion(A, B)
+    Base.@constprop :aggressive Conversion(A::Ultraspherical, B::Jacobi) = _Conversion(A, B)
+else
+    Conversion(A::Jacobi, B::Ultraspherical) = _Conversion(A, B)
+    Conversion(A::Ultraspherical, B::Jacobi) = _Conversion(A, B)
+end
 
 
 
@@ -420,10 +416,6 @@ function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,US,T},NTuple{2,UnitR
     ret
 end
 
-
-isapproxminhalf(a) = a ≈ -0.5
-isapproxminhalf(::Integer) = false
-
 function union_rule(A::Jacobi,B::Jacobi)
     if domainscompatible(A,B)
         Jacobi(min(A.b,B.b),min(A.a,B.a),domain(A))

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -133,9 +133,6 @@ end
 
 ## Conversion Operator
 
-isequalhalf(x) = x == 0.5
-isequalhalf(::Integer) = false
-
 function Conversion(A::Chebyshev, B::Ultraspherical)
     @assert domain(A) == domain(B)
     mB = order(B)
@@ -177,14 +174,16 @@ function Conversion(A::Ultraspherical,B::Ultraspherical)
         if -1 ≤ b-a ≤ 1 && (a,b) ≠ (2,1)
             return ConcreteConversion(A,B)
         elseif b-a > 1
-            r = b:-1:a+1
+            r = decreasingunitsteprange(b, a+1)
             v = [ConcreteConversion(Ultraspherical(i-1,d), Ultraspherical(i,d)) for i in r]
             if !(last(r) ≈ a+1)
                 vlast = ConcreteConversion(A, Ultraspherical(last(r)-1, d))
-                v = [v; vlast]
+                v2 = Union{eltype(v), typeof(vlast)}[v; vlast]
+            else
+                v2 = v
             end
             bwsum = isapproxinteger(b-a) ? (0, 2length(v)) : (0,ℵ₀)
-            return ConversionWrapper(TimesOperator(v, bwsum, (ℵ₀,ℵ₀), bwsum), A, B)
+            return ConversionWrapper(TimesOperator(v2, bwsum, (ℵ₀,ℵ₀), bwsum), A, B)
         end
     end
     throw(ArgumentError("please implement $A → $B"))

test/JacobiTest.jl

@@ -122,6 +122,12 @@ using Static
             @test Conversion(S1, S2) * Fun(x->x^4, S1) ≈ Fun(x->x^4, S2)
         end
 
+        C = Conversion(Jacobi(Chebyshev()), Ultraspherical(1))
+        @test C * Fun(x->x^2, Jacobi(-0.5, -0.5)) ≈ Fun(x->x^2, Ultraspherical(1))
+
+        C = Conversion(Ultraspherical(0.5), Jacobi(Chebyshev()))
+        @test C * Fun(x->x^2, Ultraspherical(0.5)) ≈ Fun(x->x^2, Jacobi(Chebyshev()))
+
         @testset "inference tests" begin
             #= Note all cases are inferred as of now,
             but as the situation eveolves in the future, more @inferred tests
@@ -165,6 +171,11 @@ using Static
 
                 CCJmix = Conversion(Chebyshev(), Jacobi(0.5,1.5))
                 @test CCJmix * Fun(Chebyshev()) ≈ Fun(Jacobi(0.5,1.5))
+
+                @inferred (P -> Conversion(P, Jacobi(1,1)))(Chebyshev());
+                if VERSION >= v"1.8"
+                    @inferred (P -> Conversion(Jacobi(1,1), P))(Ultraspherical(3));
+                end
             end
 
             @testset "Ultraspherical" begin
@@ -206,6 +217,12 @@ using Static
                 @test Legendre(Ultraspherical(0.5)) == Legendre()
                 @test_throws ArgumentError Legendre(Chebyshev())
             end
+
+            @test Jacobi(Ultraspherical(Jacobi(1,1))) === Jacobi(1,1)
+
+            @inferred (U -> Val(isinteger(Jacobi(U).a)))(Ultraspherical(half(Odd(1))))
+            @inferred (U -> Val(isinteger(Jacobi(U).a)))(Ultraspherical(1))
+            @inferred (U -> Val(isinteger(Jacobi(U).a)))(Ultraspherical(static(1)))
         end
 
         @test ApproxFunOrthogonalPolynomials.normalization(ComplexF64, Jacobi(-0.5, -0.5), 0) ≈ pi

test/UltrasphericalTest.jl

@@ -4,6 +4,8 @@ using Test
 using SpecialFunctions
 using LinearAlgebra
 using Static
+using OddEvenIntegers
+using HalfIntegers
 
 @verbose @testset "Ultraspherical" begin
     @testset "promotion" begin
@@ -40,13 +42,27 @@ using Static
             @inferred (() -> Conversion(Chebyshev(), Ultraspherical(2.5)))();
         end
 
+        Tallowed = Union{
+            typeof(Conversion(Chebyshev(), Ultraspherical(half(Odd(1))))),
+            typeof(Conversion(Chebyshev(), Ultraspherical(half(Odd(3)))))
+        }
+        @inferred Tallowed Conversion(Chebyshev(), Ultraspherical(half(Odd(1))));
+        @inferred Tallowed Conversion(Chebyshev(), Ultraspherical(half(Odd(3))));
+
         # Conversions between Ultraspherical should lead to a small union of types
         Tallowed = Union{
             typeof(Conversion(Ultraspherical(1), Ultraspherical(2))),
             typeof(Conversion(Ultraspherical(1), Ultraspherical(3)))}
         @inferred Tallowed Conversion(Ultraspherical(1), Ultraspherical(2));
         @inferred Tallowed Conversion(Ultraspherical(1), Ultraspherical(3));
 
+        Tallowed = Union{
+            typeof(Conversion(Ultraspherical(1), Ultraspherical(half(Odd(3))))),
+            typeof(Conversion(Ultraspherical(1), Ultraspherical(half(Odd(5)))))}
+
+        @inferred Tallowed Conversion(Ultraspherical(1), Ultraspherical(half(Odd(3))));
+        @inferred Tallowed Conversion(Ultraspherical(1), Ultraspherical(half(Odd(5))));
+
         @inferred Conversion(Ultraspherical(static(1)), Ultraspherical(static(4)));
         @inferred (() -> Conversion(Ultraspherical(1), Ultraspherical(4)))();
         @inferred (() -> Conversion(Ultraspherical(1.0), Ultraspherical(4.0)))();
@@ -63,6 +79,10 @@ using Static
             @test C1[1:4, 1:4] == C2[1:4, 1:4]
         end
 
+        C1 = Conversion(Chebyshev(), Ultraspherical(1.5))
+        C2 = Conversion(Chebyshev(), Ultraspherical(half(Odd(3))))
+        @test C1[1:4, 1:4] == C2[1:4, 1:4]
+
         f = Fun(x->x^2, Ultraspherical(0.5)) # Legendre
         CLC = Conversion(Ultraspherical(0.5), Chebyshev())
         @test !isdiag(CLC)