Calculating large powers of a double with modulus

I want to find large powers of (1+sqrt(3))^(n+1) where n varies from n=1 to n=1000,000,000. Can modular exponentiation by squaring work with doubles? How? I have searched a lot but no positive results yet.Any help will be highly appreciated?

The answer to the question “Can modular exponentiation by squaring work with doubles?” is no, given some assumptions I will explain.

Operations such as taking a residue modulo 1,000,000,007 are typical of various number theory calculations and other integer arithmetic. So, I will assume we want to calculate the result with integer precision. That is, we want to know the integer closest to (1+sqrt(3))^(n+1) % 1,000,000,007 for some n, 0 < n ≤ 10⁹. Consider a function xⁿ⁺¹. Its derivative with respect to x is (n+1)xⁿ. When x is 1+sqrt(3) and n is 10⁹, the derivative is about 4.65•10^436488780. This means a change in x of about 1 part in 4.65•10^436488780 will change the value of xⁿ⁺¹ by about 1, which of course also changes the residue by about 1. Therefore, to calculate the residue with the desired accuracy, we must, in effect, calculate the value of 1+sqrt(3) to over 436 million decimal places.

This is grossly beyond the precision that is provided by the double type on any modern computer. Therefore, the answer is no; no algorithm can give you the desired answer because the double type is not capable of calculating with the necessary precision.

Technically, one could use double-precision arithmetic to construct extended-precision arithmetic. (Thus, the answer above is addressed at using double-precision arithmetic in normal ways.) So, using extraordinary techniques, one can fabricate the necessary computations from double-precision operations. Discussion of the methods for doing this are beyond a Stack Overflow question. Commonly, integer operations are used for large portions of such work, although double-precision operations may play a role. And software packages are available to provide such extended-precision arithmetic, such as GMP (GNU Multiple Precision Arithmetic Library).

Calculating and retaining 436 million decimal places may strain computing resources. It may be possible to reduce the computations needed by discarding information that can be shown not to affect the final residue, so that not all 436 million decimal places are needed at the same time. However, I do not expect that such reductions would produce a computation that is feasible in double-precision arithmetic without extended-precision techniques.

And the problem is worse than that. Knowing x to about one part in 4.65•10^436488780 might help you calculate the result within an error of 1, but it will not tell you which integer is nearest. Consider that you might calculate (1+sqrt(3))ⁿ⁺¹ % 1,000,000,007 and find the result is very near 3.5. In order to determine whether the integer closest to the result is 3 or 4, you must figure out which side of 3.5 the result is on. This function, (1+sqrt(3))ⁿ⁺¹ might have some values that are extraordinarily close to such a halfway position, and you will need to calculate them even more accurately to make a proper decision.