# Summation involving floor of irrational numbers.

I need to compute the following sum for large $N$, i.e. $N = 10^{16}$.

$$S(N) = \sum_{x=2}^{x=N-1} \left(\sum_{y=x+1}^{y=min(N, \lfloor \phi \times x \rfloor)} x + y\right)$$

where $\phi = \frac{1 + \sqrt{5}}{2}$, $x$ and $y$ are positive integers.

My first step is to break the sum at $x = x_0$ such that $(x_0 \times \phi) \leq N$ and $((x_0 + 1) \times \phi) > N$. So, we now have two sums,

$$S_1(N) = \sum_{x=2}^{x=x_0} \left(\sum_{y=x+1}^{y=\lfloor \phi \times x \rfloor} x + y\right)$$ $$S_2(N) = \sum_{x=x_0+1}^{x=N-1} \left(\sum_{y=x+1}^{y=N} x + y\right)$$

Now, note that $S_2(N)$ has a closed form and can be computed very quickly. My question is, is there a way to compute $S_1(N)$ quickly, any sub-linear time algorithm would be great.

Here are some values of $(N, S(N))$ for small values of $N$, $[(2,0),(3,5),(4,12),(5,21),(6,42),(7,67),(8,109),(9,157),(10,211),(11,289),(12,375)]$.

• The reference at oeis.org/A054347 might be of some (small) help. Sep 18 '17 at 15:25
• What have you tried? In other words, $S(N)$ is an integer sequence. What are the first few values? Did you try to look it up in the OEIS? Sep 18 '17 at 18:28
• I tried looking at OEIS, but no such luck! Sep 19 '17 at 3:18
• Again, please give the first few values of $S(N)$ (perhaps ten is enough). That would help tremendously. Sep 20 '17 at 23:41
• I have added a few sample values of $(N, S(N))$. Let me know, if you need more values of $S(N)$. Sep 22 '17 at 5:50

I would like to know/learn how to compute such summations.

In this case, we can use $\phi$'s connection to Fibonacci recurrences, as well as the fact that $x+y$ is linear. Start with a plot of the summands:

Then, the problem is to find the (weighted) area. You've already noticed that the $\min(N,\cdot)$ limit is not interesting, so let's get rid of it. While we're at it, let's fill the triangle to the axis:

To find the total weight in this plot, split it recursively into smaller pieces. There is a Fibonacci-like recurrence:

(You could also derive the recurrence from the continued fraction representation $\phi = [1;1,1,\dots]$. It's not as picturesque though.)

It's not quite that simple, because these similar pieces have different weights $\sum_{x,y} x+y$. However, $x+y$ is linear, so when you shift a $k$-term piece by $(\Delta x, \Delta y)$, its weight changes by $k (\Delta x + \Delta y)$. Modulo this adjustment there are only $\mathcal{O}(\log N)$ distinct pieces, so $S(N)$ can be calculated in $\mathcal{O}(\log N)$ operations.

• Thank you very much for solving this summation for me. For background, this summation shows up in the solution to Project Euler problem 325. Sep 27 '17 at 21:14

Here's a partial solution: The best rational approximations to $\phi$ are ratios of successive Fibonacci numbers $\phi \approx F_{n+1} / F_n$. Take the first Fibonacci number $F_n >= 2 x_0$. We have $|\phi - F_{n+1} / F_n| < 1/ F_n^2 <= 1/2x_0 F_n$.

So $|\phi * x - x * F_{n+1} / F_n| < 1/2F_n$, and $\lfloor \phi * x \rfloor = \lfloor x * F_{n+1} / F_n \rfloor$.

We now have $$S_1(N) = \sum_{x=2}^{x_0} \big( \sum_{y=x+1}^{\lfloor x * F_{n+1} / F_n \rfloor } x+y \big)$$ The inner sum can be expressed as a quasi-polynomial in $x$ with period $F_n$. I suspect there is a nice closed form for sums of the form $\sum_{y=x+1}^{x + \lfloor x * p / q \rfloor } x+y$ with $p,q$ prime in terms of the Dirichlet characters mod $q$, but that's as far as I got.

• Side note : See that $\lfloor \phi x \rfloor = \lfloor \phi (x-F_n) \rfloor +F_{n+1}$ Such that $\phi \approx 1.618$ and $F_n$ is the $n$-th Fibonacci number, and $F_n <x$. Sep 25 '17 at 10:51

Thanks to Japheth Lim for providing an answer to this question. Based on his answer, I implemented a Haskell program that computes $S(10^{16})$. I needed to use an extended precision up to 50 digits for $\phi$. Here is sample output:

vamsi@vamsi-laptop:~/learn/project_euler/257_to_384/problem_325$ghci GHCi, version 8.0.2: http://www.haskell.org/ghc/ :? for help Prelude> :l problem_325_fast.hs [1 of 1] Compiling Main ( problem_325_fast.hs, interpreted ) Ok, modules loaded: Main. *Main> s 10 211 *Main> s 10000 230312207313 *Main> s (10 ^ 16) 230327668541684176515052475525037682834286696168  The code itself: import qualified Data.List as L import qualified Data.Map.Strict as M import qualified Data.Ratio as R phi :: Double phi = 0.5 * (1.0 + (sqrt 5.0)) phi_r :: Rational phi_r = 161803398874989484820458683436563811772030917980576 R.% (10 ^ 50) triangular :: Integer -> Integer triangular n | n < 0 = error "triangular: n must be >= 0." | otherwise = (n * (n + 1)) div 2 sum_of_squares :: Integer -> Integer sum_of_squares n | n < 0 = error "sum_of_squares: n must be >= 0." | otherwise = (n * (n + 1) * ((2 * n) + 1)) div 6 fibonacci :: Integer -> [Integer] fibonacci n = loop [2, 1] where loop [] = error "fibonacci: empty list found!" loop [_] = error "fibonacci: singleton list found!" loop fs@(f : g : rest) | f >= n = fs | otherwise = loop ((f + g) : fs) as_fibonacci_sum :: Integer -> [Integer] as_fibonacci_sum n = reverse$ loop fs n []
where fs = fibonacci n
loop [] x ds
| x == 0 = ds
| otherwise = error "as_fibonacci_sum: could not compute sum!"
loop fs@(f : rest) x ds
| (x == f) && (x < n) = (f : ds)
| f < x = loop rest (x - f) (f : ds)
| otherwise = loop rest x ds

as_limits :: [Integer] -> [(Integer, Integer)]
as_limits [] = error "as_limits: found empty list!"
as_limits fs@(f : rest) = reverse $loop rest [(1, f)] where loop [] ls = ls loop (f : rest) ls@((l, u) : _) = loop rest ((u + 1, u + 1 + f - 1) : ls) s_2 :: Integer -> Integer s_2 n | n < 2 = error "s_2: n must be >= 2." | otherwise = t_1 + t_2 - t_3 - t_4 where x_0 = floor ((fromIntegral n) / phi_r) ff x = (x * (x + 1) * (x + 2)) div 6 t_1 = n * ((triangular (n - 1)) - (triangular x_0)) t_2 = (triangular n) * (n - 1 - x_0) t_3 = (sum_of_squares (n - 1)) - (sum_of_squares x_0) t_4 = (ff (n - 1)) - (ff x_0) slow_F :: Integer -> Integer slow_F x_0 = sum [(x + y) | x <- [1..x_0], y <- [1..(floor (phi_r * (fromIntegral x)))]] _G :: Integer -> Integer -> Integer -> Integer _G a b h = (h * (((b * (b + 1)) div 2) - (((a - 1) * a) div 2))) + ((h * (h + 1) * (b - a + 1)) div 2) f :: Integer -> Integer f x = floor (phi_r * (fromIntegral x)) g :: Integer -> Integer g x = (x * h) + ((h * (h + 1)) div 2) where h = floor (phi_r * (fromIntegral x)) slow_k :: Integer -> Integer slow_k x = sum [f y | y <- [1..x]] type Memo = M.Map Integer Integer fast_k :: Memo -> Integer -> (Memo, Integer) fast_k m x | M.member x m = (m, m M.! x) | otherwise = let (m_1, s_1) = L.foldl fold_fn (m, 0) fd in (M.insert x s_1 m_1, s_1) where fd = as_limits$ as_fibonacci_sum x
fold_fn (m, s) x = let (m_1, s_1) = compute_k m x
in (m_1, s + s_1)
compute_k m x@(l, u)
| l == 1 = fast_k m u
| otherwise = let (m_1, k_1) = fast_k m (u - l + 1)
in (m_1, k_1 + ((u - l + 1) * ((f l) - (f 1))))

fast_F_shifted :: Memo -> Memo -> Integer -> Integer -> (Memo, Memo, Integer)
fast_F_shifted m_k m_f x y
| x <= y = (m_k_2, m_f_1, t_1 + (k * (dx + dy)) + (_G x y dy))
| otherwise = error "fast_F_shifted: x must be <= y."
where w = y - x + 1
(m_k_1, m_f_1, t_1) = fast_F m_k m_f w
(m_k_2, k) = fast_k m_k_1 w
dx = x - (f 1)
dy = (f x) - (f 1)

fast_F :: Memo -> Memo -> Integer -> (Memo, Memo, Integer)
fast_F m_k m_f x_0
| M.member x_0 m_f = (m_k, m_f, m_f M.! x_0)
| otherwise = let (m_k_1, m_f_1, s_1) = L.foldl fold_fn (m_k, m_f, 0) fd
in (m_k_1, M.insert x_0 s_1 m_f_1, s_1)
where fd = as_limits $as_fibonacci_sum x_0 -- m = M.insert 2 ((g 1) + (g 2))$ M.singleton 1 (g 1)
fold_fn (m_k, m_f, s) x = let (m_k_1, m_f_1, s_1) = compute_F m_k m_f x
in (m_k_1, m_f_1, s + s_1)
compute_F m_k m_f x@(l, u)
| l == 1 = fast_F m_k m_f u
| otherwise = fast_F_shifted m_k m_f l u

s_1 :: Integer -> Integer
s_1 n
| n < 2 = error "s_1: n must be >= 2."
| otherwise = t_0 - t_1 - t_2 - t_3
where x_0 = floor ((fromIntegral n) / phi_r)
ff x = (x * (x + 1) * (x + 2)) div 6
(_, _, t_0) = fast_F m_k m_f x_0
t_1 = 2
t_2 = (sum_of_squares x_0) - (sum_of_squares 1)
t_3 = (ff x_0) - (ff 1)
m_f = M.insert 2 ((g 1) + (g 2)) $M.singleton 1 (g 1) m_k = M.insert 2 ((f 1) + (f 2))$ M.singleton 1 (f 1)

s :: Integer -> Integer
s n
| n < 2 = error "s: n must be >= 2."
| otherwise = (s_1 n) + (s_2 n)