I recently found out that ^ only accepts integers. For example, when I did 9^2 it worked perfectly. However, when I tried 9^.5, it acted like I typed 9^0. Is this a bug or intentional?

Hmm, yeah, seems it only works with whole numbers period. For example, 9^1 works (says 9), 9^0.9 comes out as 0 (should be ~7.22), and 9^1.5 comes out as 9 (should be 27). It's basically flooring whatever value you use as the exponent. Definitely a bug.

it could have been designed that way or left for further implementation. not necessarily a bug.

i've spent several minutes researching exponentiation algorithms but it looks difficult with the limited range of math operations pico-8 provides. might be best to ask zep what the plans are.

It is correct that ^ only accepts integer exponents, and I'm fairly sure that this is the intended behavior in Pico-8. Check out this thread for a pow() function that takes decimal values: https://www.lexaloffle.com/bbs/?pid=30791#p30791

@helado: samhocevar's routine does bitwise manipulation of the exponent and accumulated n-squared or n-square-rooted versions of the base for each 1 bit, first in the whole part of the exponent, then in the fractional part of the exponent.

function pow(x,a) |

x is the base, a is the exponent.

if (a==0) return 1 if (a<0) x,a=1/x,-a |

If the exponent is 0 then the result is 1 regardless of the base. If the exponent is negative, this is equivalent to the reciprocal of the base to the absolute value of the exponent. Using a positive exponent simplifies the remaining code.

local ret,a0,xn=1,flr(a),x a-=a0 |

Initialize the return value (ret) to 1, a0 to the whole part of the exponent, and xn to the base. Remove the whole part of the exponent from a, leaving only the fractional part. Here and in a couple of other places, samhocevar is using "x,y,z=1,2,3" assignment syntax to save tokens.

while a0>=1 do if (a0%2>=1) ret*=xn xn,a0=xn*xn,shr(a0,1) end |

Iterate over the bits of the whole part of the exponent. xn tracks the value to be applied to the return value for each 1 bit. xn is squared for each bit position, but is only multiplied to the return value for each 1 bit. a0 % 2 >= 1 detects a 1 bit in the rightmost bit. a0 = shr(a0,1) slides the bits off to the right. The loop terminates when there are no more 1 bits in the exponent.

while a>0 do while a<1 do x,a=sqrt(x),a+a end ret,a=ret*x,a-1 end |

Iterate over the bits of the fractional part of the exponent. This is the same general idea as the previous loop, just using the square root to access each base value for fractional powers of 2 in the exponent, and using a+a to double the exponent at each step. The structure of the loop is a little different: the inner while loop skips all of the 0 bits, and the second line is only executed for the 1 bits.

return ret end |

Return the accumulated result.

@dddaaannn: thanks for the overview, but my question was how it was come up with in the first place, the work that went into either designing the algorithm or implementing it from an already known algorithm.

I can't help you there. It seems like a straightforward algorithm given the requirements of supporting both negative and fractional exponents with a fixed point number representation.

One nice technique used here is to take advantage of the binary representation to keep the code fast and simple. Hacker's Delight is a classic text with fun bit bashing algorithms in it, though I don't see this specific algorithm in my copy.

@dddaaannn: thanks for reminding me of hacker's delight. i checked my copy and while i also didn't find an explicit power function, i found logarithm. unfortunately, it's integer logarithm, which… i'm not so sure is helpful XD in any case my attempts failed so far, but i might fiddle with it some more later to see what i can derive.

@helado: dddaaannn gives a great breakdown of the algorithm. If you're looking for the derivation, I think it's a straightforward bit of algebra and can be understood in terms of the exponential identities, given that:

x^(y+z) = x^y * x^z

which means that you can separate the fractional component of the exponent from the whole number and compute both powers separately.

e.g. to compute x^4.3, let y=4 and z=0.3 so that

x^4.3 = x^4 * x^0.3

Obviously computing x^y where y is a whole number is trivial—we could do it with a simple "x^y" statement in Pico-8—but what you see here is an efficient implementation that groups together squares, and therefore beats the performance of a naive loop that executes "x*x*x*x*..." by doing it in ~log2(y) steps instead of y.

The tricky part is computing x^z where 0<z<1. We know that x^(1/n) = n-root(x), but in Pico-8 we only have access to a fast 2-root, or x^(1/2)—aka sqrt(). So we go back to the exponential identities:

(x^p)^q = x^(p*q)
x^p * x^q = x^(p+q)

So the clever trick here is use repeated square roots to find a close approximation to x^z such that

((x^(1/2))^(1/2))^(1/2)...) * ((x^(1/2))^(1/2))^(1/2)...) * ... =~ x^z

We do this by finding integers a, b, c... such that

0.5^a + 0.5^b + 0.5^c + ... =~ z if 0<z<1, so that

x^(0.5^a + 0.5^b + 0.5^c + ...) =~ x^z

But a,b,c... are easy to find, because z is already stored as a fractional binary number, which can be defined as the sum i*2^-1 + j*2^-2 + k*2^-3..., or i*(0.5) + j*(0.5^2) + k*(0.5^3)..., where i,j,k... are the binary digits 0 or 1. That means that a,b,c... are just the exponents of the places where binary digits are non-zero.

So, in short, the algorithm computes:

x^(y+z) = x^y * x^(0.5^a) * x^(0.5^b) * x^(0.5^c) * ...

where y is a whole number, 0<z<1, and a,b,c... are the exponents of the non-zero digits of the binary representation of z. Note that this implementation of the algorithm takes advantage of the fact that a<b<c<... and accumulates previous values of x^(0.5^...) to avoid unnecessary calls to sqrt().

Hope that helps. (also I did this quickly—so someone please fix my math if I've gotten anything wrong.)

[Please log in to post a comment]