Computations on algebraic numbers

Question

I'd like to do exact computations on regular polygons. To do so, I wrote the code you find below. But the expression cos*cos will not compile. Apparently multiplication is not defined for the algebraic number type I'm using. I guess I'll have to try some other approach. Currently there seem to be two candidates:

1. RS apparently provides some more advanced algebraic functionality, and CGAL has support for it. But I see no multiplication operators in the related headers, so I doubt that it will do multiplication the way I'd want.
2. leda::real seems to be a type for algebraic real numbers. I'll probably have to rewrite my code, but it should be possible to achieve similar results. Perhaps I could even convert the cos I computed in CGAL to such a leda::real. The LEDA header at least appears to have an operator*. LEDA is free for my use but still closed source. And leda_real.h for CGAL 4.3 looks strange: it refers to leda_real not leda::real, so perhaps it is written for an outdated version of LEDA. And it apparently includes itself, which looks pretty pointless.

Which of these alternatives would work best for construction of an exact CGAL kernel capable of describing regular n-gons for arbitrary n? Does any of these work at all? Is there another alternative I'm missing?

Since I don't have either RS or LEDA installed on my computer, I'd prefer an educated opinion before I start building them, and perhaps even writing install instructions (“ebuilds”) for my Gentoo linux.

#include <string>
#include <iostream>
#include <sstream>
#include <vector>

//define CGAL_USE_RS
#include <CGAL/Gmpz.h>
#include <CGAL/Algebraic_kernel_d_1.h>
#include <CGAL/Algebraic_kernel_rs_gmpz_d_1.h>
#include <CGAL/Homogeneous.h>
#include <CGAL/Arr_segment_traits_2.h>
#include <CGAL/Arrangement_2.h>

#define DBG(x) std::cerr << x << std::endl

typedef CGAL::Gmpz ZZ;
// typedef CGAL::Algebraic_kernel_rs_gmpz_d_1 AK;
typedef CGAL::Algebraic_kernel_d_1<ZZ> AK;
typedef AK::Polynomial_1 Polynomial;
typedef AK::Algebraic_real_1 AA;
typedef AK::Coefficient Coeff;
typedef AK::Bound Bound;
typedef AK::Multiplicity_type Multiplicity;
typedef CGAL::Homogeneous<AK> Kernel;
typedef CGAL::Arr_segment_traits_2<Kernel> Traits;
typedef Kernel::Point_2 Point;
typedef Kernel::Segment_2 Segment;
typedef CGAL::Arrangement_2<Traits> Arrangement;

static unsigned run(unsigned short n) {
  AK ak;
  AK::Construct_algebraic_real_1 to_AA = ak.construct_algebraic_real_1_object();
  AK::Solve_1 solve = ak.solve_1_object();
  Polynomial x{CGAL::shift(Polynomial(1), 1)}, twox{2*x};
  Polynomial a{1}, b{x};
  for (unsigned short i = 2; i <= n; ++i) {
    Polynomial c = twox*b - a;
    a = b;
    b = c;
  }
  std::vector<std::pair<AA, Multiplicity>> roots;
  solve(b - 1, std::back_inserter(roots));
  AA one{1}, cos{-2};
  for (auto i = roots.begin(), e = roots.end(); i != e; ++i) {
    AA cur = i->first;
    if (cur < one && cur > cos)
      cos = cur;
  }
  AA sin = CGAL::sqrt(to_AA(1) - cos*cos);
  //DBG("sin="<<CGAL::to_double(sin)<<", cos="<<CGAL::to_double(cos));
  return 0;
}

int main(int argc, char** argv) {
  for (int i = 1; i < argc; ++i) {
    unsigned short n;
    std::istringstream(argv[i]) >> n;
    std::cout << n << ": " << run(n) << std::endl;
  }
  return 0;
}

Answer 1

CGAL also comes with the CORE library, which provide the operations you need.

Here is some code (provided by the OP himself) to compute that sin and cos exactly:

#include <utility>
#include <CGAL/CORE_Expr.h>
#include <CGAL/Polynomial.h>
#include <CGAL/number_utils.h>

typedef CORE::Expr AA;
typedef CGAL::Polynomial<AA> Polynomial;

// return sin(θ) and cos(θ) for θ = 2π/n
static std::pair<AA, AA> sin_cos(unsigned short n) {
  // We actually use -x instead of x since root_of will give the k-th
  // smallest root but we want the second largest one without counting.
  Polynomial x{CGAL::shift(Polynomial(-1), 1)}, twox{2*x};
  Polynomial a{1}, b{x};
  for (unsigned short i = 2; i <= n; ++i) {
    Polynomial c = twox*b - a;
    a = b;
    b = c;
  }
  a = b - 1;
  AA cos = -CGAL::root_of(2, a.begin(), a.end());
  AA sin = CGAL::sqrt(AA(1) - cos*cos);
  return std::make_pair(sin, cos);
}