Compare commits

..

25 Commits

Author SHA1 Message Date
Andreas Fabri 543d424f58 Periodic_triangulation_3: Fix warning in demo 2025-12-03 11:36:14 +01:00
Sebastien Loriot 0fc710e95d
Fix bad 'if' condition & warning (#9126)
## Summary of Changes

```
Building CXX object test/Stream_support/CMakeFiles/issue8155.dir/issue8155.cpp.o
cd /home/cgal_tester/build/src/cmake/platforms/ArchLinux-clang-CXX20-Release/test/Stream_support && /bin/clang++ -DCGAL_DATA_DIR=\"/mnt/testsuite/data\" -DCGAL_TEST_SUITE=1 -DCGAL_USE_GMPXX=1 -I/home/cgal_tester/build/src/cmake/platforms/ArchLinux-clang-CXX20-Release/include -I/mnt/testsuite/include -isystem /usr/local/boost/include -Wall -O3 -std=c++20 -DCGAL_NDEBUG -MD -MT test/Stream_support/CMakeFiles/issue8155.dir/issue8155.cpp.o -MF CMakeFiles/issue8155.dir/issue8155.cpp.o.d -o CMakeFiles/issue8155.dir/issue8155.cpp.o -c /mnt/testsuite/test/Stream_support/issue8155.cpp
/mnt/testsuite/test/Stream_support/issue8155.cpp:27:17: warning: overlapping comparisons always evaluate to true [-Wtautological-overlap-compare]
   27 |     if (len > 0 || len != 1.0) {
      |         ~~~~~~~~^~~~~~~~~~~~~
1 warning generated.
```


https://cgal.geometryfactory.com/CGAL/testsuite/CGAL-6.2-Ic-36/Stream_support/TestReport_cgaltest_ArchLinux-clang-CXX20-Release.gz

## Release Management

* Affected package(s): `Stream_support`
* Issue(s) solved (if any): n/a
* Feature/Small Feature (if any): n/a
* License and copyright ownership: no change
2025-11-19 14:43:26 +01:00
Sebastien Loriot dee5ed8cc2
Fix edge collapse with incident non-triangular faces (#9117)
Make collapse able to handle non-triangular faces. The fix is easy as
you simply don't need to join faces in case the face won't disappear
after collapse

**TODO:**
update doc and constrained version
2025-11-19 14:42:28 +01:00
Sebastien Loriot 737bcb264e
Doc: fixing thirdparty links (#9134)
## Summary of Changes

fixed Metis/Eigen3/Qt6 links in Thirdparty
removed LEDA from doc as it is outdated

## Release Management

* Affected package(s): Documentation
2025-11-19 14:41:58 +01:00
Sebastien Loriot 2f9854b422
Fix warning about discarding return type of function with [[nodiscard]] (#9127)
## Summary of Changes

```
Building CXX object test/Generator_Demo/CMakeFiles/Generator_2.dir/Generator_2.cpp.o
cd /home/cgal_tester/build/src/cmake/platforms/ArchLinux-clang/test/Generator_Demo && /bin/clang++ -DCGAL_TEST_SUITE=1 -DCGAL_USE_GMPXX=1 -DQT_CORE_LIB -DQT_GUI_LIB -DQT_NO_KEYWORDS -DQT_OPENGLWIDGETS_LIB -DQT_OPENGL_LIB -DQT_WIDGETS_LIB -I/home/cgal_tester/build/src/cmake/platforms/ArchLinux-clang/test/Generator_Demo -I/mnt/testsuite/test/Generator_Demo -I/home/cgal_tester/build/src/cmake/platforms/ArchLinux-clang/include -I/mnt/testsuite/include -I/home/cgal_tester/build/src/cmake/platforms/ArchLinux-clang/test/AABB_tree_Demo -isystem /home/cgal_tester/build/src/cmake/platforms/ArchLinux-clang/test/Generator_Demo/Generator_2_autogen/include -isystem /usr/include/qt6/QtCore -isystem /usr/include/qt6 -isystem /usr/lib/qt6/mkspecs/linux-g++ -isystem /usr/include/qt6/QtOpenGLWidgets -isystem /usr/include/qt6/QtOpenGL -isystem /usr/include/qt6/QtGui -isystem /usr/include/qt6/QtWidgets -Wall  -fno-direct-access-external-data -MD -MT test/Generator_Demo/CMakeFiles/Generator_2.dir/Generator_2.cpp.o -MF CMakeFiles/Generator_2.dir/Generator_2.cpp.o.d -o CMakeFiles/Generator_2.dir/Generator_2.cpp.o -c /mnt/testsuite/test/Generator_Demo/Generator_2.cpp
In file included from /mnt/testsuite/test/Generator_Demo/Generator_2.cpp:27:
In file included from /mnt/testsuite/include/CGAL/Qt/DemosMainWindow.h:130:
/mnt/testsuite/include/CGAL/Qt/DemosMainWindow_impl.h:224:3: warning: ignoring return value of function declared with 'nodiscard' attribute [-Wunused-result]
  224 |   about_CGAL.open(QIODevice::ReadOnly);
      |   ^~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~
1 warning generated.
```


https://cgal.geometryfactory.com/CGAL/testsuite/CGAL-6.2-Ic-36/Generator_Demo/TestReport_gimeno_ArchLinux-clang.gz

## Release Management

* Affected package(s): `GraphicsView`
* Issue(s) solved (if any): n/a
* Feature/Small Feature (if any): n/a
* License and copyright ownership: no change
2025-11-19 14:41:07 +01:00
Sven Oesau e12a760f03 fixed Metis/Eigen3/Qt6 links in Thirdparty
removed LEDA from doc as it is outdated
2025-11-13 16:32:40 +01:00
Laurent Rineau 12445cc6e1 revert unintended modification 2025-11-13 11:52:21 +01:00
Sebastien Loriot 07347da411
Fix Ceres deprecation warning (#9121)
## Summary of Changes

Changing the used cmake target for Ceres to Ceres::ceres.

## Release Management

* Affected package(s): Installation, Polygon_mesh_processing
2025-11-13 10:02:36 +01:00
Laurent Rineau d226e504c4 fix a warning
use the correct implementation `popupAboutCGAL` from the base class
2025-11-12 16:48:59 +01:00
Laurent Rineau 097fc2c5ab fix memory leaks 2025-11-12 16:48:29 +01:00
Sébastien Loriot 066159f792 assertion should be before the removal... 2025-11-10 10:17:33 +01:00
Mael Rouxel-Labbé 31734df2ef Fix warning about discarding return type of function with [[nodiscard]] 2025-11-07 17:16:31 +01:00
Mael Rouxel-Labbé 3ef43324bb Fix bad 'if' condition & warning 2025-11-07 17:05:21 +01:00
Sébastien Loriot 3c7d507530 update doc + apply changes to the constrained version 2025-11-06 19:40:48 +01:00
Sebastien Loriot faae741666
Eigen3 5.0.0 support (#9112)
## Summary of Changes

Adding support for Eigen3 5.0.0 by updating CMakeLists.txt scripts.

## Release Management

* Issue(s) solved (if any): fix #9110
2025-11-06 10:03:20 +01:00
Sven Oesau 9e36c6744b adding support for Eigen3 5.0.0
moving the Eigen3 version check into CGAL_Eigen3_support.cmake
2025-11-06 10:01:01 +01:00
Sven Oesau 05bf3c4ffc Changing cmake target for Ceres 2025-11-04 17:04:58 +01:00
Sébastien Loriot 25005a97d8 clean up 2025-11-03 09:24:31 +01:00
Sébastien Loriot 3fe83c7ce0 handle non-triangular faces 2025-10-31 11:36:44 +01:00
Sebastien Loriot c60cc5049d
Workaround erros with recent gcc (#9105)
Should fix [those
errors](https://cgal.geometryfactory.com/CGAL/testsuite/CGAL-6.2-Ic-23/Stream_lines_2/TestReport_lrineau_Ubuntu-latest-GCC6-CXX1z.gz)
2025-10-30 20:17:22 +01:00
Sébastien Loriot 712464b690 try to please recent gcc 2025-10-21 08:26:34 +02:00
Sebastien Loriot b08f0a4aae
Handle case of identical projected points (#9100)
Do not try to insert constraint if two points are projected onto the
same 2D vertex
2025-10-17 10:06:08 +02:00
Sébastien Loriot 42068f6009 handle case of identical projected points 2025-10-15 11:03:16 +02:00
Sebastien Loriot 754e57ac3c
[OTR2] Deterministic reconstruction (#9092)
## Summary of Changes

Using timestamps to make iteration through vertices, edges and faces
deterministic

## Release Management

* Affected package(s): Optimal_transportation_reconstruction_2
* Issue(s) solved (if any): fix #8745
2025-10-13 09:33:46 +02:00
Sven Oesau 2f6e3defa7 using timestamps to make iteration through vertices, edges and faces deterministic 2025-10-01 13:42:30 +02:00
2558 changed files with 39019 additions and 219397 deletions

View File

@ -1,35 +0,0 @@
---
Language: Cpp
BasedOnStyle: LLVM
AccessModifierOffset: -2
AllowShortFunctionsOnASingleLine: true
BinPackParameters: false
BreakConstructorInitializers: BeforeComma
BreakBeforeBraces: Custom
BraceWrapping:
AfterCaseLabel: false
AfterClass: true
AfterControlStatement: MultiLine
AfterEnum: false
AfterFunction: false
AfterNamespace: false
AfterObjCDeclaration: false
AfterStruct: true
AfterUnion: false
AfterExternBlock: false
BeforeCatch: false
BeforeElse: false
BeforeLambdaBody: false
BeforeWhile: false
IndentBraces: false
SplitEmptyFunction: false
SplitEmptyRecord: false
SplitEmptyNamespace: false
ColumnLimit: 120
# Force pointers to the type for C++.
DerivePointerAlignment: false
PointerAlignment: Left
# Control the spaces around conditionals
SpacesInConditionalStatement: false
SpaceBeforeParens: false
...

2
.gitattributes vendored
View File

@ -16,7 +16,7 @@
*.py text
*.xml text
*.js text
*.html text
*.hmtl text
*.bib text
*.css text
*.ui text

View File

@ -1,11 +1,3 @@
---
name: Issue
about: Use this template for reporting issues with CGAL
title: ''
labels: ''
assignees: ''
---
_Please use the following template to help us solving your issue._
## Issue Details

View File

@ -1 +0,0 @@
blank_issues_enabled: false

View File

@ -29,7 +29,6 @@ permissions:
jobs:
pre_build_checks:
if: ${{ inputs.pr_number || startsWith(github.event.comment.body, '/build:') || startsWith(github.event.comment.body, '/force-build:') }}
runs-on: ubuntu-latest
name: Trigger the build?
outputs:

View File

@ -7,8 +7,9 @@ permissions:
jobs:
build:
if: github.repository == 'CGAL/cgal' || github.event_name != 'push'
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- name: install dependencies

View File

@ -7,8 +7,9 @@ permissions:
jobs:
cmake-testsuite:
if: github.repository == 'CGAL/cgal' || github.event_name != 'push'
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- name: install dependencies
@ -20,8 +21,9 @@ jobs:
ctest -L Installation -j $(getconf _NPROCESSORS_ONLN)
cmake-testsuite-with-qt:
if: github.repository == 'CGAL/cgal' || github.event_name != 'push'
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- name: install dependencies

View File

@ -10,18 +10,17 @@ jobs:
reuse:
runs-on: ubuntu-latest
steps:
- name: Checkout
uses: actions/checkout@v4
- name: Display reuse-tool version
uses: fsfe/reuse-action@v5
- uses: actions/checkout@v4
- name: REUSE version
uses: fsfe/reuse-action@v4
with:
args: --version
- name: REUSE Compliance Check
uses: fsfe/reuse-action@v5
- name: REUSE lint
uses: fsfe/reuse-action@v4
with:
args: --include-submodules lint
- name: REUSE SPDX SBOM
uses: fsfe/reuse-action@v5
uses: fsfe/reuse-action@v4
with:
args: spdx
- name: install dependencies
@ -30,7 +29,7 @@ jobs:
run: |
mkdir -p ./release
cmake -DDESTINATION=./release -DCGAL_VERSION=9.9 -P ./Scripts/developer_scripts/cgal_create_release_with_cmake.cmake
- name: REUSE Compliance Check of release tarball
uses: fsfe/reuse-action@v5
- name: REUSE lint release tarball
uses: fsfe/reuse-action@v4
with:
args: --root ./release/CGAL-9.9 --include-submodules lint

25
.gitignore vendored
View File

@ -1,4 +1,3 @@
build
/*build*
/*/*/*/build
/*/*/*/VC*
@ -1148,27 +1147,3 @@ Polygonal_surface_reconstruction/examples/build*
Polygonal_surface_reconstruction/test/build*
Solver_interface/examples/build*
/Mesh_3/examples/Mesh_3/indicator_0.inr.gz
/*.off
/*.xyz
/r0*
all_segments.polylines.txt
Data/data/meshes/*.*.edge
Data/data/meshes/*.*.ele
Data/data/meshes/*.*.face
Data/data/meshes/*.*.mesh
Data/data/meshes/*.*.node
Data/data/meshes/*.*.smesh
Data/data/meshes/*.*.vtk
Data/data/meshes/*.log
Data/data/meshes/*.off-cdt-output.off
dump_*.off
dump_*.txt
dump-*.binary.cgal
dump-*.polylines.txt
dump-*.xyz
dump.off.mesh
log.txt
patches_after_merge.ply
CMakeUserPresets.json
/.cache
compile_commands.json

View File

@ -1,7 +0,0 @@
{
"default": true,
"line-length": false,
"no-duplicate-heading": {
"siblings_only": true
}
}

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(AABB_traits_benchmark)
find_package(CGAL REQUIRED OPTIONAL_COMPONENTS Core)
@ -13,7 +13,7 @@ create_single_source_cgal_program("tree_construction.cpp")
find_package(benchmark QUIET)
if(benchmark_FOUND)
create_single_source_cgal_program("tree_creation.cpp")
target_link_libraries(tree_creation PRIVATE benchmark::benchmark)
target_link_libraries(tree_creation benchmark::benchmark)
else()
message(STATUS "NOTICE: The benchmark 'tree_creation.cpp' requires the Google benchmark library, and will not be compiled.")
endif()

View File

@ -1,6 +1,6 @@
# This is the CMake script for compiling the AABB tree demo.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(AABB_tree_Demo)
# Find includes in corresponding build directories
@ -14,7 +14,7 @@ find_package(Qt6 QUIET COMPONENTS Gui OpenGL)
if(CGAL_Qt6_FOUND AND Qt6_FOUND)
add_compile_definitions(QT_NO_KEYWORDS)
add_definitions(-DQT_NO_KEYWORDS)
# Instruct CMake to run moc/ui/rcc automatically when needed.
set(CMAKE_AUTOMOC ON)

View File

@ -1,4 +1,4 @@
/// \defgroup PkgAABBTreeRef Reference Manual
/// \defgroup PkgAABBTreeRef AABB Tree Reference
/// \defgroup PkgAABBTreeConcepts Concepts
/// \ingroup PkgAABBTreeRef

View File

@ -106,7 +106,7 @@ it is used only for distance computations.
\warning Having degenerate primitives in the AABB-tree is not recommended as the underlying
predicates and constructions of the traits class might not be able to handle them.
For example if one is using `CGAL::AABB_traits` with a Kernel from \cgal,
having degenerate triangles or segments in the AABB-tree will result in an undefined
having degenerate triangles or segments in the AABB-tree will results in an undefined
behavior or a crash.
\section aabb_tree_examples Examples
@ -197,9 +197,9 @@ The AABB tree example folder contains three examples of trees
constructed with custom primitives. In \ref AABB_tree/AABB_custom_example.cpp "AABB_custom_example.cpp"
the primitive contains triangles which are defined by three pointers
to custom points. In \ref AABB_tree/AABB_custom_triangle_soup_example.cpp "AABB_custom_triangle_soup_example.cpp" all input
triangles are stored into a single array as to form a triangle
soup. The primitive internally uses a `boost::iterator_adaptor`
as to provide the three functions `AABBPrimitive::id()`, `AABBPrimitive::datum()`,
triangles are stored into a single array so as to form a triangle
soup. The primitive internally uses a `boost::iterator_adaptor` so as
to provide the three functions `AABBPrimitive::id()`, `AABBPrimitive::datum()`,
and `AABBPrimitive::reference_point()`) required by the primitive concept. In
\ref AABB_tree/AABB_custom_indexed_triangle_set_example.cpp "AABB_custom_indexed_triangle_set_example.cpp" the input is an indexed
triangle set stored through two arrays: one array of points and one
@ -213,7 +213,7 @@ contains a set of polyhedron triangle facets. We measure the tree
construction time, the memory occupancy and the number of queries per
second for a variety of intersection and distance queries. The machine
used is a PC running Windows XP64 with an Intel CPU Core2 Extreme
clocked at 3.06 GHz with 4GB of RAM. By default, the kernel used is
clocked at 3.06 GHz with 4GB of RAM. By default the kernel used is
`Simple_cartesian<double>` (the fastest in our experiments). The
program has been compiled with Visual C++ 2005 compiler with the O2
option which maximizes speed.
@ -244,7 +244,7 @@ internal KD-tree). It increases to approximately 150 bytes per
primitive when constructing the internal KD-tree with one reference
point per primitive (the default mode when calling the function
`AABB_tree::accelerate_distance_queries()`). Note that the polyhedron
facet primitive stores only one facet handle as primitive id
facet primitive primitive stores only one facet handle as primitive id
and computes on the fly a 3D triangle from the facet handle stored
internally. When explicitly storing a 3D triangle in the primitive the
tree occupies approximately 140 bytes per primitive instead of 60
@ -340,7 +340,7 @@ inside the bounding box.
The experiments described above are neither exhaustive nor conclusive
as we have chosen one specific case where the input primitives are the
facets of a triangle surface polyhedron. Nevertheless, we now provide
facets of a triangle surface polyhedron. Nevertheless we now provide
some general observations and advises about how to put the AABB tree
to use with satisfactory performances. While the tree construction
times and memory occupancy do not fluctuate much in our experiments

View File

@ -27,7 +27,7 @@ void triangle_mesh(std::string fname)
typedef CGAL::AABB_tree<Traits> Tree;
TriangleMesh tmesh;
if(!CGAL::IO::read_polygon_mesh(fname, tmesh) || !CGAL::is_triangle_mesh(tmesh))
if(!CGAL::IO::read_polygon_mesh(fname, tmesh) || CGAL::is_triangle_mesh(tmesh))
{
std::cerr << "Invalid input." << std::endl;
return;

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(AABB_tree_Examples)
find_package(CGAL REQUIRED)

View File

@ -302,11 +302,12 @@ public:
typename AT::Bounding_box operator()(ConstPrimitiveIterator first,
ConstPrimitiveIterator beyond) const
{
return std::accumulate(first, beyond,
typename AT::Bounding_box{} /* empty bbox */,
[this](const typename AT::Bounding_box& bbox, const Primitive& pr) {
return bbox + m_traits.compute_bbox(pr, m_traits.bbm);
});
typename AT::Bounding_box bbox = m_traits.compute_bbox(*first,m_traits.bbm);
for(++first; first != beyond; ++first)
{
bbox = bbox + m_traits.compute_bbox(*first,m_traits.bbm);
}
return bbox;
}
};

View File

@ -199,7 +199,9 @@ namespace CGAL {
}
/// returns the axis-aligned bounding box of the whole tree.
/// \pre `!empty()`
const Bounding_box bbox() const {
CGAL_precondition(!empty());
if(size() > 1)
return root_node()->bbox();
else

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(AABB_tree_Tests)
find_package(CGAL REQUIRED)

View File

@ -49,7 +49,7 @@ We describe next the algorithm and provide examples.
\note A \ref tuto_reconstruction "detailed tutorial on surface reconstruction"
is provided with a guide to choose the most appropriate method along
with pre- and postprocessing.
with pre- and post-processing.
\section AFSR_Definitions Definitions and the Algorithm

View File

@ -1,4 +1,4 @@
/// \defgroup PkgAdvancingFrontSurfaceReconstructionRef Reference Manual
/// \defgroup PkgAdvancingFrontSurfaceReconstructionRef Advancing Front Surface Reconstruction Reference
/// \defgroup PkgAdvancingFrontSurfaceReconstructionRefConcepts Concepts
/// \ingroup PkgAdvancingFrontSurfaceReconstructionRef

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Advancing_front_surface_reconstruction_Examples)
find_package(CGAL REQUIRED)

View File

@ -382,7 +382,7 @@ namespace CGAL {
int _facet_number;
//---------------------------------------------------------------------
// For postprocessing
// For post-processing
mutable int _postprocessing_counter;
int _size_before_postprocessing;
@ -2432,7 +2432,7 @@ namespace CGAL {
std::size_t itmp, L_v_size_mem;
L_v_size_mem = L_v.size();
if ((vh_on_border_inserted != 0)&& // to postprocess only the borders
if ((vh_on_border_inserted != 0)&& // to post-process only the borders
(L_v.size() < .1 * _size_before_postprocessing))
{
{

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Advancing_front_surface_reconstruction_Tests)
find_package(CGAL REQUIRED)

View File

@ -15,7 +15,7 @@ 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 <I>number type</I>. Instead, the package distinguishes
the package avoids the term <I>number type</I>. Instead the package distinguishes
between the <I>algebraic structure</I> of a type and whether a type is embeddable on
the real axis, or <I>real embeddable</I> for short.
Moreover, the package introduces the notion of <I>interoperable</I> types which
@ -45,11 +45,11 @@ 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
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
`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
@ -80,12 +80,12 @@ 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., \link sqrt `CGAL::sqrt(x)` \endlink. This is realized via function templates using
the according functor of the traits class. For an overview see the section "Global Functions" in the
\ref PkgAlgebraicFoundationsRef.
the according functor of the traits class. For an overview see
Section \ref PkgAlgebraicFoundationsRef in the reference manual.
\subsection Algebraic_foundationsTagsinAlgebraicStructure Tags in Algebraic Structure Traits
\subsubsection Algebraic_foundationsAlgebraicCategory Algebraic Category
\subsection Algebraic_foundationsAlgebraicCategory Algebraic Category
For a type `AS`, `Algebraic_structure_traits<AS>`
provides several tags. The most important tag is the `Algebraic_category`
@ -100,7 +100,7 @@ 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`.
\subsubsection Algebraic_foundationsExactandNumericalSensitive Exact and Numerical Sensitive
\subsection Algebraic_foundationsExactandNumericalSensitive Exact and Numerical Sensitive
Moreover, `Algebraic_structure_traits<AS>` provides the tags `Is_exact`
and `Is_numerical_sensitive`, which are both `Boolean_tag`s.
@ -119,7 +119,7 @@ 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
checked for exact types only.
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
@ -134,7 +134,7 @@ category tags reflect the algebraic structure hierarchy.
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
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,
@ -190,12 +190,12 @@ 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.
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 \link Coercion_traits::Type `Coercion_traits<A,B>::Type` \endlink
for two interoperable types `A` and `B`.
Some trivial examples are `int` and `double` with coercion type double
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
@ -269,7 +269,7 @@ 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 \cgalCite{beh-eeeafcs-05},
which in the main goes back to
Lutz Kettner, Susan Hert, Arno Eigenwillig, and Michael Hemmer.
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

View File

@ -1,10 +1,11 @@
/// \defgroup PkgAlgebraicFoundationsRef Reference Manual
/// \defgroup PkgAlgebraicFoundationsRef Algebraic Foundations Reference
/// \defgroup PkgAlgebraicFoundationsAlgebraicStructuresConcepts Concepts
/// \ingroup PkgAlgebraicFoundationsRef
/*!
\addtogroup PkgAlgebraicFoundationsRef
\todo check generated documentation
\cgalPkgDescriptionBegin{Algebraic Foundations,PkgAlgebraicFoundations}
\cgalPkgPicture{Algebraic_foundations2.png}

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Algebraic_foundations_Examples)
find_package(CGAL REQUIRED)

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Algebraic_foundations_Tests)
find_package(CGAL REQUIRED COMPONENTS Core)

View File

@ -10,7 +10,7 @@ namespace CGAL {
\section Algebraic_kernel_dIntroduction Introduction
Real solving of polynomials is a fundamental problem with a wide application range.
This package provides black-box implementations of
This package is targeted at providing black-box implementations of state-of-the-art
algorithms to determine, compare, and approximate real roots of univariate polynomials
and bivariate polynomial systems. Such a black-box is called an *Algebraic Kernel*.
Since this package is aimed at providing more than one implementation, the interface of
@ -146,8 +146,8 @@ polynomial is not computed yet.
This is in particular relevant in relation to the `AlgebraicKernel_d_2`,
where `AlgebraicKernel_d_1::Algebraic_real_1` is used to represent coordinates of solutions of bivariate systems.
Hence, the design does
not allow direct access to any, seemingly obvious, members of an
`AlgebraicKernel_d_1::Algebraic_real_1`. Instead, there is, e.g.,
not allow a direct access to any, seemingly obvious, members of an
`AlgebraicKernel_d_1::Algebraic_real_1`. Instead there is, e.g.,
`AlgebraicKernel_d_1::Compute_polynomial_1` which emphasizes
that the requested polynomial may not be computed yet. Similarly,
there is no way to directly ask for the refinement of the current
@ -283,6 +283,41 @@ is not stored internally in terms of an `Algebraic_real_1` object.
Querying such a representation by calling `Compute_y_2` is a
time-consuming step, and should be avoided for efficiency reasons if possible.
\subsection Algebraic_kernel_dAlgebraicKernelsBasedon Algebraic Kernels Based on RS
The package offers two univariate algebraic kernels that are based on
the library RS \cgalCite{cgal:r-rs}, namely `Algebraic_kernel_rs_gmpz_d_1`
and `Algebraic_kernel_rs_gmpq_d_1`. As the names indicate,
the kernels are based on the library RS \cgalCite{cgal:r-rs} and support univariate
polynomials over `Gmpz` or `Gmpq`, respectively.
In general we encourage to use `Algebraic_kernel_rs_gmpz_d_1`
instead of `Algebraic_kernel_rs_gmpq_d_1`. This is caused by
the fact that the most efficient way to compute operations (such as gcd)
on polynomials with rational coefficients is to use the corresponding
implementation for polynomials with integer coefficients. That is,
the `Algebraic_kernel_rs_gmpq_d_1` is slightly slower due to
overhead caused by the necessary conversions. However, since this may
not always be a major issue, the `Algebraic_kernel_rs_gmpq_d_1`
is provided for convenience.
The core of both kernels is the implementation of the interval Descartes
algorithm \cgalCite{cgal:rz-jcam-04} of the library RS \cgalCite{cgal:r-rs},
which is used to isolate the roots of the polynomial.
The RS library restricts its attention to univariate integer
polynomials and some substantial gain of efficiency can be made by using a kernel
that does not follow the generic programming paradigm, by avoiding
interfaces between layers. Specifically, working with
only one number type allows to optimize some polynomial operations
as well as memory handling. The implementation of these kernels
make heavy use of the \mpfr \cgalCite{cgal:mt-mpfr} and \mpfi \cgalCite{cgal:r-mpfi}
libraries, and of their \cgal interfaces, `Gmpfr` and `Gmpfi`.
The algebraic numbers (roots of the polynomials) are represented
in the two RS kernels by a `Gmpfi` interval and a pointer to
the polynomial of which they are roots. See \cgalCite{cgal:lpt-wea-09}
for more details on the implementation, tests of these kernels,
comparisons with other algebraic kernels and discussions about the
efficiency.
\section Algebraic_kernel_dExamples Examples
@ -341,13 +376,12 @@ Michael Hemmer and Michael Kerber, respectively. Notwithstanding, the authors al
contribution of all authors of the \exacus project,
particularly the contribution of Arno Eigenwillig, Sebastian Limbach and Pavel Emeliyanenko.
In 2010, two univariate kernels that interface the library RS \cgalCite{cgal:r-rs} were
written by Luis Peñaranda and Sylvain Lazard \cgalCite{cgal:2009-ESA}.
The two univariate kernels that interface the library RS \cgalCite{cgal:r-rs} were
written by Luis Pe&ntilde;aranda and Sylvain Lazard.
Both models interface the library RS \cgalCite{cgal:r-rs} by Fabrice Rouillier.
The authors want to thank Fabrice Rouillier and Elias Tsigaridas for
strong support and many useful discussions that lead to the integration of RS.
Due to lack of maintenance, these kernels have been removed in 2024.
*/
} /* namespace CGAL */

View File

@ -0,0 +1,59 @@
namespace CGAL {
/*!
\ingroup PkgAlgebraicKernelDModels
\anchor Algebraic_kernel_rs_gmpq_d_1
This univariate algebraic kernel uses the Rs library to perform
rational univariate polynomial root isolation. It is a model of the
`AlgebraicKernel_d_1` concept. Due to the fact that RS can only
isolate integer polynomials, the operations of this kernel have the
overhead of converting the polynomials to integer.
\cgalModels{AlgebraicKernel_d_1}
\sa `Algebraic_kernel_rs_gmpz_d_1`
*/
class Algebraic_kernel_rs_gmpq_d_1 {
public:
/// \name Types
/// @{
/*!
It is a typedef to `CGAL::Gmpq`.
*/
typedef unspecified_type Coefficient;
/*!
It is defined as `CGAL::Polynomial<CGAL::Gmpq>`.
*/
typedef unspecified_type Polynomial_1;
/*!
Type that represents the real roots of
integer univariate polynomials, containing a pointer to the polynomial of
which the represented algebraic number is root and and a `CGAL::Gmpfi`
isolating interval.
*/
typedef unspecified_type Algebraic_real_1;
/*!
Since the isolating intervals of the roots have type
`CGAL::Gmpfi`, the bounds have type `CGAL::Gmpfr`.
*/
typedef unspecified_type Bound;
/*!
The multiplicity is an `int`.
*/
typedef unspecified_type Multiplicity_type;
/// @}
}; /* end Algebraic_kernel_rs_gmpq_d_1 */
} /* end namespace CGAL */

View File

@ -0,0 +1,57 @@
namespace CGAL {
/*!
\ingroup PkgAlgebraicKernelDModels
\anchor Algebraic_kernel_rs_gmpz_d_1
This univariate algebraic kernel uses the Rs library to perform
integer univariate polynomial root isolation. It is a model of the
`AlgebraicKernel_d_1` concept.
\cgalModels{AlgebraicKernel_d_1}
\sa `Algebraic_kernel_rs_gmpz_d_1`
*/
class Algebraic_kernel_rs_gmpz_d_1 {
public:
/// \name Types
/// @{
/*!
It is a typedef to `CGAL::Gmpz`.
*/
typedef unspecified_type Coefficient;
/*!
It is defined as `CGAL::Polynomial<CGAL::Gmpz>`.
*/
typedef unspecified_type Polynomial_1;
/*!
Type that represents the real roots of
integer univariate polynomials, containing a pointer to the polynomial of
which the represented algebraic number is root and and a `CGAL::Gmpfi`
isolating interval.
*/
typedef unspecified_type Algebraic_real_1;
/*!
Since the isolating intervals of the roots have type
`CGAL::Gmpfi`, the bounds have type `CGAL::Gmpfr`.
*/
typedef unspecified_type Bound;
/*!
The multiplicity is an `int`.
*/
typedef unspecified_type Multiplicity_type;
/// @}
}; /* end Algebraic_kernel_rs_gmpz_d_1 */
} /* end namespace CGAL */

View File

@ -9,7 +9,8 @@ algebraic functionalities on univariate polynomials of general degree \f$ d\f$.
\cgalRefines{CopyConstructible,Assignable}
\cgalHasModelsBegin
\cgalHasModels{CGAL::Algebraic_kernel_d_1}
\cgalHasModels{CGAL::Algebraic_kernel_rs_gmpz_d_1}
\cgalHasModels{CGAL::Algebraic_kernel_rs_gmpq_d_1}
\cgalHasModelsEnd
\sa `AlgebraicKernel_d_2`
@ -171,3 +172,4 @@ AlgebraicKernel_d_1::Bound_between_1 bound_between_1_object() const;
/// @}
}; /* end AlgebraicKernel_d_1 */

View File

@ -8,11 +8,6 @@ for solving and handling bivariate polynomial systems of general degree \f$ d\f$
\cgalRefines{AlgebraicKernel_d_1,CopyConstructible,Assignable}
\cgalHasModelsBegin
\cgalHasModels{CGAL::Algebraic_kernel_d_2}
\cgalHasModelsEnd
\sa `AlgebraicKernel_d_1`
*/
@ -181,3 +176,4 @@ AlgebraicKernel_d_2::Bound_between_x_2 bound_between_x_2_object() const;
/// @}
}; /* end AlgebraicKernel_d_2 */

View File

@ -1,4 +1,4 @@
/// \defgroup PkgAlgebraicKernelDRef Reference Manual
/// \defgroup PkgAlgebraicKernelDRef Algebraic Kernel Reference
/// \defgroup PkgAlgebraicKernelDConcepts Concepts
/// \ingroup PkgAlgebraicKernelDRef
@ -16,15 +16,17 @@
/*!
\addtogroup PkgAlgebraicKernelDRef
\todo check generated documentation
\cgalPkgDescriptionBegin{Algebraic Kernel,PkgAlgebraicKernelD}
\cgalPkgPicture{Algebraic_kernel_d.png}
\cgalPkgSummaryBegin
\cgalPkgAuthors{Eric Berberich, Michael Hemmer, Michael Kerber, Sylvain Lazard, Luis Peñaranda, and Monique Teillaud}
\cgalPkgDesc{Real solving of polynomials is a fundamental problem with a wide application range. This package is targeted to provide black-box implementation algorithms to determine, compare and approximate real roots of univariate polynomials and bivariate polynomial systems. Such a black-box is called an *Algebraic %Kernel*. So far the package only provides models for the univariate kernel. Nevertheless, it already defines concepts for the bivariate kernel, since this settles the interface for upcoming implementations.}
\cgalPkgDesc{Real solving of polynomials is a fundamental problem with a wide application range. This package is targeted to provide black-box implementations of state-of-the-art algorithms to determine, compare and approximate real roots of univariate polynomials and bivariate polynomial systems. Such a black-box is called an *Algebraic %Kernel*. So far the package only provides models for the univariate kernel. Nevertheless, it already defines concepts for the bivariate kernel, since this settles the interface for upcoming implementations.}
\cgalPkgManuals{Chapter_Algebraic_Kernel,PkgAlgebraicKernelDRef}
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{3.6}
\cgalPkgDependsOn{Some models depend on \ref thirdpartyRS3.}
\cgalPkgBib{cgal:bht-ak}
\cgalPkgLicense{\ref licensesLGPL "LGPL"}
\cgalPkgShortInfoEnd
@ -98,4 +100,8 @@
- `CGAL::Algebraic_kernel_d_1<Coeff>`
- `CGAL::Algebraic_kernel_d_2<Coeff>`
- `CGAL::Algebraic_kernel_rs_gmpz_d_1`
- `CGAL::Algebraic_kernel_rs_gmpq_d_1`
*/

View File

@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Algebraic_kernel_d_Examples)
find_package(CGAL REQUIRED COMPONENTS Core)

View File

@ -520,7 +520,7 @@ private:
void set_event_lines(InputIterator1 event_begin,
InputIterator1 event_end,
InputIterator2 intermediate_begin,
InputIterator2 CGAL_assertion_code(intermediate_end)) const {
InputIterator2 CGAL_precondition_code(intermediate_end)) const {
if(! this->ptr()->event_coordinates) {

View File

@ -875,7 +875,7 @@ private:
* Checks intersection with symbolic methods
*/
bool check_candidate_symbolically(Status_line_CA_1& e1,size_type ,
Status_line_CA_1& CGAL_assertion_code(e2),size_type ,
Status_line_CA_1& CGAL_precondition_code(e2),size_type ,
size_type k) const {
Polynomial_1 p = -coprincipal_subresultants(k-1);
Polynomial_1 q = principal_subresultants(k)*Coefficient(k);
@ -1173,7 +1173,7 @@ private:
/*
* \brief reduces the number of possible intersections
*
* At the position given by the event lines \c e1 and \c e2 and the slice
* At the position given by the event lins \c e1 and \c e2 and the slice
* info object \c slice, the points on the event lines are further refined
* until there are only \c n possible intersection points. The method can
* be interrupted if all possible intersection points are known to have

View File

@ -65,7 +65,7 @@ template < class PolynomialIterator,
int gen_agebraic_reals_with_mults( PolynomialIterator fac,
PolynomialIterator fac_end,
IntIterator mul,
IntIterator CGAL_assertion_code(mul_end),
IntIterator CGAL_precondition_code(mul_end),
AlgebraicRealOutputIterator oi_root,
IntOutputIterator oi_mult){

View File

@ -331,7 +331,7 @@ public:
const typename Curve_pair_analysis_2
::Curve_analysis_2& c1,
const typename Curve_pair_analysis_2
::Curve_analysis_2& CGAL_assertion_code(c2)) const
::Curve_analysis_2& CGAL_precondition_code(c2)) const
{
CGAL_precondition(0 <= j && j < number_of_events());

View File

@ -0,0 +1,114 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_ALGEBRAIC_KERNEL_RS_GMPQ_D_1
#define CGAL_ALGEBRAIC_KERNEL_RS_GMPQ_D_1
#include <CGAL/disable_warnings.h>
#include <CGAL/Gmpq.h>
#include <CGAL/Gmpfr.h>
#include <CGAL/Polynomial.h>
#include "RS/rs2_isolator_1.h"
#ifdef CGAL_USE_RS3
#include "RS/rs23_k_isolator_1.h"
#include "RS/rs3_refiner_1.h"
#include "RS/rs3_k_refiner_1.h"
#else
#include "RS/bisection_refiner_1.h"
#endif
#ifdef CGAL_RS_USE_NATIVE_GMPQ_KERNEL
#include "RS/ak_1.h"
#else
#include "RS/ak_z_1.h"
#endif
// The RS algebraic kernel for non-Gmpz types comes in two flavors. The
// native kernel converts, before each operation, the input polynomial to a
// Polynomial<Gmpz>. The z-kernel only converts to a Polynomial<Gmpz>
// before isolation and stores both polynomials in the algebraic number.
// The latter is the default behavior, but the former can be selected by
// setting the compilation flag CGAL_RS_USE_NATIVE_GMPQ_KERNEL.
namespace CGAL{
#ifdef CGAL_RS_USE_NATIVE_GMPQ_KERNEL // Use the native kernel.
// Choice of the isolator: RS default or RS-k.
#ifdef CGAL_RS_USE_K
typedef CGAL::RS23_k_isolator_1<CGAL::Polynomial<CGAL::Gmpq>,CGAL::Gmpfr>
QIsolator;
#else
typedef CGAL::RS2::RS2_isolator_1<CGAL::Polynomial<CGAL::Gmpq>,CGAL::Gmpfr>
QIsolator;
#endif
// Choice of the refiner: bisection, RS3 or RS3-k.
#ifdef CGAL_USE_RS3
#ifdef CGAL_RS_USE_K
typedef CGAL::RS3::RS3_k_refiner_1<CGAL::Polynomial<CGAL::Gmpq>,CGAL::Gmpfr>
QRefiner;
#else
typedef CGAL::RS3::RS3_refiner_1<CGAL::Polynomial<CGAL::Gmpq>,CGAL::Gmpfr>
QRefiner;
#endif
#else
typedef CGAL::Bisection_refiner_1<CGAL::Polynomial<CGAL::Gmpq>,CGAL::Gmpfr>
QRefiner;
#endif
typedef CGAL::RS_AK1::Algebraic_kernel_1<
CGAL::Polynomial<CGAL::Gmpq>,
CGAL::Gmpfr,
QIsolator,
QRefiner>
Algebraic_kernel_rs_gmpq_d_1;
#else // Use the z-kernel.
// Choice of the z-isolator: RS default or RS-k.
#ifdef CGAL_RS_USE_K
typedef CGAL::RS23_k_isolator_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
ZIsolator;
#else
typedef CGAL::RS2::RS2_isolator_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
ZIsolator;
#endif
// Choice of the z-refiner: bisection, RS3 or RS3-k.
#ifdef CGAL_USE_RS3
#ifdef CGAL_RS_USE_K
typedef CGAL::RS3::RS3_k_refiner_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
ZRefiner;
#else
typedef CGAL::RS3::RS3_refiner_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
ZRefiner;
#endif
#else
typedef CGAL::Bisection_refiner_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
ZRefiner;
#endif
typedef CGAL::RS_AK1::Algebraic_kernel_z_1<
CGAL::Polynomial<CGAL::Gmpq>,
CGAL::Polynomial<CGAL::Gmpz>,
CGAL::RS_AK1::Polynomial_converter_1<CGAL::Polynomial<CGAL::Gmpq>,
CGAL::Polynomial<CGAL::Gmpz> >,
CGAL::Gmpfr,
ZIsolator,
ZRefiner>
Algebraic_kernel_rs_gmpq_d_1;
#endif
} // namespace CGAL
#include <CGAL/enable_warnings.h>
#endif // CGAL_ALGEBRAIC_KERNEL_RS_GMPQ_D_1

View File

@ -0,0 +1,65 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_ALGEBRAIC_KERNEL_RS_GMPZ_D_1
#define CGAL_ALGEBRAIC_KERNEL_RS_GMPZ_D_1
#include <CGAL/disable_warnings.h>
#include <CGAL/Gmpz.h>
#include <CGAL/Gmpfr.h>
#include <CGAL/Polynomial.h>
#include "RS/rs2_isolator_1.h"
#ifdef CGAL_USE_RS3
#include "RS/rs23_k_isolator_1.h"
#include "RS/rs3_refiner_1.h"
#include "RS/rs3_k_refiner_1.h"
#else
#include "RS/bisection_refiner_1.h"
#endif
#include "RS/ak_1.h"
namespace CGAL{
// Choice of the isolator: RS default or RS-k.
#ifdef CGAL_RS_USE_K
typedef CGAL::RS23_k_isolator_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
Isolator;
#else
typedef CGAL::RS2::RS2_isolator_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
Isolator;
#endif
// Choice of the refiner: bisection, RS3 or RS3-k.
#ifdef CGAL_USE_RS3
#ifdef CGAL_RS_USE_K
typedef CGAL::RS3::RS3_k_refiner_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
Refiner;
#else
typedef CGAL::RS3::RS3_refiner_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
Refiner;
#endif
#else
typedef CGAL::Bisection_refiner_1<CGAL::Polynomial<CGAL::Gmpz>,CGAL::Gmpfr>
Refiner;
#endif
typedef CGAL::RS_AK1::Algebraic_kernel_1<
CGAL::Polynomial<CGAL::Gmpz>,
CGAL::Gmpfr,
Isolator,
Refiner>
Algebraic_kernel_rs_gmpz_d_1;
}
#include <CGAL/enable_warnings.h>
#endif // CGAL_ALGEBRAIC_KERNEL_RS_GMPZ_D_1

View File

@ -0,0 +1,30 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_GMPFR_MAKE_UNIQUE_H
#define CGAL_RS_GMPFR_MAKE_UNIQUE_H
#include <CGAL/Gmpfr.h>
// Make sure a number does not share references. If it does, copy it.
#ifdef CGAL_GMPFR_NO_REFCOUNT
# define CGAL_RS_GMPFR_MAKE_UNIQUE(_number,_tempvar) {};
#else // CGAL_GMPFR_NO_REFCOUNT
# define CGAL_RS_GMPFR_MAKE_UNIQUE(_number,_tempvar) \
if(!_number.is_unique()){ \
Gmpfr _tempvar(0,_number.get_precision()); \
mpfr_set(_tempvar.fr(),_number.fr(),GMP_RNDN); \
_number=_tempvar; \
CGAL_assertion_code(_tempvar=Gmpfr();) \
CGAL_assertion(_number.is_unique()); \
}
#endif // CGAL_GMPFR_NO_REFCOUNT
#endif // CGAL_RS_GMPFR_MAKE_UNIQUE_H

View File

@ -0,0 +1,172 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_AK_1_H
#define CGAL_RS_AK_1_H
#include <cstddef> // included only to define size_t
#include <CGAL/Polynomial_traits_d.h>
#include "algebraic_1.h"
#include "comparator_1.h"
#include "signat_1.h"
#include "functors_1.h"
namespace CGAL{
namespace RS_AK1{
template <class Polynomial_,
class Bound_,
class Isolator_,
class Refiner_,
class Ptraits_=CGAL::Polynomial_traits_d<Polynomial_> >
class Algebraic_kernel_1{
public:
typedef Polynomial_ Polynomial_1;
typedef typename Polynomial_1::NT Coefficient;
typedef Bound_ Bound;
private:
typedef Isolator_ Isolator;
typedef Refiner_ Refiner;
typedef Ptraits_ Ptraits;
typedef CGAL::RS_AK1::Signat_1<Polynomial_1,Bound>
Signat;
typedef CGAL::RS_AK1::Simple_comparator_1<Polynomial_1,
Bound,
Refiner,
Signat,
Ptraits>
Comparator;
public:
typedef CGAL::RS_AK1::Algebraic_1<Polynomial_1,
Bound,
Refiner,
Comparator,
Ptraits>
Algebraic_real_1;
typedef size_t size_type;
typedef unsigned Multiplicity_type;
// default constructor and destructor
public:
Algebraic_kernel_1(){};
~Algebraic_kernel_1(){};
// functors from the CGAL concept
public:
typedef CGAL::RS_AK1::Construct_algebraic_real_1<Polynomial_1,
Algebraic_real_1,
Bound,
Coefficient,
Isolator>
Construct_algebraic_real_1;
typedef CGAL::RS_AK1::Compute_polynomial_1<Polynomial_1,
Algebraic_real_1>
Compute_polynomial_1;
typedef CGAL::RS_AK1::Isolate_1<Polynomial_1,
Bound,
Algebraic_real_1,
Isolator,
Comparator,
Signat,
Ptraits>
Isolate_1;
typedef typename Ptraits::Is_square_free Is_square_free_1;
typedef typename Ptraits::Make_square_free Make_square_free_1;
typedef typename Ptraits::Square_free_factorize
Square_free_factorize_1;
typedef CGAL::RS_AK1::Is_coprime_1<Polynomial_1,Ptraits>
Is_coprime_1;
typedef CGAL::RS_AK1::Make_coprime_1<Polynomial_1,Ptraits>
Make_coprime_1;
typedef CGAL::RS_AK1::Solve_1<Polynomial_1,
Bound,
Algebraic_real_1,
Isolator,
Signat,
Ptraits>
Solve_1;
typedef CGAL::RS_AK1::Number_of_solutions_1<Polynomial_1,Isolator>
Number_of_solutions_1;
typedef CGAL::RS_AK1::Sign_at_1<Polynomial_1,
Bound,
Algebraic_real_1,
Refiner,
Signat,
Ptraits>
Sign_at_1;
typedef CGAL::RS_AK1::Is_zero_at_1<Polynomial_1,
Bound,
Algebraic_real_1,
Refiner,
Signat,
Ptraits>
Is_zero_at_1;
typedef CGAL::RS_AK1::Compare_1<Algebraic_real_1,
Bound,
Comparator>
Compare_1;
typedef CGAL::RS_AK1::Bound_between_1<Algebraic_real_1,
Bound,
Comparator>
Bound_between_1;
typedef CGAL::RS_AK1::Approximate_absolute_1<Polynomial_1,
Bound,
Algebraic_real_1,
Refiner>
Approximate_absolute_1;
typedef CGAL::RS_AK1::Approximate_relative_1<Polynomial_1,
Bound,
Algebraic_real_1,
Refiner>
Approximate_relative_1;
#define CREATE_FUNCTION_OBJECT(T,N) \
T N##_object()const{return T();}
CREATE_FUNCTION_OBJECT(Construct_algebraic_real_1,
construct_algebraic_real_1)
CREATE_FUNCTION_OBJECT(Compute_polynomial_1,
compute_polynomial_1)
CREATE_FUNCTION_OBJECT(Isolate_1,
isolate_1)
CREATE_FUNCTION_OBJECT(Is_square_free_1,
is_square_free_1)
CREATE_FUNCTION_OBJECT(Make_square_free_1,
make_square_free_1)
CREATE_FUNCTION_OBJECT(Square_free_factorize_1,
square_free_factorize_1)
CREATE_FUNCTION_OBJECT(Is_coprime_1,
is_coprime_1)
CREATE_FUNCTION_OBJECT(Make_coprime_1,
make_coprime_1)
CREATE_FUNCTION_OBJECT(Solve_1,
solve_1)
CREATE_FUNCTION_OBJECT(Number_of_solutions_1,
number_of_solutions_1)
CREATE_FUNCTION_OBJECT(Sign_at_1,
sign_at_1)
CREATE_FUNCTION_OBJECT(Is_zero_at_1,
is_zero_at_1)
CREATE_FUNCTION_OBJECT(Compare_1,
compare_1)
CREATE_FUNCTION_OBJECT(Bound_between_1,
bound_between_1)
CREATE_FUNCTION_OBJECT(Approximate_absolute_1,
approximate_absolute_1)
CREATE_FUNCTION_OBJECT(Approximate_relative_1,
approximate_relative_1)
#undef CREATE_FUNCTION_OBJECT
}; // class Algebraic_kernel_1
} // namespace RS_AK1
} // namespace CGAL
#endif // CGAL_RS_AK_1_H

View File

@ -0,0 +1,202 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_AK_Z_1_H
#define CGAL_RS_AK_Z_1_H
#include <cstddef> // included only to define size_t
#include <CGAL/Polynomial_traits_d.h>
#include "algebraic_z_1.h"
#include "comparator_1.h"
#include "signat_1.h"
#include "functors_z_1.h"
// This file defines the "Z-algebraic kernel". This kind of kernel performs
// all the internal operations using an integer polynomial type (the name
// "Z" comes from there). For this, a converter functor (passed as a
// template parameter) is used, which converts the input polynomial to the
// integer representation.
namespace CGAL{
namespace RS_AK1{
template <class ExtPolynomial_,
class IntPolynomial_,
class PolConverter_,
class Bound_,
class Isolator_,
class Refiner_,
class Ptraits_=CGAL::Polynomial_traits_d<ExtPolynomial_>,
class ZPtraits_=CGAL::Polynomial_traits_d<IntPolynomial_> >
class Algebraic_kernel_z_1{
public:
typedef ExtPolynomial_ Polynomial_1;
typedef IntPolynomial_ ZPolynomial_1;
typedef PolConverter_ PolConverter;
typedef typename Polynomial_1::NT Coefficient;
typedef Bound_ Bound;
private:
typedef Isolator_ Isolator;
typedef Refiner_ Refiner;
typedef Ptraits_ Ptraits;
typedef ZPtraits_ ZPtraits;
typedef CGAL::RS_AK1::Signat_1<ZPolynomial_1,Bound>
Signat;
typedef CGAL::RS_AK1::Simple_comparator_1<ZPolynomial_1,
Bound,
Refiner,
Signat,
ZPtraits>
Comparator;
public:
typedef CGAL::RS_AK1::Algebraic_z_1<Polynomial_1,
ZPolynomial_1,
Bound,
Refiner,
Comparator,
Ptraits,
ZPtraits>
Algebraic_real_1;
typedef size_t size_type;
typedef unsigned Multiplicity_type;
// default constructor and destructor
public:
Algebraic_kernel_z_1(){};
~Algebraic_kernel_z_1(){};
// functors from the CGAL concept
public:
typedef CGAL::RS_AK1::Construct_algebraic_real_z_1<Polynomial_1,
ZPolynomial_1,
PolConverter,
Algebraic_real_1,
Bound,
Coefficient,
Isolator>
Construct_algebraic_real_1;
typedef CGAL::RS_AK1::Compute_polynomial_z_1<Polynomial_1,
Algebraic_real_1>
Compute_polynomial_1;
typedef CGAL::RS_AK1::Isolate_z_1<Polynomial_1,
ZPolynomial_1,
PolConverter,
Bound,
Algebraic_real_1,
Isolator,
Comparator,
Signat,
Ptraits,
ZPtraits>
Isolate_1;
typedef typename Ptraits::Is_square_free Is_square_free_1;
typedef typename Ptraits::Make_square_free Make_square_free_1;
typedef typename Ptraits::Square_free_factorize Square_free_factorize_1;
typedef CGAL::RS_AK1::Is_coprime_z_1<Polynomial_1,Ptraits>
Is_coprime_1;
typedef CGAL::RS_AK1::Make_coprime_z_1<Polynomial_1,Ptraits>
Make_coprime_1;
typedef CGAL::RS_AK1::Solve_z_1<Polynomial_1,
ZPolynomial_1,
PolConverter,
Bound,
Algebraic_real_1,
Isolator,
Signat,
Ptraits,
ZPtraits>
Solve_1;
typedef CGAL::RS_AK1::Number_of_solutions_z_1<Polynomial_1,
ZPolynomial_1,
PolConverter,
Isolator>
Number_of_solutions_1;
typedef CGAL::RS_AK1::Sign_at_z_1<Polynomial_1,
ZPolynomial_1,
PolConverter,
Bound,
Algebraic_real_1,
Refiner,
Signat,
Ptraits,
ZPtraits>
Sign_at_1;
typedef CGAL::RS_AK1::Is_zero_at_z_1<Polynomial_1,
ZPolynomial_1,
PolConverter,
Bound,
Algebraic_real_1,
Refiner,
Signat,
Ptraits,
ZPtraits>
Is_zero_at_1;
typedef CGAL::RS_AK1::Compare_z_1<Algebraic_real_1,
Bound,
Comparator>
Compare_1;
typedef CGAL::RS_AK1::Bound_between_z_1<Algebraic_real_1,
Bound,
Comparator>
Bound_between_1;
typedef CGAL::RS_AK1::Approximate_absolute_z_1<Polynomial_1,
Bound,
Algebraic_real_1,
Refiner>
Approximate_absolute_1;
typedef CGAL::RS_AK1::Approximate_relative_z_1<Polynomial_1,
Bound,
Algebraic_real_1,
Refiner>
Approximate_relative_1;
#define CREATE_FUNCTION_OBJECT(T,N) \
T N##_object()const{return T();}
CREATE_FUNCTION_OBJECT(Construct_algebraic_real_1,
construct_algebraic_real_1)
CREATE_FUNCTION_OBJECT(Compute_polynomial_1,
compute_polynomial_1)
CREATE_FUNCTION_OBJECT(Isolate_1,
isolate_1)
CREATE_FUNCTION_OBJECT(Is_square_free_1,
is_square_free_1)
CREATE_FUNCTION_OBJECT(Make_square_free_1,
make_square_free_1)
CREATE_FUNCTION_OBJECT(Square_free_factorize_1,
square_free_factorize_1)
CREATE_FUNCTION_OBJECT(Is_coprime_1,
is_coprime_1)
CREATE_FUNCTION_OBJECT(Make_coprime_1,
make_coprime_1)
CREATE_FUNCTION_OBJECT(Solve_1,
solve_1)
CREATE_FUNCTION_OBJECT(Number_of_solutions_1,
number_of_solutions_1)
CREATE_FUNCTION_OBJECT(Sign_at_1,
sign_at_1)
CREATE_FUNCTION_OBJECT(Is_zero_at_1,
is_zero_at_1)
CREATE_FUNCTION_OBJECT(Compare_1,
compare_1)
CREATE_FUNCTION_OBJECT(Bound_between_1,
bound_between_1)
CREATE_FUNCTION_OBJECT(Approximate_absolute_1,
approximate_absolute_1)
CREATE_FUNCTION_OBJECT(Approximate_relative_1,
approximate_relative_1)
#undef CREATE_FUNCTION_OBJECT
}; // class Algebraic_kernel_z_1
} // namespace RS_AK1
} // namespace CGAL
#endif // CGAL_RS_AK_Z_1_H

View File

@ -0,0 +1,353 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_ALGEBRAIC_1_H
#define CGAL_RS_ALGEBRAIC_1_H
#include <boost/operators.hpp>
#include <CGAL/Real_embeddable_traits.h>
#include <CGAL/Gmpq.h>
#include <iostream>
namespace CGAL{
namespace RS_AK1{
// This class represents the simplest algebraic number one can think about.
// One algebraic number is represented by the polynomial of which it is
// root and the two endpoints of an interval that contains it, and no other
// root. Polynomial type and bound type are the first two template
// parameters.
//
// The third template parameter is a refiner, a function object that
// receives the polynomial and the bounds defining an algebraic number and
// an integer p, and modifies the two bounds until the difference between
// the two bounds is less than x*2^(-p), where x is the value of the
// represented algebraic number. The signature of a refiner must be:
// void
// Refiner_()(const Polynomial_&,Bound_&,Bound_&,int p);
//
// The fourth template argument is a comparator, a function object that
// receives the polynomials and bounds defining two algebraic numbers and
// just compares them, returning a CGAL::Comparison_result. The signature
// of a comparator must be:
// CGAL::Comparison_result
// Comparator_()(const Polynomial_&,Bound_&,Bound_&,
// const Polynomial_&,Bound_&,Bound_&);
// The comparator can modify the bounds, with the condition that the
// algebraic numbers remain consistent (one and only one root on each
// interval).
template <class Polynomial_,
class Bound_,
class Refiner_,
class Comparator_,
class Ptraits_>
class Algebraic_1:
boost::totally_ordered<Algebraic_1<Polynomial_,
Bound_,
Refiner_,
Comparator_,
Ptraits_>,
double>{
protected:
typedef Polynomial_ Polynomial;
typedef Bound_ Bound;
typedef Refiner_ Refiner;
typedef Comparator_ Comparator;
typedef Ptraits_ Ptraits;
typedef typename Ptraits::Coefficient_type Coefficient;
typedef typename Ptraits::Scale Scale;
typedef Algebraic_1<Polynomial,
Bound,
Refiner,
Comparator,
Ptraits>
Algebraic;
Polynomial pol;
mutable Bound left,right;
public:
Algebraic_1(){};
Algebraic_1(const Polynomial &p,
const Bound &l,
const Bound &r):pol(p),left(l),right(r){
CGAL_assertion(l<=r);
}
Algebraic_1(const Algebraic &a):
pol(a.pol),left(a.left),right(a.right){}
// XXX: This assumes that Gmpq is constructible from bound type and
// that polynomial coefficient type is constructible from mpz_t.
Algebraic_1(const Bound &b):left(b),right(b){
typedef typename Ptraits::Shift shift;
Gmpq q(b);
pol=Coefficient(mpq_denref(q.mpq()))*
shift()(Polynomial(1),1,0)-
Coefficient(mpq_numref(q.mpq()));
CGAL_assertion(left==right&&left==b);
}
// XXX: This implementation assumes that the bound type is Gmpfr
// and that T can be exactly converted to Gmpq. This constructor
// can handle types such as int, unsigned and long.
template <class T>
Algebraic_1(const T &t){
typedef typename Ptraits::Shift shift;
CGAL::Gmpq q(t);
pol=Coefficient(mpq_denref(q.mpq()))*
shift()(Polynomial(1),1,0)-
Coefficient(mpq_numref(q.mpq()));
left=Bound(t,std::round_toward_neg_infinity);
right=Bound(t,std::round_toward_infinity);
CGAL_assertion(left<=t&&right>=t);
}
// XXX: This constructor assumes the bound is a Gmpfr.
Algebraic_1(const CGAL::Gmpq &q){
typedef typename Ptraits::Shift shift;
pol=Coefficient(mpq_denref(q.mpq()))*
shift()(Polynomial(1),1,0)-
Coefficient(mpq_numref(q.mpq()));
left=Bound();
right=Bound();
mpfr_t b;
mpfr_init(b);
mpfr_set_q(b,q.mpq(),GMP_RNDD);
mpfr_swap(b,left.fr());
mpfr_set_q(b,q.mpq(),GMP_RNDU);
mpfr_swap(b,right.fr());
mpfr_clear(b);
CGAL_assertion(left<=q&&right>=q);
}
~Algebraic_1(){}
Algebraic_1& operator=(const Algebraic &a){
pol=a.get_pol();
left=a.get_left();
right=a.get_right();
return *this;
}
Polynomial get_pol()const{return pol;}
Bound& get_left()const{return left;}
Bound& get_right()const{return right;}
Algebraic operator-()const{
return Algebraic(Scale()(get_pol(),Coefficient(-1)),
-right,
-left);
}
#define CGAL_RS_COMPARE_ALGEBRAIC(_a) \
(Comparator()(get_pol(),get_left(),get_right(), \
(_a).get_pol(),(_a).get_left(),(_a).get_right()))
Comparison_result compare(Algebraic a)const{
return CGAL_RS_COMPARE_ALGEBRAIC(a);
};
#define CGAL_RS_COMPARE_ALGEBRAIC_TYPE(_t) \
bool operator<(_t t)const \
{Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::SMALLER;} \
bool operator>(_t t)const \
{Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::LARGER;} \
bool operator==(_t t)const \
{Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::EQUAL;}
bool operator==(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::EQUAL;}
bool operator!=(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC(a)!=CGAL::EQUAL;}
bool operator<(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::SMALLER;}
bool operator<=(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC(a)!=CGAL::LARGER;}
bool operator>(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::LARGER;}
bool operator>=(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC(a)!=CGAL::SMALLER;}
CGAL_RS_COMPARE_ALGEBRAIC_TYPE(double)
#undef CGAL_RS_COMPARE_ALGEBRAIC_TYPE
#undef CGAL_RS_COMPARE_ALGEBRAIC
#ifdef IEEE_DBL_MANT_DIG
#define CGAL_RS_DBL_PREC IEEE_DBL_MANT_DIG
#else
#define CGAL_RS_DBL_PREC 53
#endif
// This function is const because left and right are mutable.
double to_double()const{
typedef Real_embeddable_traits<Bound> RT;
typedef typename RT::To_double TD;
Refiner()(pol,left,right,CGAL_RS_DBL_PREC);
CGAL_assertion(TD()(left)==TD()(right));
return TD()(left);
}
std::pair<double,double> to_interval()const{
typedef Real_embeddable_traits<Bound> RT;
typedef typename RT::To_interval TI;
return std::make_pair(TI()(get_left()).first,
TI()(get_right()).second);
}
#undef CGAL_RS_DBL_PREC
void set_left(const Bound &l)const{
left=l;
}
void set_right(const Bound &r)const{
right=r;
}
void set_pol(const Polynomial &p){
pol=p;
}
}; // class Algebraic_1
} // namespace RS_AK1
// We define Algebraic_1 as real-embeddable
template <class Polynomial_,
class Bound_,
class Refiner_,
class Comparator_,
class Ptraits_>
class Real_embeddable_traits<RS_AK1::Algebraic_1<Polynomial_,
Bound_,
Refiner_,
Comparator_,
Ptraits_> >:
public INTERN_RET::Real_embeddable_traits_base<
RS_AK1::Algebraic_1<Polynomial_,
Bound_,
Refiner_,
Comparator_,
Ptraits_>,
CGAL::Tag_true>{
typedef Polynomial_ P;
typedef Bound_ B;
typedef Refiner_ R;
typedef Comparator_ C;
typedef Ptraits_ T;
public:
typedef RS_AK1::Algebraic_1<P,B,R,C,T> Type;
typedef CGAL::Tag_true Is_real_embeddable;
typedef bool Boolean;
typedef CGAL::Sign Sign;
typedef CGAL::Comparison_result Comparison_result;
typedef INTERN_RET::Real_embeddable_traits_base<Type,CGAL::Tag_true>
Base;
typedef typename Base::Compare Compare;
class Sgn:public CGAL::cpp98::unary_function<Type,CGAL::Sign>{
public:
CGAL::Sign operator()(const Type &a)const{
return Compare()(a,Type(0));
}
};
class To_double:public CGAL::cpp98::unary_function<Type,double>{
public:
double operator()(const Type &a)const{return a.to_double();}
};
class To_interval:
public CGAL::cpp98::unary_function<Type,std::pair<double,double> >{
public:
std::pair<double,double> operator()(const Type &a)const{
return a.to_interval();
}
};
class Is_zero:public CGAL::cpp98::unary_function<Type,Boolean>{
public:
bool operator()(const Type &a)const{
return Sgn()(a)==CGAL::ZERO;
}
};
class Is_finite:public CGAL::cpp98::unary_function<Type,Boolean>{
public:
bool operator()(const Type&)const{return true;}
};
class Abs:public CGAL::cpp98::unary_function<Type,Type>{
public:
Type operator()(const Type &a)const{
return Sgn()(a)==CGAL::NEGATIVE?-a:a;
}
};
};
template <class P,class B,class R,class C,class T>
inline
RS_AK1::Algebraic_1<P,B,R,C,T> min
BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_1<P,B,R,C,T> a,
RS_AK1::Algebraic_1<P,B,R,C,T> b){
return(a<b?a:b);
}
template <class P,class B,class R,class C,class T>
inline
RS_AK1::Algebraic_1<P,B,R,C,T> max
BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_1<P,B,R,C,T> a,
RS_AK1::Algebraic_1<P,B,R,C,T> b){
return(a>b?a:b);
}
template <class P,class B,class R,class C,class T>
inline
std::ostream& operator<<(std::ostream &o,
const RS_AK1::Algebraic_1<P,B,R,C,T> &a){
return(o<<'['<<a.get_pol()<<','<<
a.get_left()<<','<<
a.get_right()<<']');
}
// XXX: This function works, but it will be nice to rewrite it cleanly.
template <class P,class B,class R,class C,class T>
inline
std::istream& operator>>(std::istream &i,
RS_AK1::Algebraic_1<P,B,R,C,T> &a){
std::istream::int_type c;
P pol;
B lb,rb;
c=i.get();
if(c!='['){
CGAL_error_msg("error reading istream, \'[\' expected");
return i;
}
i>>pol;
c=i.get();
if(c!=','){
CGAL_error_msg("error reading istream, \',\' expected");
return i;
}
i>>lb;
c=i.get();
if(c!=','){
CGAL_error_msg("error reading istream, \',\' expected");
return i;
}
i>>rb;
c=i.get();
if(c!=']'){
CGAL_error_msg("error reading istream, \']\' expected");
return i;
}
a=RS_AK1::Algebraic_1<P,B,R,C,T>(pol,lb,rb);
return i;
}
} // namespace CGAL
#endif // CGAL_RS_ALGEBRAIC_1_H

View File

@ -0,0 +1,351 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_ALGEBRAIC_Z_1_H
#define CGAL_RS_ALGEBRAIC_Z_1_H
#include "algebraic_1.h"
namespace CGAL{
namespace RS_AK1{
// This class represents an algebraic number storing two polynomials of
// which it is root: one of them is the given polynomial; the other one is
// an integer polynomial, which is used to perform the operations such as
// refinements. Since RS works only with integer polynomials, this
// architecture permits to convert the input polynomials only once.
template <class Polynomial_,
class ZPolynomial_,
class Bound_,
class ZRefiner_,
class ZComparator_,
class Ptraits_,
class ZPtraits_>
class Algebraic_z_1:
public Algebraic_1<Polynomial_,Bound_,ZRefiner_,ZComparator_,ZPtraits_>,
boost::totally_ordered<Algebraic_z_1<Polynomial_,
ZPolynomial_,
Bound_,
ZRefiner_,
ZComparator_,
Ptraits_,
ZPtraits_>,
double>{
protected:
typedef Polynomial_ Polynomial;
typedef ZPolynomial_ ZPolynomial;
typedef Bound_ Bound;
typedef ZRefiner_ ZRefiner;
typedef ZComparator_ ZComparator;
typedef ZPtraits_ ZPtraits;
typedef typename ZPtraits::Coefficient_type ZCoefficient;
typedef typename ZPtraits::Scale ZScale;
typedef Ptraits_ Ptraits;
typedef typename Ptraits::Coefficient_type Coefficient;
typedef typename Ptraits::Scale Scale;
typedef Algebraic_1<Polynomial,Bound,ZRefiner,ZComparator,ZPtraits>
Algebraic_base;
typedef Algebraic_z_1<Polynomial,
ZPolynomial,
Bound,
ZRefiner,
ZComparator,
Ptraits,
ZPtraits>
Algebraic;
ZPolynomial zpol;
public:
Algebraic_z_1(){};
Algebraic_z_1(const Polynomial &p,
const ZPolynomial &zp,
const Bound &l,
const Bound &r):Algebraic_base(p,l,r),zpol(zp){
CGAL_assertion(l<=r);
}
Algebraic_z_1(const Algebraic &a):
Algebraic_base(a.pol,a.left,a.right),zpol(a.zpol){}
Algebraic_z_1(const Bound &b){
typedef typename Ptraits::Shift Shift;
typedef typename ZPtraits::Shift ZShift;
this->left=b;
this->right=b;
Gmpq q(b);
this->pol=Coefficient(mpq_denref(q.mpq()))*
Shift()(Polynomial(1),1,0)-
Coefficient(mpq_numref(q.mpq()));
zpol=ZCoefficient(mpq_denref(q.mpq()))*
ZShift()(ZPolynomial(1),1,0)-
ZCoefficient(mpq_numref(q.mpq()));
CGAL_assertion(this->left==this->right);
CGAL_assertion(this->left==b);
}
template <class T>
Algebraic_z_1(const T &t){
typedef typename Ptraits::Shift Shift;
typedef typename ZPtraits::Shift ZShift;
CGAL::Gmpq q(t);
this->pol=Coefficient(mpq_denref(q.mpq()))*
Shift()(Polynomial(1),1,0)-
Coefficient(mpq_numref(q.mpq()));
zpol=ZCoefficient(mpq_denref(q.mpq()))*
ZShift()(ZPolynomial(1),1,0)-
ZCoefficient(mpq_numref(q.mpq()));
this->left=Bound(t,std::round_toward_neg_infinity);
this->right=Bound(t,std::round_toward_infinity);
CGAL_assertion(this->left<=t&&this->right>=t);
}
Algebraic_z_1(const CGAL::Gmpq &q){
typedef typename Ptraits::Shift Shift;
typedef typename ZPtraits::Shift ZShift;
this->pol=Coefficient(mpq_denref(q.mpq()))*
Shift()(Polynomial(1),1,0)-
Coefficient(mpq_numref(q.mpq()));
zpol=ZCoefficient(mpq_denref(q.mpq()))*
ZShift()(ZPolynomial(1),1,0)-
ZCoefficient(mpq_numref(q.mpq()));
this->left=Bound();
this->right=Bound();
mpfr_t b;
mpfr_init(b);
mpfr_set_q(b,q.mpq(),GMP_RNDD);
mpfr_swap(b,this->left.fr());
mpfr_set_q(b,q.mpq(),GMP_RNDU);
mpfr_swap(b,this->right.fr());
mpfr_clear(b);
CGAL_assertion(this->left<=q&&this->right>=q);
}
~Algebraic_z_1(){}
ZPolynomial get_zpol()const{return zpol;}
void set_zpol(const ZPolynomial &p){zpol=p;}
Algebraic operator-()const{
return Algebraic(Scale()(this->get_pol(),Coefficient(-1)),
ZScale()(get_zpol(),ZCoefficient(-1)),
-this->right,
-this->left);
}
#define CGAL_RS_COMPARE_ALGEBRAIC_Z(_a) \
(ZComparator()(get_zpol(),this->get_left(),this->get_right(), \
(_a).get_zpol(),(_a).get_left(),(_a).get_right()))
Comparison_result compare(Algebraic a)const{
return CGAL_RS_COMPARE_ALGEBRAIC_Z(a);
};
#define CGAL_RS_COMPARE_ALGEBRAIC_Z_TYPE(_t) \
bool operator<(_t t)const \
{Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::SMALLER;} \
bool operator>(_t t)const \
{Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::LARGER;} \
bool operator==(_t t)const \
{Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::EQUAL;}
bool operator==(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::EQUAL;}
bool operator!=(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)!=CGAL::EQUAL;}
bool operator<(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::SMALLER;}
bool operator<=(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)!=CGAL::LARGER;}
bool operator>(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::LARGER;}
bool operator>=(Algebraic a)const
{return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)!=CGAL::SMALLER;}
CGAL_RS_COMPARE_ALGEBRAIC_Z_TYPE(double)
#undef CGAL_RS_COMPARE_ALGEBRAIC_Z_TYPE
#undef CGAL_RS_COMPARE_ALGEBRAIC_Z
#ifdef IEEE_DBL_MANT_DIG
#define CGAL_RS_DBL_PREC IEEE_DBL_MANT_DIG
#else
#define CGAL_RS_DBL_PREC 53
#endif
double to_double()const{
typedef Real_embeddable_traits<Bound> RT;
typedef typename RT::To_double TD;
ZRefiner()(get_zpol(),
this->get_left(),
this->get_right(),
CGAL_RS_DBL_PREC);
CGAL_assertion(TD()(this->get_left())==TD()(this->get_right()));
return TD()(this->get_left());
}
std::pair<double,double> to_interval()const{
typedef Real_embeddable_traits<Bound> RT;
typedef typename RT::To_interval TI;
return std::make_pair(TI()(this->get_left()).first,
TI()(this->get_right()).second);
}
#undef CGAL_RS_DBL_PREC
}; // class Algebraic_z_1
} // namespace RS_AK1
// We define Algebraic_z_1 as real-embeddable
template <class Polynomial_,
class ZPolynomial_,
class Bound_,
class ZRefiner_,
class ZComparator_,
class Ptraits_,
class ZPtraits_>
class Real_embeddable_traits<RS_AK1::Algebraic_z_1<Polynomial_,
ZPolynomial_,
Bound_,
ZRefiner_,
ZComparator_,
Ptraits_,
ZPtraits_> >:
public INTERN_RET::Real_embeddable_traits_base<
RS_AK1::Algebraic_z_1<Polynomial_,
ZPolynomial_,
Bound_,
ZRefiner_,
ZComparator_,
Ptraits_,
ZPtraits_>,
CGAL::Tag_true>{
typedef Polynomial_ P;
typedef ZPolynomial_ ZP;
typedef Bound_ B;
typedef ZRefiner_ R;
typedef ZComparator_ C;
typedef Ptraits_ T;
typedef ZPtraits_ ZT;
public:
typedef RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> Type;
typedef CGAL::Tag_true Is_real_embeddable;
typedef bool Boolean;
typedef CGAL::Sign Sign;
typedef CGAL::Comparison_result Comparison_result;
typedef INTERN_RET::Real_embeddable_traits_base<Type,CGAL::Tag_true>
Base;
typedef typename Base::Compare Compare;
class Sgn:public CGAL::cpp98::unary_function<Type,CGAL::Sign>{
public:
CGAL::Sign operator()(const Type &a)const{
return Compare()(a,Type(0));
}
};
class To_double:public CGAL::cpp98::unary_function<Type,double>{
public:
double operator()(const Type &a)const{return a.to_double();}
};
class To_interval:
public CGAL::cpp98::unary_function<Type,std::pair<double,double> >{
public:
std::pair<double,double> operator()(const Type &a)const{
return a.to_interval();
}
};
class Is_zero:public CGAL::cpp98::unary_function<Type,Boolean>{
public:
bool operator()(const Type &a)const{
return Sgn()(a)==CGAL::ZERO;
}
};
class Is_finite:public CGAL::cpp98::unary_function<Type,Boolean>{
public:
bool operator()(const Type&)const{return true;}
};
class Abs:public CGAL::cpp98::unary_function<Type,Type>{
public:
Type operator()(const Type &a)const{
return Sgn()(a)==CGAL::NEGATIVE?-a:a;
}
};
};
template <class P,class ZP,class B,class R,class C,class T,class ZT>
inline
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> min
BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> a,
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> b){
return(a<b?a:b);
}
template <class P,class ZP,class B,class R,class C,class T,class ZT>
inline
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> max
BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> a,
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> b){
return(a>b?a:b);
}
template <class P,class ZP,class B,class R,class C,class T,class ZT>
inline
std::ostream& operator<<(std::ostream &o,
const RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> &a){
return(o<<'['<<a.get_pol()<<','<<
a.get_zpol()<<','<<
a.get_left()<<','<<
a.get_right()<<']');
}
template <class P,class ZP,class B,class R,class C,class T,class ZT>
inline
std::istream& operator>>(std::istream &i,
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> &a){
std::istream::int_type c;
P pol;
ZP zpol;
B lb,rb;
c=i.get();
if(c!='['){
CGAL_error_msg("error reading istream, \'[\' expected");
return i;
}
i>>pol;
c=i.get();
if(c!=','){
CGAL_error_msg("error reading istream, \',\' expected");
return i;
}
i>>zpol;
c=i.get();
if(c!=','){
CGAL_error_msg("error reading istream, \',\' expected");
return i;
}
i>>lb;
c=i.get();
if(c!=','){
CGAL_error_msg("error reading istream, \',\' expected");
return i;
}
i>>rb;
c=i.get();
if(c!=']'){
CGAL_error_msg("error reading istream, \']\' expected");
return i;
}
a=RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT>(pol,zpol,lb,rb);
return i;
}
} // namespace CGAL
#endif // CGAL_RS_ALGEBRAIC_Z_1_H

View File

@ -0,0 +1,165 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
// This file contains the simplest refiner, that bisects the interval a given
// number of times.
#ifndef CGAL_RS_BISECTION_REFINER_1_H
#define CGAL_RS_BISECTION_REFINER_1_H
#include <CGAL/Polynomial_traits_d.h>
#include "signat_1.h"
#include "Gmpfr_make_unique.h"
namespace CGAL{
template <class Polynomial_,class Bound_>
struct Bisection_refiner_1{
typedef CGAL::RS_AK1::Signat_1<Polynomial_,Bound_> Signat;
void operator()(const Polynomial_&,Bound_&,Bound_&,int);
}; // class Bisection_refiner_1
// TODO: Write in a generic way, if possible (see next function).
template <class Polynomial_,class Bound_>
void
Bisection_refiner_1<Polynomial_,Bound_>::
operator()(const Polynomial_&,Bound_&,Bound_&,int){
CGAL_error_msg("bisection refiner not implemented for these types");
return;
}
// This works with any type of polynomial, but only for Gmpfr bounds.
// TODO: Beyond writing generically, optimize this function. This would
// remove the need for the next function, which is essentially the same.
template<>
void
Bisection_refiner_1<Polynomial<Gmpz>,Gmpfr>::
operator()(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec){
typedef Polynomial<Gmpz> Polynomial;
typedef Polynomial_traits_d<Polynomial> Ptraits;
typedef Ptraits::Make_square_free Sfpart;
typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
Signat;
CGAL_precondition(left<=right);
//std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
CGAL_RS_GMPFR_MAKE_UNIQUE(right,temp_right);
Polynomial sfpp=Sfpart()(pol);
Signat signof(sfpp);
CGAL::Sign sl,sc;
mp_prec_t pl,pc;
mpfr_t center;
sl=signof(left);
CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
if(sl==ZERO)
return;
pl=left.get_precision();
pc=right.get_precision();
pc=(pl>pc?pl:pc)+(mp_prec_t)prec;
mpfr_init2(center,pc);
CGAL_assertion_code(int round=)
mpfr_prec_round(left.fr(),pc,GMP_RNDN);
CGAL_assertion(!round);
CGAL_assertion_code(round=)
mpfr_prec_round(right.fr(),pc,GMP_RNDN);
CGAL_assertion(!round);
for(int i=0;i<prec;++i){
CGAL_assertion_code(round=)
mpfr_add(center,left.fr(),right.fr(),GMP_RNDN);
CGAL_assertion(!round);
CGAL_assertion_code(round=)
mpfr_div_2ui(center,center,1,GMP_RNDN);
CGAL_assertion(!round);
sc=signof(Gmpfr(center));
if(sc==ZERO){ // we have a root
CGAL_assertion_code(round=)
mpfr_set(left.fr(),center,GMP_RNDN);
CGAL_assertion(!round);
mpfr_swap(right.fr(),center);
break;
}
if(sc==sl)
mpfr_swap(left.fr(),center);
else
mpfr_swap(right.fr(),center);
}
mpfr_clear(center);
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
template<>
void
Bisection_refiner_1<Polynomial<Gmpq>,Gmpfr>::
operator()(const Polynomial<Gmpq> &pol,Gmpfr &left,Gmpfr &right,int prec){
typedef Polynomial<Gmpq> Polynomial;
typedef Polynomial_traits_d<Polynomial> Ptraits;
typedef Ptraits::Make_square_free Sfpart;
typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
Signat;
CGAL_precondition(left<=right);
//std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
CGAL_RS_GMPFR_MAKE_UNIQUE(right,temp_right);
Polynomial sfpp=Sfpart()(pol);
Signat signof(sfpp);
CGAL::Sign sl,sc;
mp_prec_t pl,pc;
mpfr_t center;
sl=signof(left);
CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
if(sl==ZERO)
return;
pl=left.get_precision();
pc=right.get_precision();
pc=(pl>pc?pl:pc)+(mp_prec_t)prec;
mpfr_init2(center,pc);
CGAL_assertion_code(int round=)
mpfr_prec_round(left.fr(),pc,GMP_RNDN);
CGAL_assertion(!round);
CGAL_assertion_code(round=)
mpfr_prec_round(right.fr(),pc,GMP_RNDN);
CGAL_assertion(!round);
for(int i=0;i<prec;++i){
CGAL_assertion_code(round=)
mpfr_add(center,left.fr(),right.fr(),GMP_RNDN);
CGAL_assertion(!round);
CGAL_assertion_code(round=)
mpfr_div_2ui(center,center,1,GMP_RNDN);
CGAL_assertion(!round);
sc=signof(Gmpfr(center));
if(sc==ZERO){ // we have a root
CGAL_assertion_code(round=)
mpfr_set(left.fr(),center,GMP_RNDN);
CGAL_assertion(!round);
mpfr_swap(right.fr(),center);
break;
}
if(sc==sl)
mpfr_swap(left.fr(),center);
else
mpfr_swap(right.fr(),center);
}
mpfr_clear(center);
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
} // namespace CGAL
#endif // CGAL_RS_BISECTION_REFINER_1_H

View File

@ -0,0 +1,81 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_COMPARATOR_1_H
#define CGAL_RS_COMPARATOR_1_H
namespace CGAL{
namespace RS_AK1{
template <class Polynomial_,
class Bound_,
class Refiner_,
class Signat_,
class Ptraits_>
struct Simple_comparator_1{
typedef Polynomial_ Polynomial;
typedef Bound_ Bound;
typedef Refiner_ Refiner;
typedef Signat_ Signat;
typedef Ptraits_ Ptraits;
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Degree Degree;
CGAL::Comparison_result
operator()(const Polynomial &p1,Bound &l1,Bound &r1,
const Polynomial &p2,Bound &l2,Bound &r2)const{
CGAL_precondition(l1<=r1&&l2<=r2);
if(l1<=l2){
if(r1<l2)
return SMALLER;
}else{
if(r2<l1)
return LARGER;
}
Polynomial G=Gcd()(p1,p2);
if(Degree()(G)==0)
return compare_unequal(p1,l1,r1,p2,l2,r2);
Signat sg(G);
CGAL::Sign sleft=sg(l1>l2?l1:l2);
if(sleft==ZERO)
return EQUAL;
CGAL::Sign sright=sg(r1<r2?r1:r2);
if(sleft!=sright)
return EQUAL;
else
return compare_unequal(p1,l1,r1,p2,l2,r2);
}
// This function compares two algebraic numbers, assuming that they
// are not equal.
CGAL::Comparison_result
compare_unequal(const Polynomial &p1,Bound &l1,Bound &r1,
const Polynomial &p2,Bound &l2,Bound &r2)const{
CGAL_precondition(l1<=r1&&l2<=r2);
int prec=CGAL::max(
CGAL::max(l1.get_precision(),
r1.get_precision()),
CGAL::max(l2.get_precision(),
r2.get_precision()));
do{
prec*=2;
Refiner()(p1,l1,r1,prec);
Refiner()(p2,l2,r2,prec);
CGAL_assertion(l1<=r1&&l2<=r2);
}while(l1<=l2?r1>=l2:r2>=l1);
return (r1<l2?SMALLER:LARGER);
}
}; // struct Simple_comparator_1
} // namespace RS_AK1
} // namespace CGAL
#endif // CGAL_RS_COMPARATOR_1_H

View File

@ -0,0 +1,430 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_DYADIC_H
#define CGAL_RS_DYADIC_H
#include <stdio.h>
#include <math.h>
#include <gmp.h>
#include <mpfr.h>
#include <CGAL/assertions.h>
// for c++, compile with -lgmpxx
#ifdef __cplusplus
#include <iostream>
#endif
#define CGALRS_dyadic_struct __mpfr_struct
#define CGALRS_dyadic_t mpfr_t
#define CGALRS_dyadic_ptr mpfr_ptr
#define CGALRS_dyadic_srcptr mpfr_srcptr
// some auxiliary defines
#define CGALRS_dyadic_set_prec(D,P) \
( mpfr_set_prec( (D), (P)>MPFR_PREC_MIN?(P):MPFR_PREC_MIN) )
#define CGALRS_dyadic_prec_round(D,P) \
( mpfr_prec_round( (D), (P)>MPFR_PREC_MIN?(P):MPFR_PREC_MIN, GMP_RNDN) )
#define CGALRS_dyadic_set_exp(D,E) \
( CGAL_assertion( (E) <= mpfr_get_emax() && \
(E) >= mpfr_get_emin() ) ,\
mpfr_set_exp(D,E) )
// init functions
#define CGALRS_dyadic_init(D) mpfr_init2(D,MPFR_PREC_MIN)
#define CGALRS_dyadic_init2(D,P) mpfr_init2(D,P)
#define CGALRS_dyadic_clear(D) mpfr_clear(D)
inline void CGALRS_dyadic_set(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op){
if(rop!=op){
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
mpfr_set(rop,op,GMP_RNDN);
}
CGAL_assertion(mpfr_equal_p(rop,op)!=0);
}
inline void CGALRS_dyadic_set_z(CGALRS_dyadic_ptr rop,mpz_srcptr z){
size_t prec;
prec=mpz_sizeinbase(z,2)-(mpz_tstbit(z,0)?0:mpz_scan1(z,0));
CGALRS_dyadic_set_prec(rop,prec);
mpfr_set_z(rop,z,GMP_RNDN);
CGAL_assertion(!mpfr_cmp_z(rop,z));
}
inline void CGALRS_dyadic_set_si(CGALRS_dyadic_ptr rop,long s){
CGALRS_dyadic_set_prec(rop,sizeof(long));
mpfr_set_si(rop,s,GMP_RNDN);
CGAL_assertion(!mpfr_cmp_si(rop,s));
}
inline void CGALRS_dyadic_set_ui(CGALRS_dyadic_ptr rop,unsigned long u){
CGALRS_dyadic_set_prec(rop,sizeof(unsigned long));
mpfr_set_ui(rop,u,GMP_RNDN);
CGAL_assertion(!mpfr_cmp_ui(rop,u));
}
inline void CGALRS_dyadic_set_fr(CGALRS_dyadic_ptr rop,mpfr_srcptr op){
if(rop!=op){
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
mpfr_set(rop,op,GMP_RNDN);
CGAL_assertion(mpfr_equal_p(rop,op)!=0);
}
}
#define CGALRS_dyadic_init_set(R,D) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set((R), (D)) )
#define CGALRS_dyadic_init_set_z(R,Z) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set_z((R), (Z)) )
#define CGALRS_dyadic_init_set_si(R,I) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set_si((R), (I)) )
#define CGALRS_dyadic_init_set_ui(R,I) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set_ui((R), (I)) )
#define CGALRS_dyadic_init_set_fr(R,F) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set_fr((R), (F)) )
#define CGALRS_dyadic_get_fr(M,D) mpfr_set(M,D,GMP_RNDN)
#define CGALRS_dyadic_get_d(D,RM) mpfr_get_d(D,RM)
inline void CGALRS_dyadic_get_exactfr(mpfr_ptr rop,CGALRS_dyadic_srcptr op){
if(rop!=op){
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
mpfr_set(rop,op,GMP_RNDN);
CGAL_assertion(mpfr_equal_p(rop,op)!=0);
}
}
#define CGALRS_dyadic_canonicalize(D) ()
// comparison functions
#define CGALRS_dyadic_sgn(D) mpfr_sgn(D)
#define CGALRS_dyadic_zero(D) mpfr_zero_p(D)
#define CGALRS_dyadic_cmp(D,E) mpfr_cmp(D,E)
// arithmetic functions
#define CGALRS_dyadic_add(R,D,E) CGALRS_dyadic_ll_add(R,D,E,0)
#define CGALRS_dyadic_sub(R,D,E) CGALRS_dyadic_ll_sub(R,D,E,0)
#define CGALRS_dyadic_mul(R,D,E) CGALRS_dyadic_ll_mul(R,D,E,0)
inline void CGALRS_dyadic_neg(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op){
if(rop!=op)
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
mpfr_neg(rop,op,GMP_RNDN);
CGAL_assertion(
rop==op||
(!mpfr_cmpabs(rop,op)&&
((CGALRS_dyadic_zero(op)&&CGALRS_dyadic_zero(rop))||
(CGALRS_dyadic_sgn(op)!=CGALRS_dyadic_sgn(rop)))));
}
// low-level addition:
// add op1 and op2 and reserve b bits for future lowlevel operations
inline void CGALRS_dyadic_ll_add(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
CGALRS_dyadic_srcptr op2,
mp_prec_t b){
mp_exp_t l,r,temp1,temp2;
mp_prec_t rop_prec;
if(mpfr_zero_p(op1)){
if(rop!=op2)
CGALRS_dyadic_set(rop,op2);
return;
}
if(mpfr_zero_p(op2)){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
l=mpfr_get_exp(op1)>mpfr_get_exp(op2)?
mpfr_get_exp(op1):
mpfr_get_exp(op2);
temp1=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1);
temp2=mpfr_get_exp(op2)-(mp_exp_t)mpfr_get_prec(op2);
r=temp1>temp2?temp2:temp1;
CGAL_assertion(l>r);
rop_prec=b+1+(mp_prec_t)(l-r);
CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
rop_prec>=mpfr_get_prec(op2));
if(rop==op1||rop==op2)
CGALRS_dyadic_prec_round(rop,rop_prec);
else
CGALRS_dyadic_set_prec(rop,rop_prec);
CGAL_assertion_code(int round=)
mpfr_add(rop,op1,op2,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_add_z(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
mpz_srcptr z){
mp_exp_t l,r;
mp_prec_t rop_prec;
if(mpfr_zero_p(op1)){
CGALRS_dyadic_set_z(rop,z);
return;
}
if(!mpz_sgn(z)){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
l=mpfr_get_exp(op1)>(mp_exp_t)mpz_sizeinbase(z,2)?
mpfr_get_exp(op1):
mpz_sizeinbase(z,2);
r=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1)<0?
mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1):
0;
CGAL_assertion(l>r);
rop_prec=1+(mp_prec_t)(l-r);
CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
rop_prec>=(mp_prec_t)mpz_sizeinbase(z,2));
if(rop==op1)
CGALRS_dyadic_prec_round(rop,rop_prec);
else
CGALRS_dyadic_set_prec(rop,rop_prec);
CGAL_assertion_code(int round=)
mpfr_add_z(rop,op1,z,GMP_RNDN);
CGAL_assertion(!round);
}
// low-level subtraction:
// subtract op2 to op1 and reserve b bits for future lowlevel operations
inline void CGALRS_dyadic_ll_sub(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
CGALRS_dyadic_srcptr op2,
mp_prec_t b){
mp_exp_t l,r,temp1,temp2;
mp_prec_t rop_prec;
if(mpfr_zero_p(op1)){
CGALRS_dyadic_neg(rop,op2);
return;
}
if(mpfr_zero_p(op2)){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
l=mpfr_get_exp(op1)>mpfr_get_exp(op2)?
mpfr_get_exp(op1):
mpfr_get_exp(op2);
temp1=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1);
temp2=mpfr_get_exp(op2)-(mp_exp_t)mpfr_get_prec(op2);
r=temp1>temp2?temp2:temp1;
CGAL_assertion(l>r);
rop_prec=b+1+(mp_prec_t)(l-r);
CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
rop_prec>=mpfr_get_prec(op2));
if(rop==op1||rop==op2)
CGALRS_dyadic_prec_round(rop,rop_prec);
else
CGALRS_dyadic_set_prec(rop,rop_prec);
CGAL_assertion_code(int round=)
mpfr_sub(rop,op1,op2,GMP_RNDN);
CGAL_assertion(!round);
}
// low-level multiplication:
// multiply op1 and op2 and reserve b bits for future lowlevel operations
inline void CGALRS_dyadic_ll_mul(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
CGALRS_dyadic_srcptr op2,
mp_prec_t b){
if(rop==op1||rop==op2)
CGALRS_dyadic_prec_round(rop,b+mpfr_get_prec(op1)+mpfr_get_prec(op2));
else
CGALRS_dyadic_set_prec(rop,b+mpfr_get_prec(op1)+mpfr_get_prec(op2));
CGAL_assertion_code(int round=)
mpfr_mul(rop,op1,op2,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_mul_z(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
mpz_srcptr z){
if(rop==op1)
CGALRS_dyadic_prec_round(
rop,
mpfr_get_prec(op1)+mpz_sizeinbase(z,2));
else
CGALRS_dyadic_set_prec(
rop,
mpfr_get_prec(op1)+mpz_sizeinbase(z,2));
CGAL_assertion_code(int round=)
mpfr_mul_z(rop,op1,z,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_mul_si(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
long s){
if(rop==op1)
CGALRS_dyadic_prec_round(rop,mpfr_get_prec(op1)+sizeof(long));
else
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op1)+sizeof(long));
CGAL_assertion_code(int round=)
mpfr_mul_si(rop,op1,s,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_mul_ui(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long u){
if(rop==op1)
CGALRS_dyadic_prec_round(
rop,
mpfr_get_prec(op1)+sizeof(unsigned long));
else
CGALRS_dyadic_set_prec(
rop,
mpfr_get_prec(op1)+sizeof(unsigned long));
CGAL_assertion_code(int round=)
mpfr_mul_ui(rop,op1,u,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_pow_ui(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long u){
if(!u){
CGAL_assertion_msg(!mpfr_zero_p(op1),"0^0");
CGALRS_dyadic_set_ui(rop,1);
return;
}
if(u==1){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
if(mpfr_zero_p(op1)){
CGAL_assertion_msg(u!=0,"0^0");
CGALRS_dyadic_set_ui(rop,0);
return;
}
if(!mpfr_cmp_ui(op1,1)){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
if(rop==op1)
CGALRS_dyadic_prec_round(rop,u*mpfr_get_prec(op1));
else
CGALRS_dyadic_set_prec(rop,u*mpfr_get_prec(op1));
CGAL_assertion_code(int round=)
mpfr_pow_ui(rop,op1,u,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_addmul(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
CGALRS_dyadic_srcptr op2){
CGALRS_dyadic_t temp;
CGALRS_dyadic_init(temp);
CGALRS_dyadic_mul(temp,op1,op2);
CGALRS_dyadic_add(rop,rop,temp);
CGALRS_dyadic_clear(temp);
}
inline void CGALRS_dyadic_addmul_si(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
long op2){
CGALRS_dyadic_t temp;
CGALRS_dyadic_init(temp);
CGALRS_dyadic_mul_si(temp,op1,op2);
CGALRS_dyadic_add(rop,rop,temp);
CGALRS_dyadic_clear(temp);
}
inline void CGALRS_dyadic_addmul_ui(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long u){
//CGALRS_dyadic_t temp;
//CGALRS_dyadic_init(temp);
//CGALRS_dyadic_mul_ui(temp,op1,u);
//CGALRS_dyadic_add(rop,rop,temp);
//CGALRS_dyadic_clear(temp);
CGALRS_dyadic_t temp;
mp_exp_t l,r,temp1,temp2;
mp_prec_t rop_prec;
if(u==0||mpfr_zero_p(op1))
return;
if(u==1){
CGALRS_dyadic_add(rop,rop,op1);
return;
}
// TODO: if(op1==1)
// calculate temp=op1*u
mpfr_init2(temp,mpfr_get_prec(op1)+sizeof(unsigned int));
CGAL_assertion_code(int round1=)
mpfr_mul_ui(temp,op1,u,GMP_RNDN);
CGAL_assertion(!round1);
// calculate the precision needed for rop
l=mpfr_get_exp(op1)>0?mpfr_get_exp(op1):0;
temp1=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1);
temp2=sizeof(unsigned long);
r=temp1>temp2?temp2:temp1;
CGAL_assertion(l>r);
rop_prec=sizeof(unsigned long)+1+(mp_prec_t)(l-r);
CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
rop_prec>=mpfr_get_prec(rop));
// set precision and add
CGALRS_dyadic_prec_round(rop,rop_prec);
CGAL_assertion_code(int round2=)
mpfr_add(rop,rop,temp,GMP_RNDN);
CGAL_assertion(!round2);
}
inline void CGALRS_dyadic_mul_2exp(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long ui){
// mpfr_mul_2ui does change the mantissa!
if(rop==op1)
CGALRS_dyadic_prec_round(
rop,
sizeof(unsigned long)+mpfr_get_prec(op1));
else
CGALRS_dyadic_set_prec(
rop,
sizeof(unsigned long)+mpfr_get_prec(op1));
CGAL_assertion_code(int round=)
mpfr_mul_2ui(rop,op1,ui,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_div_2exp(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long ui){
// mpfr_div_2ui does not change the mantissa... am I sure?
CGAL_assertion_code(int round=)
mpfr_div_2ui(rop,op1,ui,GMP_RNDN);
CGAL_assertion(!round);
}
// miscellaneous functions
#define CGALRS_dyadic_midpoint(R,D,E) \
( CGALRS_dyadic_ll_add(R,D,E,1) , mpfr_div_2ui(R,R,1,GMP_RNDN) )
#define CGALRS_dyadic_swap(D,E) mpfr_swap(D,E)
// I/O functions
#define CGALRS_dyadic_out_str(F,D) mpfr_out_str(F,10,0,D,GMP_RNDN)
#ifdef __cplusplus
inline std::ostream& operator<<(std::ostream &s,CGALRS_dyadic_srcptr op){
mp_exp_t exponent;
mpz_t mantissa;
mpz_init(mantissa);
exponent=mpfr_get_z_exp(mantissa,op);
s<<"["<<mantissa<<","<<exponent<<"]";
return s;
}
#endif
#endif // CGAL_RS_DYADIC_H

View File

@ -0,0 +1,52 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_EXACT_SIGNAT_1_H
#define CGAL_RS_EXACT_SIGNAT_1_H
#include <CGAL/Polynomial_traits_d.h>
#include "dyadic.h"
namespace CGAL{
namespace RS_AK1{
template <class Polynomial_,class Bound_>
struct ExactSignat_1{
typedef Polynomial_ Polynomial;
typedef Bound_ Bound;
typedef CGAL::Polynomial_traits_d<Polynomial> PT;
typedef typename PT::Degree Degree;
Polynomial pol;
ExactSignat_1(const Polynomial &p):pol(p){};
CGAL::Sign operator()(const Bound&)const;
}; // struct ExactSignat_1
template <>
inline CGAL::Sign
ExactSignat_1<Polynomial<Gmpz>,Gmpfr>::operator()(const Gmpfr &x)const{
int d=Degree()(pol);
if(d==0)
return pol[0].sign();
// Construct a Gmpfr containing exactly the leading coefficient.
Gmpfr h(pol[d],pol[d].bit_size());
CGAL_assertion(h==pol[d]);
// Horner's evaluation.
for(int i=1;i<d;++i){
CGALRS_dyadic_mul(h.fr(),h.fr(),x.fr());
CGALRS_dyadic_add_z(h.fr(),h.fr(),pol[d-i].mpz());
}
// TODO: We can avoid doing the last addition.
return h.sign();
}
} // namespace RS_AK1
} // namespace CGAL
#endif // CGAL_RS_EXACT_SIGNAT_1_H

View File

@ -0,0 +1,648 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_FUNCTORS_1_H
#define CGAL_RS_FUNCTORS_1_H
#include <vector>
#include <CGAL/Gmpfi.h>
namespace CGAL{
namespace RS_AK1{
template <class Polynomial_,
class Algebraic_,
class Bound_,
class Coefficient_,
class Isolator_>
struct Construct_algebraic_real_1{
typedef Polynomial_ Polynomial;
typedef Algebraic_ Algebraic;
typedef Bound_ Bound;
typedef Coefficient_ Coefficient;
typedef Isolator_ Isolator;
typedef Algebraic result_type;
template <class T>
Algebraic operator()(const T &a)const{
return Algebraic(a);
}
Algebraic operator()(const Polynomial &p,size_t i)const{
Isolator isol(p);
return Algebraic(p,isol.left_bound(i),isol.right_bound(i));
}
Algebraic operator()(const Polynomial &p,
const Bound &l,
const Bound &r)const{
return Algebraic(p,l,r);
}
}; // struct Construct_algebraic_1
template <class Polynomial_,class Algebraic_>
struct Compute_polynomial_1:
public CGAL::cpp98::unary_function<Algebraic_,Polynomial_>{
typedef Polynomial_ Polynomial;
typedef Algebraic_ Algebraic;
Polynomial operator()(const Algebraic &x)const{
return x.get_pol();
}
}; // struct Compute_polynomial_1
template <class Polynomial_,class Ptraits_>
struct Is_coprime_1:
public CGAL::cpp98::binary_function<Polynomial_,Polynomial_,bool>{
typedef Polynomial_ Polynomial;
typedef Ptraits_ Ptraits;
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Degree Degree;
inline bool operator()(const Polynomial &p1,const Polynomial &p2)const{
return Degree()(Gcd()(p1,p2))==0;
}
}; // struct Is_coprime_1
template <class Polynomial_,class Ptraits_>
struct Make_coprime_1{
typedef Polynomial_ Polynomial;
typedef Ptraits_ Ptraits;
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Degree Degree;
typedef typename Ptraits::Integral_division_up_to_constant_factor
IDiv;
bool operator()(const Polynomial &p1,
const Polynomial &p2,
Polynomial &g,
Polynomial &q1,
Polynomial &q2)const{
g=Gcd()(p1,p2);
q1=IDiv()(p1,g);
q2=IDiv()(p2,g);
return Degree()(Gcd()(p1,p2))==0;
}
}; // struct Make_coprime_1
template <class Polynomial_,
class Bound_,
class Algebraic_,
class Isolator_,
class Signat_,
class Ptraits_>
struct Solve_1{
typedef Polynomial_ Polynomial_1;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Isolator_ Isolator;
typedef Signat_ Signat;
typedef Ptraits_ Ptraits;
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Square_free_factorize_up_to_constant_factor
Sqfr;
typedef typename Ptraits::Degree Degree;
typedef typename Ptraits::Make_square_free Sfpart;
template <class OutputIterator>
OutputIterator operator()(const Polynomial_1 &p,
OutputIterator res)const{
typedef std::pair<Polynomial_1,int> polmult;
typedef std::vector<polmult> sqvec;
Polynomial_1 sfp=Sfpart()(p);
sqvec sfv;
Sqfr()(p,std::back_inserter(sfv));
Isolator isol(sfp);
int *m=(int*)calloc(isol.number_of_real_roots(),sizeof(int));
for(typename sqvec::iterator i=sfv.begin();i!=sfv.end();++i){
int k=Degree()(i->first);
Signat signof(i->first);
for(int j=0;k&&j<isol.number_of_real_roots();++j){
if(!m[j]){
CGAL::Sign sg_l=
signof(isol.left_bound(j));
CGAL::Sign sg_r=
signof(isol.right_bound(j));
if((sg_l!=sg_r)||
((sg_l==CGAL::ZERO)&&
(sg_r==CGAL::ZERO))){
m[j]=i->second;
--k;
}
}
}
}
for(int l=0;l<isol.number_of_real_roots();++l)
*res++=std::make_pair(Algebraic(p,
isol.left_bound(l),
isol.right_bound(l)),
m[l]);
free(m);
return res;
}
template <class OutputIterator>
OutputIterator operator()(const Polynomial_1 &p,
bool,
OutputIterator res)const{
Isolator isol(p);
for(int l=0;l<isol.number_of_real_roots();++l)
*res++=Algebraic(p,
isol.left_bound(l),
isol.right_bound(l));
return res;
}
template <class OutputIterator>
OutputIterator operator()(const Polynomial_1 &p,
const Bound &l,
const Bound &u,
OutputIterator res)const{
typedef std::vector<Algebraic> RV;
typedef std::pair<Polynomial_1,int> PM;
typedef std::vector<PM> PMV;
typedef typename PMV::iterator PMVI;
CGAL_precondition_msg(l<=u,
"left bound must be <= right bound");
RV roots; // all roots of the polynomial
this->operator()(p,false,std::back_inserter(roots));
size_t nb_roots=roots.size();
// indices of the first and last roots to be reported:
size_t index_l=0,index_u;
while(index_l<nb_roots&&roots[index_l]<l)
++index_l;
CGAL_assertion(index_l<=nb_roots);
if(index_l==nb_roots)
return res;
index_u=index_l;
while(index_u<nb_roots&&roots[index_u]<u)
++index_u;
CGAL_assertion(index_u<=nb_roots);
if(index_u==index_l)
return res;
// now, we have to return roots in [index_l,index_u)
PMV sfv;
Sqfr()(p,std::back_inserter(sfv)); // square-free fact. of p
// array to store the multiplicities
int *m=(int*)calloc(nb_roots,sizeof(int));
// we iterate over all the pairs <root,mult> and match the
// roots in the interval [index_l,index_u)
for(PMVI i=sfv.begin();i!=sfv.end();++i){
int k=Degree()(i->first);
Signat signof(i->first);
for(size_t j=index_l;k&&j<index_u;++j){
if(!m[j]){
CGAL::Sign sg_l=
signof(roots[j].get_left());
CGAL::Sign sg_r=
signof(roots[j].get_right());
if((sg_l!=sg_r)||
((sg_l==CGAL::ZERO)&&
(sg_r==CGAL::ZERO))){
m[j]=i->second;
--k;
}
}
}
}
for(size_t l=index_l;l<index_u;++l)
*res++=std::make_pair(roots[l],m[l]);
free(m);
return res;
}
template <class OutputIterator>
OutputIterator operator()(const Polynomial_1 &p,
bool known_to_be_square_free,
const Bound &l,
const Bound &u,
OutputIterator res)const{
typedef std::vector<Algebraic> RV;
typedef typename RV::iterator RVI;
CGAL_precondition_msg(l<=u,
"left bound must be <= right bound");
RV roots;
this->operator()(p,
known_to_be_square_free,
std::back_inserter(roots));
for(RVI it=roots.begin();it!=roots.end();it++)
if(*it>=l&&*it<=u)
*res++=*it;
return res;
}
}; // struct Solve_1
template <class Polynomial_,
class Bound_,
class Algebraic_,
class Refiner_,
class Signat_,
class Ptraits_>
class Sign_at_1:
public CGAL::cpp98::binary_function<Polynomial_,Algebraic_,CGAL::Sign>{
// This implementation will work with any polynomial type whose
// coefficient type is explicit interoperable with Gmpfi.
// TODO: Make this function generic.
public:
typedef Polynomial_ Polynomial_1;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Refiner_ Refiner;
typedef Signat_ Signat;
typedef Ptraits_ Ptraits;
private:
CGAL::Uncertain<CGAL::Sign> eval_interv(const Polynomial_1 &p,
const Bound &l,
const Bound &r)const{
typedef typename Ptraits::Substitute Subst;
std::vector<CGAL::Gmpfi> substitutions;
substitutions.push_back(CGAL::Gmpfi(l,r));
CGAL::Gmpfi eval=Subst()(p,
substitutions.begin(),
substitutions.end());
return eval.sign();
}
// This function assumes that the sign of the evaluation is not zero,
// it just refines x until having the correct sign.
CGAL::Sign refine_and_return(const Polynomial_1 &p,Algebraic x)const{
CGAL::Gmpfr xl(x.get_left());
CGAL::Gmpfr xr(x.get_right());
CGAL::Uncertain<CGAL::Sign> s;
for(;;){
Refiner()(x.get_pol(),
xl,
xr,
2*CGAL::max(xl.get_precision(),
xr.get_precision()));
s=eval_interv(p,xl,xr);
if(!s.is_same(Uncertain<CGAL::Sign>::indeterminate())){
x.set_left(xl);
x.set_right(xr);
return s;
}
}
}
public:
CGAL::Sign operator()(const Polynomial_1 &p,Algebraic x)const{
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Make_square_free Sfpart;
typedef typename Ptraits::Degree Degree;
typedef typename Ptraits::Differentiate Deriv;
CGAL::Uncertain<CGAL::Sign> unknown=
Uncertain<CGAL::Sign>::indeterminate();
CGAL::Uncertain<CGAL::Sign> s=eval_interv(p,
x.get_left(),
x.get_right());
if(!s.is_same(unknown))
return s;
// We are not sure about the sign. We calculate the gcd in
// order to know if both polynomials have common roots.
Polynomial_1 sfpp=Sfpart()(p);
Polynomial_1 gcd=Gcd()(sfpp,Sfpart()(x.get_pol()));
if(Degree()(gcd)==0)
return refine_and_return(p,x);
// At this point, gcd is not 1; we proceed as follows:
// -we refine x until having p monotonic in x's interval (to be
// sure that p has at most one root on that interval),
// -if the gcd has a root on this interval, both roots are
// equal (we return 0), otherwise, we refine until having a
// result.
// How to assure that p is monotonic in an interval: when its
// derivative is never zero in that interval.
Polynomial_1 dsfpp=Deriv()(sfpp);
CGAL::Gmpfr xl(x.get_left());
CGAL::Gmpfr xr(x.get_right());
while(eval_interv(dsfpp,xl,xr).is_same(unknown)){
Refiner()(x.get_pol(),
xl,
xr,
2*CGAL::max(xl.get_precision(),
xr.get_precision()));
}
x.set_left(xl);
x.set_right(xr);
// How to know that the gcd has a root: evaluate endpoints.
CGAL::Sign sleft,sright;
Signat sign_at_gcd(gcd);
if((sleft=sign_at_gcd(x.get_left()))==CGAL::ZERO||
(sright=sign_at_gcd(x.get_right()))==CGAL::ZERO||
(sleft!=sright))
return CGAL::ZERO;
return refine_and_return(p,x);
}
}; // struct Sign_at_1
template <class Polynomial_,
class Bound_,
class Algebraic_,
class Refiner_,
class Signat_,
class Ptraits_>
class Is_zero_at_1:
public CGAL::cpp98::binary_function<Polynomial_,Algebraic_,bool>{
// This implementation will work with any polynomial type whose
// coefficient type is explicit interoperable with Gmpfi.
// TODO: Make this function generic.
public:
typedef Polynomial_ Polynomial_1;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Refiner_ Refiner;
typedef Signat_ Signat;
typedef Ptraits_ Ptraits;
private:
CGAL::Uncertain<CGAL::Sign> eval_interv(const Polynomial_1 &p,
const Bound &l,
const Bound &r)const{
typedef typename Ptraits::Substitute Subst;
std::vector<CGAL::Gmpfi> substitutions;
substitutions.push_back(CGAL::Gmpfi(l,r));
CGAL::Gmpfi eval=Subst()(p,
substitutions.begin(),
substitutions.end());
return eval.sign();
}
public:
bool operator()(const Polynomial_1 &p,Algebraic x)const{
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Make_square_free Sfpart;
typedef typename Ptraits::Degree Degree;
typedef typename Ptraits::Differentiate Deriv;
CGAL::Uncertain<CGAL::Sign> unknown=
Uncertain<CGAL::Sign>::indeterminate();
CGAL::Uncertain<CGAL::Sign> s=eval_interv(p,
x.get_left(),
x.get_right());
if(!s.is_same(unknown))
return (s==CGAL::ZERO);
// We are not sure about the sign. We calculate the gcd in
// order to know if both polynomials have common roots.
Polynomial_1 sfpp=Sfpart()(p);
Polynomial_1 gcd=Gcd()(sfpp,Sfpart()(x.get_pol()));
if(Degree()(gcd)==0)
return false;
// At this point, gcd is not 1; we proceed as follows:
// -we refine x until having p monotonic in x's interval (to be
// sure that p has at most one root on that interval),
// -if the gcd has a root on this interval, both roots are
// equal (we return 0), otherwise, we refine until having a
// result.
// How to assure that p is monotonic in an interval: when its
// derivative is never zero in that interval.
Polynomial_1 dsfpp=Deriv()(sfpp);
CGAL::Gmpfr xl(x.get_left());
CGAL::Gmpfr xr(x.get_right());
while(eval_interv(dsfpp,xl,xr).is_same(unknown)){
Refiner()(x.get_pol(),
xl,
xr,
2*CGAL::max(xl.get_precision(),
xr.get_precision()));
}
x.set_left(xl);
x.set_right(xr);
// How to know that the gcd has a root: evaluate endpoints.
CGAL::Sign sleft,sright;
Signat sign_at_gcd(gcd);
return((sleft=sign_at_gcd(x.get_left()))==CGAL::ZERO||
(sright=sign_at_gcd(x.get_right()))==CGAL::ZERO||
(sleft!=sright));
}
}; // class Is_zero_at_1
// TODO: it says in the manual that this should return a size_type, but test
// programs assume that this is equal to int
template <class Polynomial_,class Isolator_>
struct Number_of_solutions_1:
public CGAL::cpp98::unary_function<Polynomial_,int>{
typedef Polynomial_ Polynomial_1;
typedef Isolator_ Isolator;
size_t operator()(const Polynomial_1 &p)const{
// TODO: make sure that p is square free (precondition)
Isolator isol(p);
return isol.number_of_real_roots();
}
}; // struct Number_of_solutions_1
// This functor not only compares two algebraic numbers. In case they are
// different, it refines them until they do not overlap.
template <class Algebraic_,
class Bound_,
class Comparator_>
struct Compare_1:
public CGAL::cpp98::binary_function<Algebraic_,Algebraic_,CGAL::Comparison_result>{
typedef Algebraic_ Algebraic;
typedef Bound_ Bound;
typedef Comparator_ Comparator;
CGAL::Comparison_result operator()(Algebraic a,Algebraic b)const{
Bound al=a.get_left();
Bound ar=a.get_right();
Bound bl=b.get_left();
Bound br=b.get_right();
CGAL::Comparison_result c=Comparator()(a.get_pol(),al,ar,
b.get_pol(),bl,br);
a.set_left(al);
a.set_right(ar);
b.set_left(bl);
b.set_right(br);
return c;
}
CGAL::Comparison_result operator()(Algebraic a,const Bound &b)const{
Bound al=a.get_left();
Bound ar=a.get_right();
Algebraic balg(b);
CGAL::Comparison_result c=Comparator()(a.get_pol(),al,ar,
balg.get_pol(),b,b);
a.set_left(al);
a.set_right(ar);
return c;
}
template <class T>
CGAL::Comparison_result operator()(Algebraic a,const T &b)const{
Bound al=a.get_left();
Bound ar=a.get_right();
Algebraic balg(b);
CGAL::Comparison_result c=Comparator()(a.get_pol(),
al,
ar,
balg.get_pol(),
balg.get_left(),
balg.get_right());
a.set_left(al);
a.set_right(ar);
return c;
}
}; // class Compare_1
template <class Algebraic_,
class Bound_,
class Comparator_>
struct Bound_between_1:
public CGAL::cpp98::binary_function<Algebraic_,Algebraic_,Bound_>{
typedef Algebraic_ Algebraic;
typedef Bound_ Bound;
typedef Comparator_ Comparator;
Bound operator()(Algebraic a,Algebraic b)const{
typedef Compare_1<Algebraic,Bound,Comparator> Compare;
typename Bound::Precision_type prec;
switch(Compare()(a,b)){
case CGAL::LARGER:
CGAL_assertion(b.get_right()<a.get_left());
prec=CGAL::max(b.get_right().get_precision(),
a.get_left().get_precision());
return Bound::add(b.get_right(),
a.get_left(),
1+prec)/2;
break;
case CGAL::SMALLER:
CGAL_assertion(a.get_right()<b.get_left());
prec=CGAL::max(a.get_right().get_precision(),
b.get_left().get_precision());
return Bound::add(a.get_right(),
b.get_left(),
1+prec)/2;
break;
default:
CGAL_error_msg(
"bound between two equal numbers");
}
}
}; // struct Bound_between_1
template <class Polynomial_,
class Bound_,
class Algebraic_,
class Isolator_,
class Comparator_,
class Signat_,
class Ptraits_>
struct Isolate_1:
public CGAL::cpp98::binary_function<Algebraic_,Polynomial_,std::pair<Bound_,Bound_> >{
typedef Polynomial_ Polynomial_1;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Isolator_ Isolator;
typedef Comparator_ Comparator;
typedef Signat_ Signat;
typedef Ptraits_ Ptraits;
std::pair<Bound,Bound>
operator()(const Algebraic &a,const Polynomial_1 &p)const{
std::vector<Algebraic> roots;
std::back_insert_iterator<std::vector<Algebraic> > rit(roots);
typedef Solve_1<Polynomial_1,
Bound,
Algebraic,
Isolator,
Signat,
Ptraits> Solve;
typedef Compare_1<Algebraic,Bound,Comparator> Compare;
Solve()(p,false,rit);
for(typename std::vector<Algebraic>::size_type i=0;
i<roots.size();
++i){
// we use the comparison functor, that makes both
// intervals disjoint iff the algebraic numbers they
// represent are not equal
Compare()(a,roots[i]);
}
return std::make_pair(a.left(),a.right());
}
}; // Isolate_1
template <class Polynomial_,
class Bound_,
class Algebraic_,
class Refiner_>
struct Approximate_absolute_1:
public CGAL::cpp98::binary_function<Algebraic_,int,std::pair<Bound_,Bound_> >{
typedef Polynomial_ Polynomial_1;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Refiner_ Refiner;
// This implementation assumes that Bound is Gmpfr.
// TODO: make generic.
std::pair<Bound,Bound> operator()(const Algebraic &x,const int &a)const{
Bound xl(x.get_left()),xr(x.get_right());
// refsteps=log2(xl-xr)
mpfr_t temp;
mpfr_init(temp);
mpfr_sub(temp,xr.fr(),xl.fr(),GMP_RNDU);
mpfr_log2(temp,temp,GMP_RNDU);
long refsteps=mpfr_get_si(temp,GMP_RNDU);
mpfr_clear(temp);
Refiner()(x.get_pol(),xl,xr,CGAL::abs(refsteps+a));
x.set_left(xl);
x.set_right(xr);
CGAL_assertion(a>0?
(xr-xl)*CGAL::ipower(Bound(2),a)<=Bound(1):
(xr-xl)<=CGAL::ipower(Bound(2),-a));
return std::make_pair(xl,xr);
}
}; // struct Approximate_absolute_1
template <class Polynomial_,
class Bound_,
class Algebraic_,
class Refiner_>
struct Approximate_relative_1:
public CGAL::cpp98::binary_function<Algebraic_,int,std::pair<Bound_,Bound_> >{
typedef Polynomial_ Polynomial_1;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Refiner_ Refiner;
std::pair<Bound,Bound> operator()(const Algebraic &x,const int &a)const{
if(CGAL::is_zero(x))
return std::make_pair(Bound(0),Bound(0));
Bound error=CGAL::ipower(Bound(2),CGAL::abs(a));
Bound xl(x.get_left()),xr(x.get_right());
Bound max_b=(CGAL::max)(CGAL::abs(xr),CGAL::abs(xl));
while(a>0?(xr-xl)*error>max_b:(xr-xl)>error*max_b){
Refiner()(x.get_pol(),
xl,
xr,
std::max<unsigned>(
CGAL::abs(a),
CGAL::max(xl.get_precision(),
xr.get_precision())));
max_b=(CGAL::max)(CGAL::abs(xr),CGAL::abs(xl));
}
x.set_left(xl);
x.set_right(xr);
CGAL_assertion(
a>0?
(xr-xl)*CGAL::ipower(Bound(2),a)<=max_b:
(xr-xl)<=CGAL::ipower(Bound(2),-a)*max_b);
return std::make_pair(xl,xr);
}
}; // struct Approximate_relative_1
} // namespace RS_AK1
} // namespace CGAL
#endif // CGAL_RS_FUNCTORS_1_H

View File

@ -0,0 +1,699 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_FUNCTORS_Z_1_H
#define CGAL_RS_FUNCTORS_Z_1_H
#include <vector>
#include <CGAL/Gmpfi.h>
namespace CGAL{
namespace RS_AK1{
template <class Polynomial_,
class ZPolynomial_,
class PolConverter_,
class Algebraic_,
class Bound_,
class Coefficient_,
class Isolator_>
struct Construct_algebraic_real_z_1{
typedef Polynomial_ Polynomial;
typedef ZPolynomial_ ZPolynomial;
typedef PolConverter_ PolConverter;
typedef Algebraic_ Algebraic;
typedef Bound_ Bound;
typedef Coefficient_ Coefficient;
typedef Isolator_ Isolator;
typedef Algebraic result_type;
template <class T>
Algebraic operator()(const T &a)const{
return Algebraic(a);
}
Algebraic operator()(const Polynomial &p,size_t i)const{
ZPolynomial zp=PolConverter()(p);
Isolator isol(zp);
return Algebraic(p,
zp,
isol.left_bound(i),
isol.right_bound(i));
}
Algebraic operator()(const Polynomial &p,
const Bound &l,
const Bound &r)const{
return Algebraic(p,PolConverter()(p),l,r);
}
}; // struct Construct_algebraic_real_z_1
template <class Polynomial_,class Algebraic_>
struct Compute_polynomial_z_1:
public CGAL::cpp98::unary_function<Algebraic_,Polynomial_>{
typedef Polynomial_ Polynomial;
typedef Algebraic_ Algebraic;
Polynomial operator()(const Algebraic &x)const{
return x.get_pol();
}
}; // struct Compute_polynomial_z_1
template <class Polynomial_,class Ptraits_>
struct Is_coprime_z_1:
public CGAL::cpp98::binary_function<Polynomial_,Polynomial_,bool>{
typedef Polynomial_ Polynomial;
typedef Ptraits_ Ptraits;
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Degree Degree;
inline bool operator()(const Polynomial &p1,const Polynomial &p2)const{
return Degree()(Gcd()(p1,p2))==0;
}
}; // struct Is_coprime_z_1
template <class Polynomial_,class Ptraits_>
struct Make_coprime_z_1{
typedef Polynomial_ Polynomial;
typedef Ptraits_ Ptraits;
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Degree Degree;
typedef typename Ptraits::Integral_division_up_to_constant_factor
IDiv;
bool operator()(const Polynomial &p1,
const Polynomial &p2,
Polynomial &g,
Polynomial &q1,
Polynomial &q2)const{
g=Gcd()(p1,p2);
q1=IDiv()(p1,g);
q2=IDiv()(p2,g);
return Degree()(Gcd()(p1,p2))==0;
}
}; // struct Make_coprime_z_1
template <class Polynomial_,
class ZPolynomial_,
class PolConverter_,
class Bound_,
class Algebraic_,
class Isolator_,
class Signat_,
class Ptraits_,
class ZPtraits_>
struct Solve_z_1{
typedef Polynomial_ Polynomial_1;
typedef ZPolynomial_ ZPolynomial_1;
typedef PolConverter_ PolConverter;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Isolator_ Isolator;
typedef Signat_ ZSignat;
typedef Ptraits_ Ptraits;
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Square_free_factorize_up_to_constant_factor
Sqfr;
typedef typename Ptraits::Degree Degree;
typedef typename Ptraits::Make_square_free Sfpart;
typedef ZPtraits_ ZPtraits;
typedef typename ZPtraits::Gcd_up_to_constant_factor ZGcd;
typedef typename ZPtraits::Square_free_factorize_up_to_constant_factor
ZSqfr;
typedef typename ZPtraits::Degree ZDegree;
typedef typename ZPtraits::Make_square_free ZSfpart;
template <class OutputIterator>
OutputIterator operator()(const Polynomial_1 &p,
OutputIterator res)const{
typedef std::pair<ZPolynomial_1,int> zpolmult;
typedef std::vector<zpolmult> zsqvec;
ZPolynomial_1 zp=PolConverter()(p);
Polynomial_1 sfp=Sfpart()(p);
ZPolynomial_1 zsfp=PolConverter()(sfp);
zsqvec zsfv;
ZSqfr()(zp,std::back_inserter(zsfv));
Isolator isol(zsfp);
int *m=(int*)calloc(isol.number_of_real_roots(),sizeof(int));
for(typename zsqvec::iterator i=zsfv.begin();i!=zsfv.end();++i){
int k=ZDegree()(i->first);
ZSignat signof(i->first);
for(int j=0;k&&j<isol.number_of_real_roots();++j){
if(!m[j]){
CGAL::Sign sg_l=
signof(isol.left_bound(j));
CGAL::Sign sg_r=
signof(isol.right_bound(j));
if((sg_l!=sg_r)||
((sg_l==CGAL::ZERO)&&
(sg_r==CGAL::ZERO))){
m[j]=i->second;
--k;
}
}
}
}
for(int l=0;l<isol.number_of_real_roots();++l)
*res++=std::make_pair(Algebraic(sfp,
zsfp,
isol.left_bound(l),
isol.right_bound(l)),
m[l]);
free(m);
return res;
}
template <class OutputIterator>
OutputIterator operator()(const Polynomial_1 &p,
bool,
OutputIterator res)const{
ZPolynomial_1 zp=PolConverter()(p);
Isolator isol(zp);
for(int l=0;l<isol.number_of_real_roots();++l)
*res++=Algebraic(p,
zp,
isol.left_bound(l),
isol.right_bound(l));
return res;
}
template <class OutputIterator>
OutputIterator operator()(const Polynomial_1 &p,
const Bound &l,
const Bound &u,
OutputIterator res)const{
typedef std::vector<Algebraic> RV;
typedef std::pair<Polynomial_1,int> PM;
typedef std::vector<PM> PMV;
typedef typename PMV::iterator PMVI;
CGAL_precondition_msg(l<=u,
"left bound must be <= right bound");
RV roots; // all roots of the polynomial
this->operator()(p,false,std::back_inserter(roots));
size_t nb_roots=roots.size();
// indices of the first and last roots to be reported:
size_t index_l=0,index_u;
while(index_l<nb_roots&&roots[index_l]<l)
++index_l;
CGAL_assertion(index_l<=nb_roots);
if(index_l==nb_roots)
return res;
index_u=index_l;
while(index_u<nb_roots&&roots[index_u]<u)
++index_u;
CGAL_assertion(index_u<=nb_roots);
if(index_u==index_l)
return res;
// now, we have to return roots in [index_l,index_u)
PMV sfv;
Sqfr()(p,std::back_inserter(sfv)); // square-free fact. of p
// array to store the multiplicities
int *m=(int*)calloc(nb_roots,sizeof(int));
// we iterate over all the pairs <root,mult> and match the
// roots in the interval [index_l,index_u)
for(PMVI i=sfv.begin();i!=sfv.end();++i){
ZPolynomial_1 zifirst=PolConverter()(i->first);
int k=ZDegree()(zifirst);
ZSignat signof(zifirst);
for(size_t j=index_l;k&&j<index_u;++j){
if(!m[j]){
CGAL::Sign sg_l=
signof(roots[j].get_left());
CGAL::Sign sg_r=
signof(roots[j].get_right());
if((sg_l!=sg_r)||
((sg_l==CGAL::ZERO)&&
(sg_r==CGAL::ZERO))){
m[j]=i->second;
--k;
}
}
}
}
for(size_t l=index_l;l<index_u;++l)
*res++=std::make_pair(roots[l],m[l]);
free(m);
return res;
}
template <class OutputIterator>
OutputIterator operator()(const Polynomial_1 &p,
bool known_to_be_square_free,
const Bound &l,
const Bound &u,
OutputIterator res)const{
typedef std::vector<Algebraic> RV;
typedef typename RV::iterator RVI;
CGAL_precondition_msg(l<=u,
"left bound must be <= right bound");
RV roots;
this->operator()(p,
known_to_be_square_free,
std::back_inserter(roots));
for(RVI it=roots.begin();it!=roots.end();it++)
if(*it>=l&&*it<=u)
*res++=*it;
return res;
}
}; // struct Solve_z_1
template <class Polynomial_,
class ZPolynomial_,
class PolConverter_,
class Bound_,
class Algebraic_,
class Refiner_,
class Signat_,
class Ptraits_,
class ZPtraits_>
class Sign_at_z_1:
public CGAL::cpp98::binary_function<Polynomial_,Algebraic_,CGAL::Sign>{
// This implementation will work with any polynomial type whose
// coefficient type is explicit interoperable with Gmpfi.
// TODO: Make this function generic.
public:
typedef Polynomial_ Polynomial_1;
typedef ZPolynomial_ ZPolynomial_1;
typedef PolConverter_ PolConverter;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Refiner_ Refiner;
typedef Signat_ ZSignat;
typedef Ptraits_ Ptraits;
typedef ZPtraits_ ZPtraits;
private:
CGAL::Uncertain<CGAL::Sign> eval_interv(const Polynomial_1 &p,
const Bound &l,
const Bound &r)const{
typedef typename Ptraits::Substitute Subst;
std::vector<CGAL::Gmpfi> substitutions;
substitutions.push_back(CGAL::Gmpfi(l,r));
CGAL::Gmpfi eval=Subst()(p,
substitutions.begin(),
substitutions.end());
return eval.sign();
}
// This function assumes that the sign of the evaluation is not zero,
// it just refines x until having the correct sign.
CGAL::Sign refine_and_return(const Polynomial_1 &p,Algebraic x)const{
CGAL::Gmpfr xl(x.get_left());
CGAL::Gmpfr xr(x.get_right());
CGAL::Uncertain<CGAL::Sign> s;
for(;;){
Refiner()(x.get_zpol(),
xl,
xr,
2*CGAL::max(xl.get_precision(),
xr.get_precision()));
s=eval_interv(p,xl,xr);
if(!s.is_same(Uncertain<CGAL::Sign>::indeterminate())){
x.set_left(xl);
x.set_right(xr);
return s;
}
}
}
public:
CGAL::Sign operator()(const Polynomial_1 &p,Algebraic x)const{
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Make_square_free Sfpart;
typedef typename Ptraits::Degree Degree;
typedef typename Ptraits::Differentiate Deriv;
CGAL::Uncertain<CGAL::Sign> unknown=
Uncertain<CGAL::Sign>::indeterminate();
CGAL::Uncertain<CGAL::Sign> s=eval_interv(p,
x.get_left(),
x.get_right());
if(!s.is_same(unknown))
return s;
// We are not sure about the sign. We calculate the gcd in
// order to know if both polynomials have common roots.
Polynomial_1 sfpp=Sfpart()(p);
Polynomial_1 gcd=Gcd()(sfpp,Sfpart()(x.get_pol()));
if(Degree()(gcd)==0)
return refine_and_return(p,x);
// At this point, gcd is not 1; we proceed as follows:
// -we refine x until having p monotonic in x's interval (to be
// sure that p has at most one root on that interval),
// -if the gcd has a root on this interval, both roots are
// equal (we return 0), otherwise, we refine until having a
// result.
// How to assure that p is monotonic in an interval: when its
// derivative is never zero in that interval.
Polynomial_1 dsfpp=Deriv()(sfpp);
CGAL::Gmpfr xl(x.get_left());
CGAL::Gmpfr xr(x.get_right());
while(eval_interv(dsfpp,xl,xr).is_same(unknown)){
Refiner()(x.get_zpol(),
xl,
xr,
2*CGAL::max(xl.get_precision(),
xr.get_precision()));
}
x.set_left(xl);
x.set_right(xr);
// How to know that the gcd has a root: evaluate endpoints.
CGAL::Sign sleft,sright;
ZSignat sign_at_gcd(PolConverter()(gcd));
if((sleft=sign_at_gcd(x.get_left()))==CGAL::ZERO||
(sright=sign_at_gcd(x.get_right()))==CGAL::ZERO||
(sleft!=sright))
return CGAL::ZERO;
return refine_and_return(p,x);
}
}; // struct Sign_at_z_1
template <class Polynomial_,
class ZPolynomial_,
class PolConverter_,
class Bound_,
class Algebraic_,
class Refiner_,
class Signat_,
class Ptraits_,
class ZPtraits_>
class Is_zero_at_z_1:
public CGAL::cpp98::binary_function<Polynomial_,Algebraic_,bool>{
// This implementation will work with any polynomial type whose
// coefficient type is explicit interoperable with Gmpfi.
// TODO: Make this function generic.
public:
typedef Polynomial_ Polynomial_1;
typedef ZPolynomial_ ZPolynomial_1;
typedef PolConverter_ PolConverter;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Refiner_ Refiner;
typedef Signat_ ZSignat;
typedef Ptraits_ Ptraits;
typedef ZPtraits_ ZPtraits;
private:
CGAL::Uncertain<CGAL::Sign> eval_interv(const Polynomial_1 &p,
const Bound &l,
const Bound &r)const{
typedef typename Ptraits::Substitute Subst;
std::vector<CGAL::Gmpfi> substitutions;
substitutions.push_back(CGAL::Gmpfi(l,r));
CGAL::Gmpfi eval=Subst()(p,
substitutions.begin(),
substitutions.end());
return eval.sign();
}
public:
bool operator()(const Polynomial_1 &p,Algebraic x)const{
typedef typename Ptraits::Gcd_up_to_constant_factor Gcd;
typedef typename Ptraits::Make_square_free Sfpart;
typedef typename Ptraits::Degree Degree;
typedef typename Ptraits::Differentiate Deriv;
CGAL::Uncertain<CGAL::Sign> unknown=
Uncertain<CGAL::Sign>::indeterminate();
CGAL::Uncertain<CGAL::Sign> s=eval_interv(p,
x.get_left(),
x.get_right());
if(!s.is_same(unknown))
return (s==CGAL::ZERO);
// We are not sure about the sign. We calculate the gcd in
// order to know if both polynomials have common roots.
Polynomial_1 sfpp=Sfpart()(p);
Polynomial_1 gcd=Gcd()(sfpp,Sfpart()(x.get_pol()));
if(Degree()(gcd)==0)
return false;
// At this point, gcd is not 1; we proceed as follows:
// -we refine x until having p monotonic in x's interval (to be
// sure that p has at most one root on that interval),
// -if the gcd has a root on this interval, both roots are
// equal (we return 0), otherwise, we refine until having a
// result.
// How to assure that p is monotonic in an interval: when its
// derivative is never zero in that interval.
Polynomial_1 dsfpp=Deriv()(sfpp);
CGAL::Gmpfr xl(x.get_left());
CGAL::Gmpfr xr(x.get_right());
while(eval_interv(dsfpp,xl,xr).is_same(unknown)){
Refiner()(x.get_zpol(),
xl,
xr,
2*CGAL::max(xl.get_precision(),
xr.get_precision()));
}
x.set_left(xl);
x.set_right(xr);
// How to know that the gcd has a root: evaluate endpoints.
CGAL::Sign sleft,sright;
ZSignat sign_at_gcd(PolConverter()(gcd));
return((sleft=sign_at_gcd(x.get_left()))==CGAL::ZERO||
(sright=sign_at_gcd(x.get_right()))==CGAL::ZERO||
(sleft!=sright));
}
}; // class Is_zero_at_z_1
// TODO: it says in the manual that this should return a size_type, but test
// programs assume that this is equal to int
template <class Polynomial_,
class ZPolynomial_,
class PolConverter_,
class Isolator_>
struct Number_of_solutions_z_1:
public CGAL::cpp98::unary_function<Polynomial_,int>{
typedef Polynomial_ Polynomial_1;
typedef ZPolynomial_ ZPolynomial_1;
typedef PolConverter_ PolConverter;
typedef Isolator_ Isolator;
size_t operator()(const Polynomial_1 &p)const{
// TODO: make sure that p is square free (precondition)
Isolator isol(PolConverter()(p));
return isol.number_of_real_roots();
}
}; // struct Number_of_solutions_z_1
// This functor not only compares two algebraic numbers. In case they are
// different, it refines them until they do not overlap.
template <class Algebraic_,
class Bound_,
class Comparator_>
struct Compare_z_1:
public CGAL::cpp98::binary_function<Algebraic_,Algebraic_,CGAL::Comparison_result>{
typedef Algebraic_ Algebraic;
typedef Bound_ Bound;
typedef Comparator_ Comparator;
CGAL::Comparison_result operator()(Algebraic a,Algebraic b)const{
Bound al=a.get_left();
Bound ar=a.get_right();
Bound bl=b.get_left();
Bound br=b.get_right();
CGAL::Comparison_result c=Comparator()(a.get_zpol(),al,ar,
b.get_zpol(),bl,br);
a.set_left(al);
a.set_right(ar);
b.set_left(bl);
b.set_right(br);
return c;
}
CGAL::Comparison_result operator()(Algebraic a,const Bound &b)const{
Bound al=a.get_left();
Bound ar=a.get_right();
Algebraic balg(b);
CGAL::Comparison_result c=Comparator()(a.get_zpol(),al,ar,
balg.get_zpol(),b,b);
a.set_left(al);
a.set_right(ar);
return c;
}
template <class T>
CGAL::Comparison_result operator()(Algebraic a,const T &b)const{
Bound al=a.get_left();
Bound ar=a.get_right();
Algebraic balg(b);
CGAL::Comparison_result c=Comparator()(a.get_zpol(),
al,
ar,
balg.get_zpol(),
balg.get_left(),
balg.get_right());
a.set_left(al);
a.set_right(ar);
return c;
}
}; // class Compare_z_1
template <class Algebraic_,
class Bound_,
class Comparator_>
struct Bound_between_z_1:
public CGAL::cpp98::binary_function<Algebraic_,Algebraic_,Bound_>{
typedef Algebraic_ Algebraic;
typedef Bound_ Bound;
typedef Comparator_ Comparator;
Bound operator()(Algebraic a,Algebraic b)const{
typedef Compare_z_1<Algebraic,Bound,Comparator> Compare;
typename Bound::Precision_type prec;
switch(Compare()(a,b)){
case CGAL::LARGER:
CGAL_assertion(b.get_right()<a.get_left());
prec=CGAL::max(b.get_right().get_precision(),
a.get_left().get_precision());
return Bound::add(b.get_right(),
a.get_left(),
1+prec)/2;
break;
case CGAL::SMALLER:
CGAL_assertion(a.get_right()<b.get_left());
prec=CGAL::max(a.get_right().get_precision(),
b.get_left().get_precision());
return Bound::add(a.get_right(),
b.get_left(),
1+prec)/2;
break;
default:
CGAL_error_msg(
"bound between two equal numbers");
}
}
}; // struct Bound_between_z_1
template <class Polynomial_,
class ZPolynomial_,
class PolConverter_,
class Bound_,
class Algebraic_,
class Isolator_,
class Comparator_,
class Signat_,
class Ptraits_,
class ZPtraits_>
struct Isolate_z_1:
public CGAL::cpp98::binary_function<Algebraic_,Polynomial_,std::pair<Bound_,Bound_> >{
typedef Polynomial_ Polynomial_1;
typedef ZPolynomial_ ZPolynomial_1;
typedef PolConverter_ PolConverter;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Isolator_ Isolator;
typedef Comparator_ Comparator;
typedef Signat_ Signat;
typedef Ptraits_ Ptraits;
typedef ZPtraits_ ZPtraits;
std::pair<Bound,Bound>
operator()(const Algebraic &a,const Polynomial_1 &p)const{
std::vector<Algebraic> roots;
std::back_insert_iterator<std::vector<Algebraic> > rit(roots);
typedef Solve_z_1<Polynomial_1,
ZPolynomial_1,
PolConverter,
Bound,
Algebraic,
Isolator,
Signat,
Ptraits,
ZPtraits> Solve;
typedef Compare_z_1<Algebraic,Bound,Comparator> Compare;
Solve()(p,false,rit);
for(typename std::vector<Algebraic>::size_type i=0;
i<roots.size();
++i){
// we use the comparison functor, that makes both
// intervals disjoint iff the algebraic numbers they
// represent are not equal
Compare()(a,roots[i]);
}
return std::make_pair(a.left(),a.right());
}
}; // Isolate_z_1
template <class Polynomial_,
class Bound_,
class Algebraic_,
class Refiner_>
struct Approximate_absolute_z_1:
public CGAL::cpp98::binary_function<Algebraic_,int,std::pair<Bound_,Bound_> >{
typedef Polynomial_ Polynomial_1;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Refiner_ Refiner;
// This implementation assumes that Bound is Gmpfr.
// TODO: make generic.
std::pair<Bound,Bound> operator()(const Algebraic &x,const int &a)const{
Bound xl(x.get_left()),xr(x.get_right());
// refsteps=log2(xl-xr)
mpfr_t temp;
mpfr_init(temp);
mpfr_sub(temp,xr.fr(),xl.fr(),GMP_RNDU);
mpfr_log2(temp,temp,GMP_RNDU);
long refsteps=mpfr_get_si(temp,GMP_RNDU);
mpfr_clear(temp);
Refiner()(x.get_zpol(),xl,xr,CGAL::abs(refsteps+a));
x.set_left(xl);
x.set_right(xr);
CGAL_assertion(a>0?
(xr-xl)*CGAL::ipower(Bound(2),a)<=Bound(1):
(xr-xl)<=CGAL::ipower(Bound(2),-a));
return std::make_pair(xl,xr);
}
}; // struct Approximate_absolute_z_1
template <class Polynomial_,
class Bound_,
class Algebraic_,
class Refiner_>
struct Approximate_relative_z_1:
public CGAL::cpp98::binary_function<Algebraic_,int,std::pair<Bound_,Bound_> >{
typedef Polynomial_ Polynomial_1;
typedef Bound_ Bound;
typedef Algebraic_ Algebraic;
typedef Refiner_ Refiner;
std::pair<Bound,Bound> operator()(const Algebraic &x,const int &a)const{
if(CGAL::is_zero(x))
return std::make_pair(Bound(0),Bound(0));
Bound error=CGAL::ipower(Bound(2),CGAL::abs(a));
Bound xl(x.get_left()),xr(x.get_right());
Bound max_b=(CGAL::max)(CGAL::abs(xr),CGAL::abs(xl));
while(a>0?(xr-xl)*error>max_b:(xr-xl)>error*max_b){
Refiner()(x.get_zpol(),
xl,
xr,
std::max<unsigned>(
CGAL::abs(a),
CGAL::max(xl.get_precision(),
xr.get_precision())));
max_b=(CGAL::max)(CGAL::abs(xr),CGAL::abs(xl));
}
x.set_left(xl);
x.set_right(xr);
CGAL_assertion(
a>0?
(xr-xl)*CGAL::ipower(Bound(2),a)<=max_b:
(xr-xl)<=CGAL::ipower(Bound(2),-a)*max_b);
return std::make_pair(xl,xr);
}
}; // struct Approximate_relative_z_1
} // namespace RS_AK1
} // namespace CGAL
#endif // CGAL_RS_FUNCTORS_Z_1_H

View File

@ -0,0 +1,49 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_POLYNOMIAL_CONVERTER_1_H
#define CGAL_RS_POLYNOMIAL_CONVERTER_1_H
namespace CGAL{
namespace RS_AK1{
template <class InputPolynomial_,class OutputPolynomial_>
struct Polynomial_converter_1:
public CGAL::cpp98::unary_function<InputPolynomial_,OutputPolynomial_>{
typedef InputPolynomial_ InpPolynomial_1;
typedef OutputPolynomial_ OutPolynomial_1;
OutPolynomial_1 operator()(const InpPolynomial_1&)const;
}; // class Polynomial_converter_1
template <>
Polynomial<Gmpz>
Polynomial_converter_1<Polynomial<Gmpq>,Polynomial<Gmpz> >::operator()(
const Polynomial<Gmpq> &p)const{
std::vector<Gmpz> outcoeffs;
unsigned degree=p.degree();
mpz_t lcm;
mpz_init(lcm);
mpz_lcm(lcm,mpq_denref(p[0].mpq()),mpq_denref(p[degree].mpq()));
for(unsigned i=1;i<degree;++i)
mpz_lcm(lcm,lcm,mpq_denref(p[i].mpq()));
for(unsigned i=0;i<=degree;++i){
Gmpz c;
mpz_divexact(c.mpz(),lcm,mpq_denref(p[i].mpq()));
mpz_mul(c.mpz(),c.mpz(),mpq_numref(p[i].mpq()));
outcoeffs.push_back(c);
}
mpz_clear(lcm);
return Polynomial<Gmpz>(outcoeffs.begin(),outcoeffs.end());
}
} // namespace RS_AK1
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_CONVERTER_1_H

View File

@ -0,0 +1,140 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_RS23_K_ISOLATOR_1_H
#define CGAL_RS_RS23_K_ISOLATOR_1_H
// This file includes an isolator. Its particularity is that is isolates the
// roots with RS2, and the refines them until reaching Kantorovich criterion.
// This can take long, but the later refinements will be extremely fast with
// RS3. The functor is not in RS2 neither in RS3 namespace, because it uses
// functions from both.
#include "rs2_calls.h"
#include "rs3_k_refiner_1.h"
#include <CGAL/Gmpfi.h>
#include <vector>
namespace CGAL{
template <class Polynomial_,class Bound_>
class RS23_k_isolator_1{
public:
typedef Polynomial_ Polynomial;
typedef Bound_ Bound;
private:
typedef Gmpfi Interval;
public:
RS23_k_isolator_1(const Polynomial&);
Polynomial polynomial()const;
int number_of_real_roots()const;
bool is_exact_root(int i)const;
Bound left_bound(int i)const;
Bound right_bound(int i)const;
private:
Polynomial _polynomial;
std::vector<Interval> _real_roots;
};
template <class Polynomial_,class Bound_>
RS23_k_isolator_1<Polynomial_,Bound_>::
RS23_k_isolator_1(const Polynomial_ &){
CGAL_error_msg("not implemented for these polynomial/bound types");
}
template <>
RS23_k_isolator_1<CGAL::Polynomial<CGAL::Gmpz>,Gmpfr>::
RS23_k_isolator_1(const CGAL::Polynomial<CGAL::Gmpz> &p):_polynomial(p){
typedef CGAL::Polynomial<CGAL::Gmpz> Pol;
typedef CGAL::Gmpfr Bound;
typedef CGAL::RS3::RS3_k_refiner_1<Pol,Bound> KRefiner;
int numsols;
std::vector<Gmpfi> intervals;
RS2::RS2_calls::init_solver();
RS2::RS2_calls::create_rs_upoly(p,rs_get_default_up());
set_rs_precisol(0);
set_rs_verbose(0);
rs_run_algo((char*)"UISOLE");
RS2::RS2_calls::insert_roots(std::back_inserter(intervals));
// RS2 computed the isolating intervals. Now, we use RS3 to refine each
// root until reaching Kantorovich criterion, before adding it to the
// root vector.
numsols=intervals.size();
for(int j=0;j<numsols;++j){
Gmpfr left(intervals[j].inf());
Gmpfr right(intervals[j].sup());
CGAL_assertion(left<=right);
KRefiner()(p,left,right,53);
_real_roots.push_back(intervals[j]);
}
}
template <>
RS23_k_isolator_1<CGAL::Polynomial<CGAL::Gmpq>,Gmpfr>::
RS23_k_isolator_1(const CGAL::Polynomial<CGAL::Gmpq> &qp):_polynomial(qp){
typedef CGAL::Polynomial<CGAL::Gmpz> ZPol;
typedef CGAL::Gmpfr Bound;
typedef CGAL::RS3::RS3_k_refiner_1<ZPol,Bound> ZKRefiner;
int numsols;
std::vector<Gmpfi> intervals;
CGAL::Polynomial<CGAL::Gmpz> zp=CGAL::RS_AK1::Polynomial_converter_1<
CGAL::Polynomial<Gmpq>,
CGAL::Polynomial<Gmpz> >()(qp);
RS2::RS2_calls::init_solver();
RS2::RS2_calls::create_rs_upoly(zp,rs_get_default_up());
set_rs_precisol(0);
set_rs_verbose(0);
rs_run_algo((char*)"UISOLE");
RS2::RS2_calls::insert_roots(std::back_inserter(intervals));
// RS2 computed the isolating intervals. Now, we use RS3 to refine each
// root until reaching Kantorovich criterion, before adding it to the
// root vector.
numsols=intervals.size();
for(int j=0;j<numsols;++j){
Gmpfr left(intervals[j].inf());
Gmpfr right(intervals[j].sup());
ZKRefiner()(zp,left,right,53);
_real_roots.push_back(intervals[j]);
}
}
template <class Polynomial_,class Bound_>
Polynomial_
RS23_k_isolator_1<Polynomial_,Bound_>::polynomial()const{
return _polynomial;
}
template <class Polynomial_,class Bound_>
int
RS23_k_isolator_1<Polynomial_,Bound_>::number_of_real_roots()const{
return _real_roots.size();
}
template <class Polynomial_,class Bound_>
bool
RS23_k_isolator_1<Polynomial_,Bound_>::is_exact_root(int i)const{
return _real_roots[i].inf()==_real_roots[i].sup();
}
template <class Polynomial_,class Bound_>
Bound_
RS23_k_isolator_1<Polynomial_,Bound_>::left_bound(int i)const{
return _real_roots[i].inf();
}
template <class Polynomial_,class Bound_>
Bound_
RS23_k_isolator_1<Polynomial_,Bound_>::right_bound(int i)const{
return _real_roots[i].sup();
}
} // namespace CGAL
#endif // CGAL_RS_RS23_K_ISOLATOR_1_H

View File

@ -0,0 +1,141 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_RS2_CALLS_H
#define CGAL_RS_RS2_CALLS_H
#include <CGAL/Gmpz.h>
#include <CGAL/Gmpfr.h>
#include <CGAL/Gmpfi.h>
#include <CGAL/Polynomial.h>
#include <CGAL/tss.h>
#include <rs_exports.h>
#ifdef CGAL_RS_OLD_INCLUDES
#define CGALRS_PTR(a) long int a
#else
#define CGALRS_PTR(a) void *a
#endif
// RS3 does not work with MPFR 3.1.3 to 3.1.6. In case RS3 is enabled and
// the version of MPFR is one of those buggy versions, abort the compilation
// and instruct the user to update MPFR or don't use RS3.
#ifdef CGAL_USE_RS3
static_assert(
MPFR_VERSION_MAJOR!=3 ||
MPFR_VERSION_MINOR!=1 ||
MPFR_VERSION_PATCHLEVEL<3 || MPFR_VERSION_PATCHLEVEL>6,
"RS3 does not work with MPFR versions 3.1.3 to 3.1.6. "
"Please update MPFR or disable RS3.");
#endif // CGAL_USE_RS3
namespace CGAL{
namespace RS2{
struct RS2_calls{
static void init_solver(){
CGAL_STATIC_THREAD_LOCAL_VARIABLE(bool, first,true);
if(first){
first=false;
rs_init_rs();
rs_reset_all();
}else
rs_reset_all();
}
static void create_rs_upoly(CGAL::Polynomial<CGAL::Gmpz> poly,
CGALRS_PTR(ident_pol)){
CGALRS_PTR(ident_mon);
CGALRS_PTR(ident_coeff);
rs_import_uppring((char*)"T");
for(int i=0;i<=poly.degree();++i)
if(mpz_sgn(poly[i].mpz())){ // don't add if == 0
ident_mon=rs_export_new_mon_upp_bz();
ident_coeff=rs_export_new_gmp();
rs_import_bz_gmp(ident_coeff,
TO_RSPTR_IN(&(poly[i].mpz())));
rs_dset_mon_upp_bz(ident_mon,ident_coeff,i);
rs_dappend_list_mon_upp_bz(ident_pol,
ident_mon);
}
}
static int affiche_sols_eqs(mpfi_ptr *x){
CGALRS_PTR(ident_sols_eqs);
CGALRS_PTR(ident_node);
CGALRS_PTR(ident_vect);
CGALRS_PTR(ident_elt);
int nb_elts;
ident_sols_eqs=rs_get_default_sols_eqs();
nb_elts=rs_export_list_vect_ibfr_nb(ident_sols_eqs);
ident_node=rs_export_list_vect_ibfr_firstnode(ident_sols_eqs);
mpfi_t *roots=(mpfi_t*)malloc(nb_elts*sizeof(mpfi_t));
for(int i=0;i<nb_elts;++i){
ident_vect=rs_export_list_vect_ibfr_monnode
(ident_node);
CGAL_assertion_msg(rs_export_dim_vect_ibfr
(ident_vect)==1,
"vector dimension must be 1");
ident_elt=rs_export_elt_vect_ibfr(ident_vect,0);
mpfi_ptr root_pointer=
(mpfi_ptr)rs_export_ibfr_mpfi(ident_elt);
mpfi_init2(roots[i],mpfi_get_prec(root_pointer));
mpfi_set(roots[i],root_pointer);
x[i]=roots[i];
// This doesn't work because RS relocates the
// mpfrs that form the mpfi. Nevertheless, the
// mpfi address is not changed.
//x[i]=(mpfi_ptr)rs_export_ibfr_mpfi(ident_elt);
ident_node=rs_export_list_vect_ibfr_nextnode
(ident_node);
}
return nb_elts;
}
template<class OutputIterator>
static OutputIterator insert_roots(OutputIterator x){
CGALRS_PTR(ident_sols_eqs);
CGALRS_PTR(ident_node);
CGALRS_PTR(ident_vect);
CGALRS_PTR(ident_elt);
int nb_elts;
ident_sols_eqs=rs_get_default_sols_eqs();
nb_elts=rs_export_list_vect_ibfr_nb(ident_sols_eqs);
ident_node=rs_export_list_vect_ibfr_firstnode(ident_sols_eqs);
for(int i=0;i<nb_elts;++i){
ident_vect=rs_export_list_vect_ibfr_monnode
(ident_node);
CGAL_assertion_msg(rs_export_dim_vect_ibfr
(ident_vect)==1,
"vector dimension must be 1");
ident_elt=rs_export_elt_vect_ibfr(ident_vect,0);
mpfi_ptr root_pointer=
(mpfi_ptr)rs_export_ibfr_mpfi(ident_elt);
mp_prec_t root_prec=mpfi_get_prec(root_pointer);
// Construct Gmpfr's with pointers to endpoints.
Gmpfr left(&(root_pointer->left),root_prec);
Gmpfr right(&(root_pointer->right),root_prec);
CGAL_assertion(left<=right);
// Copy them, to have the data out of RS memory.
*x++=Gmpfi(left,right,root_prec+1);
ident_node=rs_export_list_vect_ibfr_nextnode
(ident_node);
}
return x;
}
}; // struct RS2_calls
} // namespace RS2
} // namespace CGAL
#endif // CGAL_RS_RS2_CALLS_H

View File

@ -0,0 +1,114 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_RS2_ISOLATOR_1_H
#define CGAL_RS_RS2_ISOLATOR_1_H
#include "rs2_calls.h"
#include "polynomial_converter_1.h"
#include <CGAL/Gmpfi.h>
#include <vector>
namespace CGAL{
namespace RS2{
template <class Polynomial_,class Bound_>
class RS2_isolator_1{
public:
typedef Polynomial_ Polynomial;
typedef Bound_ Bound;
private:
typedef Gmpfi Interval;
public:
RS2_isolator_1(const Polynomial&);
Polynomial polynomial()const;
int number_of_real_roots()const;
bool is_exact_root(int i)const;
Bound left_bound(int i)const;
Bound right_bound(int i)const;
private:
Polynomial _polynomial;
std::vector<Interval> _real_roots;
};
template <class Polynomial_,class Bound_>
RS2_isolator_1<Polynomial_,Bound_>::
RS2_isolator_1(const Polynomial_ &){
CGAL_error_msg("not implemented for these polynomial/bound types");
}
// Isolator of Gmpz polynomials and Gmpfr bounds (the best case for RS).
template <>
RS2_isolator_1<CGAL::Polynomial<CGAL::Gmpz>,Gmpfr>::
RS2_isolator_1(const CGAL::Polynomial<CGAL::Gmpz> &p):_polynomial(p){
//mpz_t *coeffs=(mpz_t*)malloc((degree+1)*sizeof(mpz_t));
//mpfi_ptr *intervals_mpfi=(mpfi_ptr*)malloc(degree*sizeof(mpfi_ptr));
//for(unsigned i=0;i<=degree;++i)
// coeffs[i][0]=*(p[i].mpz());
RS2::RS2_calls::init_solver();
RS2::RS2_calls::create_rs_upoly(p,rs_get_default_up());
//free(coeffs);
set_rs_precisol(0);
set_rs_verbose(0);
rs_run_algo((char*)"UISOLE");
RS2::RS2_calls::insert_roots(std::back_inserter(_real_roots));
}
// Isolator of Gmpq polynomials and Gmpfr bounds. The polynomial must be
// converted to a Gmpz one with the same roots.
template <>
RS2_isolator_1<CGAL::Polynomial<CGAL::Gmpq>,Gmpfr>::
RS2_isolator_1(const CGAL::Polynomial<CGAL::Gmpq> &p):_polynomial(p){
RS2::RS2_calls::init_solver();
RS2::RS2_calls::create_rs_upoly(
CGAL::RS_AK1::Polynomial_converter_1<
CGAL::Polynomial<Gmpq>,
CGAL::Polynomial<Gmpz> >()(p),
rs_get_default_up());
set_rs_precisol(0);
set_rs_verbose(0);
rs_run_algo((char*)"UISOLE");
RS2::RS2_calls::insert_roots(std::back_inserter(_real_roots));
}
template <class Polynomial_,class Bound_>
Polynomial_
RS2_isolator_1<Polynomial_,Bound_>::polynomial()const{
return _polynomial;
}
template <class Polynomial_,class Bound_>
int
RS2_isolator_1<Polynomial_,Bound_>::number_of_real_roots()const{
return _real_roots.size();
}
template <class Polynomial_,class Bound_>
bool
RS2_isolator_1<Polynomial_,Bound_>::is_exact_root(int i)const{
return _real_roots[i].inf()==_real_roots[i].sup();
}
template <class Polynomial_,class Bound_>
Bound_
RS2_isolator_1<Polynomial_,Bound_>::left_bound(int i)const{
return _real_roots[i].inf();
}
template <class Polynomial_,class Bound_>
Bound_
RS2_isolator_1<Polynomial_,Bound_>::right_bound(int i)const{
return _real_roots[i].sup();
}
} // namespace RS2
} // namespace CGAL
#endif // CGAL_RS_RS2_ISOLATOR_1_H

View File

@ -0,0 +1,157 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_RS3_K_REFINER_1_H
#define CGAL_RS_RS3_K_REFINER_1_H
#include <CGAL/Polynomial_traits_d.h>
#include "polynomial_converter_1.h"
#include "rs2_calls.h"
#include <rs3_fncts.h>
#include "Gmpfr_make_unique.h"
// If we want assertions, we need to evaluate.
#ifndef CGAL_NO_PRECONDITIONS
#include "signat_1.h"
#endif
namespace CGAL{
namespace RS3{
template <class Polynomial_,class Bound_>
struct RS3_k_refiner_1{
void operator()(const Polynomial_&,Bound_&,Bound_&,int);
}; // class RS3_k_refiner_1
template <class Polynomial_,class Bound_>
void
RS3_k_refiner_1<Polynomial_,Bound_>::
operator()(const Polynomial_&,Bound_&,Bound_&,int){
CGAL_error_msg("RS3 k-refiner not implemented for these types");
return;
}
template<>
void
RS3_k_refiner_1<Polynomial<Gmpz>,Gmpfr>::
operator()
(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec){
typedef Polynomial<Gmpz> Polynomial;
typedef Polynomial_traits_d<Polynomial> Ptraits;
typedef Ptraits::Degree Degree;
CGAL_precondition(left<=right);
#ifndef CGAL_NO_PRECONDITIONS
typedef Ptraits::Make_square_free Sfpart;
typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
Signat;
Polynomial sfpp=Sfpart()(pol);
Signat signof(sfpp);
CGAL::Sign sl=signof(left);
CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
#endif
//std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
int deg=Degree()(pol);
mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
__mpfi_struct interval;
CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
CGAL_RS_GMPFR_MAKE_UNIQUE(right,temp_right);
interval.left=*(left.fr());
interval.right=*(right.fr());
for(int i=0;i<=deg;++i)
coefficients[i][0]=*(pol[i].mpz());
RS2::RS2_calls::init_solver();
rs3_refine_u_root(&interval,
coefficients,
deg,
prec+CGAL::max(left.get_precision(),
right.get_precision()),
1,
1);
free(coefficients);
mpfr_clear(left.fr());
mpfr_clear(right.fr());
mpfr_custom_init_set(left.fr(),
mpfr_custom_get_kind(&interval.left),
mpfr_custom_get_exp(&interval.left),
mpfr_get_prec(&interval.left),
mpfr_custom_get_mantissa(&interval.left));
mpfr_custom_init_set(right.fr(),
mpfr_custom_get_kind(&interval.right),
mpfr_custom_get_exp(&interval.right),
mpfr_get_prec(&interval.right),
mpfr_custom_get_mantissa(&interval.right));
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
template<>
void
RS3_k_refiner_1<Polynomial<Gmpq>,Gmpfr>::
operator()
(const Polynomial<Gmpq> &qpol,Gmpfr &left,Gmpfr &right,int prec){
typedef Polynomial<Gmpz> ZPolynomial;
typedef Polynomial_traits_d<ZPolynomial> ZPtraits;
typedef ZPtraits::Degree ZDegree;
CGAL_precondition(left<=right);
#ifndef CGAL_NO_PRECONDITIONS
typedef ZPtraits::Make_square_free ZSfpart;
typedef CGAL::RS_AK1::Signat_1<ZPolynomial,Gmpfr>
Signat;
#endif
//std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
Polynomial<Gmpz> zpol=CGAL::RS_AK1::Polynomial_converter_1<
CGAL::Polynomial<Gmpq>,
CGAL::Polynomial<Gmpz> >()(qpol);
#ifndef CGAL_NO_PRECONDITIONS
ZPolynomial zsfpp=ZSfpart()(zpol);
Signat signof(zsfpp);
CGAL::Sign sl=signof(left);
CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
#endif
int deg=ZDegree()(zpol);
mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
__mpfi_struct interval;
CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
CGAL_RS_GMPFR_MAKE_UNIQUE(right,temp_right);
interval.left=*(left.fr());
interval.right=*(right.fr());
for(int i=0;i<=deg;++i)
coefficients[i][0]=*(zpol[i].mpz());
RS2::RS2_calls::init_solver();
rs3_refine_u_root(&interval,
coefficients,
deg,
prec+CGAL::max(left.get_precision(),
right.get_precision()),
1,
1);
free(coefficients);
mpfr_clear(left.fr());
mpfr_clear(right.fr());
mpfr_custom_init_set(left.fr(),
mpfr_custom_get_kind(&interval.left),
mpfr_custom_get_exp(&interval.left),
mpfr_get_prec(&interval.left),
mpfr_custom_get_mantissa(&interval.left));
mpfr_custom_init_set(right.fr(),
mpfr_custom_get_kind(&interval.right),
mpfr_custom_get_exp(&interval.right),
mpfr_get_prec(&interval.right),
mpfr_custom_get_mantissa(&interval.right));
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
} // namespace RS3
} // namespace CGAL
#endif // CGAL_RS_RS3_K_REFINER_1_H

View File

@ -0,0 +1,169 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_RS3_REFINER_1_H
#define CGAL_RS_RS3_REFINER_1_H
#include <CGAL/Polynomial_traits_d.h>
#include "rs2_calls.h"
#include <rs3_fncts.h>
#include "Gmpfr_make_unique.h"
// If we want assertions, we need to evaluate.
#ifndef CGAL_NO_PRECONDITIONS
#include "signat_1.h"
#endif
namespace CGAL{
namespace RS3{
template <class Polynomial_,class Bound_>
struct RS3_refiner_1{
void operator()(const Polynomial_&,Bound_&,Bound_&,int,int=0);
}; // class RS3_refiner_1
template <class Polynomial_,class Bound_>
void
RS3_refiner_1<Polynomial_,Bound_>::
operator()(const Polynomial_&,Bound_&,Bound_&,int,int){
CGAL_error_msg("RS3 refiner not implemented for these types");
return;
}
template<>
void
RS3_refiner_1<Polynomial<Gmpz>,Gmpfr>::
operator()
(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec,int k){
typedef Polynomial<Gmpz> Polynomial;
typedef Polynomial_traits_d<Polynomial> Ptraits;
typedef Ptraits::Degree Degree;
CGAL_precondition(left<=right);
#ifndef CGAL_NO_PRECONDITIONS
typedef Ptraits::Make_square_free Sfpart;
typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
Signat;
Polynomial sfpp=Sfpart()(pol);
Signat signof(sfpp);
CGAL::Sign sl=signof(left);
CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
#endif
//std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
int deg=Degree()(pol);
mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
__mpfi_struct interval;
// Make sure the endpoints do not share references.
CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
CGAL_RS_GMPFR_MAKE_UNIQUE(right,temp_right);
// Construct the mpfi_t interval which will be refined.
interval.left=*(left.fr());
interval.right=*(right.fr());
// Construct the polynomial which will be refined (copy pointers).
for(int i=0;i<=deg;++i)
coefficients[i][0]=*(pol[i].mpz());
// Call RS.
RS2::RS2_calls::init_solver();
rs3_refine_u_root(&interval,
coefficients,
deg,
prec+CGAL::max(left.get_precision(),
right.get_precision()),
k,
k);
// Clear variables.
free(coefficients);
mpfr_clear(left.fr());
mpfr_clear(right.fr());
// Copy results back to the Gmpfr endpoints.
mpfr_custom_init_set(left.fr(),
mpfr_custom_get_kind(&interval.left),
mpfr_custom_get_exp(&interval.left),
mpfr_get_prec(&interval.left),
mpfr_custom_get_mantissa(&interval.left));
mpfr_custom_init_set(right.fr(),
mpfr_custom_get_kind(&interval.right),
mpfr_custom_get_exp(&interval.right),
mpfr_get_prec(&interval.right),
mpfr_custom_get_mantissa(&interval.right));
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
template<>
void
RS3_refiner_1<Polynomial<Gmpq>,Gmpfr>::
operator()
(const Polynomial<Gmpq> &qpol,Gmpfr &left,Gmpfr &right,int prec,int k){
typedef Polynomial<Gmpz> ZPolynomial;
typedef Polynomial_traits_d<ZPolynomial> ZPtraits;
typedef ZPtraits::Degree ZDegree;
CGAL_precondition(left<=right);
#ifndef CGAL_NO_PRECONDITIONS
typedef ZPtraits::Make_square_free ZSfpart;
typedef CGAL::RS_AK1::Signat_1<ZPolynomial,Gmpfr>
Signat;
#endif
//std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
// Construct a Gmpz polynomial from the original Gmpq polynomial.
Polynomial<Gmpz> zpol=CGAL::RS_AK1::Polynomial_converter_1<
CGAL::Polynomial<Gmpq>,
CGAL::Polynomial<Gmpz> >()(qpol);
#ifndef CGAL_NO_PRECONDITIONS
ZPolynomial zsfpp=ZSfpart()(zpol);
Signat signof(zsfpp);
CGAL::Sign sl=signof(left);
CGAL_precondition(sl!=signof(right)||(left==right&&sl==ZERO));
#endif
int deg=ZDegree()(zpol);
mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
__mpfi_struct interval;
// Make sure the endpoints do not share references.
CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
CGAL_RS_GMPFR_MAKE_UNIQUE(right,temp_right);
// Construct the mpfi_t interval which will be refined.
interval.left=*(left.fr());
interval.right=*(right.fr());
// Construct the polynomial which will be refined (copy pointers).
for(int i=0;i<=deg;++i)
coefficients[i][0]=*(zpol[i].mpz());
// Call RS.
RS2::RS2_calls::init_solver();
rs3_refine_u_root(&interval,
coefficients,
deg,
prec+CGAL::max(left.get_precision(),
right.get_precision()),
k,
k);
// Clear variables.
free(coefficients);
mpfr_clear(left.fr());
mpfr_clear(right.fr());
// Copy results back to the Gmpfr endpoints.
mpfr_custom_init_set(left.fr(),
mpfr_custom_get_kind(&interval.left),
mpfr_custom_get_exp(&interval.left),
mpfr_get_prec(&interval.left),
mpfr_custom_get_mantissa(&interval.left));
mpfr_custom_init_set(right.fr(),
mpfr_custom_get_kind(&interval.right),
mpfr_custom_get_exp(&interval.right),
mpfr_get_prec(&interval.right),
mpfr_custom_get_mantissa(&interval.right));
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
} // namespace RS3
} // namespace CGAL
#endif // CGAL_RS_RS3_REFINER_1_H

View File

@ -0,0 +1,105 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_SIGNAT_1_H
#define CGAL_RS_SIGNAT_1_H
#include <CGAL/Gmpfi.h>
#include <CGAL/Polynomial_traits_d.h>
#include "exact_signat_1.h"
#include <gmp.h>
namespace CGAL{
namespace RS_AK1{
template <class Polynomial_,class Bound_>
struct Signat_1{
typedef Polynomial_ Polynomial;
typedef Bound_ Bound;
typedef CGAL::Polynomial_traits_d<Polynomial> PT;
typedef typename PT::Degree Degree;
Polynomial pol;
Signat_1(const Polynomial &p):pol(p){};
CGAL::Sign operator()(const Bound&)const;
}; // struct Signat_1
template <class Polynomial_,class Bound_>
inline CGAL::Sign
Signat_1<Polynomial_,Bound_>::operator()(const Bound_ &x)const{
typedef Bound_ Bound;
typedef Real_embeddable_traits<Bound> REtraits;
typedef typename REtraits::Sgn BSign;
//typedef Algebraic_structure_traits<Bound> AStraits;
// This generic signat works only when Bound_ is an exact type. For
// non-exact types, an implementation must be provided.
//static_assert(std::is_same<AStraits::Is_exact,Tag_true>::value);
int d=Degree()(pol);
Bound h(pol[d]);
for(int i=1;i<=d;++i)
h=h*x+pol[d-i];
return BSign()(h);
}
template <>
inline CGAL::Sign
Signat_1<Polynomial<Gmpz>,Gmpfr>::operator()(const Gmpfr &x)const{
// In 32-bit systems, using Gmpfr arithmetic to perform exact
// evaluations can overflow. For that reason, we only use Gmpfr
// arithmetic in 64-bit systems.
#if (GMP_LIMB_BITS==64)
typedef ExactSignat_1<Polynomial,Gmpfr> Exact_sign;
#else
typedef Signat_1<Polynomial,Gmpq> Exact_sign;
#endif
// This seems to work faster for small polynomials:
// return Exact_sign(pol)(x);
int d=Degree()(pol);
if(d==0)
return pol[0].sign();
Gmpfi h(pol[d],x.get_precision()+2*d);
Uncertain<CGAL::Sign> indet=Uncertain<CGAL::Sign>::indeterminate();
if(h.sign().is_same(indet))
return Exact_sign(pol)(x);
for(int i=1;i<=d;++i){
h*=x;
h+=pol[d-i];
if(h.sign().is_same(indet))
return Exact_sign(pol)(x);
}
CGAL_assertion(!h.sign().is_same(indet));
return h.sign();
}
// This is the same code as above.
template <>
inline CGAL::Sign
Signat_1<Polynomial<Gmpq>,Gmpfr>::operator()(const Gmpfr &x)const{
typedef Signat_1<Polynomial,Gmpq> Exact_sign;
int d=Degree()(pol);
if(d==0)
return pol[0].sign();
Gmpfi h(pol[d],x.get_precision()+2*d);
Uncertain<CGAL::Sign> indet=Uncertain<CGAL::Sign>::indeterminate();
if(h.sign().is_same(indet))
return Exact_sign(pol)(x);
for(int i=1;i<=d;++i){
h*=x;
h+=pol[d-i];
if(h.sign().is_same(indet))
return Exact_sign(pol)(x);
}
CGAL_assertion(!h.sign().is_same(indet));
return h.sign();
}
} // namespace RS_AK1
} // namespace CGAL
#endif // CGAL_RS_SIGNAT_1_H

View File

@ -0,0 +1,142 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-2.1-only
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#include <CGAL/config.h>
#if defined(CGAL_USE_GMP) && defined(CGAL_USE_MPFI) && defined(CGAL_USE_RS)
#include "include/CGAL/_test_algebraic_kernel_1.h"
#ifdef CGAL_RS_TEST_LOG_TIME
#include <ctime>
#endif
#include <CGAL/Algebraic_kernel_rs_gmpq_d_1.h>
template <class AK_>
int test_ak(){
typedef AK_ AK;
typedef typename AK::Polynomial_1 Polynomial_1;
//typedef typename AK::Coefficient Coefficient;
typedef typename AK::Bound Bound;
typedef typename AK::Algebraic_real_1 Algebraic_real_1;
typedef typename AK::Multiplicity_type Multiplicity_type;
#ifdef CGAL_RS_TEST_LOG_TIME
clock_t start=clock();
#endif
AK ak; // an algebraic kernel object
CGAL::test_algebraic_kernel_1<AK>(ak); // we run standard tests first
typename AK::Solve_1 solve_1 = ak.solve_1_object();
Polynomial_1 x = CGAL::shift(Polynomial_1(1),1);
int returnvalue=0;
// 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));
if(roots.size()!=2){
returnvalue-=1;
std::cerr<<"error 1: the number of roots of x^2-2 must be 2"<<
std::endl;
}
if(-1.42>=roots[0] || -1.41<=roots[0] ||
1.41>=roots[1] || 1.42<=roots[1]){
returnvalue-=2;
std::cerr<<"error 2: the roots of x^2-2 are wrong"<<std::endl;
}
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));
if(roots.size()!=2){
returnvalue-=4;
std::cerr<<"error 3: the number of roots of (x^2-2)*(x^2-3)"<<
" between 0 and 10 must be 2"<<std::endl;
}
if(1.41>=roots[0] || 1.42<=roots[0] ||
1.73>=roots[1] || 1.74<=roots[1]){
returnvalue-=8;
std::cerr<<"error 4: the roots of (x^2-2)*(x^2-3)"<<
" between 0 and 10 are wrong"<<std::endl;
}
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));
if(mroots.size()!=2){
returnvalue-=16;
std::cerr<<"error 5: the number of roots of x^2-2 must be 2"<<
std::endl;
}
if(-1.42>=mroots[0].first || -1.41<=mroots[0].first ||
1.41>=mroots[1].first || 1.42<=mroots[1].first){
returnvalue-=32;
std::cerr<<"error 6: the roots of x^2-2 are wrong"<<std::endl;
}
if(mroots[0].second!=1 && mroots[1].second!=1){
returnvalue-=64;
std::cerr<<"error 7: the multiplicities of the"<<
" roots of x^2-2 are wrong"<<std::endl;
}
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));
if(mroots.size()!=2){
returnvalue-=128;
std::cerr<<"error 8: the number of roots of (x^2-2)*(x^2-3)"<<
" between 0 and 10 must be 2"<<std::endl;
}
if(1.41>=mroots[0].first || 1.42<=mroots[0].first ||
1.73>=mroots[1].first || 1.74<=mroots[1].first){
returnvalue-=256;
std::cerr<<"error 9: the roots of (x^2-2)*(x^2-3) are wrong"<<
std::endl;
}
if(mroots[0].second!=1 && mroots[1].second!=1){
returnvalue-=512;
std::cerr<<"error 10: the multiplicities of the roots of"<<
" (x^2-2)*(x^2-3) are wrong"<<std::endl;
}
typename AK::Number_of_solutions_1 nos_1 = ak.number_of_solutions_1_object();
if(nos_1(x*x*x-2)!=1){
returnvalue-=1024;
std::cerr<<"error 11: x^3-2 must have only one root"<<std::endl;
}
#ifdef CGAL_RS_TEST_LOG_TIME
std::cerr<<"*** test time: "<<(double)(clock()-start)/CLOCKS_PER_SEC<<
" seconds"<<std::endl;
#endif
return returnvalue;
}
int main(){
typedef CGAL::Algebraic_kernel_rs_gmpq_d_1 AK_rational;
// test and return the result
long result=0;
std::cerr<<"*** testing rational RS AK_1:";
result+=test_ak<AK_rational>();
std::cerr<<"*** result of the tests (should be 0): "<<result<<std::endl;
return result;
}
#else
int main(){
return 0;
}
#endif

View File

@ -0,0 +1,189 @@
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-2.1-only
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#include <CGAL/config.h>
#if defined(CGAL_USE_GMP) && defined(CGAL_USE_MPFI) && defined(CGAL_USE_RS)
#include "include/CGAL/_test_algebraic_kernel_1.h"
#ifdef CGAL_RS_TEST_LOG_TIME
#include <ctime>
#endif
// default RS_AK_1
#include <CGAL/Algebraic_kernel_rs_gmpz_d_1.h>
// different isolators
#include <CGAL/RS/rs2_isolator_1.h>
#ifdef CGAL_USE_RS3
#include <CGAL/RS/rs23_k_isolator_1.h>
#endif
// different refiners
#include <CGAL/RS/bisection_refiner_1.h>
#ifdef CGAL_USE_RS3
#include <CGAL/RS/rs3_refiner_1.h>
#include <CGAL/RS/rs3_k_refiner_1.h>
#endif
template <class AK_>
int test_ak(){
typedef AK_ AK;
typedef typename AK::Polynomial_1 Polynomial_1;
//typedef typename AK::Coefficient Coefficient;
typedef typename AK::Bound Bound;
typedef typename AK::Algebraic_real_1 Algebraic_real_1;
typedef typename AK::Multiplicity_type Multiplicity_type;
#ifdef CGAL_RS_TEST_LOG_TIME
clock_t start=clock();
#endif
AK ak; // an algebraic kernel object
CGAL::test_algebraic_kernel_1<AK>(ak); // we run standard tests first
typename AK::Solve_1 solve_1 = ak.solve_1_object();
Polynomial_1 x = CGAL::shift(Polynomial_1(1),1);
int returnvalue=0;
// 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));
if(roots.size()!=2){
returnvalue-=1;
std::cerr<<"error 1: the number of roots of x^2-2 must be 2"<<
std::endl;
}
if(-1.42>=roots[0] || -1.41<=roots[0] ||
1.41>=roots[1] || 1.42<=roots[1]){
returnvalue-=2;
std::cerr<<"error 2: the roots of x^2-2 are wrong"<<std::endl;
}
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));
if(roots.size()!=2){
returnvalue-=4;
std::cerr<<"error 3: the number of roots of (x^2-2)*(x^2-3)"<<
" between 0 and 10 must be 2"<<std::endl;
}
if(1.41>=roots[0] || 1.42<=roots[0] ||
1.73>=roots[1] || 1.74<=roots[1]){
returnvalue-=8;
std::cerr<<"error 4: the roots of (x^2-2)*(x^2-3)"<<
" between 0 and 10 are wrong"<<std::endl;
}
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));
if(mroots.size()!=2){
returnvalue-=16;
std::cerr<<"error 5: the number of roots of x^2-2 must be 2"<<
std::endl;
}
if(-1.42>=mroots[0].first || -1.41<=mroots[0].first ||
1.41>=mroots[1].first || 1.42<=mroots[1].first){
returnvalue-=32;
std::cerr<<"error 6: the roots of x^2-2 are wrong"<<std::endl;
}
if(mroots[0].second!=1 && mroots[1].second!=1){
returnvalue-=64;
std::cerr<<"error 7: the multiplicities of the"<<
" roots of x^2-2 are wrong"<<std::endl;
}
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));
if(mroots.size()!=2){
returnvalue-=128;
std::cerr<<"error 8: the number of roots of (x^2-2)*(x^2-3)"<<
" between 0 and 10 must be 2"<<std::endl;
}
if(1.41>=mroots[0].first || 1.42<=mroots[0].first ||
1.73>=mroots[1].first || 1.74<=mroots[1].first){
returnvalue-=256;
std::cerr<<"error 9: the roots of (x^2-2)*(x^2-3) are wrong"<<
std::endl;
}
if(mroots[0].second!=1 && mroots[1].second!=1){
returnvalue-=512;
std::cerr<<"error 10: the multiplicities of the roots of"<<
" (x^2-2)*(x^2-3) are wrong"<<std::endl;
}
typename AK::Number_of_solutions_1 nos_1 = ak.number_of_solutions_1_object();
if(nos_1(x*x*x-2)!=1){
returnvalue-=1024;
std::cerr<<"error 11: x^3-2 must have only one root"<<std::endl;
}
#ifdef CGAL_RS_TEST_LOG_TIME
std::cerr<<"*** test time: "<<(double)(clock()-start)/CLOCKS_PER_SEC<<
" seconds"<<std::endl;
#endif
return returnvalue;
}
int main(){
// We'll test three different RS-based univariate AK's:
// - the default RS one,
// - one with RS2 functions only, and
// - one with both RS2 and RS3 functions.
// the default RS kernel
typedef CGAL::Algebraic_kernel_rs_gmpz_d_1 AK_default;
typedef CGAL::Polynomial<CGAL::Gmpz> Polynomial;
typedef CGAL::Gmpfr Bound;
// the RS2-only kernel
typedef CGAL::RS2::RS2_isolator_1<Polynomial,Bound> RS2_isolator;
typedef CGAL::Bisection_refiner_1<Polynomial,Bound> B_refiner;
typedef CGAL::RS_AK1::Algebraic_kernel_1<Polynomial,
Bound,
RS2_isolator,
B_refiner> AK_RS2;
#ifdef CGAL_USE_RS3
// the RS2/RS3 kernel
typedef CGAL::RS23_k_isolator_1<Polynomial,Bound> K_isolator;
typedef CGAL::RS3::RS3_k_refiner_1<Polynomial,Bound> RS3_k_refiner;
typedef CGAL::RS_AK1::Algebraic_kernel_1<Polynomial,
Bound,
K_isolator,
RS3_k_refiner> AK_RS2_RS3;
#endif // CGAL_USE_RS3
// test all and return the result
long result=0;
std::cerr<<"*** testing default RS AK_1:";
result+=test_ak<AK_default>();
std::cerr<<"*** testing RS2 AK_1:";
result+=(2048*test_ak<AK_RS2>());
#ifdef CGAL_USE_RS3
std::cerr<<"*** testing RS2/RS3 k_AK_1:";
result+=(4096*test_ak<AK_RS2_RS3>());
#endif // CGAL_USE_RS3
std::cerr<<"*** result of the tests (should be 0): "<<result<<std::endl;
return result;
}
#else
int main(){
return 0;
}
#endif

View File

@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Algebraic_kernel_d_Tests)
# CGAL and its components
@ -10,6 +10,11 @@ if(MPFI_FOUND)
include(${MPFI_USE_FILE})
endif()
find_package(RS3 QUIET)
if(RS3_FOUND)
message(STATUS "Found RS3")
include(${RS3_USE_FILE})
endif()
# include for local directory
include_directories(BEFORE include)
@ -35,6 +40,12 @@ if(NOT CGAL_DISABLE_GMP)
create_single_source_cgal_program("Curve_pair_analysis_2.cpp")
create_single_source_cgal_program("Real_embeddable_traits_extension.cpp")
if(RS_FOUND)
create_single_source_cgal_program("Algebraic_kernel_rs_gmpq_d_1.cpp")
create_single_source_cgal_program("Algebraic_kernel_rs_gmpz_d_1.cpp")
else()
message(STATUS "NOTICE: Some tests require the RS library, and will not be compiled.")
endif()
else()
message(STATUS "NOTICE: Some tests require GMP support, and will not be compiled.")
endif()

View File

@ -46,7 +46,7 @@ class Root_for_circles_2_2 {
Root_for_circles_2_2(const Root_of_2& r1, const Root_of_2& r2)
: x_(r1), y_(r2)
{
// When it is an interval this assertion doesn't compile
// When it is an interval this assertion dont compile
//CGAL_assertion((r1.is_rational() || r2.is_rational()) ||
// (r1.gamma() == r2.gamma()));
}

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Algebraic_kernel_for_circles_Tests)
find_package(CGAL REQUIRED)

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Algebraic_kernel_for_spheres_Tests)
find_package(CGAL REQUIRED)

View File

@ -35,8 +35,8 @@ of this process in 2D (where our ice-cream spoon is simply a circle).
Alpha shapes depend on a parameter \f$ \alpha\f$ after which they are named.
In the ice-cream analogy above, \f$ \alpha\f$ is the squared radius of the
carving spoon. A very small value will allow us to eat up all the
ice-cream except the chocolate points themselves. Thus, we already see
carving spoon. A very small value will allow us to eat up all of the
ice-cream except the chocolate points themselves. Thus we already see
that the \f$ \alpha\f$-shape degenerates to the point-set \f$ S\f$ for
\f$ \alpha \rightarrow 0\f$. On the other hand, a huge value of \f$ \alpha\f$
will prevent us even from moving the spoon between two points since
@ -53,7 +53,7 @@ We distinguish two versions of alpha shapes. <I>Basic alpha shapes</I>
are based on the Delaunay triangulation. <I>Weighted alpha shapes</I>
are based on its generalization, the regular triangulation
(cf. Section \ref Section_2D_Triangulations_Regular "Regular Triangulations"),
replacing the Euclidean distance by the power to weighted points.
replacing the euclidean distance by the power to weighted points.
There is a close connection between alpha shapes and the underlying
triangulations. More precisely, the \f$ \alpha\f$-complex of \f$ S\f$ is a
@ -100,7 +100,7 @@ It provides iterators to enumerate the vertices and edges that are in
the \f$ \alpha\f$-shape, and functions that allow to classify vertices,
edges and faces with respect to the \f$ \alpha\f$-shape. They can be in
the interior of a face that belongs or does not belong to the \f$ \alpha\f$-shape.
They can be singular/regular, that is, they can be on the boundary of the \f$ \alpha\f$-shape,
They can be singular/regular, that is be on the boundary of the \f$ \alpha\f$-shape,
but not incident/incident to a triangle of the \f$ \alpha\f$-complex.
Finally, it provides a function to determine the \f$ \alpha\f$-value
@ -213,7 +213,7 @@ cube will be chosen and no optimizations will be used.
It is also recommended to switch the triangulation to 1-sheeted
covering if possible. Note that a periodic triangulation in 9-sheeted
covering space is degenerate. In this case, an exact constructions
kernel needs to be used to compute the alpha shapes. Otherwise, the
kernel needs to be used to compute the alpha shapes. Otherwise the
results will suffer from round-off problems.
\cgalExample{Alpha_shapes_2/ex_periodic_alpha_shapes_2.cpp}

View File

@ -1,4 +1,4 @@
/// \defgroup PkgAlphaShapes2Ref Reference Manual
/// \defgroup PkgAlphaShapes2Ref 2D Alpha Shapes Reference
/// \defgroup PkgAlphaShapes2Concepts Concepts
/// \ingroup PkgAlphaShapes2Ref

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Alpha_shapes_2_Examples)
find_package(CGAL REQUIRED)

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Alpha_shapes_2_Tests)
find_package(CGAL REQUIRED)

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Alpha_shapes_3_Demo)
# Find includes in corresponding build directories
@ -13,7 +13,7 @@ find_package(Qt6 QUIET COMPONENTS Widgets OpenGL)
if(CGAL_Qt6_FOUND AND Qt6_FOUND)
add_compile_definitions(QT_NO_KEYWORDS)
add_definitions(-DQT_NO_KEYWORDS)
# Instruct CMake to run moc/ui/rcc automatically when needed.
set(CMAKE_AUTOMOC ON)

View File

@ -34,8 +34,8 @@ of this process in 2D (where our ice-cream spoon is simply a circle).
Alpha shapes depend on a parameter \f$ \alpha\f$ after which they are named.
In the ice-cream analogy above, \f$ \alpha\f$ is the squared radius of the
carving spoon. A very small value will allow us to eat up all the
ice-cream except the chocolate points themselves. Thus, we already see
carving spoon. A very small value will allow us to eat up all of the
ice-cream except the chocolate points themselves. Thus we already see
that the alpha shape degenerates to the point-set \f$ S\f$ for
\f$ \alpha \rightarrow 0\f$. On the other hand, a huge value of \f$ \alpha\f$
will prevent us even from moving the spoon between two points since
@ -53,7 +53,7 @@ We distinguish two versions of alpha shapes. <I>Basic alpha shapes</I>
are based on the Delaunay triangulation. <I>Weighted alpha shapes</I>
are based on its generalization, the regular triangulation
(cf. Section \ref Triangulation3secclassRegulartriangulation "Regular Triangulations"),
replacing the Euclidean distance by the power to weighted points.
replacing the euclidean distance by the power to weighted points.
Let us consider the basic case with a Delaunay triangulation.
We first define the alpha complex of the set of points \f$ S\f$.
@ -79,7 +79,7 @@ of the alpha complex where singular faces are removed
(See \cgalFigureRef{figgenregex} for an example).
\cgalFigureBegin{figgenregex,gen-reg-ex.png}
Comparison of general and regularized alpha-shape. <B>Left:</B> Some points are taken on the surface of a torus, three points being taken relatively far from the surface of the torus; <B>Middle:</B> The general alpha-shape (for a large enough alpha value) contains the singular triangle facet of the three isolated points; <B>Right:</B> The regularized version (for the same value of alpha) does not contain any singular facet.
Comparison of general and regularized alpha-shape. <B>Left:</B> Some points are taken on the surface of a torus, three points being taken relatively far from the surface of the torus; <B>Middle:</B> The general alpha-shape (for a large enough alpha value) contains the singular triangle facet of the three isolated points; <B>Right:</B> The regularized version (for the same value of alpha) does not contains any singular facet.
\cgalFigureEnd
The alpha shapes of a set of points
@ -112,11 +112,11 @@ and radii \f$ r_1, r_2 \f$ are said to be orthogonal iff
iff \f$ C_1C_2 ^2 < r_1^2 + r_2^2\f$.
For a given value of \f$ \alpha\f$,
the weighted alpha complex is formed with the simplices of the
regular triangulation such that there is a sphere
orthogonal to the weighted points associated
regular triangulation triangulation
such that there is a sphere orthogonal to the weighted points associated
with the vertices of the simplex and suborthogonal to all the other
input weighted points. Once again the alpha shape is then defined as
the domain covered by the alpha complex and comes in general and
the domain covered by a the alpha complex and comes in general and
regularized versions.
\section Alpha_Shape_3Functionality Functionality
@ -273,7 +273,7 @@ alpha-shape, using the `Alpha_shape_3<Dt,ExactAlphaComparisonTag>` requires a ke
with exact predicates and exact constructions (or setting `ExactAlphaComparisonTag` to `Tag_true`)
while using a kernel with exact predicates is sufficient for the class `Fixed_alpha_shape_3<Dt>`.
This makes the class `Fixed_alpha_shape_3<Dt>` even more efficient in this setting.
In addition, note that the `Fixed` version is the only one of the
In addition, note that the `Fixed` version is the only of the
two that supports incremental insertion and removal of points.
We give the time spent while computing the alpha shape of a protein (considered
@ -345,7 +345,7 @@ cube will be chosen and no optimizations will be used.
It is also recommended to switch the triangulation to 1-sheeted
covering if possible. Note that a periodic triangulation in 27-sheeted
covering space is degenerate. In this case, an exact constructions
kernel needs to be used to compute the alpha shapes. Otherwise, the
kernel needs to be used to compute the alpha shapes. Otherwise the
results will suffer from round-off problems.
\cgalExample{Alpha_shapes_3/ex_periodic_alpha_shapes_3.cpp}

View File

@ -1,4 +1,4 @@
/// \defgroup PkgAlphaShapes3Ref Reference Manual
/// \defgroup PkgAlphaShapes3Ref 3D Alpha Shapes Reference
/// \defgroup PkgAlphaShapes3Concepts Concepts
/// \ingroup PkgAlphaShapes3Ref
/*!

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Alpha_shapes_3_Examples)
find_package(CGAL REQUIRED)

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Alpha_shapes_3_Tests)
find_package(CGAL REQUIRED)

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Alpha_wrap_3_Benchmark)
find_package(CGAL REQUIRED)

View File

@ -134,13 +134,13 @@ def main(argv):
avg_diff_str = str(format(abs(avg_diff_to_goal), '.2f'))
if key == "Mean_Min_Angle_(degree)" or key == "Mean_Max_Angle_(degree)":
if avg_diff_to_goal < 0 :
title += "\nIn average we lose "
title += "\nIn average we loose "
else :
title += "\nIn average we gain "
title += avg_diff_str + "° toward 60°"
elif key == "Mean_Radius_Ratio" or key == "Mean_Edge_Ratio" or key == "Mean_Aspect_Ratio" :
if avg_diff_to_goal < 0 :
title += "\nIn average we lose "
title += "\nIn average we loose "
else :
title += "\nIn average we gain "
title += avg_diff_str + " of ratio toward 1"

View File

@ -75,7 +75,7 @@ double mean_min_angle(const Mesh& mesh)
const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh);
std::array<FT, 3> angles = triangle_angles(tr);
FT min_angle = (std::min)({angles[0], angles[1], angles[2]});
FT min_angle = std::min({angles[0], angles[1], angles[2]});
min_angle = min_angle * (180.0 / CGAL_PI);
mean_min_angle += min_angle;
@ -93,7 +93,7 @@ double mean_max_angle(const Mesh& mesh)
const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh);
std::array<FT, 3> angles = triangle_angles(tr);
FT max_angle = (std::max)({angles[0], angles[1], angles[2]});
FT max_angle = std::max({angles[0], angles[1], angles[2]});
max_angle = max_angle * (180.0 / CGAL_PI);
mean_max_angle += max_angle;
@ -151,8 +151,8 @@ double mean_edge_ratio(const Mesh& mesh,
FT a = std::sqrt(CGAL::squared_distance(tr[0], tr[1]));
FT b = std::sqrt(CGAL::squared_distance(tr[1], tr[2]));
FT c = std::sqrt(CGAL::squared_distance(tr[2], tr[0]));
FT min_edge = (std::min)({a, b, c});
FT max_edge = (std::max)({a, b, c});
FT min_edge = std::min({a, b, c});
FT max_edge = std::max({a, b, c});
FT edge_ratio = max_edge / min_edge;
mean_edge_ratio += edge_ratio;
@ -181,7 +181,7 @@ double mean_aspect_ratio(const Mesh& mesh,
FT c = std::sqrt(CGAL::squared_distance(tr[2], tr[0]));
FT s = 0.5 * (a + b + c);
FT inscribed_radius = std::sqrt((s * (s - a) * (s - b) * (s - c)) / s);
FT max_edge = (std::max)({a, b, c});
FT max_edge = std::max({a, b, c});
FT aspect_ratio = max_edge / inscribed_radius;
aspect_ratio /= (2. * std::sqrt(3.)); // normalized
mean_aspect_ratio += aspect_ratio;

View File

@ -1,4 +1,4 @@
/// \defgroup PkgAlphaWrap3Ref Reference Manual
/// \defgroup PkgAlphaWrap3Ref 3D Alpha Wrapping
/// \defgroup AW3_free_functions_grp Free Functions
/// Functions to create a wrap from point clouds, triangle soups, and triangle meshes.

View File

@ -264,7 +264,7 @@ Secondly, the farther the isosurface is from the input, the more new points are
through the first criterion (i.e., through intersection with dual Voronoi edge, see Section \ref aw3_algorithm);
thus, the quality of the output improves in terms of angles of the triangle elements.
Finally, and depending on the value of the alpha parameter, a large offset can also offer defeaturing capabilities.
However, using a small offset parameter will tend to better preserve sharp features as projection
However using a small offset parameter will tend to better preserve sharp features as projection
Steiner points tend to project onto convex sharp features.
\cgalFigureAnchor{6}
@ -278,7 +278,7 @@ Impact of the offset parameter on the output.
The offset parameter is decreasing from left to right, to respectively 1/50, 1/200 and 1/1000 of the longest diagonal of the input bounding box.
The alpha parameter is equal to 1/50 of the longest diagonal of the input bounding box for all level of details.
A larger offset will produce an output less complex with better triangle quality.
However, the sharp features (red edges) are well-preserved when the offset parameter is small.
However the sharp features (red edges) are well preserved when the offset parameter is small.
\cgalFigureCaptionEnd
\cgalFigureAnchor{7}
@ -346,15 +346,15 @@ The charts below plots the computation times of the wrapping algorithm on the Th
\cgalFigureCaptionBegin{9}
Execution times and output complexity for different values of alpha on the Thingi10k data set.
Alpha increases from 1/20 to 1/200 of the length of bounding box diagonal.
The x-axis represents the complexity of the output wrap mesh in number of triangle facets.
The y-axis represents the total computation time, in seconds.
The x axis represents the complexity of the output wrap mesh in number of triangle facets.
The y axis represents the total computation time, in seconds.
The color and diameter of the dots represent the number of faces in the input triangle soup,
ranging from 10 (green) to 3154000 (blue).
\cgalFigureCaptionEnd
\section aw3_examples Examples
Here is an example with an input triangle mesh, with alpha set to 1/20 of the bounding box's longest diagonal edge length,
Here is an example with an input triangle mesh, with alpha set to 1/20 of the bounding box longest diagonal edge length,
and offset set to 1/30 of alpha (i.e., 1/600 of the bounding box diagonal edge length).
\cgalExample{Alpha_wrap_3/triangle_mesh_wrap.cpp}

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Alpha_wrap_3_Examples)
find_package(CGAL REQUIRED)

View File

@ -159,9 +159,7 @@ int main(int argc, char** argv)
// edge length of regular tetrahedron with circumradius alpha
const double l = 1.6329931618554521 * alpha; // sqrt(8/3)
CGAL::tetrahedral_isotropic_remeshing(tr, l,
CGAL::parameters::remesh_boundaries(false)
.number_of_iterations(5));
CGAL::tetrahedral_isotropic_remeshing(tr, l, CGAL::parameters::remesh_boundaries(false));
std::cout << "AFTER: " << tr.number_of_vertices() << " vertices, " << tr.number_of_cells() << " cells" << std::endl;

View File

@ -348,7 +348,7 @@ public:
#ifdef CGAL_AW3_TIMER
t.stop();
std::cout << "Manifoldness postprocessing took: " << t.time() << " s." << std::endl;
std::cout << "Manifoldness post-processing took: " << t.time() << " s." << std::endl;
t.reset();
t.start();
#endif

View File

@ -97,7 +97,7 @@ template <typename Cb>
class Cell_base_with_timestamp
: public Cb
{
std::size_t time_stamp_ = std::size_t(-2);
std::size_t time_stamp_;
public:
using Has_timestamp = CGAL::Tag_true;
@ -112,7 +112,7 @@ public:
public:
template <typename... Args>
Cell_base_with_timestamp(const Args&... args)
: Cb(args...)
: Cb(args...), time_stamp_(-1)
{ }
Cell_base_with_timestamp(const Cell_base_with_timestamp& other)

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Alpha_wrap_3_Tests)
find_package(CGAL REQUIRED)

View File

@ -17,11 +17,12 @@ concept `ApolloniusSite_2`.
\cgalHeading{I/O}
The I/O operators are defined for `std::iostream`.
The I/O operators are defined for `iostream`.
The information output in the `std::iostream` is: the point of the
The information output in the `iostream` is: the point of the
Apollonius site and its weight.
\sa `CGAL::Qt_widget`
\sa `CGAL::Apollonius_graph_traits_2<K,Method_tag>`
\sa `CGAL::Apollonius_graph_filtered_traits_2<CK,CM,EK,EM,FK,FM>`
*/
@ -49,6 +50,7 @@ Apollonius_site_2(const Apollonius_site_2<K>& other);
/*!
Inserts the
Apollonius site `s` into the stream `os`.
\note Included through `CGAL/IO/Qt_widget_Apollonius_site_2.h`.
\pre The insert operator must be defined for `Point_2` and `Weight`.
\relates Apollonius_site_2
*/
@ -57,9 +59,18 @@ std::ostream& operator<<(std::ostream& os, const Apollonius_site_2<K>& s) const;
/*!
Reads an Apollonius site from the stream `is` and assigns it
to `s`.
\note Included through `CGAL/IO/Qt_widget_Apollonius_site_2.h`.
\pre The extract operator must be defined for `Point_2` and `Weight`.
\relates Apollonius_site_2
*/
std::istream& operator>>(std::istream& is, const Apollonius_site_2<K>& s);
/*!
Inserts the Apollonius site `s` into the `Qt_widget` stream `w`.
\note Included through `CGAL/IO/Qt_widget_Apollonius_site_2.h`.
\pre The insert operator must be defined for `K::Circle_2`.
\relates Apollonius_site_2
*/
Qt_widget& operator<<(Qt_widget& w, const Apollonius_site_2<K>& s) const;
} /* end namespace CGAL */

View File

@ -1,8 +1,9 @@
/// \defgroup PkgApolloniusGraph2Ref Reference Manual
/// \defgroup PkgApolloniusGraph2Ref 2D Apollonius Graphs (Delaunay Graphs of Disks) Reference
/// \defgroup PkgApolloniusGraph2Concepts Concepts
/// \ingroup PkgApolloniusGraph2Ref
/*!
\addtogroup PkgApolloniusGraph2Ref
\todo check generated documentation
\cgalPkgDescriptionBegin{2D Apollonius Graphs (Delaunay Graphs of Disks),PkgApolloniusGraph2}
\cgalPkgPicture{CircleVoronoi.png}
\cgalPkgSummaryBegin

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Apollonius_graph_2_Examples)
find_package(CGAL REQUIRED COMPONENTS Core)

View File

@ -77,7 +77,7 @@ private:
public:
typedef Voronoi_radius_2<K> Voronoi_radius;
typedef typename K::Bounded_side Bounded_side;
typedef typename K::Bounded_side Bounded_dide;
public:
template<class Tag>

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Apollonius_graph_2_Tests)
find_package(CGAL REQUIRED)

View File

@ -1,7 +1,7 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Arithmetic_kernel_Tests)
find_package(CGAL REQUIRED COMPONENTS Core)

View File

@ -164,7 +164,7 @@
or at least issues warnings about returning a reference to temporary
variable. Can you prompt the INRIA people to fix this?
(I personally think that we should remove this example from our test-suite,
but if the INRIA people believe that its place is there, they should at
but if the INRIA people believe that it's place is there, they should at
least properly maintain it ...)
- There is a problem with the examples that use CORE on Darwin (platform #10).
Other Darwin platforms seem fine. Do we want to investigate? (perhaps it's

View File

@ -54,7 +54,7 @@ void validate(std::any & v, const std::vector<std::string> & values,
Option_parser::my_validate<Option_parser::Strategy_id>(v, values);
}
/*! constructs */
/*! Constructor */
Option_parser::Option_parser() :
m_generic_opts("Generic options"),
m_config_opts("Configuration options"),
@ -136,7 +136,7 @@ Option_parser::Option_parser() :
m_positional_opts.add("input-file", -1);
}
/*! parses the options */
/*! Parse the options */
void Option_parser::operator()(int argc, char * argv[])
{
po::store(po::command_line_parser(argc, argv).
@ -225,20 +225,20 @@ void Option_parser::operator()(int argc, char * argv[])
}
}
/*! obtains the base file-name */
/*! Obtain the base file-name */
const std::string & Option_parser::get_file_name(unsigned int i) const
{
return m_variable_map["input-file"].as<Input_path>()[i];
}
/*! obtains the full file-name */
/*! Obtain the full file-name */
const std::string & Option_parser::get_full_name(unsigned int i) const
{ return m_full_names[i]; }
/*! obtains number of type options */
/*! Obtain number of type options */
unsigned int Option_parser::get_number_opts(Type_id &)
{ return sizeof(s_type_opts) / sizeof(char *); }
/*! obtains number of strategy options */
/*! Obtain number of strategy options */
unsigned int Option_parser::get_number_opts(Strategy_id &)
{ return sizeof(s_strategy_opts) / sizeof(char *); }

View File

@ -62,17 +62,17 @@ public:
typedef Vector_strategy_id::iterator Vector_strategy_id_iter;
public:
/*! obtains number of type options */
/*! \brief obtains number of type options */
static unsigned int get_number_opts(Type_id &);
/*! obtains number of strategy options */
/*! \brief obtains number of strategy options */
static unsigned int get_number_opts(Strategy_id &);
/*! compares the i-th type option to a given option */
/*! Compare the i-th type option to a given option */
static bool compare_opt(unsigned int i, const char * opt, Type_id &)
{ return strcmp(s_type_opts[i], opt) == 0; }
/*! compares the i-th strategy option to a given option */
/*! Compare the i-th strategy option to a given option */
static bool compare_opt(unsigned int i, const char * opt, Strategy_id &)
{ return strcmp(s_strategy_opts[i], opt) == 0; }
@ -94,19 +94,19 @@ public:
Input_file_missing_error(std::string & str) : error(str) {}
};
/*! parses the options */
/*! Parse the options */
void operator()(int argc, char * argv[]);
/*! obtains the verbosity level */
/*! Obtain the verbosity level */
unsigned int get_verbose_level() const { return m_verbose_level; }
/*! obtains the number of input files */
/*! Obtain the number of input files */
unsigned int get_number_files() const { return m_number_files; }
/*! obtains the base file-name */
/*! \brief obtains the base file-name */
const std::string & get_file_name(unsigned int i) const;
/*! obtains the full file-name */
/*! \brief obtains the full file-name */
const std::string & get_full_name(unsigned int i) const;
bool get_postscript() const { return m_postscript; }
@ -117,10 +117,10 @@ public:
const char * get_strategy_name(Strategy_code id) const
{ return s_strategy_opts[id]; }
/*! obtains the window width */
/*! Obtain the window width */
unsigned int get_width() const { return m_win_width; }
/*! obtains the window height */
/*! Obtain the window height */
unsigned int get_height() const { return m_win_height; }
template <class MyId>

View File

@ -150,7 +150,7 @@
<benchArrCgalGmpqCartesianSegment print_header="1"/>
<benchArrCgalGmpqCartesianSegment/>
<benchArrCgalGmpqLazyCartesianSegment/>
<benchArrCgalGmpqLazySimpleCartesianSegment/>
<benchArrCgalGmpqLazyCartesianSegment/>
<benchArrCgalGmpqExactSegment/>
<benchArrDoubleCartesianSegment enable="false"/>
<benchArrLazyCgalGmpqCartesianSegment/>

View File

@ -232,7 +232,7 @@ protected: // methods
// Assumes that clippingRect is valid.
std::vector< X_monotone_curve_2 > visibleParts( X_monotone_curve_2 curve );
// keep only the intersection points i.e. throw out overlapping curve segments
// keep only the intersection points ie. throw out overlapping curve segments
void filterIntersectionPoints( std::vector< CGAL::Object >& res );
protected: // members

View File

@ -1,6 +1,6 @@
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Arrangement_on_surface_2_Demo)
if(NOT POLICY CMP0070 AND POLICY CMP0053)
@ -20,7 +20,7 @@ if (CGAL_Qt6_FOUND AND Qt6_FOUND)
set(CMAKE_INCLUDE_CURRENT_DIR ON)
# Arrangement package includes
add_compile_definitions(QT_NO_KEYWORDS)
add_definitions(-DQT_NO_KEYWORDS)
option(COMPILE_UTILS_INCREMENTALLY
"Compile files in Utils directory incrementally, or compile them all as a unit. \
Incremental compilation will be better for development and consume less \

View File

@ -1,6 +1,6 @@
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.12...3.31)
cmake_minimum_required(VERSION 3.12...3.29)
project(Arrangement_on_surface_2_earth_Demo)
if(NOT POLICY CMP0070 AND POLICY CMP0053)
@ -16,7 +16,7 @@ find_package(Qt6 QUIET COMPONENTS Core Gui OpenGL OpenGLWidgets Widgets Xml)
find_package(CGAL COMPONENTS Qt6)
find_package(nlohmann_json QUIET 3.9)
set(MISSING_DEPS "")
if (NOT CGAL_FOUND OR NOT CGAL_Qt6_FOUND OR NOT Qt6_FOUND OR NOT Boost_FOUND OR NOT nlohmann_json_FOUND)
if (NOT CGAL_FOUND)
set(MISSING_DEPS "the CGAL library, ${MISSING_DEPS}")
endif()
@ -26,17 +26,20 @@ endif()
if (NOT Qt6_FOUND)
set(MISSING_DEPS "the Qt6 library, ${MISSING_DEPS}")
endif()
if (NOT Boost_FOUND)
set(MISSING_DEPS "the Boost library, ${MISSING_DEPS}")
endif()
if (NOT nlohmann_json_FOUND)
set(MISSING_DEPS "JSON for Modern C++ 3.9+ (know as nlohmann_json), ${MISSING_DEPS}")
endif()
if (MISSING_DEPS)
message(STATUS "NOTICE: This project requires ${MISSING_DEPS} and will not be compiled.")
return()
endif()
add_compile_definitions(QT_NO_VERSION_TAGGING)
add_definitions(-DQT_NO_VERSION_TAGGING)
# AOS
file(GLOB source_files_aos Aos.h Aos.cpp Aos_defs.h

View File

@ -62,7 +62,7 @@ void GUI_country_pick_handler::mouse_press_event(QMouseEvent* e) {
auto sd = sqrt(d);
auto t1 = (-b - sd) / (2 * a);
auto t2 = (-b + sd) / (2 * a);
if (t1 > 0 && t2 > 0) ti = (std::min)(t1, t2);
if (t1 > 0 && t2 > 0) ti = std::min(t1, t2);
else if (t1 > 0) ti = t1;
else ti = t2;
}

View File

@ -269,7 +269,7 @@ Kml::Nodes Kml::generate_ids_approx(Placemarks& placemarks, const double eps) {
for (const auto& node : lring->nodes) {
// check if there is a node sufficiently close to the current one
auto node_index = (std::numeric_limits<std::size_t>::max)();
auto node_index = std::numeric_limits<std::size_t>::max();
for (std::size_t i = 0; i < nodes.size(); ++i) {
const auto dist = node.distance_to(nodes[i]);
if (dist < eps) {
@ -278,7 +278,7 @@ Kml::Nodes Kml::generate_ids_approx(Placemarks& placemarks, const double eps) {
}
}
if (node_index == (std::numeric_limits<std::size_t>::max)()) {
if (node_index == std::numeric_limits<std::size_t>::max()) {
// insert new node
nodes.push_back(node);
const auto node_id = nodes.size() - 1;
@ -301,7 +301,7 @@ Kml::Nodes Kml::generate_ids_approx(Placemarks& placemarks, const double eps) {
}
// find the pair of closest nodes
double min_dist = (std::numeric_limits<double>::max)();
double min_dist = std::numeric_limits<double>::max();
std::size_t ni1 = 0;
std::size_t ni2 = 0;
std::size_t num_nodes = nodes.size();

View File

@ -140,7 +140,7 @@ void Main_widget::initializeGL() {
for (auto& [country_name, triangle_points] : country_triangles_map) {
auto country_triangles = std::make_unique<Triangles>(triangle_points);
auto color = QVector4D(rndm(), rndm(), rndm(), 1);
auto m = (std::max)(color.x(), (std::max)(color.y(), color.z()));
auto m = std::max(color.x(), std::max(color.y(), color.z()));
color /= m;
color *= m_dimming_factor;
color.setW(1);

Some files were not shown because too many files have changed in this diff Show More