-
Notifications
You must be signed in to change notification settings - Fork 15
/
expanda.jl
267 lines (242 loc) · 6.84 KB
/
expanda.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
#### ExpandA
function _doexpand(x)
! isa(x,Mxpr) && return x
n = length(x)
args = margs(x)
nargs = newargs(n)
@inbounds for i in 1:n
nargs[i] = doexpand(args[i])
end
return doexpand(mxpr(mhead(x),nargs))
end
doexpand(p::Mxpr{:Power}) = do_expand_power(p,base(p),expt(p))
function doexpand(prod::Mxpr{:Times})
a = margs(prod)
len = length(a)
@inbounds for i in 1:len
a[i] = doexpand(a[i])
end
have_sum = false
j = 0
@inbounds for i in 1:len # check if we have anything to do
j += 1
if is_Mxpr(a[i],:Plus)
have_sum = true
break
end
end
! have_sum && return prod
nonsums = newargs()
@inbounds for i in 1:j-1 # none of these are sums
push!(nonsums,a[i])
end
sums = newargs()
push!(sums,a[j]) # already know its a sum
@inbounds for i in j+1:len # push more sums, if there are any
if is_Mxpr(a[i],:Plus)
push!(sums,a[i])
else
push!(nonsums,a[i])
end
end
sumres = length(sums) == 1 ? sums[1] : mulfacs(sums...)
nlen = length(nonsums)
if nlen == 0
mxout = sumres
elseif nlen == 1
mxout = mulfacs(nonsums[1],sumres)
else
mxout = mulfacs(mxpr(:Times,nonsums),sumres)
end
setfixed(mxout)
return mxout
end
doexpand(mx) = mx
do_expand_power(p,b::Mxpr{:Plus}, n::Integer) =
length(b) != 2 ? p : do_expand_binomial(p,b[1],b[2],n)
do_expand_power(p,b,ex) = p
do_expand_binomial(p,a, b, n::Integer) = expand_binomial(a,b,n)
do_expand_binomial(a,b,n) = p
function doexpand(s::Mxpr{:Plus})
args = margs(s)
@inbounds for i in 1:length(args)
args[i] = doexpand(args[i])
end
return s
end
## expand product of two sums
# a an b are the factors
# In test cases, this is fast. later, canonicalizing each term is slowest thing.
function mulfacs(a::Mxpr{:Plus},b::Mxpr{:Plus})
terms = newargs(length(a)*length(b))
i = 0
for ax in a.args
for bx in b.args
i += 1
t = flatcanon!(mxpr(:Times, ax, bx)) # TODO specialize for types of ax, bx
mergesyms(t,ax)
mergesyms(t,bx)
setfixed(t)
@inbounds terms[i] = t
end
end
mx = flatcanon!(mxpr(:Plus,terms))
for t in terms
mergesyms(mx,t)
end
setfixed(mx)
mx
end
# Not the right way to do this. We need to expand each term, as well.
mulfacs(a,b,c) = mulfacs(mulfacs(a,b),c)
mulfacs(a,b,c,xs...) = mulfacs(mulfacs(mulfacs(a,b),c),xs...)
# Should probably write this
#function mulfacs(a::Mxpr{:Plus}, b::SJSym)
#end
function _pushfacs!(facs,mx::Mxpr{:Times})
append!(facs,margs(mx))
end
function _pushfacs!(facs,b)
push!(facs,b)
end
function mulfacs(a::Mxpr{:Plus},b)
terms = newargs(length(a))
i = 0
for ax in a.args
i += 1
# if is_Mxpr(ax) && mxprtype(ax) == :Times
if isa(ax, Mxpr{:Times})
facs = copy(margs(ax))
_pushfacs!(facs,b) # I hope type inference optimizes this.
terms[i] = mxpr(:Times,facs)
else
terms[i] = doeval(mmul(b, ax))
end
canonexpr!(terms[i])
mergeargs(terms[i])
end
mx = mxpr(:Plus,terms)
mergeargs(mx)
mx
end
function mulfacs(a, b::Mxpr{:Plus})
mulfacs(b,a)
end
## construct Power. Decide whether to canonicalize according to types of args
function canonpower(base,expt::Real)
canonexpr!(mpow(base,expt))
end
function canonpower(base::SJSym,expt::Real)
canonexpr!(mpow(base,expt))
end
function canonpower(base::SJSym,expt)
mpow(base,expt)
end
function canonpower(base,expt)
mpow(base,expt)
end
##
# optimize a bit for types
function _expand_binomial_aux1(a::SJSym,b::SJSym,n)
mxpr(:Times,n,canonpower(a,(n-1)),b)
end
function _expand_binomial_aux1(a,b,n)
flatcanon!(mxpr(:Times,n,canonpower(a,(n-1)),b))
end
function _expand_binomial_aux2(a::SJSym,b::SJSym,n)
mxpr(:Times,n,a,canonpower(b,(n-1)))
end
function _expand_binomial_aux2(a,b,n)
flatcanon!(mxpr(:Times,n,a,canonpower(b,(n-1))))
end
# Be careful to construct expression in canonical form.
# Lots of ways to go wrong.
# Assume a < b in canonical order
# This is the only place we are testing meta data in Mxpr giving which
# symbols it depends on.
function expand_binomial(a,b,n::T) where T<:Integer
args = newargs(n+1)
args[1] = canonpower(a,n)
mergesyms(args[1],a)
setfixed(args[1])
args[n+1] = canonpower(b,n)
mergesyms(args[n+1],b)
setfixed(args[n+1])
if n == 2
args[2] = flatcanon!(mxpr(:Times,2,a,b)) # we have to flatcanon
else
# TODO optimize for symbols a,b, as in general case below. No flatcanon.
args[2] = _expand_binomial_aux1(a,b,n)
args[n] = _expand_binomial_aux2(a,b,n)
mergesyms(args[2],a)
mergesyms(args[2],b)
mergesyms(args[n],a)
mergesyms(args[n],b)
setfixed(args[2])
setfixed(args[n])
fac = n
k = n
l = one(n)
expand_binomial_aux(k,l,n,fac,a,b,args)
end
mx = mxprcf(:Plus,args)
mergesyms(mx,a)
mergesyms(mx,b)
apply_upvalues_to_args!(mx) # takes some time
setage(mx)
setfixed(mx)
mx
end
# Big increase in efficiency (> 10x) for both these types
# typealias ExpNoCanon Union{SJSym,Number} move to symataconstants.jl
# Expand((a+b*c)^n) is 10x slower than Expand((a+b)^n)
function _expand_mulpowers(fac,b1,e1,b2,e2)
m1 = canonpower(b1,e1) # adds 10-15% time
m2 = canonpower(b2,e2)
setfixed(m1)
setfixed(m2)
mergesyms(m1,b1)
mergesyms(m2,b2)
# return flatcanon!(mxpr(:Times, fac, m1, m2)) # flatcanon adds 10x time !, even if nothing is done
# return flatcanon!(mmul(fac,mmul(m1,m2)))
return doeval(mmul(fac,mmul(m1,m2))) # flatcannon is enough because of nested call.
end
function _expand_mulpowers(fac,b1::T,e1,b2::V,e2) where {T<:ExpNoCanon, V<:ExpNoCanon}
m1 = b1^e1
m2 = b2^e2
setfixed(m1)
setfixed(m2)
mergesyms(m1,b1)
mergesyms(m2,b2)
return mxpr(:Times, fac, m1, m2)
end
function _expand_mulpowers(fac,b1::T,e1,b2,e2) where T<:ExpNoCanon
m1 = b1^e1
m2 = canonpower(b2,e2)
setfixed(m1)
setfixed(m2)
mergesyms(m1,b1)
mergesyms(m2,b2)
return flatcanon!(mxpr(:Times, fac, m1, m2))
end
function _expand_mulpowers(fac,b1,e1,b2::T,e2) where T<:ExpNoCanon
m1 = canonpower(b1,e1)
m2 = b2^e2
setfixed(m1)
setfixed(m2)
mergesyms(m1,b1)
mergesyms(m2,b2)
return flatcanon!(mxpr(:Times, fac, m1, m2))
end
function expand_binomial_aux(k,l,n,fac,a,b,args)
@inbounds for j in 2:n-2
k = k - 1
l = l + 1
fac = mmul(fac,k)
fac = div(fac,l)
args[j+1] = _expand_mulpowers(fac,a,n-j,b,j)
mergesyms(args[j+1],a)
mergesyms(args[j+1],b)
setfixed(args[j+1])
end
end