tl;dr: Use the right data structure for the job. I implemented a Fenwick tree for Fortran90 with a couple neat features, available here.
Our fictitious story begins with an eccentric billionaire and Christmas aficionado who is holding a contest for the whole USA: they are giving away $1 billion over the span of 1 hour, all people have to do to compete is try and put as many ornaments on a Christmas tree as possible and submit their total to the contest website! If they manage to fit even more on, or if some fall off they need to update their total.
You are the intern in the team tasked by the billionaire to handle the drawing of winners. The rest of the team has already set up the website, database, etc. The team reserved this simplest of tasks for the intern, surely you can at least do what they ask!
Your boss tells you they expect around 100 million contestants, they will award $250,000 once a second for a full hour. In order to reward contestants with more ornaments on their tree with a better chance of winning they choose winners in the following way: there is an array with an entry for each contestant that contains their total, perform a cumulative sum over this array and pick a random number in the range from 1 to the total number of ornaments, whichever entry in the cumulative sum is closest to the random number is a winner! All you need to do is pick a winner once a second and update any contestants entries if they changed since the second before.
You think a bit about how they are drawing winners, is that fair? After some thought you agree; people with bigger entries are proportionally more likely to be chosen, while those with smaller entries are proportionately less likely. “This implementation should be dirt simple!” you say, and type up a prototype in Fortran90:
use binary_search
use contest_stuff !provides contest related functions by team
integer, dimension(num_contestants) :: contest_entry, cumulative_sum
real :: rand_real
integer :: i,t,rand_int
do t=1,3600 !once a second, for an hour
call wait_and_update(contest_entry) !pulls in from database etc.
cumulative_sum(1) = contest_entry(1)
do i=2, num_contestants
cumulative_sum(i) = cumulative_sum(i-1) + contest_entry(i)
end do
call random_number(rand_real)
rand_int = int(rand_real * cumulative_sum(num_contestants+1))
winner = binary_search(cumulative_sum, rand_int) !returns winner index
call pick_winner(winner) !sends winner out to database etc.
end doSlick! You even used a binary search over the cumulative sum (you noticed that it will be monotonically increasing since all entries are positive, so you can do this), so it’s cost will be 
You re-run and profile the code, and are surprised to see that the inner do loop takes 99% of running time! In retrospect this isn’t too surprising, after all the cost of the cumulative sum is 
do loop itself is an appreciable fraction of the total loop cost!). But what room for improvement is there in practically the simplest loop possible?
This is a silly story, but it proves a point. It’s clear that part of the problem is you are reading/writing a single array entry at a time; you could loop unroll and process a vector stride at a time but you’ll find that compiling with -O3 already does this for you. The loop unrolling gave a 3x speedup on my system, but that still leaves the asymptotic complexity untouched.
The Fenwick tree to the rescue! Also known as a binary indexed tree, it gives a fast way to compute cumulative sums and update them based on changes in the original parent array. With our example code we needed to compute the cumulative sum each second, with a cost scaling of 
The Fenwick tree is also quite compact to implement as well, partly because it never actually implements a tree and instead stores it implicitly in a flat array. Here’s some example code to build a tree from a parent array:
pure function build_fen(arr) result(fen)
integer, dimension(:), intent(in) :: arr
integer, dimension(size(arr)) :: fen
integer :: i,j,n
fen = 0
n = size(arr)
do i=1,n
j = i
do while (j<=n)
fen(j) = fen(j) + arr(i)
j = j + iand(j,-j)
end do
end do
end function build_fenNote: This function uses assumed-shape arrays, as such if you want to use it you either need to give it an explicit interface or stick it in a module. Check links at the top and bottom of this post to my F90 module that implements all this for you.
“But wait!” you say, aren’t there 
It isn’t usually mentioned, but you can actually search the tree in
There’s another refinement to implementation performance to be had as well. It is often the case that contiguous chunks of updates are made to the parent array. Naive updates would then have a O(M log N) cost for tree size N and M updates. Instead we can break these updates into two regions: inside the contiguous chunk and everything after. As a result we have reduced the cost to O(log N + M log M). For “large” M and N>>M this can translate to a significant savings, but even for M=2 or 3 it is still faster than 2 or 3 naive sequential updates.
When I came across this problem I wasn’t able to find a Fortran90 implementation, so I made one (along with the improvements to chunk updates and binary search I’ve mentioned). You can find it here on my Github!
work (where 
is some constant, in this example c~=300 FLOPs). Memory latency and bandwidth has not scaled at nearly the same rate as CPU throughput, and so in cases of weak data-reuse it tends to dominate the behavior. Compare this to say a naive n-body calculation where computation scales as 
despite still using just ![[a,b]](https://s0.wp.com/latex.php?latex=%5Ba%2Cb%5D&bg=ffffff&fg=000000&s=0 "[a,b]")
as:
space the inner product is essentially the dot product of the two vectors. Vectors that are orthogonal to each other have a dot product of 0. The concept of orthogonality can be abstracted and extended to vectors in any vector space, so that vectors that have a zero inner product are orthogonal.
on the domain ![I_k =[x_{k-^1\!/_2},x_{k+^1\!/_2}]](https://s0.wp.com/latex.php?latex=I_k+%3D%5Bx_%7Bk-%5E1%5C%21%2F_2%7D%2Cx_%7Bk%2B%5E1%5C%21%2F_2%7D%5D&bg=ffffff&fg=000000&s=0 "I_k =[x_{k-^1\!/_2},x_{k+^1\!/_2}]")
may not permit an analytical representation and so most be approximated in some way. An easy choice is a polynomial approximation of order M. A natural choice for the polynomial basis is the monomials, 
. These are satisfactory for low orders, but we encounter a problem at higher orders.
norm 
we arrive at 
. If the inner product of two basis functions is 0, they are orthogonal and therefore linearly independent. If many of our basis functions end up linearly dependent it is difficult to reconstruct the solution approximation.
the inner product is close to 1, indicating near linear dependence between basis functions.
which are orthogonal on the interval ![[-1,1]](https://s0.wp.com/latex.php?latex=%5B-1%2C1%5D+&bg=ffffff&fg=000000&s=0 "[-1,1]")
(these are defined as the solutions to the 
). They can be defined using a recurrence formula or the explicit Rodrigue’s formula.
![P_m(x) = \frac{1}{2^m m!}\frac{\mathrm{d}^2}{\mathrm{d}x^m}\left[(x^2-1)^m\right]](https://s0.wp.com/latex.php?latex=P_m%28x%29+%3D+%5Cfrac%7B1%7D%7B2%5Em+m%21%7D%5Cfrac%7B%5Cmathrm%7Bd%7D%5E2%7D%7B%5Cmathrm%7Bd%7Dx%5Em%7D%5Cleft%5B%28x%5E2-1%29%5Em%5Cright%5D+&bg=ffffff&fg=000000&s=0 "P_m(x) = \frac{1}{2^m m!}\frac{\mathrm{d}^2}{\mathrm{d}x^m}\left[(x^2-1)^m\right]")
and 
. In the orthonormal case we would need to explicitly evaluate 
at the endpoints for every 
(it is important to be clear with what one means by normalization in this case: orthonormal with respect to the 
norm or normalized with respect to the maximal value of 
on 
to transform 
to 
. The newly defined shifted Legendre polynomials are now orthogonal over 
and the inner product becomes
for arbitrary order M using basis functions 
and basis weights 
order spatial modal ![\sum_{m=0}^M \left[ \frac{\mathrm{d}\overset{m}{a}}{\mathrm{d} t} \int_{I_k}\! \overset{m}{\psi} \, \overset{n}{\phi} \,\mathrm{d}x \right] + [g(v^-\!,v^+) \; \overset{n}{\phi}] \Big\rvert_{x_{k-^1\!/_2}}^{x_{k+^1\!/_2}} \!-\! \int_{I_k}\! f(v) \,\frac{\mathrm{d} \overset{n}{\phi}}{\mathrm{d} x} \,\mathrm{d}x = 0\;\;for \; 0 \leq n \leq N](https://s0.wp.com/latex.php?latex=%5Csum_%7Bm%3D0%7D%5EM+%5Cleft%5B+%5Cfrac%7B%5Cmathrm%7Bd%7D%5Coverset%7Bm%7D%7Ba%7D%7D%7B%5Cmathrm%7Bd%7D+t%7D+%5Cint_%7BI_k%7D%5C%21+%5Coverset%7Bm%7D%7B%5Cpsi%7D+%5C%2C+%5Coverset%7Bn%7D%7B%5Cphi%7D+%5C%2C%5Cmathrm%7Bd%7Dx+%5Cright%5D+%2B+%5Bg%28v%5E-%5C%21%2Cv%5E%2B%29+%5C%3B+%5Coverset%7Bn%7D%7B%5Cphi%7D%5D+%5CBig%5Crvert_%7Bx_%7Bk-%5E1%5C%21%2F_2%7D%7D%5E%7Bx_%7Bk%2B%5E1%5C%21%2F_2%7D%7D+%5C%21-%5C%21+%5Cint_%7BI_k%7D%5C%21+f%28v%29+%5C%2C%5Cfrac%7B%5Cmathrm%7Bd%7D+%5Coverset%7Bn%7D%7B%5Cphi%7D%7D%7B%5Cmathrm%7Bd%7D+x%7D+%5C%2C%5Cmathrm%7Bd%7Dx+%3D+0%5C%3B%5C%3Bfor+%5C%3B+0+%5Cleq+n+%5Cleq+N+&bg=ffffff&fg=000000&s=0 "\sum_{m=0}^M \left[ \frac{\mathrm{d}\overset{m}{a}}{\mathrm{d} t} \int_{I_k}\! \overset{m}{\psi} \, \overset{n}{\phi} \,\mathrm{d}x \right] + [g(v^-\!,v^+) \; \overset{n}{\phi}] \Big\rvert_{x_{k-^1\!/_2}}^{x_{k+^1\!/_2}} \!-\! \int_{I_k}\! f(v) \,\frac{\mathrm{d} \overset{n}{\phi}}{\mathrm{d} x} \,\mathrm{d}x = 0\;\;for \; 0 \leq n \leq N")
![\frac{\mathrm{d}\overset{n}{a}}{\mathrm{d} t}\left(\frac{\Delta x}{2n+1}\right)= \!\! \sum_{\underset{by \: 2}{m=n-1}}^0 \overset{m}{a}_k-\sum_{m=0}^M \left[ \overset{m}{a}_k - \overset{m}{a}_{k-1}(-1)^m \right]\;\;\; for \; 0 \leq n \leq N](https://s0.wp.com/latex.php?latex=%5Cfrac%7B%5Cmathrm%7Bd%7D%5Coverset%7Bn%7D%7Ba%7D%7D%7B%5Cmathrm%7Bd%7D+t%7D%5Cleft%28%5Cfrac%7B%5CDelta+x%7D%7B2n%2B1%7D%5Cright%29%3D+%5C%21%5C%21+%5Csum_%7B%5Cunderset%7Bby+%5C%3A+2%7D%7Bm%3Dn-1%7D%7D%5E0+%5Coverset%7Bm%7D%7Ba%7D_k-%5Csum_%7Bm%3D0%7D%5EM+%5Cleft%5B+%5Coverset%7Bm%7D%7Ba%7D_k+-+%5Coverset%7Bm%7D%7Ba%7D_%7Bk-1%7D%28-1%29%5Em+%5Cright%5D%5C%3B%5C%3B%5C%3B+for+%5C%3B+0+%5Cleq+n+%5Cleq+N+&bg=ffffff&fg=000000&s=0 "\frac{\mathrm{d}\overset{n}{a}}{\mathrm{d} t}\left(\frac{\Delta x}{2n+1}\right)= \!\! \sum_{\underset{by \: 2}{m=n-1}}^0 \overset{m}{a}_k-\sum_{m=0}^M \left[ \overset{m}{a}_k - \overset{m}{a}_{k-1}(-1)^m \right]\;\;\; for \; 0 \leq n \leq N")
where 
are constants. Determine the constants so that these functions are mutually orthogonal in (-1,1) and thus obtain the functions.