CGAL 4.7 - Algebraic Foundations: User Manual

\( \newcommand{\E}{\mathrm{E}} \) \( \newcommand{\A}{\mathrm{A}} \) \( \newcommand{\R}{\mathrm{R}} \) \( \newcommand{\N}{\mathrm{N}} \) \( \newcommand{\Q}{\mathrm{Q}} \) \( \newcommand{\Z}{\mathrm{Z}} \) \( \def\ccSum #1#2#3{ \sum_{#1}^{#2}{#3} } \def\ccProd #1#2#3{ \sum_{#1}^{#2}{#3} }\)

CGAL 4.7 - Algebraic Foundations

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Groups Pages
User Manual

Author
  Michael Hemmer

# Introduction

CGAL is targeting towards exact computation with non-linear objects, in particular objects defined on algebraic curves and surfaces. As a consequence types representing polynomials, algebraic extensions and finite fields play a more important role in related implementations. This package has been introduced to stay abreast of these changes. Since in particular polynomials must be supported by the introduced framework the package avoids the term number type. Instead the package distinguishes between the algebraic structure of a type and whether a type is embeddable on the real axis, or real embeddable for short. Moreover, the package introduces the notion of interoperable types which allows an explicit handling of mixed operations.

# Algebraic Structures

The algebraic structure concepts introduced within this section are motivated by their well known counterparts in traditional algebra, but we also had to pay tribute to existing types and their restrictions. To keep the interface minimal, it was not desirable to cover all known algebraic structures, e.g., we did not introduce concepts for such basic structures as groups or exceptional structures as skew fields.

Figure 1.1 Concept Hierarchy of Algebraic Structures

Figure 1.1 shows the refinement relationship of the algebraic structure concepts. `IntegralDomain`, `UniqueFactorizationDomain`, `EuclideanRing` and `Field` correspond to the algebraic structures with the same name. `FieldWithSqrt`, `FieldWithKthRoot` and `FieldWithRootOf` are fields that in addition are closed under the operations 'sqrt', 'k-th root' and 'real root of a polynomial', respectively. The concept `IntegralDomainWithoutDivision` also corresponds to integral domains in the algebraic sense, the distinction results from the fact that some implementations of integral domains lack the (algebraically always well defined) integral division. Note that `Field` refines `IntegralDomain`. This is because most ring-theoretic notions like greatest common divisors become trivial for `Field`s. Hence we see `Field` as a refinement of `IntegralDomain` and not as a refinement of one of the more advanced ring concepts. If an algorithm wants to rely on gcd or remainder computation, it is trying to do things it should not do with a `Field` in the first place.

The main properties of an algebraic structure are collected in the class `Algebraic_structure_traits`. In particular the (most refined) concept each concrete model `AS` fulfills is encoded in the tag `Algebraic_structure_traits<AS>::Algebraic_category` . An algebraic structure is at least `Assignable`, `CopyConstructible`, `DefaultConstructible` and `EqualityComparable`. Moreover, we require that it is constructible from `int`. For ease of use and since their semantic is sufficiently standard to presume their existence, the usual arithmetic and comparison operators are required to be realized via C++ operator overloading. The division operator is reserved for division in fields. All other unary (e.g., sqrt) and binary functions (e.g., gcd, div) must be models of the well known STL-concepts `AdaptableUnaryFunction` or `AdaptableBinaryFunction` concept and local to the traits class (e.g., `Algebraic_structure_traits<AS>::Sqrt()(x)` ). This design allows us to profit from all parts in the STL and its programming style and avoids the name-lookup and two-pass template compilation problems experienced with the old design using overloaded functions. However, for ease of use and backward compatibility all functionality is also accessible through global functions defined within namespace `CGAL`, e.g., `CGAL::sqrt(x)` . This is realized via function templates using the according functor of the traits class. For an overview see Section Algebraic Foundations Reference in the reference manual.

## Tags in Algebraic Structure Traits

## Algebraic Category

For a type `AS`, `Algebraic_structure_traits<AS>` provides several tags. The most important tag is the `Algebraic_category` tag, which indicates the most refined algebraic concept the type `AS` fulfills. The tag is one of; `Integral_domain_without_division_tag`, `Integral_domain_tag`, `Field_tag`, `Field_with_sqrt_tag`, `Field_with_kth_root_tag`, `Field_with_root_of_tag`, `Unique_factorization_domain_tag`, `Euclidean_ring_tag`, or even `Null_tag` in case the type is not a model of an algebraic structure concept. The tags are derived from each other such that they reflect the hierarchy of the algebraic structure concept, e.g., `Field_with_sqrt_tag` is derived from `Field_tag`.

## Exact and Numerical Sensitive

Moreover, `Algebraic_structure_traits<AS>` provides the tags `Is_exact` and `Is_numerical_sensitive`, which are both `Boolean_tag`s.

An algebraic structure is considered exact, if all operations required by its concept are computed such that a comparison of two algebraic expressions is always correct.

An algebraic structure is considered as numerically sensitive, if the performance of the type is sensitive to the condition number of an algorithm. Note that there is really a difference among these two notions, e.g., the fundamental type `int` is not numerical sensitive but considered inexact due to overflow. Conversely, types as `leda_real` or `CORE::Expr` are exact but sensitive to numerical issues due to the internal use of multi precision floating point arithmetic. We expect that `Is_numerical_sensitive` is used for dispatching of algorithms, while `Is_exact` is useful to enable assertions that can be check for exact types only.

Tags are very useful to dispatch between alternative implementations. The following example illustrates a dispatch for `Field`s using overloaded functions. The example only needs two overloads since the algebraic category tags reflect the algebraic structure hierarchy.

File Algebraic_foundations/algebraic_structure_dispatch.cpp

#include <CGAL/basic.h>
#include <CGAL/IO/io.h>
#include <CGAL/Algebraic_structure_traits.h>
template< typename NT > NT unit_part(const NT& x);
template< typename NT >
NT unit_part_(const NT& x, CGAL::Field_tag);
template< typename NT >
NT unit_part_(const NT& x, CGAL::Integral_domain_without_division_tag);
template< typename NT >
NT unit_part(const NT& x){
// the unit part of 0 is defined as 1.
if (x == 0 ) return NT(1);
typedef CGAL::Algebraic_structure_traits<NT> AST;
typedef typename AST::Algebraic_category Algebraic_category;
return unit_part_(x,Algebraic_category());
}
template< typename NT >
NT unit_part_(const NT& x, CGAL::Integral_domain_without_division_tag){
// For many other types the only units are just -1 and +1.
return NT(int(CGAL::sign(x)));
}
template< typename NT >
NT unit_part_(const NT& x, CGAL::Field_tag){
// For Fields every x != 0 is a unit.
// Therefore, every x != 0 is its own unit part.
return x;
}
int main(){
// Function call for a model of EuclideanRing, i.e. int.
std::cout<< "int: unit_part(-3 ): " << unit_part(-3 ) << std::endl;
// Function call for a model of FieldWithSqrt, i.e. double
std::cout<< "double: unit_part(-3.0): " << unit_part(-3.0) << std::endl;
return 0;
}
// Note that this is just an example
// This implementation for unit part won't work for some types, e.g.,
// types that are not RealEmbeddable or types representing structures that have
// more units than just -1 and +1. (e.g. MP_Float representing Z[1/2])
// From there Algebraic_structure_traits provides the functor Unit_part.

# Real Embeddable

Most number types represent some subset of the real numbers. From those types we expect functionality to compute the sign, absolute value or double approximations. In particular we can expect an order on such a type that reflects the order along the real axis. All these properties are gathered in the concept `RealEmbeddable`. The concept is orthogonal to the algebraic structure concepts, i.e., it is possible that a type is a model of `RealEmbeddable` only, since the type may just represent values on the real axis but does not provide any arithmetic operations.

As for algebraic structures this concept is also traits class oriented. The main functionality related to `RealEmbeddable` is gathered in the class `Real_embeddable_traits`. In particular, it porivdes the boolean tag `Is_real_embeddable` indicating whether a type is a model of `RealEmbeddable`. The comparison operators are required to be realized via C++ operator overloading. All unary functions (e.g. sign, to_double) and binary functions (e.g. compare ) are models of the STL-concepts `AdaptableUnaryFunction` and `AdaptableBinaryFunction` and are local to `Real_embeddable_traits`.

In case a type is a model of `IntegralDomainWithoutDivision` and `RealEmbeddable` the number represented by an object of this type is the same for arithmetic and comparison. It follows that the ring represented by this type is a superset of the integers and a subset of the real numbers and hence has characteristic zero.

In case the type is a model of `Field` and `RealEmbeddable` it is a superset of the rational numbers.

# Real Number Types

Every CGAL `Kernel` comes with two real number types (number types embeddable into the real numbers). One of them is a `FieldNumberType`, and the other a `RingNumberType`. The coordinates of the basic kernel objects (points, vectors, etc.) come from one of these types (the `FieldNumberType` in case of Cartesian kernels, and the `RingNumberType` for Homogeneous kernels).

The concept `FieldNumberType` combines the requirements of the concepts `Field` and `RealEmbeddable`, while `RingNumberType` combines `IntegralDomainWithoutDivision` and `RealEmbeddable`. Algebraically, the real number types do not form distinct structures and are therefore not listed in the concept hierarchy of Figure 1.1.

# Interoperability

This section introduces two concepts for interoperability of types, namely `ImplicitInteroperable` and `ExplicitInteroperable`. While `ExplicitInteroperable` is the base concept, we start with `ImplicitInteroperable` since it is the more intuitive one.

In general mixed operations are provided by overloaded operators and functions or just via implicit constructor calls. This level of interoperability is reflected by the concept `ImplicitInteroperable`. However, within template code the result type, or so called coercion type, of a mixed arithmetic operation may be unclear. Therefore, the package introduces `Coercion_traits` giving access to the coercion type via `Coercion_traits<A,B>::Type` for two interoperable types `A` and `B`.

Some trivial example are `int` and `double` with coercion type double or `Gmpz` and `Gmpq` with coercion type `Gmpq`. However, the coercion type is not necessarily one of the input types, e.g. the coercion type of a polynomial with integer coefficients that is multiplied by a rational type is supposed to be a polynomial with rational coefficients.

`Coercion_traits` is also required to provide a functor `Coercion_traits<A,B>::Cast()` , that converts from an input type into the coercion type. This is in fact the core of the more basic concept `ExplicitInteroperable`. `ExplicitInteroperable` has been introduced to cover more complex cases for which it is hard or impossible to guarantee implicit interoperability. Note that this functor can be useful for `ImplicitInteroperable` types as well, since it can be used to void redundant type conversions.

In case two types `A` and `B` are `ExplicitInteroperable` with coercion type `C` they are valid argument types for all binary functors provided by `Algebraic_structure_traits` and `Real_embeddable_traits` of `C`. This is also true for the according global functions.

## Examples

The following example illustrates how two write code for `ExplicitInteroperable` types.

File Algebraic_foundations/interoperable.cpp

#include <CGAL/basic.h>
#include <CGAL/Coercion_traits.h>
#include <CGAL/IO/io.h>
// this is an implementation for ExplicitInteroperable types
// the result type is determined via Coercion_traits<A,B>
template <typename A, typename B>
typename CGAL::Coercion_traits<A,B>::Type
binary_func(const A& a , const B& b){
typedef CGAL::Coercion_traits<A,B> CT;
// check for explicit interoperability
CGAL_static_assertion((CT::Are_explicit_interoperable::value));
// CT::Cast is used to to convert both types into the coercion type
typename CT::Cast cast;
// all operations are performed in the coercion type
return cast(a)*cast(b);
}
int main(){
// Function call for the interoperable types
std::cout<< binary_func(double(3), int(5)) << std::endl;
// Note that Coercion_traits is symmetric
std::cout<< binary_func(int(3), double(5)) << std::endl;
return 0;
}

The following example illustrates a dispatch for `ImplicitInteroperable` and `ExplicitInteroperable` types. The binary function (that just multiplies its two arguments) is supposed to take two `ExplicitInteroperable` arguments. For `ImplicitInteroperable` types a variant that avoids the explicit cast is selected.

File Algebraic_foundations/implicit_interoperable_dispatch.cpp

#include <CGAL/basic.h>
#include <CGAL/Coercion_traits.h>
#include <CGAL/Quotient.h>
#include <CGAL/Sqrt_extension.h>
#include <CGAL/IO/io.h>
// this is the implementation for ExplicitInteroperable types
template <typename A, typename B>
typename CGAL::Coercion_traits<A,B>::Type
binary_function_(const A& a , const B& b, CGAL::Tag_false){
std::cout << "Call for ExplicitInteroperable types: " << std::endl;
typedef CGAL::Coercion_traits<A,B> CT;
typename CT::Cast cast;
return cast(a)*cast(b);
}
// this is the implementation for ImplicitInteroperable types
template <typename A, typename B>
typename CGAL::Coercion_traits<A,B>::Type
binary_function_(const A& a , const B& b, CGAL::Tag_true){
std::cout << "Call for ImpicitInteroperable types: " << std::endl;
return a*b;
}
// this function selects the correct implementation
template <typename A, typename B>
typename CGAL::Coercion_traits<A,B>::Type
binary_func(const A& a , const B& b){
typedef CGAL::Coercion_traits<A,B> CT;
typedef typename CT::Are_implicit_interoperable Are_implicit_interoperable;
return binary_function_(a,b,Are_implicit_interoperable());
}
int main(){
CGAL::set_pretty_mode(std::cout);
// Function call for ImplicitInteroperable types
std::cout<< binary_func(double(3), int(5)) << std::endl;
// Function call for ExplicitInteroperable types
CGAL::Quotient<int> rational(1,3); // == 1/3
CGAL::Sqrt_extension<int,int> extension(1,2,3); // == 1+2*sqrt(3)
CGAL::Sqrt_extension<CGAL::Quotient<int>,int> result = binary_func(rational, extension);
std::cout<< result << std::endl;
return 0;
}

# Fractions

Beyond the need for performing algebraic operations on objects as a whole, there are also number types which one would like to decompose into numerator and denominator. This does not only hold for rational numbers as `Quotient`, `Gmpq`, `mpq_class` or `leda_rational`, but also for compound objects as `Sqrt_extension` or `Polynomial` which may decompose into a (scalar) denominator and a compound numerator with a simpler coefficient type (e.g. integer instead of rational). Often operations can be performed faster on these denominator-free multiples. In case a type is a `Fraction` the relevant functionality as well as the numerator and denominator type are provided by `Fraction_traits`. In particular `Fraction_traits` provides a tag `Is_fraction` that can be used for dispatching.

A related class is `Rational_traits` which has been kept for backward compatibility reasons. However, we recommend to use `Fraction_traits` since it is more general and offers dispatching functionality.

## Examples

The following example show a simple use of `Fraction_traits`:
File Algebraic_foundations/fraction_traits.cpp

#include <CGAL/basic.h>
#include <CGAL/Fraction_traits.h>
#include <CGAL/IO/io.h>
#ifdef CGAL_USE_GMP
#include <CGAL/Gmpz.h>
#include <CGAL/Gmpq.h>
int main(){
typedef CGAL::Fraction_traits<CGAL::Gmpq> FT;
typedef FT::Numerator_type Numerator_type;
typedef FT::Denominator_type Denominator_type;
CGAL_static_assertion((boost::is_same<Numerator_type,CGAL::Gmpz>::value));
CGAL_static_assertion((boost::is_same<Denominator_type,CGAL::Gmpz>::value));
Numerator_type numerator;
Denominator_type denominator;
CGAL::Gmpq fraction(4,5);
FT::Decompose()(fraction,numerator,denominator);
CGAL::set_pretty_mode(std::cout);
std::cout << "decompose fraction: "<< std::endl;
std::cout << "fraction : " << fraction << std::endl;
std::cout << "numerator : " << numerator<< std::endl;
std::cout << "denominator: " << denominator << std::endl;
std::cout << "re-compose fraction: "<< std::endl;
fraction = FT::Compose()(numerator,denominator);
std::cout << "fraction : " << fraction << std::endl;
}
#else
int main(){ std::cout << "This examples needs GMP" << std::endl; }
#endif

The following example illustrates the integralization of a vector, i.e., the coefficient vector of a polynomial. Note that for minimizing coefficient growth `Fraction_traits<Type>::Common_factor` is used to compute the least common multiple of the denominators.

File Algebraic_foundations/integralize.cpp

#include <CGAL/basic.h>
#include <CGAL/Fraction_traits.h>
#include <CGAL/IO/io.h>
#include <vector>
template <class Fraction>
std::vector<typename CGAL::Fraction_traits<Fraction>::Numerator_type >
integralize(
const std::vector<Fraction>& vec,
typename CGAL::Fraction_traits<Fraction>::Denominator_type& d
) {
typedef CGAL::Fraction_traits<Fraction> FT;
typedef typename FT::Numerator_type Numerator_type;
typedef typename FT::Denominator_type Denominator_type;
typename FT::Decompose decompose;
std::vector<Numerator_type> num(vec.size());
std::vector<Denominator_type> den(vec.size());
// decompose each coefficient into integral part and denominator
for (unsigned int i = 0; i < vec.size(); i++) {
decompose(vec[i], num[i], den[i]);
}
// compute 'least' common multiple of all denominator
// We would like to use gcd, so let's think of Common_factor as gcd.
typename FT::Common_factor gcd;
d = 1;
for (unsigned int i = 0; i < vec.size(); i++) {
d *= CGAL::integral_division(den[i], gcd(d, den[i]));
}
// expand each (numerator, denominator) pair to common denominator
for (unsigned int i = 0; i < vec.size(); i++) {
// For simplicity ImplicitInteroperability is expected in this example
num[i] *= CGAL::integral_division(d, den[i]);
}
return num;
}
#ifdef CGAL_USE_GMP
#include <CGAL/Gmpz.h>
#include <CGAL/Gmpq.h>
int main(){
std::vector<CGAL::Gmpq> vec(3);
vec[0]=CGAL::Gmpq(1,4);
vec[1]=CGAL::Gmpq(1,6);
vec[2]=CGAL::Gmpq(1,10);
std::cout<< "compute an integralized vector" << std::endl;
std::cout<<"input vector: ["
<< vec[0] << "," << vec[1] << "," << vec[2] << "]" << std::endl;
CGAL::Gmpz d;
std::vector<CGAL::Gmpz> integral_vec = integralize(vec,d);
std::cout<<"output vector: ["
<< integral_vec[0] << ","
<< integral_vec[1] << ","
<< integral_vec[2] << "]" << std::endl;
std::cout<<"denominator : "<< d <<std::endl;
}
#else
int main(){ std::cout << "This examples needs GMP" << std::endl; }
#endif

# Design and Implementation History

The package is part of CGAL since release 3.3. Of course the package is based on the former Number type support of CGAL. This goes back to Stefan Schirra and Andreas Fabri. But on the other hand the package is to a large extend influenced by the experience with the number type support in Exacus [1], which in the main goes back to Lutz Kettner, Susan Hert, Arno Eigenwillig and Michael Hemmer. However, the package abstracts from the pure support for number types that are embedded on the real axis which allows the support of polynomials, finite fields, and algebraic extensions as well. See also related subsequent chapters.