Problem 235: An Arithmetic Geometric Sequence

Detailed mathematical approach

Problem Summary

The sequence \(u(k)=(900-3k)\,r^{k-1}\) multiplies the arithmetic progression \(897, 894, 891, \dots\), which falls by \(3\) at every step, by the geometric progression \(1, r, r^2, \dots\); such a product is called an arithmetic-geometric sequence. With \(s(n)=u(1)+u(2)+\dots+u(n)\), we must find the value of \(r\) for which \(s(5000)=-600\,000\,000\,000\), rounded to \(12\) places after the decimal point. ^[1] The statement gives no checkpoint values. We write \(a_k=900-3k\) for the coefficient of the \(k\)-th term, \(S(r)\) for \(s(5000)\) as a function of \(r\), and \(T=-6\cdot10^{11}\) for the target.

The whole solution in one sentence: the coefficients \(900-3k\) change sign only once, so after division by a suitable power of \(r\) the difference \(S(r)-T\) falls steadily as \(r\) grows and the equation has exactly one positive solution, which lies between \(1\) and \(1.2\); halving this interval again and again, with the sum added up term by term in floating-point arithmetic, pins it down far more precisely than \(12\) decimals.

Mathematical Approach

We find a closed form, prove that there is exactly one solution and that it lies in the interval searched by the programs, and explain why bisection in floating-point arithmetic is precise enough; a five-term model serves as the worked example.

1. The coefficients and the value at \(r=1\)

The coefficients start with \(a_1=897\), \(a_2=894\), \(a_3=891\), fall by \(3\) per step to \(a_{299}=3\), vanish at \(k=300\) and are negative afterwards, down to \(a_{5000}=-14100\). So \(299\) terms are positive, one is zero and \(4700\) are negative, and the negative ones carry the high powers of \(r\). At \(r=1\) every power equals \(1\), and \(S(1)\) is the sum of an arithmetic progression, the number of terms times the mean of the first and the last coefficient:

$$\begin{aligned}S(1)&=5000\cdot\frac{897-14100}{2}\\&=-33\,007\,500.\end{aligned}$$

The positive coefficients add up to \(134\,550\) and the negative ones to \(-33\,142\,050\). The target is about \(18\,000\) times larger in magnitude, so the negative terms must be amplified, and the factor \(r^{k-1}\) amplifies the late terms most when \(r\gt1\).

2. A five-term model

With \(9\) in place of \(900\) and five terms we get the small sum

$$\begin{aligned}q(r)&=\sum_{k=1}^{5}(9-3k)\,r^{k-1}\\&=6+3r-3r^3-6r^4,\end{aligned}$$

whose coefficients \(6, 3, 0, -3, -6\) are again positive, zero and negative in this order, and we solve \(q(r)=-6\). The value \(q(1)=0\) lies above the target, and so does \(q(1.1)=-3.4776\), while \(q(1.2)=-8.0256\) lies below it, so the model equation has a solution between \(1.1\) and \(1.2\).

3. A closed form

The geometric series \(1+r+\dots+r^{m-1}\)\({}=(1-r^m)/(1-r)\), for \(r\ne1\), is found by subtracting \(r\) times the sum from the sum, which cancels all terms but two. ^[3] The same trick works here, because neighbouring coefficients always differ by \(a_k-a_{k-1}=-3\). Write \(s=s(n)\). The sum \(r\,s\) has the coefficient \(a_{k-1}\) at \(r^{k-1}\), so

$$s-r\,s=897-3\sum_{k=1}^{n-1}r^k-a_n\,r^n.$$

Summing the geometric series, multiplying by \(1-r\) and collecting terms gives, for \(r\ne1\), ^[2]

$$\begin{aligned}&(1-r)^2\,s(n)=897-900r\\&\quad+r^n\bigl(3n-897-(3n-900)\,r\bigr).\end{aligned}$$

For \(n=2\) and \(r=2\) the right-hand side is \(897-1800+4\cdot(-891+1788)\)\({}=2685\), and indeed \((1-2)^2\,s(2)=897+894\cdot2=2685\). For the problem itself, \(n=5000\), the formula reads

$$\begin{aligned}&(1-r)^2\,S(r)=897-900r\\&\quad+r^{5000}\,(14103-14100r).\end{aligned}$$

At \(r=1\) both sides vanish. The programs add up the terms directly instead, which needs no special case at \(r=1\) and avoids the loss of digits the formula suffers near \(r=1\), where it subtracts nearly equal numbers; here the formula serves to rule out negative solutions.

4. Exactly one solution

For \(r\lt0\) both parts of the right-hand side are positive: \(897-900r\gt0\), and \(r^{5000}\gt0\) multiplies \(14103-14100r\gt0\). So \(S(r)\gt0\) for \(r\lt0\), and \(S(0)=a_1=897\); no solution has \(r\le0\). For \(r\gt0\) consider

$$g(r)=\frac{S(r)-T}{r^{299}},$$

which has the same sign as \(S(r)-T\). Since \(a_{300}=0\), dividing every term of \(S(r)-T\) by \(r^{299}\) gives

$$\begin{aligned}g(r)&=\frac{897-T}{r^{299}}+\sum_{k=2}^{299}\frac{a_k}{r^{300-k}}\\&\quad+\sum_{k=301}^{5000}a_k\,r^{k-300}.\end{aligned}$$

Each term strictly decreases as \(r\) grows: in the first two parts a positive number is divided by a growing power of \(r\), where \(897-T\gt0\) because \(T\lt0\), and in the last part a negative number is multiplied by a growing power. So \(g\) is strictly decreasing for \(r\gt0\), and \(S(r)=T\) has at most one solution. Moreover \(g(1)=S(1)-T\gt0\), so a solution must exceed \(1\). This is the simplest case of Descartes' rule of signs: the nonzero coefficients of the polynomial \(S(r)-T\) change sign exactly once, so it has exactly one positive root. ^[4] The same argument with \(0\) in place of \(T\) shows that \(S(r)/r^{299}\) is negative and decreasing for \(r\ge1\); hence \(S\) itself is negative and strictly decreasing there.

For the model, dividing \(q(r)+6\) by \(r^2\) gives \(12/r^2+3/r-3r-6r^2\), which is strictly decreasing for the same reason, so its solution between \(1.1\) and \(1.2\) is its only positive one; with an odd number of terms it also has a negative one, as \(q(-1)=0\) and \(q(-2)=-72\), which the even exponent \(5000\) rules out for \(S(r)\).

5. The search interval

The programs search between \(1\) and \(1.2\). At \(r=1\) we have \(S(1)\gt T\). At \(r=1.2\) each positive term \(a_k\,1.2^{k-1}\), \(k\le299\), is at most \(a_k\,1.2^{298}\), so together they are at most \(134\,550\cdot1.2^{298}\approx5.3\cdot10^{28}\), while the last term alone is \(-14100\cdot1.2^{4999}\approx-9.5\cdot10^{399}\) and all remaining terms are negative or zero. Hence \(S(1.2)\lt T\) by an enormous margin; in fact \(S(1.2)\approx-5.7\cdot10^{400}\). A polynomial is continuous, so by the intermediate value theorem \(S\) takes the value \(T\) between \(1\) and \(1.2\), ^[5] and by Section 4 this is the only solution. We call it \(r^*\).

6. Bisection

Bisection keeps an interval \([\ell,h]\) with \(S(\ell)\gt T\ge S(h)\), which by Section 4 contains \(r^*\). It evaluates the sum at the midpoint \(m=(\ell+h)/2\) and keeps the half that still has this property: \(\ell\) becomes \(m\) if \(S(m)\gt T\), and otherwise \(h\) becomes \(m\). ^[5] Every step halves the width. For the model equation \(q(r)=-6\), starting from \([1,1.2]\):

• \(m=1.1\): \(q(m)=-3.4776\gt-6\), new interval \([1.1,1.2]\);
• \(m=1.15\): \(q(m)\approx-5.6067\gt-6\), new interval \([1.15,1.2]\);
• \(m=1.175\): \(q(m)\approx-6.7785\lt-6\), new interval \([1.15,1.175]\);
• \(m=1.1625\): \(q(m)\approx-6.1833\lt-6\), new interval \([1.15,1.1625]\);
• \(m=1.15625\): \(q(m)\approx-5.8927\gt-6\), new interval \([1.15625,1.1625]\).

After five steps the width is \(0.2/2^5=0.00625\), and further steps close in on the model solution, approximately \(1.158569\). In general \(i\) steps leave the width \(0.2/2^i\), which is below \(10^{-13}\) after \(41\) steps.

7. Floating-point numbers

A double-precision number has a \(53\)-bit significand; between \(1\) and \(2\) neighbouring doubles are \(2^{-52}\approx2.2\cdot10^{-16}\) apart, and the largest finite double is about \(1.8\cdot10^{308}\). ^[6] With the default rounding to nearest, every addition, subtraction and multiplication returns the exact result rounded to the nearest double, so for results in the normal range its relative error is at most \(u=2^{-53}\approx1.1\cdot10^{-16}\). ^[7] Two consequences matter. First, the interval cannot shrink below one spacing: in the actual run its ends become neighbouring doubles after \(50\) steps, as \(0.2/2^{50}\lt2^{-52}\) suggests, and the remaining \(150\) steps change nothing. Second, \(1.2^{4999}\approx6.7\cdot10^{395}\) would overflow, but the programs never evaluate the sum at \(1.2\); at the largest point they evaluate, the first midpoint \(1.1\), even \(1.1^{5000}\approx9.2\cdot10^{206}\) is finite.

8. How accurate the computed sum is

The programs form the powers by repeated multiplication, \(p_1=1\) and \(p_{k+1}=p_k\,r\), and add the terms \(a_k\,p_k\) one after another; the coefficients are small integers and exact. The computed \(p_k\) has gone through \(k-1\) roundings, the product \(a_k\,p_k\) adds one, and the running total is rounded once in each of the at most \(n-k+1\) additions from the one that brings this term in to the last one. So every term enters the computed sum \(\widehat S(r)\) multiplied by at most \(n+1\) factors \(1+\delta\) with \(|\delta|\le u\), and with \(\gamma=(1+u)^{n+1}-1\approx(n+1)\,u\),

$$\bigl|\widehat S(r)-S(r)\bigr|\le\gamma\sum_{k=1}^{n}|a_k|\,r^{k-1}.$$

The contribution of the additions is the classical error estimate for summing numbers one by one. ^[7] At \(r^*\) the sum of absolute values is about \(6.0\cdot10^{11}\), almost all of it from the negative terms, and \(\gamma\approx5.6\cdot10^{-13}\), so the computed sum is off by less than \(0.34\). In practice rounding errors partly cancel: at \(2001\) points near \(r^*\) the double-precision sum is never off by more than \(0.0062\).

9. From the sum to the twelve decimals

An error \(\varepsilon\) in the computed sum can reverse a comparison only at points \(m\) with \(|S(m)-T|\le\varepsilon\). By Section 4 the sum is strictly decreasing for \(r\ge1\), and near \(r^*\) it falls steeply, \(S'(r)\approx-2.76\cdot10^{15}\). With \(\varepsilon\lt0.34\), such points lie within about \(0.34/(2.76\cdot10^{15})\approx1.2\cdot10^{-16}\) of \(r^*\), half a spacing of doubles, so the final pair of neighbouring doubles lies within about one and a half spacings of \(r^*\); in the actual run the returned value is \(1.8\cdot10^{-16}\) away from \(r^*\). The rounded value jumps only at the numbers \((j+\tfrac12)\cdot10^{-12}\) with an integer \(j\), and \(r^*\) is more than \(10^{-13}\) away from all of them, so double precision is ample. This is an estimate, but an exact proof is easy too: if \(c\) is the printed value, the numbers \(c\mp\tfrac12\cdot10^{-12}\) are rational, and evaluating \((1-r)^2\,(S(r)-T)\) at them with exact integer arithmetic gives a positive value at the smaller one and a negative value at the larger one. By Sections 4 and 5, \(r^*\) lies strictly between them, so \(c\) is its correct rounding.

10. Other ways

Newton's method, \(r\leftarrow r-(S(r)-T)/S'(r)\), converges quadratically once it is close to a simple root, but the start matters. ^[8] From \(r=1\), where \(S'(1)=-113\,752\,245\,000\), the first step jumps to \(r\approx6.27\), where the sum overflows in double precision. Even in exact arithmetic the method then crawls back: there \(S\) behaves roughly like a multiple of \(r^{4999}\), each step lowers \(r\) by only about \(r/4999\), and \(9168\) steps are needed to get within \(10^{-15}\) of \(r^*\). Applied to \(\ln(-S(r))-\ln(-T)\), defined for \(r\ge1\) by Section 4, Newton's method from \(r=1\) needs only \(4\) steps, because the slope of \(\ln(-S(r))\) stays between \(3446\) and \(4861\) on \([1,1.1]\). Bisection needs neither a derivative nor a good start.

How the Code Works

The C++ program ^[10] keeps \(n=5000\), the target as a long double and a checkpoint flag in an Options record. parse_arguments recognises --n=, --target= and --skip-checkpoints, and main returns exit code \(1\) if parsing fails or \(n\lt1\); the parsers neither guard against overflow of \(n\) nor reject trailing characters or NaN in the target. Unless it is skipped, run_checkpoints requires series_sum(5000, 1) to lie within \(10^{-6}\) of \(-33\,007\,500\), the value of Section 1, and otherwise main returns exit code \(2\); the loop meets this value exactly, because at \(r=1\) every partial sum is an integer of absolute value far below \(2^{53}\).

series_sum(n, r) is the loop of Section 8: power starts at \(1\), and the \(k\)-th pass adds \(900-3k\) times power to sum and then multiplies power by \(r\). solve(n, target) is the bisection of Section 6 with lo \(=1\), hi \(=1.2\) and a fixed \(200\) passes; it moves lo to the midpoint if the sum there exceeds the target and hi otherwise, and returns the midpoint of the final interval, which main casts to double and prints with std::fixed and std::setprecision(12). Neither the bracket nor the direction of the comparison is checked for other parameters: if the solution lies outside \([1,1.2]\), or if the sum increases with \(r\), as for \(n\le300\), the loop silently returns an end of the interval, and --target=0 prints 1.000000000000.

The precision of long double depends on the platform: typical GCC and Clang builds for x86-64 Linux use the 80-bit extended format with a \(64\)-bit significand, for which an emulation stops changing the interval after \(61\) steps instead of \(50\), while with Microsoft's compiler and on 64-bit ARM under macOS or Windows it has the format of double. ^[9] Either way \(200\) steps are more than enough. The Python and Java programs repeat the two loops in double precision; Python's solve takes \(n\) and the target as optional arguments, Java's solve() takes none, both return the formatted result, and neither program has command-line options or checkpoints. Java formats with String.format(Locale.US, "%.12f", ans), which guarantees a decimal point instead of, say, the decimal comma of a German locale.

Complexity and Verification

One evaluation of the sum makes \(n=5000\) passes of a few floating-point operations and the bisection makes \(200\) evaluations, so a run performs \(10^6\) term updates, plus \(5000\) for the C++ checkpoint, in \(O(1)\) memory. In general the cost is \(O(nI)\) for \(I\) steps, and in double precision about \(\log_2(0.2/2^{-52})\approx50\) steps suffice. On the test machine all three programs finish in well under a tenth of a second. The only built-in check is the C++ checkpoint at \(r=1\).

For this article the following checks were run separately; none of them is part of the published solutions.

• Exact rational arithmetic confirms the closed form of Section 3 against the direct sum for every \(n\le40\) at eight rational values of \(r\) between \(-2\) and \(5\), and all values of the model.
• A 60-digit bisection on the closed form, rounded to \(12\) decimals, agrees with the output of the three programs, which print identical strings, and the exact sign test of Section 9 confirms this rounding.
• The double-precision bisection stops changing after \(50\) steps and ends \(1.8\cdot10^{-16}\) from \(r^*\); the emulation with a \(64\)-bit significand stops after \(61\) steps, \(2.1\cdot10^{-20}\) from \(r^*\). The other figures of Sections 8 to 10 come from computations with \(30\) to \(60\) digits.

Footnotes and References

The references below supply the background facts used above; the uniqueness argument, the error bounds and the exact rounding test are derived in this article.

1. Project Euler 235 — An Arithmetic Geometric Sequence. The official statement defines u(k) = (900 − 3k)·r^(k−1) and s(n) = u(1) + … + u(n) and asks for the r with s(5000) = −600 000 000 000, rounded to 12 decimal places.
2. Wikipedia — Arithmetico-geometric sequence. The first n terms of (a + (k − 1)d)·b·r^(k−1) add up to [ab − (a + nd)·b·rⁿ]/(1 − r) + d·b·r·(1 − rⁿ)/(1 − r)², found by subtracting r times the sum from the sum; a = 897, d = −3, b = 1 give Section 3.
3. Wikipedia — Geometric series. Subtracting r times a finite geometric sum from the sum leaves its first and last terms only, which gives the closed form for r ≠ 1.
4. Wikipedia — Descartes' rule of signs. With its nonzero terms ordered by exponent, a real polynomial has as many positive roots, counted with multiplicity, as there are sign changes between consecutive coefficients, or fewer by an even number; with at most one sign change the two numbers are equal.
5. Wikipedia — Bisection method. For a continuous function with values of opposite sign at two points, the intermediate value theorem gives a root between them; each step keeps the half with the sign change.
6. Wikipedia — Double-precision floating-point format. 53 bits of precision, 52 of them stored; numbers between 1 and 2 are 2⁻⁵² apart, the largest is about 1.8 × 10³⁰⁸, and rounding to the nearest representable number has a relative error of at most 2⁻⁵³.
7. David Goldberg — What Every Computer Scientist Should Know About Floating-Point Arithmetic. Computing Surveys, 1991. The IEEE standard requires the basic operations to be exactly rounded, by default to the nearest floating-point number; a sum added up term by term equals the exact sum of perturbed terms, the first perturbed by about nε and the last by ε.
8. Wikipedia — Newton's method. For a continuously differentiable f with a nonzero derivative and a second derivative at the root, x − f(x)/f′(x) converges at least quadratically from a close enough start, and it can fail from a poor one.
9. Wikipedia — long double. Most compilers on x86, GCC among them, implement it as the 80-bit extended format, whose 64-bit significand is described in Wikipedia — Extended precision; with Visual C++, and on 64-bit ARM under Windows and macOS, it has the same format as double.
10. C++ — Euler235.cpp. The linked immutable revision contains series_sum, solve and the checkpoint; it is the implementation record for this article.