On Wednesday, July 16, 2014 11:28:12 AM UTC-4, Dahua Lin wrote:
> sumabs2(x)
>
That's cheating! Neal, I think what you wrote will work just fine, and
sumabs2 will work even better. It's a little tough to see what all is
being done in its implementation because the base library is trying to eek
every last ounce of performance out of the code, and do so generically for
many operations. Some of the optimizations may include:
1. Using an explicit loop. Unlike other languages (Matlab/Python/R),
Julia's for loops are faster than its vectorized alternatives. That's
because it needs to allocate an entire copy of an array for every
operation. Something like this will get you most of the way there:
s = zero(T) # Use a real zero of the same type as the complex number so it
doesn't change type
for xi in x
s += real(xi)*real(xi)+imag(xi)*imag(xi)
end
2. But for arrays, we can do better. Iteratively summing like that has
major floating point issues; if you're summing lots of small numbers, `s`
slowly becomes so large that the rest of the small numbers get lost within
the precision. The base implementation of sum for arrays is smarter — it
sums things "pairwise." Instead of summing everything into one large pot,
it sums things into small pots, then sums those together, and so on. This
requires an array as not all iterables support random access.
I think those two things will get you most of the way to being competitive
with the base algorithm.
