Summation involving floor of irrational numbers.

Question

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)]$.

Answer 1

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.

Answer 2

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 \begin{equation} S_1(N) = \sum_{x=2}^{x_0} \big( \sum_{y=x+1}^{\lfloor x * F_{n+1} / F_n \rfloor } x+y \big) \end{equation} 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.

Answer 3

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)