Struggling to find loop invariant in power function

I am struggling to find a good loop invariant for the following function, which returns a^b where a is a real number and b is a natural number:

power <- function(a, b){
     c <- 1
     while(b > 0){
          if(b %% 2 == 1){
               c <- c * a
          }
          b <- floor(b / 2)
          a <- a * a
     }
     return c
}

I've ran through the loop with a couple of examples, and I see that it has 2 kinds of cases; when b is even or odd. I also understand that on the kth iteration, a = a_0^(2^k), but I am struggling to find a proper invariant as there is no real iterating variable to use.

For the invariant to be useful, it will have to have c = a_0^b_0 as a special case after the while loop terminates, which occurs when b = 0.

For the invariant to be true, we have to get a_0^b_0 on the left hand side before the first iteration. We already know the left hand side has a c in it, and before the first iteration c = 1, so multiplication seems like a good idea.

Whatever we multiply by must end up being 1 after the loop terminates, which (as before) occurs when b = 0. Getting something to equal 1 when b is 0 suggests we want b to be an exponent, and the desired a_0^b_0 on the left hand side also suggests b should be an exponent.

Putting this all together, the invariant will either be c * a_0^b = a_0^b_0 or c * a^b = a_0^b_0. I will leave it to you to determine which one of those is correct, and to prove that it is in fact invariant.