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

Question

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.

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

Here's the entry: http://oeis.org/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

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'

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)

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

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 :-)