Thread: Conceptual question

Hello everyone ,

I'm not really asking a coding question, but more of a concept question for my homework that is baffling me right now.

I will explain it below, but a direct link to my homework is here too if you want to skim it...

http://ranger.uta.edu/~weems/NOTES23.../lab1sum13.pdf

Basically we are given a random amount of cartesian coordinates and the first step is to sort these points via mergesort by the angle they have to the x-axis. My first idea was to change the cartesian coordinates into polar coordinates as it is suggested in the homework. The only problem is that my professor wants us to only use integers (no doubles, floats, etc.) in the code.

This is where I run into my problem of figuring out how to find their angles when dividing two integers can give rounding error.

Any ideas on how to find the angles would be greatly appreciated. No code necessary, just trying to figure out how to convert these accurately at the moment.

Thanks!

From the pdf file

Since this is inherently slow (e.g. computing the arctangent with atan2()) ...

On my system atan2() takes about 55 nano-seconds, I'm not sure if that qualifies as inherently slow, but then it's never advisable to argue with your professor.

Note the problem states that the value of the inputs will be limited to -32767 to +32767, and assuming your program runs on a 32 bit machine then integers can range from -2147483647 (- (2^(31) - 1)) to + 2147483647 (+ (2^(31) - 1)). See if you can figure out a way to take advantage of this to get more accuracy from your integer based math.

Thanks for the advice rcgldr. However, I forgot to mention that in class he said he forgot to take that off and does not want us to use arctangent or atan2(). I appreciate the idea though. If all else fails, I'll go see one of my T.A.'s tomorrow and see if I can get a better idea on it.

Originally Posted by mgracecar

My first idea was to change the cartesian coordinates into polar coordinates as it is suggested in the homework. The only problem is that my professor wants us to only use integers (no doubles, floats, etc.) in the code.

OK, so write down the math you would use to convert cartesian coordinates into polar coordinates, and see if you can figure out a way to do this with integers.

One issue is how accurate do the calculated angles have to be? Is rounding to the nearest integer degree acceptable? or perhaps to the nearest 1/1000th of pi?

Note the problem states that the value of the inputs will be limited to -32767 to +32767, and assuming your program runs on a 32 bit machine then integers can range from -2147483647 (- (2^(31) - 1)) to + 2147483647

Unless of course you're using a 16 bit compiler. What compiler are you using?

Jim

Originally Posted by jimblumberg

Unless of course you're using a 16 bit compiler. What compiler are you using?

Jim

I'm using Dev C personally. I believe the professor is going to be testing it on SSH.

And @rcgldr, thank you for that last post. It didn't even cross my mind on how accurate the degree needed to be. I will have to get this specified tomorrow with my teacher. For now, I'm going to work with the formula and see what algebra I can do. I will keep yall posted in the next few days for curiosities sake.

As always I appreciate the help and input,

Thanks

Originally Posted by mgracecar

the first step is to sort these points via mergesort by the angle they have to the x-axis.

If that were so, it'd been very cleverly phrased!

The trick is that

|tan(y/x)|

sorts in the same order as

|y| / sqrt(x*x + y*y)

within precision limits of course; and that can be squared while keeping the ordering the same, so the key would have been

(y*y) / (x*x + y*y)

To realize why this works, consider the (x,y) as vectors, and extend them to the unit circle (circle at origin, radius 1). The point where the (extended) vector intersects the unit circle is (x', y'):

x' = x / sqrt(x*x + y*y)
y' = y / sqrt(x*x + y*y)

Now, if you look at the angles (0 to 90 degrees), you can see that the angle to the x axis increases as the magnitude of the y coordinate on the unit circle increases.

---

However, the exercise states that the points are to be sorted based on their angle around the origin. This is the same as angle to the positive x axis, plus or minus some constant angle as the exercise does not state the direction of zero degrees.

Fear not. You can extend my initial suggestion further to cover this case too, it just takes some more effort.

Considering the signs of x and y, you have four quadrants: ++, +-, -+, and --.
Split each quadrant into two, one with (|x|>=|y|), the other with (|x|<|y|). Now, you have eight octants.

If you first determine the octant, you can use either (x*x)/(x*x+y*y) (if (|x|<|y)), or (y*y)/(x*x+y*y) (if (|x|>=|y)), to determine the fractional "angle" within the octant. The result is always between 0 and 0.5, and "angle" calculated from the y or x axis, respectively, so you just need to add or substract the "angle" from a value depending on the quadrant.

The overall angle obtained this way is exactly what you'd get if, if instead of a circle, the angles were numbered at constant intervals around an unit square:

Code:

  6/16  5/16  4/16  3/16  2/16     0.3750 0.3125 0.2500 0.1876 0.1250
  7/16                    1/16     0.4375                      0.0625
  8/16                    0/16     0.5000                      0.0000
  9/16                   15/16     0.5625                      0.9375
 10/16 11/16 12/16 13/16 14/16     0.6250 0.6875 0.7500 0.8125 0.8750

Because the order of any two angles is the same in both square angles (as pictured above two compasses) and normal angles, angles measured either way sort in the same order.

My suggestion is you use an unsigned integer type variable, perhaps unsigned long, for the angles. If you #include <limits.h>, UINT_MAX tells you the largest possible unsigned int, and ULONG_MAX the largest possible unsigned long. If you reserve a quarter of the range for each quadrant, you can use

((unsigned long)(y * y) * ULONG_MAX / (unsigned long)4) / ((unsigned long)(x * x) + (unsigned long)(y * y))
((unsigned long)(y * y) * ULONG_MAX / (unsigned long)4) / ((unsigned long)(x * x) + (unsigned long)(y * y))

to obtain the offset within each quadrant -- you just need to pick the correct one, the offset based on the quadrant, and whether to add or substract from the offset based on octant (half of the quadrant).

Questions?

But don't forget this requirement:

you must do all computations on integers without using division.

Highlight mine.

Jim

I only noticed the integer limitation, not the do-not-divide one.

There still is a simple solution: implement arctan2 using CORDIC. A simple web search on CORDIC and atan2 should provide a number of interesting results. Note that this version returns only the angle, not the distance, since the distance is not needed.

For example:

Code:

/* Array of atan2(1, 1<<i)*2147483648/(2*Pi).
 * This yields 2147483648 = 360 degrees.
*/
const int iangle[] = {
    268435456, /* 45.00000000000000000000 degrees */
    158466703, /* 26.56505117707799001892 degrees */
     83729454, /* 14.03624346792647692439 degrees */
     42502378, /* 7.12501634890179769144 degrees */
     21333666, /* 3.57633437499735151732 degrees */
     10677233, /* 1.78991060824606940116 degrees */
      5339919, /* 0.89517371021107439155 degrees */
      2670123, /* 0.44761417086055305115 degrees */
      1335082, /* 0.22381050036853808449 degrees */
       667543, /* 0.11190567706620689614 degrees */
       333772, /* 0.05595289189380366762 degrees */
       166886, /* 0.02797645261700367619 degrees */
        83443, /* 0.01398822714226501639 degrees */
        41722, /* 0.00699411367535291827 degrees */
        20861, /* 0.00349705685070401126 degrees */
        10430, /* 0.00174852842698044954 degrees */
         5215, /* 0.00087426421369378026 degrees */
         2608, /* 0.00043713210687233457 degrees */
         1304, /* 0.00021856605343934784 degrees */
          652, /* 0.00010928302672007150 degrees */
          326, /* 0.00005464151336008544 degrees */
          163, /* 0.00002732075668004893 degrees */
           81, /* 0.00001366037834002524 degrees */
           41, /* 0.00000683018917001272 degrees */
           20, /* 0.00000341509458500637 degrees */
           10, /* 0.00000170754729250319 degrees */
            5, /* 0.00000085377364625159 degrees */
            3, /* 0.00000042688682312580 degrees */
            1, /* 0.00000021344341156290 degrees */
            1, /* 0.00000010672170578145 degrees */
            0
};

static int iatan2(int y, int x)
{
    size_t           i = 0;
    int              angle, direction;

    /* Which half-plane? */
    if (y < 0) {
        x = -x;
        y = -y;
        angle = -1073741824;
    } else
        angle = 0;

    /* Which quadrant? */
    if (x < 0) {
        x = -x;
        angle += 1073741824;
        direction = -1;
    } else
        direction = +1;

    /* Which octant? */
    if (y > x) {
        const int oldy = y;
        y = x;
        x = oldy;
        angle += direction * 536870912;
        direction = -direction;
    }

    /* Along axis? */
    if (y == 0)
        return angle;

    /* Diagonal? */
    if (y == x)
        return angle + direction * 268435456;

    /* Scale, for optimum results. x > y > 0.
     * If coordinates are < 32768 in magnitude,
     * you could just multiply by 32768.
    */
    while (x < 268435456) {
        x *= 2;
        y *= 2;
    }

    /* Adjust direction. */
    if (direction < 0)
        y = -y;

    /* CORDIC for angle, 0 < y < x. */
    while (iangle[i]) {
        const int oldx = x;
        const int oldy = y;

        if (y > 0) {
            x += (oldy >> i);
            y -= (oldx >> i);
            angle += iangle[i++];
        } else
        if (y < 0) {
            x -= (oldy >> i);
            y += (oldx >> i);
            angle -= iangle[i++];
        } else
            break;
    }

    return angle;
}

iatan2(y,x) yields angles [-536870912, 536870912], analogously to atan2(y,x) returning angles [-Pi, Pi], so each unit is 360/1073741824 degrees; one degree is about 2982616 units.

For x=-32768..32768, y=-32768..32768, the absolute error in angles is -13 to 13 units; i.e. the error of iatan2(y,x) is less than 0.00000436 degrees.

By the way, this is NOT optimized code. This, too, is in public domain, so you can do whatever you want with it.

Originally Posted by rcgldr

Using arctan2 (via a table) isn't required for this program. [...] (hint: think slope)

Slope requires division, which is forbidden. If you sort by e.g. octant (say in bits 16 to 18), and max(|x|, |y|) in bits 0 to 15, binary search will fail to find a point in the same direction but further from the origin.

Ah, do you mean a loop using bit shifts to scale (|x|,|y|) so that max(|x|, |y|) == 32768; in which case you can use the quadrant, and the coordinate which does not have absolute magnitude 32768, to get the "square angle" coordinates?

Ahhh excellent points by all of you. My T.A. hasn't gotten back to me and I have been at work all day so I am just now getting back to this. I will continue to work on it tonight. I see what you guys are saying. Let me see if I can write down/mess around with some formulas and draw up some code. Thank you for everyone's help so far. It's really appreciated.