Writing a python function for raising irrational numbers to higher powers

I want to write a function in python called compute(a,b,p,n) where a,b,n are integers and n>0 and p is prime, that outputs the tuple (x,y) (x,y being integers) such that

Note that the constraints of the problem state that you’re not allowed to use any external libraries.

Now, I set about to solve this problem by writing a recursive function. To that end, let us assume that

But, we also have

Therefore, we have

Equating coefficients now yields

With

Now, making use of the equations in (1) and (2), we can write a simple and straightforward recursive python function as

def compute(a,b,p,n):
    
     if n==1:
        return  (a,b)

     if n>1:
        return (a*compute(a,b,p,n-1)[0]+b*p*compute(a,b,p,n-1)[1], b*compute(a,b,p,n-1)[0]+a*compute(a,b,p,n-1)[1])

This code works fine for small values of n, say upto 10, but for large values of n it is horribly inefficient. Is there some alternate way to approach this problem which works faster, or perhaps some way to optimize this code so that it runs faster?

Answer

The biggest improvements in speed often come from choosing a different algorithm. That’s the approach I’m taking here. Your current implementation works by reducing the exponent n by 1 on each iteration, which is (at best) linear in n. If you can find a way to cut n in half on each iteration, the complexity is reduced to O(log n). This is not only asymptotically faster, it greatly reduces the recursive stack requirements—you can handle much much larger exponents without causing stack overflows.

You can achieve a repeated halving solution for raising base B to the n^th power by leveraging the mathematical identities Bⁿ = (B²)^n/2 for n even, or B*(B²)^n/2 for n odd. As stated above, this results in an algorithm with O(log n) computational complexity, yielding a significant improvement for large values of n.

For your particular problem, B = (a + b*sqrt(p)). B² evaluates to (a²+b²*p + 2ab*sqrt(p)), while multiplying (a + b*sqrt(p))(c + d*sqrt(p)) yields ((ac+bdp) + (ad+bc)sqrt(p)). Using these results for the n even and n odd cases described above leads the following implementation in python:

def compute(a,b,p,n):
    if n == 1: return (a,b)
    c, d = compute(a*a + b*b*p, 2*a*b, p, n // 2)
    if n % 2 == 0:     # n is even
        return (c, d)
    else:              # n is odd
        return (a*c + b*d*p, a*d + b*c)

In coding this I noticed that your base case is incorrect. The problem specification says that n is strictly greater than 0, so the base case should be (a,b) for n == 1 since (a + b*sqrt(p))¹ = (a + b*sqrt(p)). If you want to include zero as a possible power, it would become (1,0) for n == 0. You can confirm that this is correct by making that change to your code and comparing the results to a hand-calculation of compute(1, 2, 9, 2), i.e., a simple square of (1 + 2sqrt(9))² = 49 should be (37 + 4sqrt(9)) = 49. Your code currently produces (19 + 3sqrt(9)) = 28 but yields the correct answer when the base case is corrected as above.