FFT Slide Rule

[formerdeathcorps]

The PSX is built on an integer processor.  For those of you not versed in Computer Science, this means every value stored in the game is an negative or positive whole number.  And yet, as secondadvent and I found out, the archer's arcing attacks requires the calculation of a parabola's trajectory to determine whether a target in range can actually be hit.  This requires FFT to solve the quadratic formula, leaving us both wondering how FFT calculates the square root.  After some digging by SA (http://ffhacktics.com/paste/7hJtayK4), we found that Square coded an entire functioning slide rule / logarithm table into FFT starting at SCUS location 0x1c0b8 (0x2b8b8 in RAM).  (For those of you who are math illiterate, PM me and I'll explain more.)  Furthermore, previous hex values before the given table suggests that Square also hard-coded similar tables for evaluating trig functions (necessary for map and unit movement graphics).
The usual Squarecode wastage aside, Square's table is inaccurate for large inputs.  Small numbers like sqrt(2) or sqrt(6) have errors of no more than 1 / 5000, but a 31-bit signed word I picked at random has an error of 22%, or a 1000x increase in error.

Seeing this mess, I decided to fix it.  My solution (http://ffhacktics.com/paste/QSTmeHZ3) will take a bit more runtime, but should be perfectly accurate (to the nearest integer) and cuts out the need for tables altogether.  Given 32-bit inputs, it need not iterate more than 5 times.

[Pride]

This is... really interesting and awesome = o

[formerdeathcorps]

My method has its own rounding error; if a positive number is of the form n(n+1), it will always round up instead of down.  In other words, my routine yields:
Round(sqrt(2)) = 2 instead of 1
Round(sqrt(6)) = 3 instead of 2
Round(sqrt(12)) = 4 instead of 3
...
but as N approaches infinity, the % error drops to 0%, so I'm not too concerned.

[Choto]

That's all pretty crazy! I'm wondering what aspects are covered by this aside from arrow trajectories. Have you tried playing with your version and noticed anything significant?

[Glain]

Isn't using the nearest integer already inaccurate though?  Shouldn't the square roots be fixed point numbers?

[formerdeathcorps]

EDIT: Sorry for misreading your post.  I read floating point at first.

Square effectively did that in their original routine by outputting outputting sqrt{r4} * 64.  If you're concerned about accuracy, just multiply your inputs by 4096 (sll r4, r4, 0x0A) to get the same result Square's original routine did.  That will also save on the tedium of deleting sra rY, rX, 0x06 from every original function call involving the square root.
Fixed point is nice, but the necessary adjusting of the decimal point erases any space gain I get over Square's original routine.  Nor do I really appreciate having an input range smaller than vanilla FFT.  Since we are human and not computer ASM compilers and thus understand how the same mathematical output from the same routine can have different meanings due implied context, I'd rather the coders using my routine have the flexibility of creating an implied fixed point by multiplying inputs by an appropriate factor for accuracy BUT also retain the option of running it as a full 32-bit value as well.  (This is why I didn't hard-code the multiplier of 4096 into my routine.)

As for accuracy:
Nearest integer has a maximum error of ~42.4% for the sqrt{2}.
sqrt{6} has an error of ~22%.
sqrt{12} has an error of ~15%.
sqrt{20} has an error of ~11.8%.
....
IN general, the maximum error function is bounded by 0.6 / sqrt{N^2 + N}.  When N^2 + N ~ 100, we have a 6% error, ~1000 yields 2%, and ~5000 (assuming you all take my advice and multiply inputs by 4096) yields < 1% error.  I think that's fine.

Also, I misstated Square's accuracy.  Their table is indeed very accurate and if Square had used floating point (or fixed point), then their error for sqrt{X} for small X would be < 1/5000, BUT SINCE the final output is [sqrt{r4} * 64], their actual error is around 1 / 200, which is exactly the same as my routine.
E.g.
A scientific calculator that performs FFT's routine for the sqrt{2} would yield 90.5, but that rounds down to 90 in FFT.  90.5 / 64 has an error of 1 / 6000 but 90 / 64 has an error of 1 / 200.

[Glain]

Oh, so the reason it's sqrt * 64 is so that it can essentially return fractions in the form of n/64.  I see.

I'm kinda worried about sqrt(2).  Returning 1 or 2 for that, if it's used as a multiplier somewhere, is going to mess calculations up.  (Same to a lesser extent for other small square roots).

Edit: I don't understand why multiplying inputs by 4096 would work.  Wouldn't that just multiply the (already inaccurate) result by 64, instead of keeping precision? (Also wouldn't you want to shift by 12, not 10? 2^10 = 1024)

[formerdeathcorps]

Glain, the input has no error, so multiplying it by 4096 has no effect on ultimate accuracy.  However, multiplying the input by 4096 yields [sqrt{r4} * 64], not 64 * [sqrt{r4}].  When r4 is small, the difference between these two is not negligible.

For example:
If your input is 2, you will get sqrt{2} = 2 as an answer.  Error = 42%
If your input is 2 * 4096 = 8192, you will get sqrt{2} * 64 = 91 (not 64 * 2 = 128) as an output.  Error = 0.5%

[Glain]

Oh... I believe I understand.  For this approximate square root function, increasing the input means decreasing the error, so for the most accurate result, multiply the argument by the highest number you can without causing it to overflow that is both a power of 2 and a perfect square, then divide by its square root at the end? (Making use of the fact that mathematically, sqrt(x * y) = sqrt(x) * sqrt(y))

approx_sqrt(2) = 2

approx_sqrt(2 * 16) = approx_sqrt(32) = 6
6 / 4 = 1.5

approx_sqrt(2 * 4096) = approx_sqrt(8192) = 91
91 / 64 = ~1.422

approx_sqrt(2 * 65536) = approx_sqrt(131072) = 363
363 / 256 = ~1.418

(But [approx_sqrt(r4) * 64] and 64 * [approx_sqrt(r4)] are totally the same thing, unless you're suggesting multiplication isn't commutative.  Perhaps what we mean to say is that approx_sqrt(r4 * 64^2) isn't the same as approx_sqrt(r4) * 64, because of the error of approx_sqrt() being decreased with large inputs?)

[formerdeathcorps]

But [approx_sqrt(r4) * 64] and 64 * [approx_sqrt(r4)] are totally the same thing, unless you're suggesting multiplication isn't commutative.  Perhaps what we mean to say is that approx_sqrt(r4 * 64^2) isn't the same as approx_sqrt(r4) * 64, because of the error of approx_sqrt() being decreased with large inputs?

Yes, that's what I meant.  Since [] = Greatest Integer Function, I meant that mathematically, [sqrt{2}] * 64 =/= [sqrt{2} * 64], which is equivalent to saying that approx_sqrt(r4 * 64^2) =/= approx_sqrt(r4) * 64, for the same reasons.

[Glain]

Ah, I see now!  Yeah, that makes sense. 

I've never seen brackets used to denote the ceiling function before; typically, in my experience, brackets are interchangable with parentheses. The symbols for the ceiling function are like brackets but with the bottom parts removed... much like for the floor function, the top parts are removed. (If you use brackets for ceiling(), how do you denote floor()?)  I'm not sure there's a good way to represent those with computers so I think I've only ever seen it as ceiling() or ceil().