When I do ternary search for floating points, I normally do it like this:

while (R - L > EPS) {
  double M1 = (2*L + R) / 3;
  double M2 = (L + 2*R) / 3;
  if (f(M1) <= f(M2)) R = M2;
  else L = M1;
}
return f(L);

It works perfectly, but when I try to do it with bitonic arrays (with integer indices), I sometimes get the wrong value, either because I'm off by 1 or because some edge cases don't work.

while (L < R) {
  int M1 = (2*L + R) / 3;
  int M2 = (L + 2*R) / 3;
  if (a[M1] <= a[M2]) R = M2;
  else L = M1;
}
return a[L]; // doesn't work :(

I try to address this by looping at the end when the search space is small enough. This addresses the small edge cases and off by 1 errors.

while (L + 3 < R) {
  int M1 = (2*L + R) / 3;
  int M2 = (L + 2*R) / 3;
  if (a[M1] <= a[M2]) R = M2;
  else L = M1;
}
int mn = a[L++];
while (L <= R) mn = min(mn, a[L++]);
return mn; // works, but loop is so, yeah

Is there maybe a cleaner way of doing ternary search on arrays without the brute force loop at the end? Like, do I need to tweak the formula for M1 and M2 or something similar?

Now that we know how to sieve properly, let's look at hacks to the Eratosthenes sieve to get some other interesting sieves.

Introducing a problem

Let's say we want to count the number of divisors of a number. One way is to check all numbers up to √n and check if n divides that number. The other way is to get its prime factorization and get the product of (exponent + 1) through combinatorics. Either method is O(√n) on average, thus O(n√n) if done for all numbers up to n.

But what if a problem asks us to print the number of divisors of all numbers from 1 to 10^7 under 3 seconds? An O(n√n) algorithm will be too slow! When I tried counting divisors using the square root method, it ran for about 61 seconds on my computer. That definitely won't run in time.

Fortunately, we can tweak the Eratosthenes sieve to count the number of divisors more efficiently and elegantly. And you will see that this technique works not only for number of divisors, but also for generating sum of divisors, totient function, biggest prime divisor, basically all functions that have to do with divisors!

Divisor Sieve O(n log n)

int divisors[n + 1];
for (int i = 1; i <= n; ++i)
 for (int j = i; j <= n; j += i)
  ++divisors[j];

OK, so what's up with this rather short code? This short code generates the number of divisors of all numbers less than or equal to n. Oh wow, just solved our problem! But wait, aren't we too rushed? Does this code even work? And how does it even work in the first place? How did we come up with O(n log n)? Don't worry, we'll answer those questions one by one.

Correctness

We want to count the number of divisors of a number. On another perspective, we can instead start from the divisor then increment the count of all its multiples. Do that for all divisors, we now have all divisor counts of all numbers up to n.Hurray.

Complexity

Now that we have proved that the algorithm is correct, how are we sure that the complexity is O(n log n) when it looks like O(n^2)? The answer is because we are summing a harmonic series. The inner loop runs ⌊n / i⌋ iterations, therefore, the total number of iterations is:

If you know calculus, the count just approximates to the Riemann sum of the harmonic series, which we can integrate to get n log n. Mathemagic at its finest.

In general, O(n log n) for n = 10^7 runs for approximately 1.700s on a normal computer, and even faster on online judges that do cloud computing. This already solves our problem, and heck, with really short code!

Sum of Divisors Sieve O(n log n)

int sumdiv[n + 1];
for (int i = 1; i <= n; ++i)
 for (int j = i; j <= n; j += i)
  sumdiv[j] += i;

We can also use this technique to get sum of divisors. Just increment by the divisor instead of just incrementing by 1.

Euler Totient Sieve O(n log log n)

int totient[n + 1];
for (int i = 1; i <= n; ++i) totient[i] = i;
for (int i = 2; i <= n; ++i)
 if (totient[i] == i)
  for (int j = i; j <= n; j += i)
   totient[j] -= totient[j] / i;

This technique could also generate the Euler totient function, where totient[a] is the number of positive integers less than a that is relatively prime to a. It's O(n log log n) because it does the inner loop only if the number is prime (see prime harmonic series).

You might ask why it's totient[j] -= totient[j] / i. This is due to the nature of the Euler totient function, which needs some background of number theory to prove. This wonderful blog entry by PraveenDhinwa provides a good explanation of it if you want the extensive proof.

Biggest Prime Divisor Sieve O(n log log n)

int big[n + 1] = {1, 1};
for (int i = 1; i <= n; ++i)
 if (big[i] == 1)
  for (int j = i; j <= n; j += i)
   big[j] = i;

We can tabulate the biggest prime divisor per number. This is useful for pruned prime checking (when big[p] == p) and easier prime factorization. You don't need to iterate through all the primes to prime factorize anymore, you just need to a single while loop, something like while (n > 1) { factors.push_back(big[n]); n /= big[n]; }.

There are many more extensions of the wonderful Eratosthenes sieve. If you know some interesting ones, please let me know as well so I can add it (and credit you) in this blog.

Disclaimer: I was the one who named the "sieves" so they're not official names or anything. Moreover, they're not technically sieves anymore (sieves are filters, but the functions above are number generators), but I thought the sieve label was cool so I put it, if you don't mind.

Recently, I created a blog entry about using C++ std::complex as an alternative to points for computational geometry. We will now apply this technique for more geometry for C++!

Triangle centers are important for they create a connection between triangles, circles, and angles in geometric construction. But the classical way to determine them is a little hassle to either derive, or code.

For example, we can get the circumcenter by constructing two perpendicular bisectors and intersecting them. Math for dummies provides a brief explanation of some triangle center formulae if you want to know what I'm talking about.

But can we generalize the way to get ALL kinds of triangle centers?

Barycentric Coordinates

Barycentric coordinates uses a kind of coordinate system through three vertices of the triangle as basis. Basically, a coordinate has three real numbers (a, b, c) determining the "weight" of the vertex at vertices (A, B, C) respectively. This paper provides a well-formed definition of the Barycentric coordinates.

Barycentric points can be determined using the formula (A*a + B*b + C*c) / (a + b + c). You might relate this formula to the center of mass or the weighted average of three objects in space in physics. In addition, you can observe that (a, b, c) == (k*a, k*b, k*c), meaning, a coordinate is not unique when mapped from a point.

The Bary Function

We need to have a function that converts Barycentric coordinates to Cartesian points. Well, the formula is rather straightforward so it's rather easy. Since std::complex allows us vector addition and scalar multiplication, it's a 1-liner. Remember, we're using complex numbers so don't forget to typedef std::complex<double> point.

point bary(point A, point B, point C, double a, double b, double c) {
    return (A*a + B*b + C*c) / (a + b + c);
}

Triangle Centers

How can we get the triangle centers from Barycentric coordinates? Here's the cheat sheet.

point centroid(point A, point B, point C) {
    // geometric center of mass
    return bary(A, B, C, 1, 1, 1);
}

point circumcenter(point A, point B, point C) {
    // intersection of perpendicular bisectors
    double a = norm(B - C), b = norm(C - A), c = norm(A - B);
    return bary(A, B, C, a*(b+c-a), b*(c+a-b), c*(a+b-c));
}

point incenter(point A, point B, point C) {
    // intersection of internal angle bisectors
    return bary(A, B, C, abs(B-C), abs(A-C), abs(A-B));
}

point orthocenter(point A, point B, point C) {
    // intersection of altitudes
    double a = norm(B - C), b = norm(C - A), c = norm(A - B);
    return bary(A, B, C, (a+b-c)*(c+a-b), (b+c-a)*(a+b-c), (c+a-b)*(b+c-a));
}

point excenter(point A, point B, point C) {
    // intersection of two external angle bisectors
    double a = abs(B - C), b = abs(A - C), c = abs(A - B);
    return bary(A, B, C, -a, b, c);

    //// NOTE: there are three excenters
    // return bary(A, B, C, a, -b, c);
    // return bary(A, B, C, a, b, -c);
}

Getting diabetes from the syntax sugar? I hope you did. If you want to know how it works, Wolfram-alpha provides a brief summary and some equations for each triangle center. The site also provides other centers I did not mention here, such as the Symmedian point.

For the proofs, some of them are straightforward. For example, for the centroid, you just need to show that if a = b = c = 1, then (A*a + B*b + C*c) / (a + b + c) == (A + B + C) / 3. The other triangle centers, however, have rather extensive proofs (e.g. orthocenter proof), so I suggest that you just look up papers online on your own.

You might also want to see Silvester's book "Geometry — Ancient & Modern" for he also showed useful theorems built around Barycentric coordinates, though it's not available online so you should buy it in a bookstore. If you also found some research papers that look useful, I would also like to know.

UPD: These formulas work in 3D and higher dimensions as well! For example, the circumcenter code will get the center of the circumsphere in 3D space.

UPD2: Pardon for the confusion, I'm still talking about centers of 2D triangles on higher dimensional spaces. The circumcenter code is for the circumcircle, not the circumsphere because there can be many circumspheres if the figure is determined by just 3 points.

It's always a hassle to define our 2D Geometry library during a contest. Is there a way to make our computational geometry lives easier in any way? Fortunately for us, there is, at least in C++, using complex numbers.

Complex numbers are of the form a + bi, where a is the real part and b imaginary. Thus, we can let a be the x-coordinate and b be the y-coordinate. Whelp, complex numbers can be represented as 2D vectors! Therefore, we can use complex numbers to define a point instead of defining the class ourselves. You can look at std::complex reference here.

Defining our point class

We can define our point class by typing typedef complex<double> point; at the start of our program. To access our x- and y-coordinates, we can macro the real() and imag() functions by using #define. Of course, don't forget to #include <complex> before anything.

#include <iostream>
#include <complex>
using namespace std;

// define x, y as real(), imag()
typedef complex<double> point;
#define x real()
#define y imag()

// sample program
int main() {
	point a = 2;
	point b(3, 7);
	cout << a << ' ' << b << endl; // (2, 0) (3, 7)
	cout << a + b << endl;         // (5, 7)
}

Oh goodie! We can use std:cout for debugging! We can also add points as vectors without having to define operator+. Nifty. And apparently, we can overall add points, subtract points, do scalar multiplication without defining any operator. Very nifty indeed.

Example

point a(3, 2), b(2, -7);

// vector addition and subtraction
cout << a + b << endl;   // (5,-5)
cout << a - b << endl;   // (1,9)

// scalar multiplication
cout << 3.0 * a << endl; // (9,6)
cout << a / 5.0 << endl; // (0.6,0.4)

Functions using std::complex

What else can we do with complex numbers? Well, there's a lot that is really easy to code.

1. Vector addition: a + b

2. Scalar multiplication: r * a

3. Dot product: (conj(a) * b).x

4. Cross product: (conj(a) * b).y
   
   notice: conj(a) * b = (ax*bx + ay*by) + i (ax*by — ay*bx)

5. Squared distance: norm(a - b)

6. Euclidean distance: abs(a - b)

7. Angle of elevation: arg(b - a)

8. Slope of line (a, b): tan(arg(b - a))

9. Polar to cartesian: polar(r, theta)

10. Cartesian to polar: point(abs(p), arg(p))

11. Rotation about the origin: a * polar(1.0, theta)

12. Rotation about pivot p: (a-p) * polar(1.0, theta) + p
    
    UPD: added more useful functions

13. Angle ABC: abs(remainder(arg(a-b) - arg(c-b), 2.0 * M_PI))
    
    remainder normalizes the angle to be between [-PI, PI]. Thus, we can get the positive non-reflex angle by taking its abs value.

14. Project p onto vector v: v * dot(p, v) / norm(v);

15. Project p onto line (a, b): a + (b - a) * dot(p - a, b - a) / norm(b - a)

16. Reflect p across line (a, b): a + conj((p - a) / (b - a)) * (b - a)

17. Intersection of line (a, b) and (p, q):

point intersection(point a, point b, point p, point q) {
  double c1 = cross(p - a, b - a), c2 = cross(q - a, b - a);
  return (c1 * q - c2 * p) / (c1 - c2); // undefined if parallel
}

Drawbacks

Using std::complex is very advantageous, but it has one disadvantage: you can't use std::cin or scanf. Also, if we macro x and y, we can't use them as variables. But that's rather minor, don't you think?

EDIT: Credits to Zlobober for pointing out that std::complex has issues with integral data types. The library will work for simple arithmetic like vector addition and such, but not for polar or abs. It will compile but there will be some errors in correctness! The tip then is to rely on the library only if you're using floating point data all throughout.