mmarco/SIROCCO2
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Merge pull request #3 from embray/icomplex-polynomial-eval

Browse Source 
Fix complex interval polynomial evaluation
This commit is contained in:
Miguel Marco 2020-02-11 14:25:08 +01:00 committed by GitHub
commit a42ea6032f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
10 changed files with 57 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
.gitignore vendored
View File

@ -17,5 +17,10 @@ lib/*.la
libtool
m4/
tests/.deps/
tests/.libs/
tests/*.log
tests/*.o
tests/*.trs
tests/*_test
tests/Makefile
tests/Makefile.in

17
.travis.yml Normal file
View File

@ -0,0 +1,17 @@
language: cpp
os: linux
dist: bionic
compiler:
- gcc
- clang
addons:
apt:
packages:
- libmpfr-dev
script:
- autoreconf -f -i
- ./configure
- make
- make check
after_failure:
- tail -n +1 tests/*.log 2>/dev/null

1
configure.ac
View File

@ -1,7 +1,6 @@
AC_INIT([libsirocco], [2.0], [mmarco@unizar.es])
AC_CONFIG_AUX_DIR([build-aux])
AC_CONFIG_MACRO_DIR([m4])
: ${CXXFLAGS=-O0 -g}
AM_INIT_AUTOMAKE([foreign -Wall])
AM_PROG_AR
AC_PROG_CXX

20
include/polynomial.hpp
View File

@ -73,7 +73,7 @@ template <class T>
std::ostream & operator<< (std::ostream &output, const Polynomial<T> &op) {
int allZero = 1;
for (int i=0; i<=op.degree; ++i)
for (int j=0; j<=i; ++j){
for (int j=0; j<=i; ++j) {
if (op.coef[i*(i+1)/2 + j] == 0.0) continue;
allZero = 0;
if (op.coef[i*(i+1)/2 + j] == -1.0)
@ -85,26 +85,28 @@ std::ostream & operator<< (std::ostream &output, const Polynomial<T> &op) {
if (i == 0)
if (op.coef[0] == T(1))
output << "1";
if (i != 0)
if (i-j == 0)
if (i != 0) {
if (i-j == 0) {
if (j == 1)
output << "y";
else
output << "y^" << j;
else if (i-j == 1)
} else if (i-j == 1) {
if (j == 0)
output << "x";
else if (j == 1)
output << "x*y";
else
output << "x*y^" << j;
else
} else {
if (j == 0)
output << "x^" << i;
else if (j == 1)
output << "x^" << (i-j) << "*y";
else
output << "x^" << (i-j) << "*y^" << j;
}
}
}
if (allZero) output << "0";
return output;
@ -230,6 +232,10 @@ T Polynomial<T>::operator() (const T &x, const T &y) const {
return this->evalPolClassic (x,y);
}
// Defined in lib/polynomial.cpp
template <> IComplex Polynomial<IComplex>::operator() (const IComplex &x, const IComplex &y) const;
template <> MPIComplex Polynomial<MPIComplex>::operator () (const MPIComplex &x, const MPIComplex &y) const;
template <class T>
T Polynomial<T>::diffX (const T &x, const T &y) const {
// coef[(i*(i+1))/2 + j] is coeficient of monomial of degree 'i',
@ -311,6 +317,10 @@ T Polynomial<T>::diffY (const T &x, const T &y) const {
return this->evalPolYClassic (x,y);
}
// Defined in lib/polynomial.cpp
template <> IComplex Polynomial<IComplex>::diffY (const IComplex &x, const IComplex &y) const;
template <> MPIComplex Polynomial<MPIComplex>::diffY (const MPIComplex &x, const MPIComplex &y) const;
template <class T>
T Polynomial<T>::diffXX (const T &x, const T &y) const {

4
include/sirocco.h
View File

@ -5,7 +5,9 @@
#include <mpfr.h>
// ACTIVATE ROUND MODE CONTROL
#pragma STDC FENV_ACCESS on
#ifndef __clang__
#pragma STDC FENV_ACCESS ON
#endif
double * homotopyPath (int degree, double *_coef, double _y0R, double _y0I);

4
include/sirocco.hpp
View File

@ -8,7 +8,9 @@
// ACTIVATE ROUND MODE CONTROL
#pragma STDC FENV_ACCESS on
#ifndef __clang__
#pragma STDC FENV_ACCESS ON
#endif
typedef std::complex<double> Complex;

6
lib/icomplex.cpp
View File

@ -5,7 +5,9 @@
#define MIN(a,b) (((a) < (b))? (a) : (b))
#define MAX(a,b) (((a) > (b))? (a) : (b))
#pragma STDC FENV_ACCESS on
#ifndef __clang__
#pragma STDC FENV_ACCESS ON
#endif
// CONSTRUCTORS
@ -386,4 +388,4 @@ double f (double x, double c) {
}
double g (double x, double c) {
return -c / (x*x + c*c);
}
}

6
lib/interval.cpp
View File

@ -5,7 +5,9 @@
// ACTIVATE ROUND MODE CONTROL
#pragma STDC FENV_ACCESS on
#ifndef __clang__
#pragma STDC FENV_ACCESS ON
#endif
// MACROS FOR MIN AND MAX
#define MIN(a,b) (((a) < (b))? (a) : (b))
@ -365,4 +367,4 @@ Interval operator/ (double op1, const Interval &op2) {
// RESTORE OLD ROUNDING MODE
fesetround (oldRound);
return rop;
}
}

4
lib/mp_interval.cpp
View File

@ -80,21 +80,25 @@ MPInterval & MPInterval::operator= (const MPInterval &op) {
// UPDATE VALUES
mpfr_set (this->a, op.a, MPFR_RNDD);
mpfr_set (this->b, op.b, MPFR_RNDU);
return *this;
}
MPInterval & MPInterval::operator= (const Interval &op) {
// UPDATE VALUES
mpfr_set_d (this->a, op.a, MPFR_RNDD);
mpfr_set_d (this->b, op.b, MPFR_RNDU);
return *this;
}
MPInterval & MPInterval::operator= (const mpfr_t op) {
// UPDATE VALUES
mpfr_set (this->a, op, MPFR_RNDD);
mpfr_set (this->b, op, MPFR_RNDU);
return *this;
}
MPInterval & MPInterval::operator= (double op) {
// UPDATE VALUES
mpfr_set_d (this->a, op, MPFR_RNDD);
mpfr_set_d (this->b, op, MPFR_RNDU);
return *this;
}

2
lib/sirocco.cpp
View File

@ -17,7 +17,9 @@ using namespace std;
// z2 = t(1+0.01i)
******************************************/
#ifndef __clang__
#pragma STDC FENV_ACCESS ON
#endif
// PERFORM A MINIMUM OF 5 ITERATIONS OF NEWTON METHOD TO CORRECT THE APPROXIMATION OF THE ROOT (SUPPOSED TO BE GOOD ENOUGH)
template <>