The problem.

the_problem.cpp

// This function finds a pair of approximations
// [c1 -s1]       [c2 -s2]
// [s1  c1]  and  [s2  c2]
// to given planar rotations subject to the following constraints:
// 1. c1, s1, c2, s2 are all integer
// 2. c1*c1 + s1*s1 == c2*c2 + s2*s2, i.e. the approximations have the same norm.
// This is done while simultaneously trying to minimize:
// 1. The relative errors of the approximations
// 2. The size of the integers in question (smaller is better)
// 3. The cost of computing this approximation when no integer multiplies are available.
//    Computing the minimum cost of an integer multiply is fairly tricky in and of itself.
//    Here, an added complication is that this is still a (scaled) rotation matrix, and
//    the rotation
//      x' = c*x - s*y
//      y' = s*x + c*y
//    can be equivalently factored as either
//      t = c*(x+y)
//      x' = t - (c+s)*y
//      y' = t + (s-c)*x
//    or
//      t = s*(x-y)
//      x' = t + (c-s)*x
//      y' = t + (c+s)*y
//    which has less integer multiplies and might be cheaper to compute, depending on how the
//    factors come out. Oh, and there might be a way to share computations between the products,
//    but this function doesn't consider that.