-
Notifications
You must be signed in to change notification settings - Fork 38
/
FormalPowerSeries.jl
314 lines (270 loc) · 8.45 KB
/
FormalPowerSeries.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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
# Excursions in formal power series (fps)
#
# Jiahao Chen <[email protected]> 2013
#
# The primary reference is
# [H] Peter Henrici, "Applied and Computational Complex Analysis", Volume I:
# Power Series---Integration---Conformal Mapping---Location of Zeros,
# Wiley-Interscience: New York, 1974
import Base: zero, one, inv, length, ==, +, -, *, ^
import Base.Broadcast: broadcasted
export FormalPowerSeries, fps, tovector, trim, isunit, isnonunit,
MatrixForm, reciprocal, derivative, isconstant, compose,
isalmostunit, FormalLaurentSeries
#############################
# Formal power series (fps) #
#############################
# Definition [H, p.9] - we limit this to finitely many coefficients
mutable struct FormalPowerSeries{Coefficients}
c :: Dict{BigInt, Coefficients}
FormalPowerSeries{Coefficients}(c::Dict{BigInt, Coefficients}) where Coefficients = new(c)
function FormalPowerSeries{Coefficients}(v::Vector{Coefficients}) where Coefficients
c=Dict{BigInt, Coefficients}()
for i in eachindex(v)
if v[i] != 0
c[i-1] = v[i] #Note this off by 1 - allows constant term c[0] to be set
end
end
FormalPowerSeries{Coefficients}(c)
end
end
FormalPowerSeries(v::Vector{T}) where T = FormalPowerSeries{T}(v)
#Convenient abbreviation for floats
fps = FormalPowerSeries{Float64}
zero(::Type{FormalPowerSeries{T}}) where {T} = FormalPowerSeries(T[])
one(::Type{FormalPowerSeries{T}}) where {T} = FormalPowerSeries(T[1])
zero(::FormalPowerSeries{T}) where {T} = zero(FormalPowerSeries{T})
one(::FormalPowerSeries{T}) where {T} = one(FormalPowerSeries{T})
#Return truncated vector with c[i] = P[n[i]]
function tovector(P::FormalPowerSeries{T}, n::AbstractVector{Index}) where {T,Index<:Integer}
nn = length(n)
c = zeros(nn)
for (k,v) in P.c, i in eachindex(n)
n[i]==k && (c[i]=v)
end
c
end
tovector(P::FormalPowerSeries, n::Integer)=tovector(P,1:n)
# Basic housekeeping and properties
# Remove extraneous zeros
function trim(P::FormalPowerSeries{T}) where {T}
for (k,v) in P.c
if v==0
delete!(P.c, k)
end
end
return P
end
length(P::FormalPowerSeries{T}) where {T} = maximum([k for (k,v) in P.c])
function ==(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) where {T}
for (k,v) in P.c
if v==0 #ignore explicit zeros
continue
elseif !haskey(Q.c, k)
return false
elseif Q.c[k] != v
return false
end
end
for (k,v) in Q.c
if v==0 #ignore explicit zeros
continue
elseif !haskey(P.c, k)
return false
elseif P.c[k] != v
return false
end
end
return true
end
# Basic arithmetic [H, p.10]
function +(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) where {T}
c = Dict{BigInt, T}()
for (k,v) in P.c
haskey(c,k) ? (c[k]+=v) : (c[k]=v)
end
for (k,v) in Q.c
haskey(c,k) ? (c[k]+=v) : (c[k]=v)
end
FormalPowerSeries{T}(c)
end
function -(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) where {T}
c = Dict{BigInt, T}()
for (k,v) in P.c
c[k] = get(c,k,zero(T)) + v
end
for (k,v) in Q.c
c[k] = get(c,k,zero(T)) - v
end
FormalPowerSeries{T}(c)
end
#negation
function -(P::FormalPowerSeries{T}) where {T}
c = Dict{BigInt, T}()
for (k,v) in P.c
c[k] = -v
end
FormalPowerSeries{T}(c)
end
#multiplication by scalar
function *(P::FormalPowerSeries{T}, n::Number) where {T}
c = Dict{BigInt, T}()
for (k,v) in P.c
c[k] = n*v
end
FormalPowerSeries{T}(c)
end
*(n::Number, P::FormalPowerSeries{T}) where {T} = *(P, n)
#Cauchy product [H, p.10]
*(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) where {T} = CauchyProduct(P, Q)
function CauchyProduct(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) where {T}
c = Dict{BigInt, T}()
for (k1, v1) in P.c, (k2, v2) in Q.c
c[k1+k2]=get(c, k1+k2, zero(T))+v1*v2
end
FormalPowerSeries{T}(c)
end
#Hadamard product [H, p.10] - the elementwise product
broadcasted(::typeof(*), P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) where T =
HadamardProduct(P, Q)
function HadamardProduct(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) where {T}
c = Dict{BigInt, T}()
for (k,v) in P.c
if v!=0 && haskey(Q.c,k) && Q.c[k] !=0
c[k] = v * Q.c[k]
end
end
FormalPowerSeries{T}(c)
end
#The identity element over the complex numbers
#replacement for previous function eye(P::FormalPowerSeries{T})
function FormalPowerSeries{T}(s::UniformScaling) where {T}
c = Dict{BigInt, T}()
c[0] = 1
return FormalPowerSeries{T}(c)
end
isunit(P::FormalPowerSeries{T}) where {T <: Number} = P==FormalPowerSeries{Float64}(I)
# [H, p.12]
isnonunit(P::FormalPowerSeries{T}) where {T} = (!haskey(P.c, 0) || P.c[0]==0) && !isunit(P)
#Constructs the top left m x m block of the (infinite) semicirculant matrix
#associated with the fps [H, Sec.1.3, p.14]
#[H] calls it the semicirculant, but in contemporary nomenclature this is an
#upper triangular Toeplitz matrix
#This constructs the dense matrix - Toeplitz matrices don't exist in Julia yet
function MatrixForm(P::FormalPowerSeries{T}, m :: Integer) where {T}
m<0 && error("Invalid matrix dimension $m requested")
M=zeros(T, m, m)
for (k,v) in P.c
if k < m
for i=1:m-k
M[i, i+k]=v
end
end
end
M
end
#Reciprocal of power series [H, Sec.1.3, p.15]
#
#This algorithm is NOT that in [H], since it uses determinantal inversion
#formulae for the inverse. It is much quicker to find the inverse of the
#associated infinite triangular Toeplitz matrix!
#
#Every fps is isomorphic to such a triangular Toeplitz matrix [H. Sec.1.3, p.15]
#
#As most reciprocals of finite series are actually infinite,
#allow the inverse to have a finite truncation
#
#TODO implement the FFT-based algorithm of one of the following
# doi:10.1109/TAC.1984.1103499
# https://cs.uwaterloo.ca/~glabahn/Papers/tamir-george-1.pdf
function reciprocal(P::FormalPowerSeries, n::Integer)
n<0 && error(sprintf("Invalid inverse truncation length %d requested", n))
a = tovector(P, 0:n-1) #Extract original coefficients in vector
a[1]==0 && error("Inverse does not exist")
inv_a1 = inv(a[1])
b = zeros(n) #Inverse
b[1] = inv_a1
for k=2:n
b[k] = -inv_a1 * sum([b[k+1-i] * a[i] for i=2:k])
end
FormalPowerSeries(b)
end
#Derivative of fps [H. Sec.1.4, p.18]
function derivative(P::FormalPowerSeries{T}) where T
c = Dict{BigInt, T}()
for (k,v) in P.c
if k != 0 && v != 0
c[k-1] = k*v
end
end
FormalPowerSeries{T}(c)
end
#[H, Sec.1.4, p.18]
function isconstant(P::FormalPowerSeries)
for (k,v) in P.c
if k!=0 && v!=0
return false
end
end
return true
end
# Power of fps [H, Sec.1.5, p.35]
function ^(P::FormalPowerSeries{T}, n::Integer) where T
if n == 0
return FormalPowerSeries{T}([one(T)])
elseif n > 1
#The straightforward recursive way
return P^(n-1) * P
#TODO implement the non-recursive formula
elseif n==1
return P
else
error("I don't know what to do for power n = $n")
end
end
# Composition of fps [H, Sec.1.5, p.35]
# This is the quick and dirty version
function compose(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) where {T}
sum([v * Q^k for (k, v) in P.c])
end
#[H, p.45]
function isalmostunit(P::FormalPowerSeries{T}) where {T}
(haskey(P.c, 0) && P.c[0]!=0) ? (return false) :
(haskey(P.c, 1) && P.c[1]!=0) ? (return true) : (return false)
end
# Reversion of a series (inverse with respect to composition)
# P^[-1]
# [H. Sec.1.7, p.47, but more succinctly stated on p.55]
# Constructs the upper left nxn subblock of the matrix representation
# and inverts it
function reversion(P::FormalPowerSeries{T}, n :: Integer) where {T}
n>0 ? error("Need non-negative dimension") : nothing
Q = P
A = zeros(n, n)
#Construct the matrix representation (1.9-1), p.55
for i=1:n
Q = P
ai = tovector(Q, n) #Extract coefficients P[1]...P[n]
A[i,:] = ai
i<n && (Q *= P)
end
#TODO I just need the first row of the inverse
B = inv(A)
FormalPowerSeries{T}(B[1, :]'[:,1])
end
inv(P::FormalPowerSeries, n::Integer) = reversion(P, n)
###############################
# Formal Laurent series (fLs) #
# [H., Sec. 1.8, p. 51] #
###############################
mutable struct FormalLaurentSeries{Coefficients}
c :: Dict{BigInt, Coefficients}
FormalLaurentSeries{Coefficients}(c::Dict{BigInt, Coefficients}) where Coefficients = new(c)
function FormalLaurentSeries{Coefficients}(v::Vector{Coefficients}) where Coefficients
c = new(c)
for i=1:length(v)
c[i] = v[i]
end
c
end
end