1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
// Copyright (c) 2005 Stanford University (USA).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of the License,
// or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/next/Kinetic_data_structures/src/CGAL/JAMA_numeric_solver.cpp $
// $Id: JAMA_numeric_solver.cpp 67093 2012-01-13 11:22:39Z lrineau $
//
//
// Author(s) : Daniel Russel <drussel@alumni.princeton.edu>
#include <CGAL/Polynomial/internal/numeric_solvers_support.h>
#include <CGAL/Polynomial/internal/numeric_solvers.h>
#ifdef CGAL_HAVE_TNT
#include <TNT/tnt_array2d.h>
#include <TNT/tnt_array1d.h>
#include <TNT/jama_eig.h>
#endif
#include <algorithm>
#include <functional>
#include <iterator>
//#include <iomanip>
namespace CGAL { namespace POLYNOMIAL { namespace internal {
#if CGAL_HAVE_TNT
//static const double max_error_value =0.00005;
template <bool CLEAN, class NT>
static void jama_compute_roots(const NT *begin, const NT *end, NT lb,
NT ub, std::vector<NT> &roots)
{
int degree= end-begin-1;
TNT::Array2D<NT> arr(degree, degree, 0.0);
for (int i=0; i< degree; ++i) {
arr[0][i]=-begin[degree-i-1]/begin[degree];
}
for (int i=0; i+1< degree; ++i) {
arr[i+1][i]=1;
}
JAMA::Eigenvalue<NT> ev(arr);
TNT::Array1D<NT> real, imag;
ev.getImagEigenvalues(imag);
ev.getRealEigenvalues(real);
CGAL_Polynomial_assertion(imag.dim1()== real.dim1());
/*NT tol;
if (CLEAN) tol=.00005;
else tol=0;*/
for (int i=0; i< real.dim1(); ++i) {
if (root_is_good(real[i], imag[i], lb-tol, ub)) {
roots.push_back(real[i]/*polish_root(begin, end, real[i])*/);
} else {
}
}
std::sort(roots.begin(), roots.end(), std::greater<NT>());
if (CLEAN) filter_roots(begin, end, lb, roots);
}
#endif
void jama_polynomial_compute_roots(const double *begin, const double *end,
double lb, double ub,
std::vector<double> &roots)
{
std::size_t degree= end-begin-1;
switch( degree) {
case -1:
case 0:
break;
case 1:
compute_linear_roots(begin,end, lb, ub, roots);
break;
case 2:
compute_quadratic_roots(begin, end, lb, ub, roots);
break;
default:
#ifdef CGAL_HAVE_TNT
jama_compute_roots<false>(begin, end, lb, ub, roots);
#else
CGAL_error();
#endif
//jama_compute_roots<false>(begin, end, lb, ub, roots);
}
}
void jama_polynomial_compute_cleaned_roots(const double *begin, const double *end,
double lb, double ub,
std::vector<double> &roots)
{
std::size_t degree= end-begin-1;
switch( degree) {
case -1:
case 0:
break;
case 1:
compute_linear_cleaned_roots(begin,end, lb, ub, roots);
break;
case 2:
compute_quadratic_cleaned_roots(begin, end, lb, ub, roots);
break;
default:
#ifdef CGAL_HAVE_TNT
jama_compute_roots<true>(begin, end, lb, ub, roots);
#else
CGAL_error();
#endif
}
}
} } } //namespace CGAL::POLYNOMIAL::internal
|
