Since your function is actually $f(k) = (2^{2k-1} + 2^{k-1})\pmod{10^9+7}$, it suffices to tabulate values of $2^k$.
$10^9+7$ is prime, so by Fermat's theorem, we know that $2^{10^9+6}\equiv 1\pmod{10^9+7}$.
So suppose you want to calculate $f(k)$. Calculating powers is expensive, so we will calculate $z=2^{k-1}\pmod{10^9+7}$ and then $f(k) = z\cdot(2z+1) \pmod{10^9+7}$.
First calculate $k' = k-1\pmod{10^9+6}$. Then calculate $2^{k-1} \equiv 2^{k'}\pmod{10^9+7}$ using the usual fast recursive exponentiation algorithm:
pow2(k): result = 1 power = 1 while k > 0: power = (power * 2) mod 1000000007 if k is odd, result = (result * power) mod 1000000007 k = int(k/2) return result
(In C, use k >>= 1
instead of k=int(k/2)
, power <<= 1
, etc.)
We have $k < 10^9 + 7$, so this will iterate at most about 30 times. If this is too slow, you can cache the results for $0 ≤ k ≤ 10^6$ and then it will iterate at most ten times per call.
Let $z = \mathrm{pow2}(k-1\pmod{10^9+6})$, and then your answer is $(2z+1)\cdot z\pmod{(10^9+6}$.
Wolfram|Alpha appears to have gone temporarily insane.