Trigonometric basis

src/MultivariateBases.jl

@@ -47,7 +47,7 @@ include("scaled.jl")
 export AbstractMultipleOrthogonal,
     ProbabilistsHermite, PhysicistsHermite, Laguerre
 export AbstractGegenbauer,
-    Legendre, Chebyshev, ChebyshevFirstKind, ChebyshevSecondKind
+    Legendre, Chebyshev, ChebyshevFirstKind, ChebyshevSecondKind, Trigonometric
 export algebra_element,
     sparse_coefficients,
     univariate_orthogonal_basis,
@@ -62,6 +62,7 @@ include("hermite.jl")
 include("laguerre.jl")
 include("legendre.jl")
 include("chebyshev.jl")
+include("trigonometric.jl")
 include("lagrange.jl")
 include("quotient.jl")
 

src/chebyshev.jl

@@ -24,11 +24,23 @@ const Chebyshev = ChebyshevFirstKind
 
 # https://en.wikipedia.org/wiki/Chebyshev_polynomials#Properties
 # `T_n * T_m = T_{n + m} / 2 + T_{|n - m|} / 2`
-function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial)
+function univariate_mul!(::Mul{Chebyshev}, terms, var, a, b)
+    I = eachindex(terms)
+    for i in I
+        mono = MP.monomial(terms[i]) * var^(a + b)
+        terms[i] = MA.mul!!(terms[i], var^abs(a - b))
+        terms[i] = MA.operate!!(/, terms[i], 2)
+        α = MA.copy_if_mutable(MP.coefficient(terms[i]))
+        push!(terms, MP.term(α, mono))
+    end
+    return
+end
+
+function (mul::Mul{B})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) where {B<:AbstractMonomialIndexed}
     terms = [MP.term(1 // 1, MP.constant_monomial(a * b))]
-    vars_a = MP.variables(a)
+    vars_a = MP.effective_variables(a)
     var_state_a = iterate(vars_a)
-    vars_b = MP.variables(b)
+    vars_b = MP.effective_variables(b)
     var_state_b = iterate(vars_b)
     while !isnothing(var_state_a) || !isnothing(var_state_b)
         if isnothing(var_state_a) ||
@@ -48,16 +60,7 @@ function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial)
         else
             var_a, state_a = var_state_a
             var_b, state_b = var_state_b
-            d_a = MP.degree(a, var_a)
-            d_b = MP.degree(b, var_b)
-            I = eachindex(terms)
-            for i in I
-                mono = MP.monomial(terms[i]) * var_a^(d_a + d_b)
-                terms[i] = MA.mul!!(terms[i], var_a^abs(d_a - d_b))
-                terms[i] = MA.operate!!(/, terms[i], 2)
-                α = MA.copy_if_mutable(MP.coefficient(terms[i]))
-                push!(terms, MP.term(α, mono))
-            end
+            univariate_mul!(mul, terms, var_a, MP.degree(a, var_a), MP.degree(b, var_b))
             var_state_a = iterate(vars_a, state_a)
             var_state_b = iterate(vars_b, state_b)
         end

src/trigonometric.jl

@@ -0,0 +1,96 @@
+"""
+    struct Trigonometric <: AbstractMonomialIndexed end
+
+Univariate trigonometric basis is
+```
+a0 + a1 cos(ωt) + a2 sin(ωt) + a3 cos(2ωt) + a4 sin(2ωt)
+```
+"""
+struct Trigonometric <: AbstractMonomialIndexed end
+
+_cos_id(d) = iszero(d) ? 0 : 2d - 1
+_sin_id(d) = 2d
+_is_cos(d) = isodd(d)
+_is_sin(d) = isodd(d)
+_id(d) = div(d + 1, 2)
+
+# https://en.wikipedia.org/wiki/Chebyshev_polynomials#Properties
+# Using
+# sin(a + b) = sin(a) cos(b) + cos(a) sin(b)
+# sin(a - b) = sin(a) cos(b) - cos(a) sin(b)
+# If a > b
+# sin(a) cos(b) = sin(a + b) + sin(a - b)
+# sin(a) cos(b) = sin(a + b) + sin(a - b)
+function univariate_mul!(::Mul{Trigonometric}, terms, var, a, b)
+    @assert !iszero(a)
+    @assert !iszero(b)
+    I = eachindex(terms)
+    da = _id(a)
+    db = _id(b)
+    for i in I
+        if _is_cos(a) == _is_cos(b)
+            # Chebyshev first kind
+            mono = MP.monomial(terms[i]) * var^(_cos_id(da + db))
+            terms[i] = MA.mul!!(terms[i], var^_cos_id(abs(da - db)))
+            terms[i] = MA.operate!!(/, terms[i], 2)
+            α = MA.copy_if_mutable(MP.coefficient(terms[i]))
+            push!(terms, MP.term(α, mono))
+            # cos(a + b) = cos(a) cos(b) - sin(a) sin(b)
+            # cos(a - b) = cos(a) cos(b) + sin(a) sin(b)
+            # cos(a - b) - cos(a + b)
+            if _is_sin(a)
+                terms[end] = MA.operate!!(*, terms[end], -1)
+            end
+        else
+            if _is_cos(a)
+                da, db = db, da
+            end
+            # sin(da) * cos(db)
+            if da == db
+                # sin(da) * cos(da) = sin(2da) / 2
+                terms[i] = MA.mul!!(terms[i], var^_cos_id(da + db))
+                terms[i] = MA.operate!!(/, terms[i], 2)
+            else
+                # Using
+                # sin(a + b) = sin(a) cos(b) + cos(a) sin(b)
+                # sin(a - b) = sin(a) cos(b) - cos(a) sin(b)
+                # If a > b
+                # sin(a) cos(b) = (sin(a + b) + sin(a - b)) / 2
+                # If a < b
+                # sin(a) cos(b) = (sin(b + a) - sin(b - a)) / 2
+                mono = MP.monomial(terms[i]) * var^(_sin_id(da + db))
+                terms[i] = MA.mul!!(terms[i], var^_sin_id(abs(da - db)))
+                terms[i] = MA.operate!!(/, terms[i], 2)
+                α = MA.copy_if_mutable(MP.coefficient(terms[i]))
+                push!(terms, MP.term(α, mono))
+                if da < db
+                    terms[i] = MA.operate!!(*, terms[i], -1)
+                end
+            end
+        end
+    end
+    return
+end
+
+function degree_one_univariate_polynomial(::Type{Trigonometric}, variable)
+    MA.@rewrite(variable + 0)
+end
+
+function recurrence_eval(
+    ::Type{Trigonometric},
+    previous::AbstractVector,
+    value,
+    degree,
+)
+    d = _id(degree)
+    if _is_cos(degree)
+        # Chebyshev first order
+        return 2 * value * previous[_cos_id(d - 1) + 1] - previous[_cos_id(d - 2) + 1]
+    else
+        return sqrt(1 - previous[degree]^2)
+    end
+end
+
+function _promote_coef(::Type{T}, ::Type{Trigonometric}) where {T}
+    return _promote_div(T)
+end