Commit 28bd87af authored by David Bommes's avatar David Bommes
Browse files

added Polynomials class and minor improvement of TinyAD interface

parent ecaafa49
Pipeline #18025 failed with stages
in 5 minutes and 10 seconds
......@@ -75,7 +75,7 @@ public:
auto f = ElementT::eval_f(x_ad,_c);
auto H = f.Hess;
auto H = f.Hess.template selfadjointView<Eigen::Lower>();
_triplets.clear();
for(unsigned int i=0; i<H.rows(); ++i)
......
//=============================================================================
//
// OpenFlipper
// Copyright (C) 2008 by Computer Graphics Group, RWTH Aachen
// www.openflipper.org
//
//-----------------------------------------------------------------------------
//
// License
//
// OpenFlipper is free software: 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.
//
// OpenFlipper is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenFlipper. If not, see <http://www.gnu.org/licenses/>.
//
//-----------------------------------------------------------------------------
//
// $Author: David Bommes $
// $Date: 2021-06-21 $
//
//=============================================================================
//=============================================================================
//
// CLASS Polynomials
//
//=============================================================================
#ifndef COMISO_POLYNOMIALS_HH_INCLUDED
#define COMISO_POLYNOMIALS_HH_INCLUDED
//=============================================================================
namespace COMISO
{
class Polynomials
{
public:
// robustly find **positive** roots of f(x) = _a x^2 + _b x + _c
// return {-1, 0, 1, 2} which corresponds to number of roots, or -1 which means infinitely many roots for a=b=c=0
// the potential roots are returned as _x0 and _x1 with _x0 < _x1
static int positive_roots_of_quadratic(const double _a, const double _b, const double _c, double& _x0, double& _x1, const double _eps=1e-12)
{
// get all roots
int n = roots_of_quadratic(_a, _b, _c, _x0, _x1, _eps);
// remove negative ones
if (n == 2 && _x0 < 0.0) {
std::swap(_x0, _x1);
--n;
}
// remove negative ones
if (n == 1 && _x0 < 0.0)
--n;
return n;
}
// robustly find roots of f(x) = _a x^2 + _b x + _c
// return {-1, 0, 1, 2} which corresponds to number of roots, or -1 which means infinitely many roots for a=b=c=0
// the potential roots are returned as _x0 and _x1 with _x0 < _x1
static int roots_of_quadratic(const double _a, const double _b, const double _c, double& _x0, double& _x1, const double _eps=1e-12)
{
double abs_max = std::max(std::max(std::abs(_a), std::abs(_b)), std::abs(_c));
// infinitely many roots?
// a==b==c==0
if(abs_max < _eps)
return -1;
else
{
double a = _a/abs_max;
double b = _b/abs_max;
double c = _c/abs_max;
if( std::abs(a) < _eps)
// a==0
return roots_of_linear(_b, _c, _x0, _eps);
else
{
// no constant part?
if(std::abs(c) < _eps)
{
if(std::abs(b) < _eps)
{
// a!=0, b==0, c==0
_x0 = 0.0;
return 1;
}
else
{
// a!=0, b!=0, c==0
_x0 = 0.0;
_x1 = -_b / _a;
if( _x0 > _x1)
std::swap(_x0,_x1);
return 2;
}
}
else // a!=0 && c!=0
{
// calculate discriminant
double d = _b*_b - 4.0*_a*_c;
if( d < -_eps*_eps) // d<-eps
return 0;
else
if( d<0.0) // -eps < d < 0
{
_x0 = -0.5*_b/_a;
return 1;
}
else
{
double sd = std::sqrt(d);
_x0 = 0.5*(-_b + sd)/_a;
_x1 = 0.5*(-_b - sd)/_a;
if(_x0 > _x1) std::swap(_x0,_x1);
return 2;
}
}
}
}
}
// robustly find roots of f(x) = _a x + _b
// return {-1, 0, 1} which corresponds to number of roots, or -1 which means infinitely many roots for a=b=0
// the potential roots are returned as _x0 and _x1 with _x0 < _x1
static int roots_of_linear(const double _a, const double _b, double& _x0, const double _eps=1e-12)
{
double abs_max = std::max(std::abs(_a), std::abs(_b));
// infinitely many roots?
if(abs_max < _eps)
return -1;
if(std::abs(_a/abs_max) < _eps)
return 0;
else
{
_x0 = -_b/_a;
return 1;
}
}
};
}
#endif // COMISO_POLYNOMIALS_HH_INCLUDED
//=============================================================================
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment