CS494 Lecture Notes - The TrianglesContainOrigin Topcoder problem

Acos.cpp and acos_2.cpp

Before getting into the topcoder problem, I want to show you two programs that should make you feel uneasy about roundoff error and trigometric implementations on various computers.

The first is src/acos.cpp, and I'll just show the header comment so you can see what it does:

/* This program demonstrates fun with roundoff error.  What I do is 
   use some general code to calculate the angle of point (x,x) for x = 1, 2, 3, ..., 9.

   You'll note that this angle should be the same.  When I calculate it in a reasonable
   way, the angle differs from point to point.

   I'm calculating it by calculating the hypotenuse of the triangle (0,0), (x,0), (x,x).
   I then divide x by this hypotenuse and call acos().  As you can see below, these
   values do not match for successive values of x.
 */

When we run it, you'll see the angles don't match, even though one feels like they should:

UNIX> bin/acos
Point: (1,1).  Angle: 0.785398
Point: (2,2).  Angle: 0.785398.  Equal to last angle? true
Point: (3,3).  Angle: 0.785398.  Equal to last angle? false
Point: (4,4).  Angle: 0.785398.  Equal to last angle? false
Point: (5,5).  Angle: 0.785398.  Equal to last angle? true
Point: (6,6).  Angle: 0.785398.  Equal to last angle? false
Point: (7,7).  Angle: 0.785398.  Equal to last angle? true
Point: (8,8).  Angle: 0.785398.  Equal to last angle? false
Point: (9,9).  Angle: 0.785398.  Equal to last angle? false
UNIX> 
The second program, src/acos_2.cpp, generates 10,000 random values betwen -1 and 1, and calls acos() on them. I have run this on my mac, a hydra, an old laptop from my neighbor Larry, a pi5 and a pi3. What you see is that the mac gets different values. To show that it's the acos(), and not the random number generation, the program drand48() just shows the random values:
UNIX> openssl md5 txt/acos* txt/drand*
MD5(txt/acos_2_hydra.txt)= b269e01f2c4732b77eaf811e52fb2b36
MD5(txt/acos_2_larry.txt)= b269e01f2c4732b77eaf811e52fb2b36
MD5(txt/acos_2_mac--.txt)= 3db7942856e69a2944509ba0c4fc12c4
MD5(txt/acos_2_pi3--.txt)= b269e01f2c4732b77eaf811e52fb2b36
MD5(txt/acos_2_pi5--.txt)= b269e01f2c4732b77eaf811e52fb2b36
MD5(txt/drand48_hydra.txt)= 67f52a2c405c1246437733009e4870a7
MD5(txt/drand48_larry.txt)= 67f52a2c405c1246437733009e4870a7
MD5(txt/drand48_mac--.txt)= 67f52a2c405c1246437733009e4870a7
MD5(txt/drand48_pi3--.txt)= 67f52a2c405c1246437733009e4870a7
MD5(txt/drand48_pi5--.txt)= 67f52a2c405c1246437733009e4870a7
UNIX> 
Of the 10,000 values, 369 of them do not match:
UNIX> diff -y txt/acos_2_mac--.txt txt/acos_2_hydra.txt | grep '|' | wc
     369    1107   17712
UNIX> 
I can tell you from experience that this can be a real problem if you have applications that need trigonometric functions.

TrianglesContainOrigin

I have the problem description, examples and an input generator at http://web.eecs.utk.edu/~jplank/topcoder-writeups/2014/TrianglesContainOrigin/.

In this directory, I have my presentation in TrianglesContainOrigin.odp, converted to PDF in TrianglesContainOrigin.pdf.

Below is the commented code (src/make_jgraph.cpp) which creates the jgraph in the presentation, and also solves the problem in the fastest way:

/* TrianglesContainOrigin/src/make_jgraph.cpp
   James S. Plank
   September, 2020.

   This program takes input to the TrianglesContainOrigin Topcoder problem on 
   standard input, and then draws the jgraph which can help you visualize a solution.
   It then solves the problem and prints it as a comment on the last line in the 
   jgraph file.
 */

#include <string>
#include <vector>
#include <cmath>
#include <map>
#include <iostream>
#include <cstdio>
#include <cstdlib>
using namespace std;

/* Bundling up all of my information in this data structure makes the problem easier. */

class Point {
  public:
    double x;          // The point's x value
    double y;          // The point's x value
    double angle;      // The point's angle from the origin (radians)
    int index;         // The point's number in the Points vector, which is sorted by angle.
};

int main()
{
  double hyp, a;                          // Hypotenuse, angle
  double pi;                              // Convenient to have pi in a variable
  map <double, Point *> m;                // Map of points sorted by angle
  vector <Point *> points;                // Vector of points sorted by angle
  vector <string> colors;                 // This is to draw the colored regions.
  long long num;                          // The number of triangles

  int i;                                  // Temporary variable
  int p1, p2, p3, p4;                     // These match the presentation
  double dx, dy;                          
  Point *p;
  map <double, Point *>::iterator mit;
  const char *c;

  /* Set pi, and put two colors (pink, blue) into colors */

  pi = acos(-1.0);
  colors.push_back("1 .7 .7");
  colors.push_back(".7 1 1");

  /* Read the points, calculate their angles from the origin and put them into the map */

  while (cin >> dx >> dy) {
    p = new Point;
    p->x = dx;
    p->y = dy;
    hyp = sqrt(dx*dx + dy*dy);
    p->angle = acos(dx/hyp);
    if (dy < 0) p->angle = 2 * pi - p->angle;
    m[p->angle] = p;
  }

  /* Now calculate their indices, and make the points array, sorted by angle. */

  for (mit = m.begin(); mit != m.end(); mit++) {
    p = mit->second;
    p->index = points.size();
    points.push_back(p);
  }

  /* Draw the jgraph axes */

  printf("newgraph\n");
  printf("xaxis min -10000 max 10000 size 6 nodraw\n");
  printf("yaxis min -10000 max 10000 size 6 nodraw\n");
  printf("clip\n");

  /* Draw the colored regions opposite each adjacent pair of points.
  If there are an odd number of points, make the first region light
  green, and alternate blue/pink for the others (I do this because
  otherwise the first and last regions are the same color).  */

  printf("\n(*Draw the colored regions *)\n\n");  // This is handy when you want to look
                                                   // at the jgraph.

  for (i = 0; i < (int) points.size(); i++) {
    p = points[i];   
    a = p->angle + pi;
    if (i == 0 && points.size()%2 == 1) {
      c = ".5 1 .5";
    } else {
      c = colors[i%colors.size()].c_str();
    }

    /* Multiplying by 100,000 means that points of the triangle will be outside of
       the axes.  By specifying "clip" above, anything outside of the axes won't
       be drawn. */

    printf("newline color %s poly pcfill %s pts %lg %lg 0 0 ", 
           c, c, 
           100000 * cos(a),
           100000 * sin(a));
    p = points[(i+1)%points.size()];   
    a = p->angle + pi;
    printf("%lg %lg\n", 100000 * cos(a), 100000 * sin(a));
  }

  printf("\n(*Draw the lines from the origin to each point *)\n\n");
  
  for (i = 0; i < (int) points.size(); i++) {
    printf("newline pts 0 0 %lg %lg\n", points[i]->x, points[i]->y);
  }

  printf("\n(*Draw the points *)\n\n");

  printf("newcurve marktype circle pts\n");
  for (i = 0; i < (int) points.size(); i++) {
    printf("  %lg %lg\n", points[i]->x, points[i]->y);
  }

  printf("\n(*Draw the origin *)\n\n");

  printf("newcurve marktype box color 1 0 0 pts 0 0\n");

  printf("\n(*Draw the labels if there are <= 26 points *)\n\n");

  if (points.size() <= 26) {
    for (i = 0; i < (int) points.size(); i++) {
      p = points[i];
      printf("newstring fontsize 14 hjc vjc x %lg y %lg : %c\n", 
             p->x + cos(p->angle)*500,
             p->y + sin(p->angle)*500,
             'A'+i);
    }
  }

  /* Do the calculation in the fastest way */

  /* Put a big ol sentinel onto the array. */

  p = new Point;
  p->x = 0;
  p->y = 0;
  p->angle = 4*pi; /* Actually, 3*pi will do */
  p->index = points.size();
  points.push_back(p);

  /* Now, enumerate p1 and p2, and incrementally calculate p3 and p4, and keep
     adding the number of points between p3 and p4 to the answer. */

  num = 0;
  p3 = 0;
  for (p1 = 0; p1 < (int) points.size()-1; p1++) {
    while (points[p3]->angle <= points[p1]->angle + pi) p3++;
    p4 = p3;
    for (p2 = p1+1; p2 < p3; p2++) {
      while (points[p4]->angle <= points[p2]->angle + pi) p4++;
      num += (points[p4]->index - points[p3]->index);
    }
  }
  
  printf("\n(*Number of triangles: %lld. *)\n", num);

  return 0;
}