Trigonometric basis (#49)

* Trigonometric basis * Fix format * Add tests
JuliaAlgebra · Jul 13, 2024 · 37d9cf3 · 37d9cf3
1 parent e4737e2
commit 37d9cf3
Show file tree

Hide file tree

Showing 5 changed files with 143 additions and 14 deletions.
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -38,4 +38,5 @@ AbstractGegenbauer
 Legendre
 ChebyshevFirstKind
 ChebyshevSecondKind
+Trigonometric
 ```
diff --git a/src/MultivariateBases.jl b/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")
 

diff --git a/src/chebyshev.jl b/src/chebyshev.jl
@@ -24,11 +24,26 @@ 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 +63,13 @@ 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

diff --git a/src/trigonometric.jl b/src/trigonometric.jl
@@ -0,0 +1,100 @@
+"""
+    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 <: AbstractMultipleOrthogonal end
+
+_cos_id(d) = iszero(d) ? 0 : 2d - 1
+_sin_id(d) = 2d
+_is_cos(d) = isodd(d)
+_is_sin(d) = d > 0 && iseven(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
+
+# FIXME The cos part is, like Chebysev, maybe the sin part too ? We should do better here, this is just a stopgap
+even_odd_separated(::Type{Trigonometric}) = false
diff --git a/test/trigonometric.jl b/test/trigonometric.jl
@@ -0,0 +1,15 @@
+using Test
+using MultivariateBases
+using DynamicPolynomials
+
+@testset "StarAlgebras" begin
+    @polyvar x
+    a = MB.Polynomial{MB.Trigonometric}(x)
+    b = a * a
+    @test b.coeffs == MB.sparse_coefficients(1 // 2 + 1 // 2 * x^3)
+    c = b * b
+    @test c.coeffs ==
+          MB.sparse_coefficients(3 // 8 + 1 // 2 * x^3 + 1 // 8 * x^7)
+    @test a * MB.Polynomial{MB.Trigonometric}(constant_monomial(typeof(x))) ==
+          a * MB.Polynomial{MB.Trigonometric}(x^0)
+end