From 37d9cf3872fc4be5d6835d8e0aab2135215ce570 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 13 Jul 2024 21:31:59 +0200 Subject: [PATCH] Trigonometric basis (#49) * Trigonometric basis * Fix format * Add tests --- docs/src/index.md | 1 + src/MultivariateBases.jl | 3 +- src/chebyshev.jl | 38 ++++++++++----- src/trigonometric.jl | 100 +++++++++++++++++++++++++++++++++++++++ test/trigonometric.jl | 15 ++++++ 5 files changed, 143 insertions(+), 14 deletions(-) create mode 100644 src/trigonometric.jl create mode 100644 test/trigonometric.jl diff --git a/docs/src/index.md b/docs/src/index.md index 9200511..bcc9504 100644 --- 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 index 8e35e58..bd23b67 100644 --- 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 index 35619c1..aeb35d4 100644 --- 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 new file mode 100644 index 0000000..7b6ce57 --- /dev/null +++ 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 new file mode 100644 index 0000000..1c1fda7 --- /dev/null +++ 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