3
$\begingroup$

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

$\endgroup$
7
  • $\begingroup$ The reference at oeis.org/A054347 might be of some (small) help. $\endgroup$ Sep 18 '17 at 15:25
  • $\begingroup$ 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? $\endgroup$
    – Somos
    Sep 18 '17 at 18:28
  • $\begingroup$ I tried looking at OEIS, but no such luck! $\endgroup$
    – vamsikal
    Sep 19 '17 at 3:18
  • $\begingroup$ Again, please give the first few values of $S(N)$ (perhaps ten is enough). That would help tremendously. $\endgroup$
    – Somos
    Sep 20 '17 at 23:41
  • $\begingroup$ I have added a few sample values of $(N, S(N))$. Let me know, if you need more values of $S(N)$. $\endgroup$
    – vamsikal
    Sep 22 '17 at 5:50
2
$\begingroup$

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:

plot for N=20

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:

plot for N=20, filled

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

hierarchical decomposition

(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.

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

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.

$\endgroup$
1
  • $\begingroup$ 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$. $\endgroup$
    – Ahmad
    Sep 25 '17 at 10:51
0
$\begingroup$

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)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.