For loop computing recurrence relation takes very long
<pre><code>Q(x)=[Q(xâ1)+Q(xâ2)]^2
Q(0)=0, Q(1)=1
</code></pre>

<p>I need to find Q(29). I wrote a code in python but it is taking too long. How to get the output (any language would be fine)?</p>

<p>Here is the code I wrote:</p>

<pre><code>a=0
b=1
for i in range(28):
    c=(a+b)*(a+b)
    a=b
    b=c
print(b)
</code></pre>

<p>I don't think this is a tractable problem with programming. The reason why your code is slow is that the numbers within grow <em>very</em> rapidly, and python uses infinite-precision integers, so it takes its time computing the result.</p>

<p>Try your code with double-precision floats:</p>

<pre><code>a=0.0
b=1.0
for i in range(28):
    c=(a+b)*(a+b)
    a=b
    b=c
print(b)
</code></pre>

<p>The answer is <code>inf</code>. This is because the answer is much much larger than the largest representable double-precision number, which is <a href="http://stackoverflow.com/a/3477332/5067311">rougly 10^308</a>. You could try using finite-precision integers, but those will have an even smaller representable maximum. Note that using doubles will lead to loss of precision, but surely you don't want to know every single digit of your huuuge number (side note: <a href="http://math.stackexchange.com/questions/1948492/sum-of-digits-of-p29-the-sum-of-two-recurrence-series">I happen to know that you do</a>, making your job even harder).</p>

<p>So here's some math background for my skepticism: Your recurrence relation goes</p>

<pre><code>Q[k] = (Q[k-2] + Q[k-1])^2
</code></pre>

<p>You can formulate a more tractable sequence from the square root of this sequence:</p>

<pre><code>P[k] = sqrt(Q[k])
P[k] = P[k-2]^2 + P[k-1]^2
</code></pre>

<p>If you can solve for <code>P</code>, you'll know <code>Q = P^2</code>.</p>

<p>Now, consider this sequence:</p>

<pre><code>R[k] = R[k-1]^2
</code></pre>

<p>Starting from the same initial values, this will always be smaller than <code>P[k]</code>, since</p>

<pre><code>P[k] = P[k-2]^2 + P[k-1]^2 &gt;= P[k-1]^2
</code></pre>

<p>(but this will be a "pretty close" lower bound as the first term will always be insignificant compared to the second). We can construct this sequence:</p>

<pre><code>R[k] = R[k-1]^2 = R[k-2]^4 = R[k-3]^6 = R[k-m]^(2^m) = R[0]^(2^k)
</code></pre>

<p>Since <code>P[1 give or take]</code> starts with value 2, we should consider</p>

<pre><code>R[k] = 2^(2^k)
</code></pre>

<p>as a lower bound for <code>P[k]</code>, give or take a few exponents of 2. For <code>k=28</code> this is</p>

<pre><code>P[28] &gt; 2^(2^28) = 2^(268435456) = 10^(log10(2)*2^28) ~ 10^80807124
</code></pre>

<p>That's at least <code>80807124</code> digits for the final value of <code>P</code>, which is the square root of the number you're looking for. That makes <code>Q[28]</code> larger than <code>10^1.6e8</code>. If you printed that number into a text file, it would take more than 150 megabytes.</p>

<p>If you imagine you're trying to handle these integers exactly, you'll see why it takes so long, and why you should reconsider your approach. What if you could compute that huge number? What would you do with it? How long would it take python to print that number on your screen? None of this is trivial, so I suggest that you try to solve your problem on paper, or find a way around it.</p>

<hr>

<p>Note that you can use a symbolic math package such as <code>sympy</code> in python to get a feeling of how hard your problem is:</p>

<pre><code>import sympy as sym
a,b,c,b0 = sym.symbols('a,b,c,b0')
a = 0
b = b0
for k in range(28):
    c = (a+b)**2
    a = b
    b = c
print(c)
</code></pre>

<p>This will take a while, but it will fill your screen with the explicit expression for <code>Q[k]</code> with only <code>b0</code> as parameter. You would "only" have to substitute your values into that monster to obtain the exact result. You could also try <code>sym.simplify</code> on the expression, but I couldn't wait for that to return anything meaningful.</p>

<hr>

<p>During lunch time I let your loop run, and it finished. The result has</p>

<pre><code>&gt;&gt;&gt; import math
&gt;&gt;&gt; print(math.log10(c))
49287457.71120789
</code></pre>

<p>So my lower bound for <code>k=28</code> is a bit large, probably due to off-by-one errors in the exponent. The memory needed to store this integer is</p>

<pre><code>&gt;&gt;&gt; import sys
&gt;&gt;&gt; sys.getsizeof(c)
21830612
</code></pre>

<p>that is roughly 20 MB.</p>

python loops recurrence