How do I find the millionth number in the series: 2 3 4 6 9 13 19 28 42 63 ... ?

Question

It takes about minute to achieve 3000 in my comp but I need to know the millionth number in the series. The definition is recursive so I cannot see any shortcuts except to calculate everything before the millionth number. How can you fast calculate millionth number in the series?

Series Def

n_{i+1} = \floor{ 3/2 * n_{i} } and n_{0}=2.

Interestingly, only one site list the series according to Google: this one.

Too slow Bash code

#!/bin/bash

function series 
{
        n=$( echo "3/2*$n" | bc -l | tr '\n' ' ' | sed -e 's@\\@@g' -e 's@ @@g' );
                                        # bc gives \ at very large numbers, sed-tr for it
        n=$( echo $n/1 | bc )           #DUMMY FLOOR func
}

n=2
nth=1

while [ true ]; #$nth -lt 500 ];
do
        series $n                        # n gets new value in the function through global value
        echo $nth $n
        nth=$( echo $nth + 1 | bc )     #n++
done

Answer 1

Hmm, bash is not what I'd be using for high speed numerical processing. Get yourself a copy of GMP and put together a C program to do it.

There may well be a mathematical formula to give it to you quickly but, in the time it takes you to figure it out, GMP could probably throw you the result :-)

Answer 2

You almost found it. Next time, check out the Online Encyclopedia of Integer Series.

Here's the entry:

http://www.research.att.com/~njas/sequences/A061418

     FORMULA    

a(n) = ceiling[K*(3/2)^n] where K=1.08151366859...

The constant K is 2/3*K(3) (see A083286). - Ralf Stephan, May 29, 2003 

That said:

>>> def f(n):
...     K = 1.08151366859
...     return int(ceil(K*(1.5)**n))

Sanity test:

>>> for x in range(1, 10):
...     print f(x)
...     
2
3
4
6
9
13
19
28
42

Awesome! Now how about 1000000:

>>> f(1000000)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "<input>", line 3, in f
OverflowError: (34, 'Result too large')

Well, I tried. :] But you get the idea.

Edit again: Solution is found! See Timo or Lasse V. Karlsen's answers.

Edit: Using Timo's bit-shifting idea:

import gmpy
n=gmpy.mpz(2)
for _ in xrange(10**6-1):
    n+=n>>1
print(n)

yields

1963756763...226123087 (176092 digits)

% time test.py > test.out

real    0m21.735s
user    0m21.709s
sys 0m0.024s

Answer 3

This is almost a first order recurrence relation, except for the floor, which messes things up. If you didn't want the floor,

http://en.wikipedia.org/wiki/Recurrence_relation

Also, don't use bash.

Answer 4

This is identified as sequence A061418 in the sequences site (AKA "The On-Line Encyclopedia of Integer Sequences"); per the relevant page,

FORMULA a(n) =A061419(n)+1 = ceiling[K*(3/2)^n] where K=1.08151366859... The constant K is 2/3*K(3) (see A083286).

and with a suitable high-precision library (GMP as already suggested, or MPIR, and maybe a wrapper on top like my baby gmpy is for Python) you can used the closed-form formula for much speedier computation of "the millionth item in the series" and the like.

It's often possible to put recursively specified recurrences into closed formulas. For an extensive beginner's introduction to the subject, Concrete Mathematics (by Graham, Knuth and Patashnik) is really hard to beat.

Answer 5

I used the following Java program:

import java.math.*;

public class Series {
    public static void main(String[] args) {
        BigInteger num = BigInteger.valueOf(2);
        final int N = 1000000;

        long t = System.currentTimeMillis();
        for (int i = 1; i < N; i++) {
            num = num.shiftLeft(1).add(num).shiftRight(1);
        }
        System.out.println(System.currentTimeMillis() - t);
        System.out.println(num);
    }
}

The output, cropped: (full output on pastebin)

516380 (milliseconds)
196375676351034182442....29226123087

So it took about 8.5 minutes on my modest machine. I used -Xmx128M, but not sure if that was really necessary.

There are probably better algorithms out there, but this took a total of 10 minutes, including writing the naive implementation and running the program.

Sample runs

• N=500 (agrees with OEIS 061418)
• N=100K (1445232038814....522773508)

Answer 6

The recursive formulation is going to take quite a long time under most curcumstances because it has to maintain the machine stack. Why not use dynamic programming instead?

i.e. (pseudocode)

bignum x = 2
for (int i = 1; i < 1000000; i++) {
    x = floor(3.0/2*x)
}

Of course, for a meaningful result you'll need a high precision number library.

Answer 7

The reason your script is so slow is that it's spawning bc three times, tr once and sed once in a loop.

Rewrite the whole thing in bc and do the sed at the end. My test shows that the bc-only version is over 600 times faster. It took just under 16 minutes on an old slow system for the bc version to find the 100,000th value (only printing the last).

Also, note that your "floor" function is actually "int".

#!/usr/bin/bc -lq
define int(n) {
    auto oscale
    oscale = scale
    scale = 0
    n = n/1
    scale = oscale
    return n
}

define series(n) {
    return int(3/2*n)
}

n = 2
end = 1000
for (count = 1; count < end; count++ ) {
    n = series(n)
}
print count, "\t", n, "\n"
quit

Note that print is an extension and some versions of bc may not have it. If so just reference the variable by itself and its value should be output.

Now you can do chmod +x series.bc and call it like this:

./series.bc | tr -d '\n' | sed 's.\\..g'

Answer 8

You can probably get a bit closer by using a more suitable language, e.g., Scheme:

(define (series n) (if (= n 0) 2
                       (quotient (* 3 (series (- n 1))) 2)))

This computes the 17610 digits of (series 100000) in about 8 seconds on my machine. Unfortunately, (series 1000000) still takes far too long for even a much newer/faster machine to have any hope of finishing in a minute.

Switching to C++ with a big-integer library (NTL, in this case):

#include <NTL/zz.h>
#include <fstream>
#include <time.h>
#include <iostream>

int main(int argc, char **argv) {
    NTL::ZZ sum;

    clock_t start = clock();
    for (int i=0; i<1000000; i++) 
        sum = (sum*3)/2;
    clock_t finish = clock();
    std::cout << "computation took: " << double(finish-start)/CLOCKS_PER_SEC << " seconds.\n";

    std::ofstream out("series_out.txt");
    out << sum;

    return 0;
}

This computes the series for 1000000 in 4 minutes, 35 seconds on my machine. That's fast enough that I can almost believe a really fast, new machine might at least come close to finishing in a minute (and yes, I checked what happened when I used shifts instead of multiplication/division -- it was slower).

Unfortunately, the closed-form computation others have suggested seems to be of little help. To use it you need to compute the constant K to sufficient precision. I don't see a closed-form computation of K, so this really just shifts the iteration to computing K, and it looks like computing K to sufficient precision is little (if any) faster that computing the original series.

Answer 9

Here's a Python version that on my 10-year old laptop takes about 220 seconds to run:

import math;
import timeit;

def code():
  n = 2
  nth = 1

  while nth < 1000000:
    n = (n * 3) >> 1
    nth = nth + 1

  print(n);

t = timeit.Timer(setup="from __main__ import code", stmt="code()");
print(t.timeit(1));

It produces the same result that this answer has on the pastebin (that is, I verified the start and end of it, not everything.)

Answer 10

You can easily solve this by thinking about the problem in binary. Floor(3/2*i) is basically shift right, truncate and add. In pseudo code:

0b_n[0]   = 10              // the start is 2
0b_n[1]   = 10 + 1(0) = 11  // shift right, chop one bit off and add 
0b_n[i+1] = 0b_n[i] + Truncate(ShiftRight(0b_n[i]))

This should be quite fast to implement in any form.

I just made an implementation of this in Mathematica and it seems that the BitShiftRight operation also chops off the bit past the unit position, so that gets taken care of automatically. Here is the one liner:

In[1] := Timing[num = Nest[(BitShiftRight[#] + #) &, 2, 999999];]
Out[2] = {16.6022, Null}

16 seconds, the number prints just fine, it is quite long though:

In[2] := IntegerLength[num]
Out[2] = 176092

In[3] := num
Out[3] = 1963756...123087

Full result here.

Answer 11

I converted Timo's ideas to elisp. It fails with 100, giving negative number. FAIL, please, see no BigNums!

(progn
  (let ((a 2)
        (dummy 0))
    (while (< dummy 100)
      (setq a (+ a (lsh a -1)))
      (setq dummy (1+ dummy)))
    (message "%d" a)))
-211190189 #WRONG evalution

Answer 12

It's very easy to do in Pari:

n=2;for(k=1,10^6,n+=n>>1)

This takes 14 seconds on my machine. There are faster ways, of course -- GMP comes to mind -- but why bother? You won't be able to shave more than 10 seconds off the runtime, and development time would be on the order of minutes. :)

Minor point: It's ambiguous in the original formulation whether the millionth term n₉₉₉₉₉₉ is desired or n_1000000, the number with index one million; I give the latter since I see the former is already calculated above.