I have problems with the series that appear evaluating this integral $$\int_0^n \{x^2\}\,\text{d}x$$ where $\{\cdot\}$ is the fractional part function. Now, I'm pretty sure that $$I=\sum_{k=0}^{n^2-1}\int_\sqrt k^\sqrt{k+1} (x^2-k)\,\text{d}x$$ since $\{x^2\}=x^2-\lfloor x^2\rfloor=x^2-k$ for $x\in[\sqrt k,\sqrt{k+1})$. Now, if that is correct, I proceed integrating and I obtain the series $$I=\frac13\sum_{k=0}^{n^2-1} ((k+1)^\frac32-k^\frac32)-\sum_{k=0}^{n^2-1} k(\sqrt{k+1}-\sqrt k)$$ The first telescopes to $\frac{n^3}3$ while for the second I thought to apply the formula of summation by parts, that is $\sum_{n=k}^Na_nb_n=S_Nb_N-S_{k-1}b_k-\sum_{n=k}^{N-1}S_n(b_{n+1}-b_n)$ where $S_N=\sum_{n=1}^Na_n$. Choosing $a_n=\sqrt{k+1}-\sqrt k$, since its $(n^2-1)$-th partial sum telescopes to $n-1$, the formula gives $$\sum_{k=0}^{n^2-1} k(\sqrt{k+1}-\sqrt k)=(n-1)(n^2-1)-\sum_{k=0}^{n^2-2} (k-1)$$ that seems to be false because $\sum_{k=0}^{n^2-2} (k-1)>(n-1)(n^2-1)$ while $\sum_{k=0}^{n^2-1} k(\sqrt{k+1}-\sqrt k)>0$. I sincerely think I have done something wrong but I don't understand what… Could someone explain or find another solution (the answer given by the book is $-\frac{2n^3}3+\sum_{k=1}^{n^2} \sqrt k$)?
Evaluating $\int_0^n \{x^2\}\,\text{d}x$
definite integralsfractional-partintegrationsummation-by-partstelescopic-series
Related Solutions
If you are using IEEE 754 compliant floating point numbers, it may be that square root operations are faster than you might suppose, as the square root operation is directly supported in the standard and any compliant hardware is guaranteed to round correctly (unlike for sine and cosine) and is often using a hardware implementation of the Newton-Raphson method Alijah described.
That said, the algorithm you linked to only uses the fractional part of the square root to calculate pixel opacity, and consequently the final opacity value ranges only from 0 to 255. Because of the small range, floating point numbers may be overkill and a fixed-point integer representation might work better. If the range is truly only a byte and the maximum size of your radii aren't too large, you can use a look-up table with fixed-point input to skip the expensive square root calculation. A 16-bit fixed-point number would give you a 64KB look-up table, which isn't too bad.
You might also be able to avoid a division and square root operation for calculating the $45^\circ$ value (ffd
in the algorithm) by using the Fast Inverse Square Root hack.
Now for the question of whether there is a method to calculate the fractional part of a square root knowing only the integer part, there is one approach that iteratively calculates a square root and minimizes divisions: Continued Fractions. The simplest approach I know of for what you are wanting to do (fast convergence) is detailed here and works as follows:
$a_0 = x - \lfloor\sqrt x\rfloor^2$
$b_0 = 2\lfloor\sqrt x\rfloor$
$a_n = a_0b_{n-1}$
$b_n = b_0b_{n-1} + a_{n-1}$
which gives you quadratically better and better approximations of the fractional part of $\sqrt x$, and you divide $\frac{a_n}{b_n}$ when you've done enough iterations to ensure accuracy. If you are ultimately needing only byte sized opacity values, it should only take 3 or 4 iterations or so, and we save the division until the end, which is the only significant difference between this and the Newton-Raphson method other than the fact that it gives you the fractional part directly.
If you really want to pursue incrementally calculating variables as far as it can go, you can use Gosper's continued fraction algorithms (see especially the section on square roots) and calculate all the variables involved as continued fractions one term at a time. This allows you to avoid square root calculations and divisions other than bit shifts, as well as abort as soon as you know the correct pixel opacity (or whatever else you're calculating) without wasting time calculating digits you don't need, but it involves a serious overhaul to the algorithm you linked, and if I went into detail this answer would turn into a book.
So essentially, if you have memory to spare and the max length of your radii isn't huge and your output size is small, go with a look-up table with fixed-point numbers. If you want simplicity of implementation, go with floating point or fixed-point numbers. If you want to absolutely avoid square root calculation without a look-up table go with Newton-Raphson or the continued fraction variant on it. If you want to absolutely minimize wasted computation at the expense of some up-front overhead, go with Gosper's continued fraction algorithms.
Here's half a solution—an asymptotic formula for the terms in the sum with $a$ even. Hopefully the method can be carried through with minor modifications for $a$ odd.
We start with \begin{align*} \sum_{\substack{a = \lceil\sqrt{N + 2}\rceil \\ a\text{ even}}}^{\lfloor 2N/3 \rfloor} \biggl\{ \frac{{a}^{2}-N}{2a} \biggr\} &= \sum_{\substack{a=1 \\ a\text{ even}}}^{\lfloor 2N/3 \rfloor} \biggl\{ \frac{{a}^{2}-N}{2a} \biggr\} + O(\sqrt N) \\ &= \sum_{\substack{1\le a\le 2N/3 \\ a\text{ even}}} \biggl\{ \frac{-N}{2a} \biggr\} + O(\sqrt N) = \sum_{1\le b\le N/3} \biggl\{ \frac{-N}{b} \biggr\} + O(\sqrt N) \end{align*} after setting $b=2a$. Note that $\{-y\}=1-\{y\}$ unless $y$ is an integer, while the number of $b$ such that $\frac{-N}b$ is an integer is at most $\tau(N)$, the number of divisors of $N$, which is certainly $O(\sqrt N)$. Therefore \begin{align*} \sum_{b\le N/3} \biggl\{ \frac{-N}{b} \biggr\} &= \sum_{b\le N/3} \biggl( 1 - \biggl\{ \frac{N}{b} \biggr\} \biggr) + O(\sqrt N) \\ &= \frac N3 - \sum_{b\le N/3} \biggl( \frac{N}{b} - \biggl\lfloor \frac{N}{b} \biggr\rfloor \biggr) + O(\sqrt N) \\ &= \frac N3 - N\biggl( \log \frac N3 + \gamma + O\biggl( \frac1N \biggr) \biggr) + \sum_{b\le N/3} \biggl\lfloor \frac{N}{b} \biggr\rfloor + O(\sqrt N) \\ &= -N\log N + N(\tfrac13 + \log3 - \gamma) + \sum_{b\le N/3} \biggl\lfloor \frac{N}{b} \biggr\rfloor + O(\sqrt N), \end{align*} where $\gamma\approx 0.577216$ is Euler's constant. On the other hand, \begin{align*} \sum_{b\le N/3} \biggl\lfloor \frac{N}{b} \biggr\rfloor &= \sum_{b\le N/3} \sum_{\substack{d\le N \\ b\mid d}} 1 = \sum_{d\le N} \sum_{\substack{b\le N/3 \\ b\mid d}} 1 \\ &= \sum_{d\le N} \bigl( \tau(d) - [1\text{ if } d>\tfrac N3] - [1\text{ if } d>\tfrac{2N}3\text{ and $d$ is even}] \bigr) \\ &= \sum_{d\le N} \tau(d) - \frac{2N}3 - \frac N6 + O(1) \\ &= N\log N + (2\gamma-1)N - \frac{2N}3 - \frac N6 + O(\sqrt N) \\ &= N\log N + N\bigl( 2\gamma-\tfrac{11}6 \bigr) + O(\sqrt N). \end{align*} Consequently, \begin{align*} \sum_{\substack{a = \lceil\sqrt{N + 2}\rceil \\ a\text{ even}}}^{\lfloor 2N/3 \rfloor} \biggl\{ \frac{{a}^{2}-N}{2a} \biggr\} &= -N\log N + N(\tfrac13 + \log3 - \gamma) \\ &\qquad{}+ N\log N + N\bigl( 2\gamma-\tfrac{11}6 \bigr) + O(\sqrt N) \\ &= N(\log3 + \gamma - \tfrac32) + O(\sqrt N). \end{align*} Here the constant $\log3 + \gamma - \tfrac32 \approx 0.175828$ is borne out by numerical experiments.
Best Answer
Ok, I just understood my mistake and it's very silly. We have that $S_k=\sqrt{k+1}-1$ and instead I used $S_{k^2-1}$ as $S_k$. Then the answer comes easily expanding, re-elaborating and re-indexing the sums. So stupid...