ES - Editorial (Unofficial)



Author: Gaoyuan Chen
Editorialist: Tiny Wong




Beatty sequence,
Continued fractions.


Given n \leq 10^{4000}, calculate the sum

\sum_{k=1}^n \left\lfloor k \mathrm{e}\right\rfloor,

where \mathrm{e} = 2.718\dots is Euler’s number.


Let S(n, \alpha) = \sum\limits_{k=1}^n \left\lfloor k \alpha \right\rfloor. If 1 < \alpha < 2, we can get that

S(n, \alpha) = \frac 1 2 (n+n') (n+n'+1) - S(n', \beta),

where n' = \left\lfloor (\alpha-1)n \right\rfloor and \beta = \dfrac \alpha {\alpha-1}.


#1. Universal Solution

Let \alpha > 0 be an irrational number. Then the sequence B_\alpha = \{ \left\rfloor \alpha \right\rfloor, \left\rfloor 2\alpha \right\rfloor, \dots, \left\rfloor n\alpha \right\rfloor, \dots \} is a Beatty sequence, whose prefix sum is denoted as

S(n, \alpha) = \sum_{k=1}^n \left\rfloor k \alpha\right\rfloor.

If \alpha > 2, let \beta = \alpha-1, then

S(n, \alpha) = \sum_{k=1}^n \left\lfloor{k \alpha}\right\rfloor = \sum_{k=1}^n \left\lfloor{k \beta + k}\right\rfloor = \sum_{k=1}^n \left\lfloor{k \beta}\right\rfloor + \sum_{k=1}^n k = S(n, \beta) + \frac 1 2 k (k+1).

If \alpha < 1, let \beta = \alpha+1, then

S(n, \alpha) = S(n, \beta) - \frac 1 2 k (k+1).

Then we can assume that 1 < \alpha < 2.

Here is an important theorem for solving this problem.

Theorem (Rayleigh). Let \alpha and \beta be two irrational numbers with \dfrac 1 \alpha + \dfrac 1 \beta = 1, then B_\alpha and B_\beta partition the set of positive numbers, i.e. B_\alpha \cup B_\beta = \mathbb{N}^+ and B_\alpha \cap B_\beta = \emptyset.

By Rayleigh Theorem, if 1 < \alpha < 2, let \beta = \dfrac \alpha {\alpha - 1}, m = \left\lfloor{n \alpha}\right\rfloor and n' = \left\lfloor{\dfrac m \beta}\right\rfloor, then

S(n, \alpha) + S(n', \beta) = \sum_{i=1}^m i = \frac 1 2 m ( m+1 ).

Note that

n' = \left\lfloor{\frac m \beta}\right\rfloor = \left\lfloor{\frac {m(\alpha-1)} \alpha}\right\rfloor = m - \left\lceil{\frac m \alpha}\right\rceil = m-n = \left\lfloor{(\alpha-1)n}\right\rfloor,


S(n, \alpha) = \frac 1 2 m ( m+1 ) - S(n', \beta) = \frac 1 2 (n+n') (n+n'+1) - S(n', \beta). (*)

If 1 < \alpha < 2, n' = \left\lfloor{(\alpha-1)n}\right\rfloor < (\alpha-1)n < n, we can calculate S(n, \alpha) by (*).

#2. Continued Fractions

A positive irrational number \alpha can be written as a continued fraction

\alpha = a_0 + \dfrac 1 {a_1 + \dfrac 1 {a_2 + \dfrac 1 {a_3 + \ddots}}},

where \{ a_0, a_1, a_2, \dots \} is an infinite sequence and a_1, a_2, \dots \in \mathbb{N}^+. We can use the shorthand

\alpha = [a_0; a_1, a_2, a_3, \dots].

Note that the integer part of a is [\alpha] = a_0, and the fraction part is

\{ \alpha \} = [0; a_1, a_2, a_3, \dots].


\frac 1 {\{ \alpha \}} = [a_1; a_2, a_3, \dots].

In the former analysis, we assume that 1 < \alpha < 2, as a result a_0 = 1, and

\beta = 1 + \frac 1 {\alpha - 1} = 1 + \frac 1 {\{ \alpha \}} = [a_1+1; a_2, a_3, \dots].
n' = \left\lfloor{(\alpha-1)n}\right\rfloor = \left\lfloor{\{\alpha\} n}\right\rfloor < \frac n {a_1}.

#3. Complexity Analysis

Now we describe the algorithm formally. Let \alpha^{(0)} = \alpha and n^{(0)} = n. Now for each k, consider
the calculation of S(n^{(k)}, \alpha^{(k)}). Let \alpha^{(k)} = [a_0^{(k)}; a_1^{(k)}, a_2^{(k)}, a_3^{(k)}, \dots].

Let $n^{(k+1)} = \left\lfloor{ {\alpha^{(k)}} n^{(k)}}\right\rfloor$,$\beta^{(k)} = 1+\dfrac 1 {{\alpha^{(k)}}}, by (*)$,

S(n^{(k)}, \alpha^{(k)}) = \frac 1 2 (n^{(k)}+n^{(k+1)}) (n^{(k)}+n^{(k+1)}+1) - S(n^{(k+1)}, \beta^{(k)}).

Notice that

\beta^{(k)} = [a_1^{(k)}+1; a_2^{(k)}, a_3^{(k)}, \dots].

And then

S(n^{(k+1)}, \beta^{(k)}) = S(n^{(k+1)}, \beta^{(k)}-a_1^{(k)}) + \frac 1 2 a_1^{(k)} n^{(k+1)} (n^{(k+1)}+1).

If we let \alpha^{(k+1)} = \beta^{(k)}-a_1^{(k)} = [1; a_2^{(k)}, a_3^{(k)}, \dots], then a_1^{(k+1)} = a_2^{(k)} = \dots = a_{k+1}^{(0)}, and we get that

S(n^{(k)}, \alpha^{(k)}) = \frac 1 2 (n^{(k)}+n^{(k+1)}) (n^{(k)}+n^{(k+1)}+1) - \frac 1 2 a_{k+1}^{(0)} n^{(k+1)} (n^{(k+1)}+1) - S(n^{(k+1)}, \alpha^{(k+1)}).

Note that

\alpha^{(k)} = [1; a_k^{(0)}, a_{k+1}^{(0)}, \dots].

and that

\begin{aligned} n^{(k+1)} & = \left\lfloor{\{\alpha^{(k)}\} n^{(k)}}\right\rfloor = \dots \\ & = \left\lfloor{\{\alpha^{(k)}\}}\right\rfloor \left\lfloor{\{\alpha^{(k-1)}\}}\right\rfloor \dots \left\lfloor{\{\alpha^{(0)}\}}\right\rfloor n^{(0)} \\ & < \{\alpha^{(k)}\} \{\alpha^{(k-1)}\} \dots \{\alpha^{(0)}\} n^{(0)}. \end{aligned}

Let \gamma^{(k)} = \dfrac 1 {\{\alpha^{(k)}\}} = [a_1^{(k)}; a_2^{(k)}, a_3^{(k)}, \dots] = [a_{k+1}^{(0)}; a_{k+2}^{(0)}, a_{k+3}^{(0)}, \dots], then

n^{(k+1)} = \frac {n^{(0)}} {\gamma^{(k)} \gamma^{(k-1)} \dots \gamma^{(0)}}.

Note that

\gamma^{(k-1)} = a_{k}^{(0)} + \frac 1 {\gamma^{(k)}},


\gamma^{(k)} \gamma^{(k-1)} = a_{k}^{(0)} \gamma^{(k)} +1.

Here we have \gamma^{(1)} \gamma^{(0)} = a_1^{(0)} \gamma^{(1)} + 1. We may assume that \gamma^{(k)} \gamma^{(k-1)} \dots \gamma^{(0)} = p^{(k)} \gamma^{(k)} + q^{(k)}, and then

\begin{aligned} & \quad \gamma^{(k+1)} \gamma^{(k)} \gamma^{(k-1)} \dots \gamma^{(0)} \\ & = \gamma^{(k+1)} (p^{(k)} \gamma^{(k)} + q^{(k)}) \\ & = p^{(k)} \gamma^{(k+1)} \gamma^{(k)} + q^{(k)} \\ & = p^{(k)} \left( a_{k+1}^{(0)} \gamma^{(k+1)} +1 \right) + q^{(k)} \gamma^{(k+1)} \\ & = \left(p^{(k)}a_{k+1}^{(0)}+q^{(k)}\right) \gamma^{(k+1)} + p^{(k)} \\ & = p^{(k+1)} \gamma^{(k+1)} + q^{(k+1)}. \end{aligned}

As a result, p^{(k+1)} = p^{(k)}a_{k+1}^{(0)}+q^{(k)} and q^{(k+1)} = p^{(k)}, and then

p^{(k+1)} = p^{(k)} a_{k+1}^{(0)}+p^{(k-1)} \geq p^{(k)}+p^{(k-1)}.

Note that p^{(0)} = 1 and p^{(1)} = a_1^{(0)} \geq 1, then p^{(k)} \geq f_k, where \{ f_k \} is the Fibonacci sequence, which is f_0 = f_1 = 1, and f_k = f_{k-1}+f_{k-2} for all $k \geq 2$。

Again we have

\begin{aligned} & \quad \gamma^{(k)} \gamma^{(k-1)} \dots \gamma^{(0)} \\ & = \frac {\gamma^{(k+1)} \gamma^{(k)} \gamma^{(k-1)} \dots \gamma^{(0)}} {\gamma^{(k+1)}} \\ & = \frac {p^{(k+1)} \gamma^{(k+1)} + q^{(k+1)}} {\gamma^{(k+1)}} \\ & = p^{(k+1)} + \frac {p^{(k)}} {\gamma^{(k+1)}} \\ & \geq p^{(k+1)}+p^{(k)} \\ & \geq f_{k+1}+f_{k} = f_{k+2}. \end{aligned}


n^{k+1} = \frac {n^{(0)}} {\gamma^{(k)} \gamma^{(k-1)} \dots \gamma^{(0)}} \leq \frac {n^{(0)}} {f_{k+2}} = \frac {n} {O(\phi^{k})},

where \phi = \dfrac {1+\sqrt5} 2.

As a result, the number of iterations of (*) is $O(\log n)$。

#4. Details

In practice, we need only two values, which are n^{(k)} and \{ \alpha^{(k)} \} and can be calculate by the followings.

\begin{aligned} & n^{(k+1)} = \left\lfloor{ \{ \alpha^{(k)} \} n^{(k)} }\right\rfloor, \\ & \{ \alpha^{k+1} \} = \left\{ \frac 1 {\{\alpha^{(k)}\}} \right\}. \end{aligned}

Algorithm 1. Calculate all floating numbers after dealing the errors.

Algorithm 2. Use continued fractions instead of painful floatings.

An irration number \alpha can be approximated by a sequence of continued fractions

\{ [a_0], [a_0; a_1], [a_0; a_1, a_2], \dots, [a_0; a_1, a_2, \dots, a_n], \dots \}.

whose each term is an approximate representaion of \alpha and the precision is increasing.

Suppose the k-th term of the above sequence is \dfrac {P_k} {Q_k}, then

\frac {P_k} {Q_k} = \frac {a_k P_{k-1} + P_{k-2}} {a_k Q_{k-1} + Q_{k-2}}.

What we are facing is to calculate n^{(k+1)} = \left\lfloor{ \{ \alpha^{(k)} \} n^{(k)} }\right\rfloor. If \{ \alpha^{(k)} \} \approx \dfrac P Q with Q > n, then we can just calculate it forward

n^{(k+1)} = \left\lfloor{ \frac {P n^{(k)}} Q }\right\rfloor.

No matter which algorithm you choose, the overall time complexity is O(\log^2 n \log \log n). That is because O(\log n) for iterations, O(\log n \log \log n) for big integer multiplications and divisions (More precisely, multiplication is O(\log n) and division is O(\log n \log \log n)).


There may be other ways to solve the problem, please share!


Editorialist’s solution can be found here.


My preview is correct. I don’t know why it looks awful T_T!


Indeed, the forum latex preview is inconsistent with the real thing. You can try changing \floor{} to \lfloor\rfloor

Any C++ soln of this problem??

@hikarico It was a hard time for me. Now it looks correctly.

1 Like

@vijju123 There are no C++ solutions because of the limit of the length of source code.

1 Like

Can someone explain why this problem cannot have a CPP solution?

because of the limit of the length of source code.

Only a suggestion, but can we have lesser problems of these types that need intense research and have less reusability in other contests. Being one of the hardest problems the research and efforts do pay off but I personally would like to see problems with rarer paradigms/algorithms that can be used (at least partially) in short contests for hard problems like reduction to multiple usual algorithms and they fit together in some complex way, etc

Like SANDWICH problem last time was good, I actually learned Wilson’s Theorem, Lucas Theorem and bunch of other mathematical properties - they weren’t enough to get full but still useful for other contests. To get full I had to get into generalized Lucas that paid off with 100 points.

Anyways just a suggestion, the problem set was good and challenging for my level.

On a side note, Is it just me that haven’t seen a medium-hard DP problem in a top 7 - 8 most solved problems in Long recently ? :smiley:


@ayushgupta321 it has to do with the precision of e. In order to compute this Beatty sequence with n \leqslant 10^{4000}, e must have a 4000-digit precision at least. I think this is too much for cpp.

I actually used the formula from Quick Explanation section but alas got only 50 pts. Decimal is too slow and generating e upto 4000 digits takes me 4.5s and N as large as 10^4000 takes about <= 2000 recursion depth to decay to zero

from sys import stdin, stdout
from decimal import Decimal as D, getcontext as gc, ROUND_FLOOR

def s(a, n):
	if n == 0:
		return 0
	_n = ((a - D(1)) * n).to_integral_exact()
	b = a / (a - D('1.0'))
	k = b.to_integral_exact() - D(1)
	b = b - k
	n_n = n + _n
	tmp = (n_n * n_n)
	tmp = tmp + n_n
	tmp = tmp - (k * (_n) * (_n + D(1)))
	tmp = tmp / D(2)
	tmp = tmp - s(b, _n)
	return tmp

m = stdin.readline()
gc().prec = min(len(m) * 2 + 3, 4003)
gc().rounding = ROUND_FLOOR

m = D(m)
e = D('1.5')
f = D('0.5')
for i in range(3, 1506):
	f = f / D(i)
	e = e + f

ans = ((D(m) * D(m + 1))/D(2)) + s(e, m)

My solution uses this awesome paper :
T.C. Brown and P.J.-S. Shiue, Sums of fractional parts of integer multiples of an irrational, J.
Number Theory 50 (1995), 181–192.

It ran in 0.38s on python :


I think there are some ideas from ES that can be used in the short contests. I ever met a similar problem in short contests, which also calculates the sum but this time e becomes a simple fraction, i.e. calculates the grid point in the triangle, and to solve this problem also needs such ideas.

After all, these problems are “projecteuler”-style problems.

And the same as you, I would like some medium DP problems in long challenge too~

1 Like

It looks cool! Could you please give us a brief idea about it?

This can’t have a solution for cpp…wht’s the use then? -_-

1 Like

I also complain that there are no cpp solutions. However, I think the author’s making the source code limit under 1000 bytes is to tell us there are no complex codes but cute ideas?

1 Like

I suggest you use just big integers instead of big decimals, which would cause a
significant speed up.

I am basically using lemma 4.
In the first part, we are calculating the p/q form for e till q < Number.
using continued fraction approach and storing then in num , denom array.

in the solver function :

S(N) = solution for N,

since N > qn , we write :

N = b * qn + k,
from lemma 4:

S(b * qi + c) = S(b * qi) + S(k) + k * b * pi

S(b * qi) = (1/2)* b*(term) where term = b * pi * qi - qi + pi + (-1)^i
here qi is ith denominator, similarly for pi.

we find qi for K, and using same above forumula, solve for S(K).
this will surely end, since q0 and q1 are 1.
recursion wont work.


Just one (lame, maybe) question, how is any codes’ the size calculated?

the count of ASCII chars, including ‘\n’, ’ ', etc.