CGALの計算基盤

CGALの幾何計算の基盤には、以下のような群論の本に出てくる得体の知れない概念がある。

  • IntegralDomainWithoutDivision:割り算のない整域
  • IntegralDomain(整域):整数の環を一般化したもの
  • UniqueFactorizationDomain(一意分解環):因数分解が一意。複素数は異なる。
  • Euclidean Ring(ユークリッド環):
  • PrincipalIdeal(主イデアル):環の、和と差に関して閉じていて、なおかつかけ算にも閉じているような部分集合(イデアル)のうち、元が1つのもの。
  • Field(体, Körper):四則演算に制限がない有理数のような集合。
  • FieldWithSqrt/FieldWithKthRoot/FieldWithRootOf:sqrt、Kルート、多項式ルートに関して閉じているような体。

これらの対象に関して、CGALのフレームワークは、比較可能かとか絶対値などの、普段は意識しない数的操作の要件をAlgebraic_structure_traitsとして指定している。そのうち、最も基本的な条件は、ExactとNumericalSensitiveの概念だ。Exactとは、すべての演算によって2値の比較は常に正しいということで、NumericalSensitiveとは、アルゴリズムの条件に対して敏感かどうかという意味。
例えば、int型は、NumericalSensitiveではなく、オーバーフローで精度が落ちるのでexactではない。一方、leda_realのような高精度型は、オーバーフローなどの精度悪化はないが、内部に複雑なデータを抱えているのでNumericalSensitiveである。

#include <CGAL/basic.h>
#include <CGAL/IO/io.h>
#include <CGAL/Algebraic_structure_traits.h>
// 特殊な体(field)を指定する
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;
}

実数軸埋め込み可能性 :
コンピュータで扱う数の型は本質的には実数のサブセットでしかないが、我々は当然、その符号、絶対値、ダブル精度の概数値が使えることを期待する。とりわけ、実数軸上で順番に並べられることは重要で、それらの属性を「実数軸埋め込み可能性」として定義している。(これは代数演算の概念とは無関係)。当然、比較演算子が使える対象として実装されている。
CGALでは体型と環型の二つの実数型がある。体型はデカルト座標で、環型は斉次座標で使われる。

相互演算可能性(Interoperability):
型の異なる数の演算ではどの型に合わせられるかがはっきりしない。CGALはこれをCoercion_trait::Typeとして把握する。最も簡単な例は、整数型と浮動小数点型の演算結果で、浮動小数点になる。複雑な例では、整数係数の多項式型と有理数のかけ算では、有理数係数の多項式型になる。

明示的演算可能型の例。Coercion_traitsはA,Bの型の強制変換を司る。

#include <CGAL/basic.h>
#include <CGAL/Coercion_traits.h>
#include <CGAL/IO/io.h>
// 掛け算の関数
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;
 // 明示的に強制変換が設定されているか
 CGAL_static_assertion((CT::Are_explicit_interoperable::value));
 // 強制変換にCT::Castが使われる
 typename CT::Cast cast;
 return cast(a)*cast(b);
}
int main(){
 std::cout<< binary_func(double(3), int(5)) << std::endl;
 std::cout<< binary_func(int(3), double(5)) << std::endl;
 return 0;
}

暗黙演算可能型と明示演算可能型でのバイナリー関数の例

#include <CGAL/basic.h>
#include <CGAL/Coercion_traits.h>
#include <CGAL/Quotient.h>
#include <CGAL/Sqrt_extension.h>
#include <CGAL/IO/io.h>
// 明示的掛け算の実装
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);
}
// 暗黙掛け算の実装
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;
}
// 呼び出す関数を選択(最後に_がないことに注意)
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;
}

分数型(有理数型):数を分数型のまま扱う型で、Quotient,Gmpq,leda_rationalなどがある。

#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 << "分数 : " << fraction << std::endl;
 std::cout << "分子 : " << numerator<< std::endl;
 std::cout << "分母: " << denominator << std::endl;
 fraction = FT::Compose()(numerator,denominator);
 std::cout << "再構成された分数 : " << fraction << std::endl;
}
#else
int main(){ std::cout << "GMPがCGALに組み込まれていない" << std::endl; }
#endif

Integralization of a vector, i.e., the coefficient vector of a polynomialの例。分母の最小公倍数計算にFraction_traits<Type>::Common_factor が使われている。

#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());
 for (unsigned int i = 0; i < vec.size(); i++) {
  decompose(vec[i], num[i], den[i]);
 }
 // 全分母の最小公倍数を求める
 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]));
 }
 // 最小公倍数を分母としたときの各分子を計算する
 for (unsigned int i = 0; i < vec.size(); i++) {
  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<<"分数のベクトル: [" << vec[0] << "," << vec[1] << "," << vec[2] << "]" << std::endl;
 CGAL::Gmpz d;
 std::vector<CGAL::Gmpz> integral_vec = integralize(vec,d);
 std::cout<<"分子のベクトル: [" << integral_vec[0] << "," << integral_vec[1] << "," << integral_vec[2] << "]" << std::endl;
 std::cout<<"共通分母: "<< d <<std::endl;
}
#else
int main(){ std::cout << "GMPがCGALに組み込まれていない" << std::endl; }
#endif

CGALで使う数型


  • float,double,long double:高速だが本質的には概数。2線の交点がその線上にあるという当たり前の判定に失敗する。
  • CGAL型:Quotient(商型)は整数による分数表現。MP_Float型は2の累乗による多精度浮動小数点。Lazy_exact_ntは高速化のためまず浮動小数点比較を試みる。
  • GMP型(GNU Multiple Preision):CGAL/Gmpz.hで導入する。
  • その他、LEDA型、CORE型もある。
  • INTERVAL_NT型:ある式の符号を求めるとき、2点で挟み込んでその間隔を極小まで変化させながら正負をチェックするという手法などで使う。
  • Residue型:素数(静的に定義されている)で割った余りを扱う。中国の剰余定理(3で割ると2余り、5で割ると3余り、7で割ると2余る数は何か)のような問題を扱うためにモジュラー算術が定義されている。


多項式の最大公約数を高速化する例

#include <CGAL/basic.h>
#ifdef CGAL_USE_GMP
#include <CGAL/Gmpz.h> 
#include <CGAL/Polynomial.h>
// 多項式の関数定義
template< typename Polynomial > bool may_have_common_factor(const Polynomial& p1, const Polynomial& p2, CGAL::Tag_true){
 std::cout<< "数型はモジュラー可能か" << std::endl; 
 // IEEE double精度で丸める
 CGAL::Protect_FPU_rounding<true> pfr(CGAL_FE_TONEAREST);
 typedef CGAL::Modular_traits<Polynomial> MT;
 typedef typename MT::Residue_type MPolynomial;
 typedef typename MT::Modular_image Modular_image;

 MPolynomial mp1 = Modular_image()(p1);
 MPolynomial mp2 = Modular_image()(p2);
 // 度数をチェックする
 typename CGAL::Polynomial_traits_d<Polynomial>::Degree degree;
 typename CGAL::Polynomial_traits_d<MPolynomial>::Degree mdegree;
 if ( degree(p1) != mdegree(mp1)) return true; 
 if ( degree(p2) != mdegree(mp2)) return true; 
 // 最大公約数gcdを計算する
 MPolynomial mg = CGAL::gcd(mp1,mp2);
 if ( mdegree(mg) > 0 ){
  std::cout << "gcdは自明でない" << std::endl;
  return true;
 }else{
  std::cout << "gcdは自明" << std::endl;
  return false; 
 }
}
template< typename Polynomial > bool may_have_common_factor(const Polynomial&, const Polynomial&, CGAL::Tag_false){
 std::cout<< "モジュラー可能でない" << std::endl; 
 return true; 
}
template< typename Polynomial > Polynomial modular_filtered_gcd(const Polynomial& p1, const Polynomial& p2){
 typedef CGAL::Modular_traits<Polynomial> MT;
 typedef typename MT::Is_modularizable Is_modularizable;
 // よけいな計算をしないために
 if (may_have_common_factor(p1,p2, Is_modularizable())){
  return CGAL::gcd(p1,p2);
 }else{
  typename CGAL::Polynomial_traits_d<Polynomial>::Univariate_content content;
  typename CGAL::Polynomial_traits_d<Polynomial>::Construct_polynomial construct;
  return construct(CGAL::gcd(content(p1),content(p2)));
 }
}
int main(){
 CGAL::set_pretty_mode(std::cout);
 typedef CGAL::Gmpz NT; 
 typedef CGAL::Polynomial<NT> Poly; 
 CGAL::Polynomial_traits_d<Poly>::Construct_polynomial construct;

 Poly f1=construct(NT(2), NT(6), NT(4)); 
 Poly f2=construct(NT(12), NT(4), NT(8));
 Poly f3=construct(NT(3), NT(4));
 std::cout << "f1 : " << f1 << std::endl;
 std::cout << "f2 : " << f2 << std::endl;

 Poly g1 = modular_filtered_gcd(f1,f2);
 std::cout << "gcd(f1,f2): " << g1 << std::endl;
 Poly p1 = f1*f3;
 Poly p2 = f2*f3;
 std::cout << "f3 : " << f3 << std::endl;
 std::cout << "p1=f1*f3 : " << p1 << std::endl;
 std::cout << "p2=f2*f3 : " << p2 << std::endl;

 Poly g2 = modular_filtered_gcd(p1,p2);
 std::cout << "gcd(p1,p2): " << g2 << std::endl;
}
#else
int main(){ std::cout << "GMPがCGALに組み込まれていない" << std::endl; }
#endif 

多項式型:
一般的な多項式の構築は、iteratorで係数をPolynomial_traits_d::Coefficient_typeとして与える方法と、std::pair<Exponent_vector, Polynomial_traits_d::Innermost_coefficient_type>として与える方法がある。すこし分かりにくいのは、shiftによるエレガントな多項式構築。

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
int main(){
 CGAL::set_pretty_mode(std::cout);
 typedef CGAL::Polynomial_type_generator<int,2>::Type Poly_2;
 typedef CGAL::Polynomial_traits_d<Poly_2> PT_2;
 typedef PT_2::Coefficient_type Poly_1;
 typedef PT_2::Innermost_coefficient_type Integer; 
 PT_2::Construct_polynomial construct_polynomial;
//  intによる多項式の構築
 Poly_2 two(2); // 2のこと 
 std::cout << "A constant polynomial: " << two << std::endl;
// iteratorによる多項式構築
 std::list<Poly_1> univariate_coeffs;
 univariate_coeffs.push_back(Poly_1(3));
 univariate_coeffs.push_back(Poly_1(0));
 univariate_coeffs.push_back(Poly_1(5));
 Poly_2 F = // 5*y^2 + 3のこと
 construct_polynomial(univariate_coeffs.begin(),univariate_coeffs.end());
 std::cout << "The bivariate polynomial F: " << F << std::endl;
// 単項式のiteratorによる構築
 std::list<std::pair<CGAL::Exponent_vector, Integer> > innermost_coeffs;
 innermost_coeffs.push_back(std::make_pair(CGAL::Exponent_vector(0,0),-2));
 innermost_coeffs.push_back(std::make_pair(CGAL::Exponent_vector(3,5),2));
 Poly_2 G = // (2*x^3)*y^5 + (-2)のこと
 construct_polynomial(innermost_coeffs.begin(),innermost_coeffs.end());
 std::cout << "The bivariate polynomial G: " << G << std::endl;
// shiftによる構築
 PT_2::Shift shift;
 Poly_2 x = shift(Poly_2(1),1,0); // Poly_2(1) * x0 ^ 1
 Poly_2 y = shift(Poly_2(1),1,1); // Poly_2(1) * x1 ^ 1
 Poly_2 H = 5 * x * y + 3 * y * y; // = 3*y^2 + (5*x)*yのこと
 std::cout << "The bivariate polynomial H: " << H << std::endl;
}

共役多項式を求める例

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
int main(){
 CGAL::set_pretty_mode(std::cout);
 typedef CGAL::Polynomial_type_generator<int,1>::Type Poly_1;
 typedef CGAL::Polynomial_traits_d<Poly_1> PT_1;
 PT_1::Shift shift;
 PT_1::Gcd gcd;
 PT_1::Gcd_up_to_constant_factor gcd_utcf;
 PT_1::Multivariate_content mcontent;
 PT_1::Canonicalize canonicalize;
 //construction using shift 
 Poly_1 x = shift(Poly_1(1),1,0); // x^1
 // common factor 7 * (x^2-2) 
 Poly_1 F = 21*(x-5)*(x*x-2); // = 21*x^3 + (-105)*x^2 + (-42)*x + 210
 Poly_1 G = 14*(x-3)*(x*x-2); // = 14*x^3 + (-42)*x^2 + (-28)*x + 84
 std::cout << "The univariate polynomial F: " << F << std::endl;
 std::cout << "The univariate polynomial G: " << G << std::endl;
 std::cout << "Common multivariate content: " 
  << CGAL::gcd(mcontent(F),mcontent(G)) << std::endl;
 std::cout << "The gcd of F and G: " << gcd(F,G) // = 7*x^2 + (-14)
  << std::endl;
 std::cout << "The gcd up to constant factor of F and G: " 
  << gcd_utcf(F,G) // = x^2 + (-2)
  << std::endl;
 std::cout << "Same as canonicalized gcd of F and G: " 
  << canonicalize(gcd_utcf(F,G)) // = x^2 + (-2)
  << std::endl;
}

代入と評価:要するに具体的な数値を代入した場合。

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
int main(){
 CGAL::set_pretty_mode(std::cout);
 typedef CGAL::Polynomial_type_generator<int,2>::Type Poly_2;
 typedef CGAL::Polynomial_traits_d<Poly_2> PT_2;
 //construction using shift 
 Poly_2 x = PT_2::Shift()(Poly_2(1),1,0); // x^1
 Poly_2 y = PT_2::Shift()(Poly_2(1),1,1); // y^1
 Poly_2 F = 2*x*y + 3*CGAL::ipower(y,3);
 std::cout << "The bivariate polynomial F: " << F // = 3*y^3 + (2*x)*y
  << std::endl << std::endl;
 PT_2::Evaluate evaluate;
 PT_2::Evaluate_homogeneous hevaluate;
 // Evaluation considers a polynomials as univariate:
 std::cout << "F(5): " << evaluate(F,5) // = 10*x + 375
  << std::endl; 
 // Evaluate_homogeneous considers F as a homogeneous polynomial in 
 // the outermost variable only, that is, F is interpreted as 
 // F(u,v) = 2*x*u*v^2 + 3 * u^3 
 std::cout << "F(5,7): " << hevaluate(F,5,7) // = 490*x + 375
  << std::endl << std::endl;
 PT_2::Substitute substitute;
 PT_2::Substitute_homogeneous hsubstitute;
 // Substitute considers a polynomials as multivariate, that is, the 
 // new values for the variables are given by an iterator range
 // Note that the value type only has to be interoperable with the innermost 
 // coefficient
 std::list<Poly_2> replacements;
 replacements.push_back(x-1); // replace x by x-1
 replacements.push_back(y); // replace y by y, i.e., do nothing
 std::cout << "The bivariate polynomial F: " << F // = 3*y^3 + (2*x)*y
  << std::endl;
 std::cout << "F(x-1,y): " // = 3*y^3 + (2*x + (-2))*y
  << substitute(F,replacements.begin(),replacements.end())
  << std::endl;
 // Substitute_homogeneous considers F as a homogeneous polynomial in 
 // all variable, that is, F is interpreted as 
 // F(x,y,w) = 2*x*y*w + 3 * y^3 
 replacements.push_back(y); // replace z by y 
 std::cout << "F(x-1,y,y): " // = 3*y^3 + (2*x + (-2))*y^2
  << hsubstitute(F,replacements.begin(),replacements.end())
  << std::endl;
}

終結式 (resultant)、主部分終結式(subresultant)など:終結式、ましてやSturm-Habicht sequence がなにかなんて知らないが、根の数を数えることができるらしい。

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
#include <CGAL/Gmpz.h>
int main(){
CGAL::set_pretty_mode(std::cout);
typedef CGAL::Gmpz Int;
typedef CGAL::Polynomial_type_generator<Int,1>::Type Poly_1;
typedef CGAL::Polynomial_traits_d<Poly_1> PT_1;
//construction using shift 
Poly_1 x = PT_1::Shift()(Poly_1(1),1); // x^1
Poly_1 F // = (x+1)^2*(x-1)*(2x-1)=2x^4+x^3-3x^2-x+1
= 2 * CGAL::ipower(x,4) + 1 * CGAL::ipower(x,3)
- 3 * CGAL::ipower(x,2) - 1 * CGAL::ipower(x,1)
+ 1 * CGAL::ipower(x,0);
std::cout << "F=" << F << std::endl;
Poly_1 G // = (x+1)*(x+3)=x^2+4*x+3
= 1 * CGAL::ipower(x,2) + 4 * CGAL::ipower(x,1) + 3 * CGAL::ipower(x,0);
std::cout << "G=" << G << std::endl;
// Resultant computation:
PT_1::Resultant resultant;
std::cout << "The resultant of F and G is: " << resultant(F,G) << std::endl;
// It is zero, because F and G have a common factor
// Real root counting:
PT_1::Principal_sturm_habicht_sequence stha;
std::vector<Int> psc;
stha(F,std::back_inserter(psc));
int roots = CGAL::number_of_real_roots(psc.begin(),psc.end());
std::cout << "The number of real roots of F is: " << roots << std::endl; // 3
roots = CGAL::number_of_real_roots(G);
std::cout << "The number of real roots of G is: " << roots << std::endl; // 2
return 0;
}

代数カーネル:以下の例は、CGALをMPFIで構築しないと動かない。そもそも、何をやっているか言葉が難しい。

代数の実数の構築の例

#include <CGAL/basic.h>
#ifdef CGAL_USE_MPFI 
#include <CGAL/Algebraic_kernel_d_1.h>
#include <CGAL/Gmpz.h>
#include <vector>
#include <iostream>
typedef CGAL::Algebraic_kernel_d_1<CGAL::Gmpz> AK;
typedef AK::Polynomial_1 Polynomial_1;
typedef AK::Algebraic_real_1 Algebraic_real_1;
typedef AK::Coefficient Coefficient;
typedef AK::Bound Bound;
typedef AK::Multiplicity_type Multiplicity_type;
int main(){
 AK ak; // an object of 
 AK::Construct_algebraic_real_1 construct_algreal_1 = ak.construct_algebraic_real_1_object();
 std::cout << "Construct from int : " << construct_algreal_1(int(2)) << "\n";
 std::cout << "Construct from Coefficient : " << construct_algreal_1(Coefficient(2)) << "\n";
 std::cout << "Construct from Bound : " << construct_algreal_1(Bound(2)) << "\n\n";
 Polynomial_1 x = CGAL::shift(AK::Polynomial_1(1),1); // the monomial x
 std::cout << "Construct by index : " << construct_algreal_1(x*x-2,1) << "\n"
  << to_double(construct_algreal_1(x*x-2,1)) << "\n";
 std::cout << "Construct by isolating interval : "
  << construct_algreal_1(x*x-2,Bound(0),Bound(2)) << "\n"
  << to_double(construct_algreal_1(x*x-2,Bound(0),Bound(2))) << "\n\n";
 return 0;
}
#else
int main(){std::cout << "CGALはMPFIライブラリで構築しなければならない" << std::endl;return 0;}
#endif

多項式の解の例

#include <CGAL/basic.h>
#ifdef CGAL_USE_MPFI 
#include <CGAL/Algebraic_kernel_d_1.h>
#include <CGAL/Gmpz.h>
#include <vector>
typedef CGAL::Algebraic_kernel_d_1<CGAL::Gmpz> AK;
typedef AK::Polynomial_1 Polynomial_1;
typedef AK::Algebraic_real_1 Algebraic_real_1;
typedef AK::Bound Bound;
typedef AK::Multiplicity_type Multiplicity_type;
int main(){
 AK ak; // an object of 
 AK::Solve_1 solve_1 = ak.solve_1_object();
 Polynomial_1 x = CGAL::shift(AK::Polynomial_1(1),1); // the monomial x
 // variant using a bool indicating a square free polynomial, multiplicities are not computed
 std::vector<Algebraic_real_1> roots;
 solve_1(x*x-2,true, std::back_inserter(roots));
 std::cout << "Number of roots is : " << roots.size() << "\n";
 std::cout << "First root should be -sqrt(2): " << CGAL::to_double(roots[0]) << "\n";
 std::cout << "Second root should be sqrt(2): " << CGAL::to_double(roots[1]) << "\n";
 roots.clear();
 
 // variant for roots in a given range of a square free polynomial
 solve_1((x*x-2)*(x*x-3),true, Bound(0),Bound(10),std::back_inserter(roots));
 std::cout << "Number of roots is : " << roots.size() << "\n";
 std::cout << "First root should be sqrt(2): " << CGAL::to_double(roots[0]) << "\n";
 std::cout << "Second root should be sqrt(3): " << CGAL::to_double(roots[1]) << "\n\n";
 roots.clear();

 // variant computing all roots with multiplicities
 std::vector<std::pair<Algebraic_real_1,Multiplicity_type> > mroots;
 solve_1((x*x-2), std::back_inserter(mroots));
 std::cout << "Number of roots is : " << mroots.size() << "\n";
 std::cout << "First root should be -sqrt(2): " << CGAL::to_double(mroots[0].first) << ""
  << " with multiplicity " << mroots[0].second << "\n";
 std::cout << "Second root should be sqrt(2): " << CGAL::to_double(mroots[1].first) << ""
  << " with multiplicity " << mroots[1].second << "\n\n";
 mroots.clear();

 // variant computing roots with multiplicities for a range
 solve_1((x*x-2)*(x*x-3),Bound(0),Bound(10),std::back_inserter(mroots));
 std::cout << "Number of roots is : " << mroots.size() << "\n";
 std::cout << "First root should be sqrt(2): " << CGAL::to_double(mroots[0].first) << ""
  << " with multiplicity " << mroots[0].second << "\n";
 std::cout << "Second root should be sqrt(3): " << CGAL::to_double(mroots[1].first) << ""
  << " with multiplicity " << mroots[1].second << "\n\n";
 return 0;
}
#else
int main(){std::cout << "CGALはMPFIライブラリで構築しなければならない" << std::endl;return 0;}
#endif


比較の例

#include <CGAL/basic.h>
#ifdef CGAL_USE_MPFI 
#include <CGAL/Algebraic_kernel_d_1.h>
#include <CGAL/Gmpz.h>
#include <vector>
typedef CGAL::Algebraic_kernel_d_1<CGAL::Gmpz> AK;
typedef AK::Coefficient Coefficient;
typedef AK::Polynomial_1 Polynomial_1;
typedef AK::Algebraic_real_1 Algebraic_real_1;
typedef AK::Bound Bound;
typedef std::pair<Bound,Bound> Interval;
int main(){
 AK ak;
 AK::Construct_algebraic_real_1 construct_algebraic_real_1 = ak.construct_algebraic_real_1_object();
 Polynomial_1 x = CGAL::shift(AK::Polynomial_1(1),1); // the monomial x
 Algebraic_real_1 a = construct_algebraic_real_1(x*x-2,1); // sqrt(2)
 Algebraic_real_1 b = construct_algebraic_real_1(x*x-3,1); // sqrt(3)
 // Algebraic_real_1 is RealEmbeddable (just some functions:)
 std::cout << "sign of a is : " << CGAL::sign(a) << "\n";
 std::cout << "double approximation of a is : " << CGAL::to_double(a) << "\n";
 std::cout << "double approximation of b is : " << CGAL::to_double(b) << "\n";
 std::cout << "double lower bound of a : " << CGAL::to_interval(a).first << "\n";
 std::cout << "double upper bound of a : " << CGAL::to_interval(a).second << "\n";
 std::cout << "LessThanComparable (a<b) : " << (a<b) << "\n\n";
 // use compare_1 with int, Bound, Coefficient, Algebraic_real_1
 AK::Compare_1 compare_1 = ak.compare_1_object();
 std::cout << " compare with an int : " << compare_1(a ,int(2)) << "\n";
 std::cout << " compare with an Coefficient : " << compare_1(a ,Coefficient(2)) << "\n";
 std::cout << " compare with an Bound : " << compare_1(a ,Bound(2)) << "\n";
 std::cout << " compare with another Algebraic_real_1: " << compare_1(a ,b) << "\n\n";
 // get a value between two roots
 AK::Bound_between_1 bound_between_1 = ak.bound_between_1_object();
 std::cout << " value between sqrt(2) and sqrt(3) " << bound_between_1(a,b) << "\n";
 std::cout << " is larger than sqrt(2) " << compare_1(bound_between_1(a,b),a) << "\n";
 std::cout << " is less than sqrt(3) " << compare_1(bound_between_1(a,b),b) << "\n\n";
 // approximate with relative precision
 AK::Approximate_relative_1 approx_r = ak.approximate_relative_1_object();
 std::cout << " lower bound of a with at least 100 bits: "<< approx_r(a,100).first << "\n";
 std::cout << " upper bound of a with at least 100 bits: "<< approx_r(a,100).second << "\n\n";
 // approximate with absolute error
 AK::Approximate_absolute_1 approx_a = ak.approximate_absolute_1_object();
 std::cout << " lower bound of b with error less than 2^-100: "<< approx_a(b,100).first << "\n";
 std::cout << " upper bound of b with error less than 2^-100: "<< approx_a(b,100).second << "\n\n";
 return 0;
}
#else
int main(){std::cout << "CGALはMPFIライブラリで構築しなければならない" << std::endl;return 0;}
#endif


4.4 Isolation of Algebraic Real Numbers with respect to roots of other polynomials

#include <CGAL/basic.h>
#ifdef CGAL_USE_MPFI 
#include <CGAL/Algebraic_kernel_d_1.h>
#include <CGAL/Gmpz.h>
#include <vector>
typedef CGAL::Algebraic_kernel_d_1<CGAL::Gmpz> AK;
typedef AK::Polynomial_1 Polynomial_1;
typedef AK::Algebraic_real_1 Algebraic_real_1;
typedef AK::Coefficient Coefficient;
typedef AK::Bound Bound;
typedef AK::Multiplicity_type Multiplicity_type;
int main(){
 AK ak; // an object of 
 AK::Construct_algebraic_real_1 construct_algreal_1 = ak.construct_algebraic_real_1_object();
 AK::Isolate_1 isolate_1 = ak.isolate_1_object();
 AK::Compute_polynomial_1 compute_polynomial_1 = ak.compute_polynomial_1_object();
 // construct an algebraic number from an integer
 Algebraic_real_1 frominteger=construct_algreal_1(int(2));
 std::cout << "Construct from int: " << frominteger << "\n";
 // the constructed algebraic number is root of a polynomial
 Polynomial_1 pol=compute_polynomial_1(frominteger);
 std::cout << "The constructed number is root of: " << pol << "\n";
 // construct an algebraic number from a polynomial and an isolating interval
 Polynomial_1 x = CGAL::shift(AK::Polynomial_1(1),1); // the monomial x
 Algebraic_real_1 frominterval=construct_algreal_1(x*x-2,Bound(0),Bound(2));
 std::cout << "Construct from isolating interval: " << frominterval << "\n";
 // isolate the second algebraic number from the first: this is to say,
 // isolating the second algebraic number with respect to the polynomial
 // of which the first constructed number is root
 std::pair<Bound,Bound> isolation1 = isolate_1(frominterval,pol);
 std::cout << "Isolating the second algebraic number gives: ["
  << isolation1.first << "," << isolation1.second << "]\n";
 // isolate again the same algebraic number, this time with respect to
 // the polynomial 10*x-14 (which has root 1.4, close to this algebraic
 // number)
 std::pair<Bound,Bound> isolation2 = isolate_1(frominterval,10*x-14);
 std::cout << "Isolating again the second algebraic number gives: ["
  << isolation2.first << "," << isolation2.second << "]\n";
 return 0;
}
#else
int main(){std::cout << "CGALはMPFIライブラリで構築しなければならない" << std::endl;return 0;}
#endif

別の例

#include <CGAL/basic.h>
#ifdef CGAL_USE_MPFI 
#include <CGAL/Algebraic_kernel_d_1.h>
#include <CGAL/Gmpz.h>
#include <vector>
typedef CGAL::Algebraic_kernel_d_1<CGAL::Gmpz> AK;
typedef AK::Polynomial_1 Polynomial_1;
typedef AK::Algebraic_real_1 Algebraic_real_1;
typedef AK::Coefficient Coefficient;
typedef AK::Bound Bound;
typedef AK::Multiplicity_type Multiplicity_type;
int main(){
 AK ak;
 AK::Solve_1 solve_1 = ak.solve_1_object();
 AK::Sign_at_1 sign_at_1 = ak.sign_at_1_object();
 AK::Is_zero_at_1 is_zero_at_1 = ak.is_zero_at_1_object();
 // 多項式p=x^2-5、q=x-2の構築
 Polynomial_1 x = CGAL::shift(AK::Polynomial_1(1),1); // the monomial x
 Polynomial_1 p = x*x-5;
 std::cout << "Polynomial p: " << p << "\n";
 Polynomial_1 q = x-2;
 std::cout << "Polynomial q: " << q << "\n";
 // 根を求める
 std::vector<Algebraic_real_1> roots_p,roots_q;
 solve_1(p,true, std::back_inserter(roots_p));
 solve_1(q,true, std::back_inserter(roots_q));
 // evaluate the second root of p in q
 std::cout << "qの中のpの根2における値の符号: " << sign_at_1(q,roots_p[1]) << "\n";
 // evaluate the root of q in p
 std::cout << "qの中のpの根1における値の符号: " << sign_at_1(p,roots_q[0]) << "\n";
 // check whether the evaluation of the first root of p in p is zero
 std::cout << "pの中のpの根1の値はゼロか: " << is_zero_at_1(p,roots_p[0]) << "\n";
 return 0;
}
#else
int main(){std::cout << "CGALはMPFIライブラリで構築しなければならない" << std::endl;return 0;}
#endif


行列の高速化走査


#include <CGAL/Random.h>
#include <CGAL/Cartesian_matrix.h>
#include <CGAL/sorted_matrix_search.h>
#include <vector>
#include <algorithm>
#include <iterator>
#include <functional>
typedef int                                     Value;
typedef std::vector<Value>                      Vector;
typedef Vector::iterator                        Value_iterator;
typedef std::vector<Vector>                     Vector_cont;
typedef CGAL::Cartesian_matrix<std::plus<int>,Value_iterator,Value_iterator>  Matrix;
int main(){
 // ベクトルを乱数で埋める:
 Vector a;
 const int n = 5;
 for (int i = 0; i < n; ++i) a.push_back(CGAL::default_random(100));
 std::sort(a.begin(), a.end());
 std::cout << "a = ( ";
 std::copy(a.begin(), a.end(), std::ostream_iterator<int>(std::cout," "));
 std::cout << ")\n";

 // ベクトルからデカルト行列を作る
 Matrix M(a.begin(), a.end(), a.begin(), a.end());
 for(int i = 0;i<n;++i){
  std::cout << M(i,0) << ","  << M(i,1) << ","  << M(i,2) << ","  << M(i,3) << ","  << M(i,4) << std::endl;
 }
 // aの最大値の上界(最大値より大きい最小の数)
 Value bound = a[n-1];
 Value upper_bound = CGAL::sorted_matrix_search(
   &M, &M + 1,
   CGAL::sorted_matrix_search_traits_adaptor(std::bind2nd(std::greater_equal<Value>(), bound), M)
 );
 std::cout << "Upper bound for " << bound << " is " << upper_bound << "." << std::endl;
 return 0;
}

線形計画法



// 二次計画法の例
#include <iostream>
#include <cassert>
#include <CGAL/basic.h>
#include <CGAL/QP_models.h>
#include <CGAL/QP_functions.h>
#ifdef CGAL_USE_GMP
#include <CGAL/Gmpz.h>
typedef CGAL::Gmpz ET;
#else
#include <CGAL/MP_Float.h>
typedef CGAL::MP_Float ET;
#endif
typedef CGAL::Quadratic_program<int> Program;
typedef CGAL::Quadratic_program_solution<ET> Solution;

int main() {
 // 非負QP with Ax <= b
 Program qp (CGAL::SMALLER, true, 0, false, 0);
 const int X = 0;
 const int Y = 1;
 qp.set_a(X, 0,  1); qp.set_a(Y, 0, 1); qp.set_b(0, 7);  //  x + y  <= 7のこと
 qp.set_a(X, 1, -1); qp.set_a(Y, 1, 2); qp.set_b(1, 4);  // -x + 2y <= 4のこと
 qp.set_u(Y, true, 4);                                   //       y <= 4のこと
 qp.set_d(X, X, 2); qp.set_d (Y, Y, 8); // !!specify 2D!!    x^2 + 4 y^2のこと
 qp.set_c(Y, -32);                                       // -32yのこと
 qp.set_c0(64);                                          // +64のこと
 Solution s = CGAL::solve_quadratic_program(qp, ET());
 assert (s.solves_quadratic_program(qp));
 std::cout << s;
 return 0;
}