How can I rotate a 16-bit signed integer vector?

Question

I am trying to rotate a vector made of two signed 16-bit coordinates (-32768 to 32767). Using a LUT for sin/cos, I can use the standard method for rotating a vector:

newX = x * cos(theta) - y * sin(theta);
newY = x * sin(theta) + y * cos(theta);

However, since this requires multiplying the coordinates by sin/cos (in my case, ranging from -127 to 127), large vectors can easily overflow in the middle of the process before being scaled back down. This means either the length of the vector or the accuracy of the sin/cos is severely limited, despite the original and resulting rotated vector being safely in range.

Is there another method for rotating integer vectors?

Edit: The numbers cannot be cast to floats or higher-bit types.

Answer 1

Use 32-bit integers for the intermediate calculations.

They'll usually be divided down into the 16-bit range at the end. If they're not - for example, rotating \$(30000,30000)\$ by an eighth-turn should give \$(0,42426)\$ which doesn't fit into 16 bits - check it and produce an error at this point.

Answer 2

There are two issues to address here, each of which has a couple of solutions.

Size of vectors

Your first problem is how you do any processing at all without overflows. A vector (32767, 0) is clearly going to rotate nicely to anywhere. However a vector (32767, 32767) is going to become (49339, 0) when it's rotated 45 degrees, and that'll overflow. Even (-32768, 0) is going to overflow if you rotate it by 180 degrees.

Limit vector magnitude

You might solve this by having something limit the vectors entered so that their magnitude is small enough to fit. If data is being entered from a PC, you can use that to only give you valid data. Calculating a magnitude needs a square root though, which is going to be tricky on a constrained device.

Or only store 15-bit values ####8

If these checks are too much of a problem, you could just limit the range of values for vector entry to signed 15-bit. Processing can then take place in signed 16-bit arithmetic, knowing that the result cannot overflow. This would be my preferred solution, not just for simplicity, but because you have engineered out any possibility of failure.

Calculation of rotation

And then we come to how you do the rotation calculation without overflow. The sin/cos clearly is a fraction. I would strongly recommend having an extra bit in the sin/cos lookup value so that +1 can be exactly represented (e.g. +/-128 in a 9-bit value; or +/-64 if constrained to 8-bit). But from there, how to avoid overflow?

Use your processor's multiply instruction registers

On pretty much all processors with a multiply instruction, you have two registers in and two registers out. If you start with two 16-bit values, you get one 32-bit value held in two registers. If you're using +/-64 to represent +/-1, you'll need to divide the result by 64 - which conveniently just means discarding the lower 7 bits. With some simple bit-shifting and masking, it's easy to get a single word containing the result you want.

Break the multiply into lots of conditional adds

If you've got a sin/cos value of 64 (representing 1.0) then the input value passes through unchanged. Each lower bit after that represents a fraction of the input - 32 represents a half, 16 represents a quarter, and so on. You can add all these fractions together, knowing that the result can never overflow. Of course this takes a bunch of operations, but you're optimising for storage over speed. (You might also get some rounding errors in the lsbs, but again, this is prioritising storage over all else.)

But more generally...?

I'd have to wonder about where these limitations came from. If you're on an FPGA, it's unlikely you're so constrained by storage as to be worried by this. And if you're using a micro, it's rare that you'd be this limited these days.

Answer 3

If you don't have a way to get the high half of a 16bit multiplication then you can instead scale the vector back before multiplying. (s and c are your byte-sized sin and cos)

x = (x>>8) * c - (y>>8) * s;
y = (x>>8) * s + (y>>8) * c;

That will end up canceling out the upscale of the sin and cos.

If you don't mind a bit more operations per rotation you can also compute the bit you cut off:

x = (x>>8 * c  + ((x<0?-1:1)*(x&x0ff)*c)>>8 ) - (y>>8*s + ((y<0?-1:1)*(y&0xff)*s)>>8 );
y = (x>>8 * s  + ((x<0?-1:1)*(x&x0ff)*s)>>8 ) + (y>>8*c + ((y<0?-1:1)*(y&0xff)*c)>>8 );

Answer 4

I think the trick here is to use the fact that the hardware is usually more capable than the compiler exposes. In particular, when you multiply a 16 bit value by a 16 bit value, you will get a 32 bit value, but the compiler will immediately throw away the upper half (unless you are looking at its internal library to do 32 bit math).

So what you want to do is get your sin(),cos() values in the -32768 ... 32767 range (by multiplying by 32768 and special casing the overflow), and then do the signed multiply, left shift the resultant register pair one bit left (typically rotate-through-carry-left or RCL), and then return the upper half of the result. You might also benefit from using the unsigned multiply operations, but then you need to track the signs yourself.

Answer 5

So...

• you do not want to use 32 bits per integer;
• you do not want to use float (or double);
• you do not want to lose precision.

That cannot be done for the full range of numbers, but the precision can be dealt with, respecting the first 2 conditions 100%.

---

The trick is that:

• x and y are both max 15 bits + 1 bit sign;
• sine and cos are both max 7 bits + 1 bit sign;
• the result is max 15 bits + 1 bit sign.

To speed things up a bit, you need to keep a parallel table with the LUT for sin and cos, to know the number of relevant bits of each value (leading zeros are irrelevant).

And somehow you find the number of relevant bits for x and y (separately).

An now the magic: If the result has 15 bits (max), and the s / c has k bits, and x / y has q bits, you need to scale the number^(*) by (k+q-15) if (k+q) is bigger than 15, perform the multiplication, and then reverse scale the result (if needed / desired).

I hope that my explanation makes sense.

^(*) I recommend to scale the bigger number, to keep the precision of the smaller number. Just compare e.g., x to sin, and scale whichever is bigger. Otherwise you risk to transform a non-zero number to a "very" zero, which is more undesirable than some loss of precision.

---

It might be the case that the resulting coordinate will be out of screen / out of visible area. In this case, just by counting bits you might find out if the resulting coordinates are worth calculating (with any precision) or not at all.

---

Disclaimer: you will lose some speed, though... It is in your responsibility to fond the sweet spot of the optimization. Trial and error, I guess...

Answer 6

First of all, as Kromster commented, you need to use immediate variables in order to not accidentally use the new x value when calculating the new y value, unless you're actually using a vectorized version. Second, when you say your sin/cos is encoded in the range -127..127 (shouldn't one of them be 128?) your formula needs to include the rescaling and should not call the encoded values sin/cos since those must be the proper floats in the -1..1 range:

x_new = (x*ct - y*st)/TRIG_SCALE
y_new = (x*st + y*ct)/TRIG_SCALE

where TRIG_SCALE=127, ct = TRIG_SCALE*cos(theta) and st=TRIG_SCALE*sin(theta).

Next, you can split x and y into to 8-bit parts each and their sign, such that x = sign_x *(x_high * 256 + x_low) (and analogously for y). Then

x_new = (sign_x*(x_high*256 + x_low)*ct - sign_y*(y_high*256 + y_low)*st) / TRIG_SCALE
      =   (sign_x*sign_ct*x_high*abs(ct) - sign_y*sign_st*y_high*abs(st))*256/TRIG_SCALE
        + (sign_x*sign_ct*x_low*abs(ct)  - sign_y*sign_st*y_low*abs(st))/TRIG_SCALE

and similar for y_new. Both parentheses in the second and third line are limited to 14/15 bit plus sign (7/8 bit from high/low bits plus 7 bit from abs(cos/sin)) so they can be handled individually via two 16 bit variables per coordinate which can then be individually scaled down by TRIG_SCALE before the addition. Note that since 256/TRIG_SCALE=2 the high part actually just gets doubled, thus yielding a 15 bit number as well. For the low part you'll loose precision due to the integer division, but that's to be expected here. Due to the mathematical nature of sin/cos those two 15 bit numbers shouldn't overflow when added, but I'll leave proving that as an exercise to the reader ;)