-
Notifications
You must be signed in to change notification settings - Fork 5
/
trigonometric.jl
100 lines (93 loc) · 3.19 KB
/
trigonometric.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
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