Implement Binary Exponentiation

Question

Background

In programming, there is a recursive algorithm called binary exponentiation, which allows for large integer powers to be calculated in a faster way. Given a non-zero base \$x\$ and a non-negative exponent \$n\$, the algorithm goes something like this (based on the example code from Wikipedia):

Function exp_by_squaring(x, n)
    if n = 0  then return  1;
    else if n is even  then return exp_by_squaring(x * x,  n / 2);
    else if n is odd  then return x * exp_by_squaring(x * x, (n - 1) / 2);

Basically, the code "reduces" the exponent term by first checking whether the current exponent is odd or even, then if it is even, just square root it; otherwise, divide by the base then square root it. Then repeat until the exponent reaches \$0\$. If the exponent is initially \$0\$, then just return \$1\$ directly.

This allows for an exponent term to be calculated faster than just multiplying the base by itself one at a time.

Example

Here is an example of the algorithm being applied on \$x=3,n=21\$.

1. \$n=21\$ is odd and non-zero, so we divide by the base, then square root. In this case, the number reduces to \$\sqrt{\frac{3^{21}}3}=3^{10}\$.
2. \$n=10\$ is even and non-zero, so we simply take the square root. \$\sqrt{3^{10}}=3^5\$.
3. \$n=5\$ is odd and non-zero, so we divide by the base, then square root. \$\sqrt{\frac{3^5}3}=3^2\$.
4. Continuing the process, we get \$3^1\$ then \$3^0\$, after which the recursion stops.

Task

Notice how at each step in the example above, we have a value which resulted from reducing the original number. Your task is to return a list of these numbers, given a non-zero integer base \$x\ne-1\$ and a non-negative integer exponent \$n\$. The list should always contain the initial value as the first element. The list can be returned in reverse order if you want to.

Test Cases

   x, n     -> Output
   3, 21    -> [10460353203, 59049, 243, 9, 3, 1]
1000, 0     -> [1]
   2, 15    -> [32768, 128, 8, 2, 1]
   1, 40000 -> [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  21, 3     -> [9261, 21, 1]
   2, 30    -> [1073741824, 32768, 128, 8, 2, 1]

---

This is code-golf, so shortest code in bytes wins!

Answer 1

Vyxal, 5 bytes

‡½⌊↔e

Try it Online!

Takes inputs in reverse order.

Jelly, 6 bytes

HḞ$Ƭ*@

Try it online!

Takes inputs in reverse order.

Explanation

HḞ$Ƭ*@  Dyadic link f(n, x)
   Ƭ    Repeatedly apply on n until the results are no longer unique
  $     Last two links as a monad:
H        Halve
 Ḟ       Floor
    *@  x to the power of each

Answer 2

Jelly, 5 bytes

:Ƭ2*@

Try it online!

Takes n on the left, and x on the right. Based on Steffan's answer, so be sure to give them an upvote as well.

How it works

:Ƭ2*@ - Main link. Takes n on the left and x on the right
 Ƭ    - Repeatedly apply until a duplicate value is found:
: 2   -   Floor divide n by 2
   *@ - Raise x to each power

Answer 3

05AB1E, 7 bytes

bη0šÙCm

Try it online!

b    convert the exponent to binary
η    takes all prefixes
0š   append the prefix "0", to include 1 in the list
Ù    uniquify, since if the exponent was 0 then "0" was already a prefix
C    convert each prefix from binary to a number
m    exponentiation

05AB1E, 7 bytes

·.Γ;ï}m

Try it online!

A port of the Vyxal and Jelly answer

Answer 4

JavaScript (Node.js), 30 bytes

x=>g=n=>[x**n,...n?g(n>>1):[]]

Try it online!

Answer 5

PARI/GP, 34 bytes

f(x,n)=[x^m\=2|i<-binary(m=2*n+1)]

Attempt This Online!

---

PARI/GP, 29 bytes

f(a,n)=if(n,f(a,n\2)*x+a^n,1)

Attempt This Online!

This returns a polynomial whose coefficients are the results.

Answer 6

C (GCC), 66 53 bytes

-13 bytes thanks to @att

r;f(x,n){printf("%d ",r=n?f(x,n/2),(n%2?x:1)*r*r:1);}

Attempt This Online!

Answer 7

Desmos, 42 bytes

f(x,n)=x^{floor(n/2^{[0...ceil(log_2n)]})}

Try it on Desmos!

Answer 8

Python 2, 37 bytes (-1 @Steffan, -2 @dingledooper)

f=lambda x,n:[1][n:]or[x**n]+f(x,n/2)

Attempt This Online!

Old Python, 40 bytes

f=lambda x,n:[x**n]+(n*[n]and f(x,n//2))

Attempt This Online!

Doesn't actually implement the algorithm. But produces correct output.

Answer 9

Java, 59 bytes

x->n->{for(n*=2;n>0;)System.out.println(Math.pow(x,n/=2));}

Try it online!

Answer 10

BQN, 16 bytes

⊣⋆{×𝕩?𝕩∾𝕊⌊𝕩÷2;0}

Try it at BQN REPL

  {          }    # Recursive function acting on right argument 𝕩 :
     ?            # if
   ×𝕩             # sign of 𝕩 is non-zero
      𝕩∾          # join 𝕩 to
        𝕊         # the result(s) of a recursive call to this function
                  # with argument:
         ⌊        # floor of
          𝕩÷2     # 𝕩 divided by 2.
                  # Now, get
⊣                 # the left argument
 ⋆                # raised to these powers

Answer 11

Prolog (SWI), 45 bytes

_*0*[1].
X*N*[P|R]:-P is X^N,T is N//2,X*T*R.

Try it online!

Answer 12

x86-64 machine code, 24 bytes

Or 1 less for 32-bit mode, with 1-byte dec ecx. If anyone has any other ideas for golfing further, please comment; I expect there's room for improvement. A straightforward binary exponentiation that doesn't have to produce the exact same proof-of-"work" array output is only 19 bytes.

Actually implements a binary-exponentiation algorithm, not doing a separate pow for every output element. Only one or two other answers do this, and I think this is the only iterative (non-recursive) one of those.

This is a function that takes an output pointer in RDI, to an array of DWORDs. x in ESI, n in EBX. (EDX is clobbered by mul. This can be signed if you want; mul and imul are the same binary operation on EAX, differing only in the EDX high half. They're the same size.)

Outputs: Array of steps in the memory the caller passed, RDI points to one-past-end. EAX = x**n.

trimmed NASM listing: address, machine-code (the answer), source

                   expint: ; uint32_t *out (rdi),  x = esi   n = ebx
 00 6A01               push 1
 02 58                 pop  rax        ;    mov eax,1
 03 AB                 stosd           ;  *rdi++ = eax
 04 0FBDCB             bsr  ecx, ebx   ; bit-scan-reverse.  we need to loop bits high to low; a set bit cares about (gets squared by) how many bits are *below* it, not above.
                                ; ECX = bit-index of highest set bit
 07 740E               jz   .n0        ; exponent=0 is a special case
                   .loop:
 09 F7E0               mul  eax     ; square total
 0B 0FB3CB             btr  ebx, ecx          ; bit-test-reset.  Or just BT, clearing n as we go wasn't helpful.
 0E 7302               jnc .was_even
 10 F7E6                mul  esi     ; total *= x;   (destroys EDX with the high-half result, but this saves a byte vs. imul eax, esi)
                   .was_even:
 12 AB                 stosd
 13 FFC9               dec  ecx
 15 7DF2               jge .loop     ; }while(--bitpos >= 0);  // loop runs for ECX=0, can't use LOOP unless we left-shifted EBX or something
                   .n0:
 17 C3                 ret
                   
size: 0x18 = 24 bytes

It feels like there should be room to golf more, but it's tricky to cover the corner cases (like n==0 is costing 2 bytes) and the off-by-ones (e.g. having to run the loop body for ECX=0 means I can't use loop).

Starting with 1 isn't as much of a burden as I'd thought it might be, mostly just requiring an extra stosd to output before the first multiply and after the last. We need to start with the MSB, not LSB, so we always have a 1 bit so the first iteration can do 1*1 * x to get us started with total=x. That does let us fall into the loop with total=1 after all. If not for having to store the n, though, I'd likely start with EAX=x and adjust the loop somehow.

(The test-case I asked for with a non-trivial even power turned out not to be relevant for this; I hadn't realized you needed to start from the high end for an iterative version. But it makes sense if you think about it: a 1 in the lowest bit of n results in a *x at the end. But if you do that at the start, a bunch of squaring happens so that odd power becomes an even power. If you just shr ebx, 1 / jnc to divide by 2 and check for odd/even, you'd effectively be bit-reversing the exponent, or something like that. Which works for exponents like 15 that have no significant zeros, but not others.)

Update: you don't have to start from the MSB, I just had a brain fart: the normal algorithm would do x *= x; every iteration, only multiplying it into the total for set bits in n. The algorithm is identical to multiplication by shifting-and-adding, in case it helps to think about how adding partial products works. (Stepanov's lectures on algorithms uses binary exponentiation as one of the main examples, pointing out that replacing + with * makes exponentiation instead of multiplication.) This question's presentation in terms of square-rooting the total in reverse got me stuck on squaring the total instead of x, and if you do that then you have to go from MSB downward.

No CPU I'm aware of has an integer sqrt instruction, or a single-instruction pow to give a starting point. (x87 has some exponential stuff, but takes multiple instructions to do pow. Still maybe worth considering.

Hopefully we can still recover all the same intermediate values this question wants, while doing shr ecx,1 / jnc. But that would make starting with 1 inconvenient.

---

BMI2 mulx is a 5-byte instruction, so even if we could reuse FLAGS to save a dec or test, it wouldn't pay for itself.

Try it online! with _start test caller (doesn't print: use GDB display *(int[8]*)($rsp+8) to display the output array state as you single-step, given that the caller passes RDI=RSP as the output.)

19 bytes, different array from a more normal algorithm:

expint_v2: ; uint32_t *out (rdi),  x = esi   n = ebx
    push 1
    pop  rax      ; mov eax,1
    stosd
.loop:              ;do{
    shr  ebx, 1  ; n /= 2,  CF = old low bit
    jnc .was_even
    mul  esi          ; total *= x;   (destroys EDX with the high-half result, but this saves a byte vs. imul eax, esi)
.was_even:
    stosd             ; *out++ = eax
    imul esi, esi     ; square x
    test ebx, ebx
    jnz  .loop      ;}while(n!=0);
    ret

Same inputs as before, same result in EAX and the highest array element (since the last bit shifted out before it becomes zero is a 1, unless n=0 in which case we leave the loop without having done mul esi.)

e.g. for 3^21 it produces these steps, where 1870418611 is 3^21 mod 2^32.

{1, 3, 3, 243, 243, 1870418611}

The x*=x steps where n doesn't have a set bit don't change the running product. (Like for multiply by addition, where you're not adding a partial product at that shift.)

I don't think it's possible to recover the same proof-of-work temporaries from this (even with a mov/imul into another temp), because we're working from the low bit up, not from the high bit down.

Starting with 1 is actually helpful, since we actually need to not start with total = x if n is even, but after we square x this iteration we no longer have any odd power of the original x left to multiply by, so the total would remain an odd power.

Starting with total=1 also makes the n=0 case Just Work for free.

Anyway, this order of temporaries is perhaps less interesting, but could probably also be recovered from a recursive implementation. I'm curious whether another challenge, perhaps titled iterative binary exponentiation, would get different answers in other languages. Probably still easier to use ** or ^ built-ins in languages that have them, after generating the right array of what, prefix-sum of bit-indices of non-zero bits or something?

Answer 13

Retina, 69 bytes

^
1,
\d+$
*
+`(.+),(.+),\b(_*)\3(_?)
$1¶$.($($.4*$2,1)***),$2,$3
,.+

Try it online! Link includes test cases. Outputs the list in ascending order. Explanation:

^
1,

Prefix a 1 which is the working value for the result.

\d+$
*

Convert n to unary.

(.+),(.+),\b(_*)\3(_?)

Match the previous result, x, n/2 and n%2.

$1¶$.($($.4*$2,1)***),$2,$3

Keep a copy of the previous result, then square it and multiply it by x or 1 depending on whether n was odd. (Retina only uses the first integer it finds in a string, so multiplying 3,1, 59049,3,_ and 59049,3,_ is equivalent to multiplying 3, 59049 and 59049.) Append x and n/2 for the next pass through the loop.

+`

Repeat the above until n=0.

,.+

Delete x.

Retina doesn't have exponentiation. Retina 1 does however have both eval and multiplication, so exponentiation can be (slowly) performed in 30 bytes:

~`(.+),(.+)
.+¶$$.($2*$($1$*)_

Try it online! Link includes examples.

Answer 14

Husk, 6 bytes

M^U¡÷2

Try it online!

Another port of Steffan's answers: upvote those.

   ¡      # Apply function repeatedly, collect values in infinite list:
          #  (using arg1 as starting value)
    ÷2    #  Integer-divide by 2
  U       # Longest prefix with all unique elements
M         # Map over this:
 ^        #  arg2 to the power of this

Answer 15

APL (Dyalog Unicode), 20 bytes ^SBCS

{⍵=0:1⋄(⍺*⍵),⍺∇⌊⍵÷2}

Try it on APLgolf!

-2 bytes thanks to Vadim Tukaev

Answer 16

Charcoal, 18 bytes

ＩＸθＥ⊕Ｌ↨η²↨²✂↨粦⁰ι

Try it online! Link is to verbose version of code. Outputs the list in ascending order (+1 byte to output in the order given in the question). Explanation:

  θ                 First input `x`
 Ｘ                  Vectorised raised to power
       η            Second input `n`
      ↨ ²           Converted to base 2
     Ｌ              Take the length
    ⊕               Incremented
   Ｅ                Map over implicit range
             η      Second input `n`
            ↨ ²     Converted to base 2
           ✂    ⁰   Sliced to length
                 ι  Current value
         ↨²         Converted from base 2
Ｉ                   Cast to string
                    Implicitly print

Although Charcoal's BaseString() function would convert 0 to "0", its Base() function returns an array rather than a string, and this array is empty when n=0, so it doesn't have the problem that, say, the 05AB1E answer has. Instead, it has the problem that CycleChop([], 0) throws a ZeroDivisionError, so I have to use Slice() instead.