Add functions to treat the factors of the polynomial separatedly.
This commit is contained in:
parent
8dffa82a98
commit
dbdc52aeb3
GPG Key ID:
DD0114EDB7410592
32
README.md
32
README.md
@ -10,7 +10,7 @@ Just install it as any other optional Sage package:
|
||||
sage -i sirocco
|
||||
```
|
||||
|
||||
and you can use it from Sage as explained in the Sage documentation for
|
||||
and you can use it from Sage as explained in the Sage documentation for
|
||||
[affine](https://doc.sagemath.org/html/en/reference/curves/sage/schemes/curves/affine_curve.html#sage.schemes.curves.affine_curve.AffinePlaneCurve.fundamental_group)
|
||||
and [projective](https://doc.sagemath.org/html/en/reference/curves/sage/schemes/curves/projective_curve.html#sage.schemes.curves.projective_curve.ProjectivePlaneCurve.fundamental_group) curves.
|
||||
|
||||
@ -277,8 +277,32 @@ make install
|
||||
|
||||
as before.
|
||||
|
||||
## To do:
|
||||
## Multiple components:
|
||||
|
||||
- include an option for passing the factors of the polynomial, for the cases
|
||||
where the original polynomial is reducible. This might speed up the computations
|
||||
In the case where the polynomial is reducible, its roots can be tracked for each
|
||||
factor separatedly. In that case, each root is tracked for its corresponding factor,
|
||||
making sure that there is no root of the other factors in the corresponding bounding
|
||||
box.
|
||||
|
||||
The functions to do so are `homotopyPath_comps` and `homotopyPath_mp_comps`.
|
||||
Their interface is similar to the ones of `homotopyPath` and `homotopyPath_mp`, but
|
||||
including the other factors. That is:
|
||||
|
||||
- `int degree` : the degree of the factor whose root is being tracked
|
||||
- `double *_coef` : a malloc'ed array with the coefficients of that factor
|
||||
- `double _y0R` : the real part of the starting root
|
||||
- `double _y0I` : the imaginary part of the starting root
|
||||
- `int nothercomps` : the number of the other factors of the polynomial
|
||||
- `int *degreescomps` : an array with the degrees of the other factors
|
||||
- `double *_coefscomps` : an array with the coefficients of the other factors, all concatenated in just one array.
|
||||
|
||||
The multiple precission version works similarly, but changing the doubles by mpfr numbers, and having an extra parameter with the used precission:
|
||||
|
||||
- `int degree` : the degree of the factor whose root is being tracked
|
||||
- `mpfr_t *_coef` : a malloc'ed array with the coefficients of that factor
|
||||
- `mpfr_t _y0R` : the real part of the starting root
|
||||
- `mpfr_t _y0I` : the imaginary part of the starting root
|
||||
- `int prec` : the used precission
|
||||
- `int nothercomps` : the number of the other factors of the polynomial
|
||||
- `int *degreescomps` : an array with the degrees of the other factors
|
||||
- `mpfr_t *_coefscomps` : an array with the coefficients of the other factors, all concatenated in just one array.
|
||||
|
||||
@ -13,4 +13,8 @@
|
||||
double * homotopyPath (int degree, double *_coef, double _y0R, double _y0I);
|
||||
mpfr_t * homotopyPath_mp (int degree, mpfr_t *_coef, mpfr_t _y0R, mpfr_t _y0I, int prec);
|
||||
|
||||
double * homotopyPath_comps (int degree, double *_coef, double _y0R, double _y0I,
|
||||
int nothercomps, int *degreescomps, double *_coefscomps);
|
||||
mpfr_t * homotopyPath_mp_comps (int degree, mpfr_t *_coef, mpfr_t _y0R, mpfr_t _y0I, int prec, int nothercomps, int *degreescomps, mpfr_t *_coefscomps);
|
||||
|
||||
#endif
|
||||
|
||||
@ -5,6 +5,7 @@
|
||||
#include <fenv.h>
|
||||
#include <mp_complex.hpp>
|
||||
#include <polynomial.hpp>
|
||||
#include <stdio.h>
|
||||
|
||||
|
||||
// ACTIVATE ROUND MODE CONTROL
|
||||
@ -17,6 +18,9 @@ typedef std::complex<double> Complex;
|
||||
double * homotopyPath (int degree, double *_coef, double _y0R, double _y0I);
|
||||
mpfr_t * homotopyPath_mp (int degree, mpfr_t *_coef, mpfr_t _y0R, mpfr_t _y0I, int prec);
|
||||
|
||||
double * homotopyPath_comps (int degree, double *_coef, double _y0R, double _y0I, int nothercomps, int *degreescomps, double *_coefscomps);
|
||||
mpfr_t * homotopyPath_mp_comps (int degree, mpfr_t *_coef, mpfr_t _y0R, mpfr_t _y0I, int prec, int nothercomps, int *degreescomps, mpfr_t *_coefscomps);
|
||||
|
||||
|
||||
template <class T> // T -> Complex, IComplex, MPComplex or MPIComplex
|
||||
void correctRoot (const Polynomial<T>, const T &x0, T &y0) {}
|
||||
@ -72,4 +76,35 @@ int validate (const Polynomial<T> &If, const T &Ix0, const T &Iy0, const T &IY0)
|
||||
|
||||
return 1;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
int validate (const Polynomial<T> &If, const T &Ix0, const T &Iy0, const T &IY0, int nothercomps, const Polynomial<T> othercomps[]) {
|
||||
|
||||
|
||||
T IY1;
|
||||
if (!newtonTest<T> (If, Ix0, Iy0, IY0)) return 0;
|
||||
|
||||
|
||||
// COMPUTES AN ENCLOSURE THREE TIMES BIGGER FOR N2
|
||||
// CREATE A 3-TIMES BIGGER ENCLOSURE (ROUNDING TO BIGGER)
|
||||
IY1 = IY0 + IY0;
|
||||
IY1 = IY1 - IY0;
|
||||
|
||||
if (!newtonTest<T> (If, Ix0, Iy0, IY1)) return 0;
|
||||
|
||||
T evalcomp;
|
||||
for (int i=0; i<nothercomps; ++i) {
|
||||
|
||||
evalcomp = othercomps[i](Ix0, IY1);
|
||||
|
||||
|
||||
if ( evalcomp.containsZero() ) {
|
||||
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
return 1;
|
||||
}
|
||||
#endif
|
||||
|
||||
566
lib/sirocco.cpp
566
lib/sirocco.cpp
@ -4,6 +4,7 @@
|
||||
#include <sirocco.hpp>
|
||||
#include <list.hpp>
|
||||
|
||||
|
||||
using namespace std;
|
||||
|
||||
// unsigned int *_Polynomial::combinatorial;
|
||||
@ -69,6 +70,284 @@ void correctRoot<MPComplex> (const Polynomial<MPComplex> f, const MPComplex &x0,
|
||||
return;
|
||||
}
|
||||
|
||||
double * homotopyPath_comps (int degree, double *_coef, double _y0R, double _y0I, int nothercomps, int *degreescomps, double *_coefscomps) {
|
||||
|
||||
int tdegree = degree;
|
||||
|
||||
|
||||
for (int i=0; i< nothercomps; ++i) {
|
||||
|
||||
tdegree = max(tdegree, degreescomps[i]);
|
||||
}
|
||||
|
||||
_Polynomial::startCombinatorialNumbers (tdegree);
|
||||
|
||||
// /************************************************/
|
||||
// /* VARIABLE DECLARATION AND DATA INITIALIZATION */
|
||||
// /************************************************/
|
||||
Complex aux; // Complex aux variable.
|
||||
|
||||
double stepsize = 0.0001; // step size
|
||||
double eps; // radius of the box
|
||||
|
||||
|
||||
int nCoef = ((degree+1)*(degree+2)) / 2; // number of monomials
|
||||
Complex coef[nCoef]; // coefficient list for floating point representation of
|
||||
// the polynomial.
|
||||
IComplex Icoef[nCoef]; // coefficient list for Interval representation of
|
||||
// the polynomial.
|
||||
|
||||
std::vector<Cell<double> > outputList;
|
||||
|
||||
|
||||
|
||||
// INITIALIZE COEFICIENT DATA FROM ARGUMENT INTO APPROPIATE DATA STRUCTURES
|
||||
for (int i=0; i<nCoef; ++i) {
|
||||
Icoef[i] = IComplex (_coef[4*i] ,_coef[4*i+1] ,_coef[4*i+2] ,_coef[4*i+3]);
|
||||
coef[i] = Complex ((_coef[4*i] + _coef[4*i+1])/2.0 ,
|
||||
(_coef[4*i+2] + _coef[4*i+3])/2.0);
|
||||
}
|
||||
|
||||
|
||||
Polynomial<Complex> f (degree, coef); // complex polynomials (with fp arithmetic)
|
||||
// to be used in Newton method correction
|
||||
Polynomial<IComplex> If (degree, Icoef); // complex polynomials (with Interval arithmetic)
|
||||
Polynomial<IComplex> Ig (degree); // to be used in validated homotopy
|
||||
Polynomial<IComplex> Ih (degree);
|
||||
|
||||
|
||||
Polynomial<IComplex> othercomponents[nothercomps];
|
||||
Polynomial<IComplex> Igs[nothercomps];
|
||||
Polynomial<IComplex> Ihs[nothercomps];
|
||||
|
||||
int acumcoefscomps = 0;
|
||||
|
||||
for (int i=0; i< nothercomps; ++i) {
|
||||
int cdegree = degreescomps[i];
|
||||
int nCoefcomp = ((cdegree+1)*(cdegree+2)) / 2;
|
||||
IComplex Icoefcomp[nCoefcomp];
|
||||
for (int j=0; j<nCoefcomp; ++j) {
|
||||
Icoefcomp[j] = IComplex (_coefscomps[4*acumcoefscomps] ,_coefscomps[4*acumcoefscomps+1] ,_coefscomps[4*acumcoefscomps+2] ,_coefscomps[acumcoefscomps*4+3]);
|
||||
acumcoefscomps++;
|
||||
}
|
||||
Polynomial<IComplex> Ifcomp (cdegree, Icoefcomp);
|
||||
othercomponents[i] = Ifcomp;
|
||||
Polynomial<IComplex > Igcomp (cdegree);
|
||||
Polynomial<IComplex > Ihcomp (cdegree);
|
||||
Igs[i] = Igcomp;
|
||||
Ihs[i] = Ihcomp;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
// INITIALIZE FLOATING POINT COMPLEX APPROXIMATION OF THE ROOT
|
||||
Complex x0 (0.0); // x0 of the root (x0.imag = 0)
|
||||
Complex y0 (_y0R, _y0I); // approximated y0 of the root
|
||||
IComplex Ix0; // Interval representation for x0
|
||||
IComplex Iy0; // Interval representation for y0
|
||||
IComplex IY0; // Interval enclosure for y0 for Newton Method
|
||||
IComplex Iy1, IY1; // Interval representation of zeros for rotated polynomial
|
||||
IComplex Ix1;
|
||||
|
||||
Complex a; // y'(x0) as Complex
|
||||
IComplex Ia; // y'(x0) as IComplex
|
||||
|
||||
|
||||
#ifdef PROFILING
|
||||
int nsteps = 0; // number of steps performed
|
||||
int nrejectedEps = 0; // number of rejected validations of boxes
|
||||
int nrejectedSteps = 0; // number of rejected validations of tubes
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
outputList.push_back (Cell<double> (0.0, y0.real (), y0.imag ()));
|
||||
|
||||
/*******************/
|
||||
/* CORRECT ROOT */
|
||||
/*******************/
|
||||
|
||||
|
||||
correctRoot<Complex> (f, x0, y0);
|
||||
|
||||
|
||||
/*******************/
|
||||
/** MAIN LOOP */
|
||||
/*******************/
|
||||
while (x0.real () < 1.0) {
|
||||
|
||||
#ifdef PROFILING
|
||||
nsteps++; // used for profiling only
|
||||
#endif
|
||||
|
||||
/***********************************************************/
|
||||
/*** BOX SIZE ESTIMATION. F RESTRICTED TO X0 IS INJECTIVE. */
|
||||
/***********************************************************/
|
||||
aux = f.diffYY (x0, y0);
|
||||
if (abs (aux) < 1.0e-10) {
|
||||
eps = 0.5;
|
||||
} else {
|
||||
eps = (abs (f.diffY(x0,y0) / aux)) / 8.0;
|
||||
}
|
||||
|
||||
eps /= 2.0;
|
||||
// cout << "estima eps = " << eps << endl;
|
||||
|
||||
/**************************/
|
||||
/*** VALIDATE y0 AS POINT.*/
|
||||
/**************************/
|
||||
Ix0 = x0;
|
||||
Iy0 = y0;
|
||||
IY0 = y0 + IComplex (-eps, eps, -eps, eps);
|
||||
|
||||
|
||||
while (!validate<IComplex> (If, Ix0, Iy0, IY0, nothercomps, othercomponents)) {
|
||||
#ifdef PROFILING
|
||||
nrejectedEps++; // used for profiling
|
||||
#endif
|
||||
if (eps < 1.0e-13) {
|
||||
printf ("Singularity detected. Abort Sirocco\n");
|
||||
/*******************/
|
||||
/* CLEAN VARIABLES */
|
||||
/*******************/
|
||||
|
||||
_Polynomial::releaseCombinatorialNumbers ();
|
||||
return NULL;
|
||||
}
|
||||
eps *= 0.5;
|
||||
IY0 = Iy0;
|
||||
IY0 = y0 + IComplex (-eps, eps, -eps, eps);
|
||||
}
|
||||
// cout << "acepta eps = " << eps << endl;
|
||||
|
||||
// COMPUTE a = -fx(x0,y0) / fy(x0.y0)
|
||||
// both fp and Interval
|
||||
a = getA<Complex> (f, x0, y0);
|
||||
Ia = a;
|
||||
|
||||
|
||||
/************************/
|
||||
/* ESTIMATE STEP SIZE ***/
|
||||
/************************/
|
||||
// we have: aux = fyy
|
||||
// store all in aux2
|
||||
|
||||
aux = a*a*aux + 2.0*a*f.diffXY(x0,y0) + f.diffXX (x0, y0);
|
||||
|
||||
// DETECT INFLEXION POINT
|
||||
if (abs(aux) < 1.0e-10)
|
||||
stepsize=1.0;
|
||||
else
|
||||
stepsize = sqrt(abs((eps * f.diffY (x0, y0) )/ aux));
|
||||
if (stepsize + x0.real () > 1.0) stepsize = 1.000001 - x0.real ();
|
||||
|
||||
|
||||
/*******************/
|
||||
/** VALIDATE STEP */
|
||||
/*******************/
|
||||
Ix1 = Interval(0,stepsize);
|
||||
|
||||
|
||||
// PERFORM TRANSLATION AND ROTATION OF THE POLYNOMIAL
|
||||
Ig = If.IlinearVarChange2 (Ix0);
|
||||
Ih = Ig.IlinearVarChange (Ia);
|
||||
|
||||
for (int jj = 0; jj<nothercomps; ++jj) {
|
||||
|
||||
Igs[jj] = othercomponents[jj].IlinearVarChange2(Ix0);
|
||||
Ihs[jj] = Igs[jj].IlinearVarChange (Ia);
|
||||
}
|
||||
|
||||
|
||||
while (!validate<IComplex> (Ih, Ix1, Iy0, IY0, nothercomps, Ihs)) {
|
||||
#ifdef PROFILING
|
||||
nrejectedSteps++; // used for profiling
|
||||
#endif
|
||||
stepsize *= 0.5;
|
||||
eps*=0.95;
|
||||
IY0 = y0 + IComplex (-eps, eps, -eps, eps);
|
||||
Ix1 = Interval (0, stepsize);
|
||||
|
||||
if (stepsize < 1.0e-13) {
|
||||
/*******************/
|
||||
/* CLEAN VARIABLES */
|
||||
/*******************/
|
||||
_Polynomial::releaseCombinatorialNumbers ();
|
||||
return NULL;
|
||||
}
|
||||
|
||||
}
|
||||
#ifdef DEVELOPER
|
||||
std::cout << "x0 = " << x0.real () << " step = " << stepsize << std::endl;
|
||||
#endif
|
||||
|
||||
|
||||
// /**************************/
|
||||
// /** VALIDATE FINAL STEP */
|
||||
// /**************************/
|
||||
if (x0.real () + stepsize > 1.0) {
|
||||
stepsize = 1.0 - x0.real ();
|
||||
Ix1.r.b = stepsize;
|
||||
|
||||
#ifdef DEVELOPER
|
||||
std::cout << "final step = " << stepsize << std::endl;
|
||||
#endif
|
||||
|
||||
while (!validate<IComplex> (Ih, Ix1, Iy0, IY0, nothercomps, Ihs));
|
||||
}
|
||||
|
||||
|
||||
// /********************/
|
||||
// /* UPDATE VARIABLES */
|
||||
// /********************/
|
||||
x0 = Complex (Ix0.r.a + stepsize, 0.0);
|
||||
y0 = Complex (y0.real () + Ia.r.a * stepsize,
|
||||
y0.imag () + Ia.i.a * stepsize);
|
||||
|
||||
|
||||
|
||||
|
||||
// /*******************/
|
||||
// /* CORRECT ROOT */
|
||||
// /*******************/
|
||||
correctRoot<Complex> (f, x0, y0);
|
||||
|
||||
/**********************************************/
|
||||
/* CHECK WE HAVE NOT JUMPED TO ANOTHER ROOT */
|
||||
/**********************************************/
|
||||
if (!(IY0+Ia*stepsize).contains (y0)) {
|
||||
printf ("error! Jumped to other thread!\n");
|
||||
_Polynomial::releaseCombinatorialNumbers ();
|
||||
return NULL;
|
||||
}
|
||||
|
||||
|
||||
outputList.push_back (Cell<double> (x0.real (), y0.real (), y0.imag ()));
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// /***********************************************/
|
||||
/* PREPARE OUTPUT AS A BINARY ARRAY OF DOUBLES */
|
||||
/***********************************************/
|
||||
|
||||
double *rop = new double [((3 * outputList.size () + 1))];
|
||||
rop[0] = outputList.size ();
|
||||
for (int i=0; i<outputList.size (); i++) {
|
||||
rop[3*i+1] = outputList[i].x[0];
|
||||
rop[3*i+2] = outputList[i].x[1];
|
||||
rop[3*i+3] = outputList[i].x[2];
|
||||
}
|
||||
_Polynomial::releaseCombinatorialNumbers ();
|
||||
return rop;
|
||||
|
||||
}
|
||||
|
||||
|
||||
double * homotopyPath (int degree, double *_coef, double _y0R, double _y0I) {
|
||||
|
||||
_Polynomial::startCombinatorialNumbers (degree);
|
||||
@ -556,3 +835,290 @@ mpfr_set_default_prec(oldprec);
|
||||
_Polynomial::releaseCombinatorialNumbers ();
|
||||
return rop;
|
||||
}
|
||||
|
||||
|
||||
mpfr_t * homotopyPath_mp_comps (int degree, mpfr_t *_coef, mpfr_t _y0R, mpfr_t _y0I, int prec, int nothercomps, int *degreescomps, mpfr_t *_coefscomps) {
|
||||
mpfr_prec_t oldprec = mpfr_get_default_prec();
|
||||
mpfr_set_default_prec(prec);
|
||||
|
||||
|
||||
int tdegree = degree;
|
||||
for (int i=0; i< nothercomps; ++i) {
|
||||
|
||||
tdegree = max(tdegree, degreescomps[i]);
|
||||
}
|
||||
|
||||
_Polynomial::startCombinatorialNumbers (tdegree);
|
||||
|
||||
|
||||
// // /************************************************/
|
||||
// // /* VARIABLE DECLARATION AND DATA INITIALIZATION */
|
||||
// // /************************************************/
|
||||
MPComplex aux; // Complex aux variable.
|
||||
|
||||
mpfr_t stepsize , eps, calc, tol;
|
||||
mpfr_inits (stepsize, eps, calc, tol, NULL); // step size and radius of the box
|
||||
mpfr_set_d (stepsize, 0.0001, MPFR_RNDN);
|
||||
mpfr_set_d (tol, mpfr_get_default_prec(), MPFR_RNDN);
|
||||
mpfr_sub_si (tol, tol, 10, MPFR_RNDN);
|
||||
mpfr_neg (tol, tol, MPFR_RNDN);
|
||||
mpfr_exp2 (tol, tol, MPFR_RNDN);
|
||||
|
||||
|
||||
int nCoef = ((degree+1)*(degree+2)) / 2; // number of monomials
|
||||
MPComplex coef[nCoef]; // coefficient list for floating point representation of
|
||||
// the polynomial.
|
||||
MPIComplex Icoef[nCoef]; // coefficient list for Interval representation of
|
||||
// the polynomial.
|
||||
|
||||
std::vector<Cell<mpfr_t> > outputList;
|
||||
|
||||
|
||||
|
||||
// INITIALIZE COEFICIENT DATA FROM ARGUMENT INTO APPROPRIATE DATA STRUCTURES
|
||||
for (int i=0; i<nCoef; ++i) {
|
||||
Icoef[i] = MPIComplex (_coef[4*i] ,_coef[4*i+1] ,_coef[4*i+2] ,_coef[4*i+3]);
|
||||
coef[i] = 0.5*(MPComplex (_coef[4*i] , _coef[4*i+2]) +
|
||||
MPComplex (_coef[4*i+1] , _coef[4*i+3]));
|
||||
}
|
||||
|
||||
|
||||
Polynomial<MPComplex> f (degree, coef); // complex polynomials (with fp arithmetic)
|
||||
// to be used in Newton method correction
|
||||
Polynomial<MPIComplex> If (degree, Icoef); // complex polynomials (with Interval arithmetic)
|
||||
Polynomial<MPIComplex> Ig (degree); // to be used in validated homotopy
|
||||
Polynomial<MPIComplex> Ih (degree);
|
||||
|
||||
Polynomial<MPIComplex> othercomponents[nothercomps];
|
||||
Polynomial<MPIComplex> Igs[nothercomps];
|
||||
Polynomial<MPIComplex> Ihs[nothercomps];
|
||||
|
||||
int acumcoefscomps = 0;
|
||||
|
||||
for (int i=0; i< nothercomps; ++i) {
|
||||
int cdegree = degreescomps[i];
|
||||
int nCoefcomp = ((cdegree+1)*(cdegree+2)) / 2;
|
||||
MPIComplex MPIcoefcomp[nCoefcomp];
|
||||
for (int j=0; j<nCoefcomp; ++j) {
|
||||
MPIcoefcomp[j] = MPIComplex (_coefscomps[4*acumcoefscomps] ,_coefscomps[4*acumcoefscomps+1] ,_coefscomps[4*acumcoefscomps+2] ,_coefscomps[acumcoefscomps*4+3]);
|
||||
acumcoefscomps++;
|
||||
}
|
||||
Polynomial<MPIComplex> Ifcomp (cdegree, MPIcoefcomp);
|
||||
othercomponents[i] = Ifcomp;
|
||||
Polynomial<MPIComplex > Igcomp (cdegree);
|
||||
Polynomial<MPIComplex > Ihcomp (cdegree);
|
||||
Igs[i] = Igcomp;
|
||||
Ihs[i] = Ihcomp;
|
||||
}
|
||||
|
||||
|
||||
|
||||
// // INITIALIZE FLOATING POINT COMPLEX APPROXIMATION OF THE ROOT
|
||||
MPComplex x0 (0.0); // x0 of the root (x0.imag = 0)
|
||||
MPComplex y0 (_y0R, _y0I); // approximated y0 of the root
|
||||
MPIComplex Ix0; // Interval representation for x0
|
||||
MPIComplex Iy0; // Interval representation for y0
|
||||
MPIComplex IY0; // Interval enclosure for y0 for Newton Method
|
||||
MPIComplex Iy1, IY1; // Interval representation of zeros for rotated polynomial
|
||||
MPIComplex Ix1;
|
||||
|
||||
MPComplex a; // y'(x0) as Complex
|
||||
MPIComplex Ia; // y'(x0) as MPIComplex
|
||||
|
||||
|
||||
#ifdef PROFILING
|
||||
int nsteps = 0; // number of steps performed
|
||||
int nrejectedEps = 0; // number of rejected validations of boxes
|
||||
int nrejectedSteps = 0; // number of rejected validations of tubes
|
||||
#endif
|
||||
|
||||
|
||||
outputList.push_back (Cell<mpfr_t> (x0.r, y0.r, y0.i));
|
||||
|
||||
// /*******************/
|
||||
// /* CORRECT ROOT */
|
||||
// /*******************/
|
||||
|
||||
|
||||
correctRoot<MPComplex> (f, x0, y0);
|
||||
|
||||
|
||||
// /*******************/
|
||||
// /** MAIN LOOP */
|
||||
// /*******************/
|
||||
while (mpfr_cmp_si(x0.r, 1.0) < 0) {
|
||||
|
||||
#ifdef PROFILING
|
||||
nsteps++; // used for profiling only
|
||||
#endif
|
||||
|
||||
/***********************************************************/
|
||||
/*** BOX SIZE ESTIMATION. F RESTRICTED TO X0 IS INJECTIVE. */
|
||||
/***********************************************************/
|
||||
aux = f.diffYY (x0, y0);
|
||||
abs (calc, aux);
|
||||
if (mpfr_cmp_d (calc, 1.0e-10) < 0) {
|
||||
mpfr_set_d (eps, 0.5, MPFR_RNDN);
|
||||
} else {
|
||||
abs (eps, (f.diffY(x0,y0) / aux));
|
||||
mpfr_div_si (eps, eps, 8, MPFR_RNDN);
|
||||
}
|
||||
|
||||
/**************************/
|
||||
/*** VALIDATE y0 AS POINT.*/
|
||||
/**************************/
|
||||
Ix0 = x0;
|
||||
Iy0 = y0;
|
||||
IY0 = MPIComplex(y0) + eps*MPIComplex (-1.0, 1.0, -1.0, 1.0);
|
||||
|
||||
|
||||
while (!validate<MPIComplex> (If, Ix0, Iy0, IY0, nothercomps, othercomponents)) {
|
||||
#ifdef PROFILING
|
||||
nrejectedEps++; // used for profiling
|
||||
#endif
|
||||
if (mpfr_cmp (eps ,tol) < 0) {
|
||||
printf ("Singularity detected. Abort Sirocco\n");
|
||||
/*******************/
|
||||
/* CLEAN VARIABLES */
|
||||
/*******************/
|
||||
mpfr_clears (stepsize, eps, calc, tol, NULL);
|
||||
|
||||
mpfr_set_default_prec(oldprec);
|
||||
|
||||
_Polynomial::releaseCombinatorialNumbers ();
|
||||
return NULL;
|
||||
}
|
||||
mpfr_div_si (eps, eps, 2, MPFR_RNDN);
|
||||
IY0 = Iy0;
|
||||
IY0 = MPIComplex(y0) + eps*MPIComplex (-1.0, 1.0, -1.0, 1.0);
|
||||
}
|
||||
|
||||
// COMPUTE a = -fx(x0,y0) / fy(x0.y0)
|
||||
// both fp and Interval
|
||||
a = getA<MPComplex> (f, x0, y0);
|
||||
Ia = a;
|
||||
|
||||
|
||||
/************************/
|
||||
/* ESTIMATE STEP SIZE ***/
|
||||
/************************/
|
||||
// we have: aux = fyy
|
||||
// store all in aux2
|
||||
|
||||
aux = a*a*aux + 2.0*a*f.diffXY(x0,y0) + f.diffXX (x0, y0);
|
||||
|
||||
// DETECT INFLEXION POINT
|
||||
abs (calc, aux);
|
||||
if (mpfr_cmp_d(calc, 1.0e-10) < 0) {
|
||||
mpfr_set_si (stepsize, 1, MPFR_RNDN);
|
||||
} else {
|
||||
abs (calc, (eps * f.diffY (x0, y0) )/ aux);
|
||||
mpfr_sqrt (stepsize, calc, MPFR_RNDN);
|
||||
}
|
||||
|
||||
// if (stepsize + x0.real () > 1.0) stepsize = 1.000001 - x0.real ();
|
||||
mpfr_add (calc, x0.r, stepsize, MPFR_RNDU);
|
||||
if (mpfr_cmp_si (calc, 1) > 0) {
|
||||
mpfr_si_sub (stepsize, 1, x0.r, MPFR_RNDU);
|
||||
}
|
||||
|
||||
|
||||
/*******************/
|
||||
/** VALIDATE STEP */
|
||||
/*******************/
|
||||
mpfr_set_si (calc, 0, MPFR_RNDN);
|
||||
Ix1 = MPInterval(calc,stepsize);
|
||||
|
||||
|
||||
// PERFORM TRANSLATION AND ROTATION OF THE POLYNOMIAL
|
||||
Ig = If.IlinearVarChange2 (Ix0);
|
||||
Ih = Ig.IlinearVarChange (Ia);
|
||||
|
||||
|
||||
for (int jj = 0; jj<nothercomps; ++jj) {
|
||||
|
||||
Igs[jj] = othercomponents[jj].IlinearVarChange2(Ix0);
|
||||
Ihs[jj] = Igs[jj].IlinearVarChange (Ia);
|
||||
}
|
||||
|
||||
|
||||
while (!validate<MPIComplex> (Ih, Ix1, Iy0, IY0, nothercomps, Ihs)) {
|
||||
#ifdef PROFILING
|
||||
nrejectedSteps++; // used for profiling
|
||||
#endif
|
||||
mpfr_div_si (stepsize, stepsize, 2, MPFR_RNDN);
|
||||
mpfr_mul_d (eps, eps, 0.95, MPFR_RNDN);
|
||||
|
||||
IY0 = MPIComplex(y0) + eps*MPIComplex (-1.0, 1.0, -1.0, 1.0);
|
||||
Ix1 = MPInterval (calc, stepsize);
|
||||
|
||||
if (mpfr_cmp_d (stepsize, 1.0e-13) < 0) {
|
||||
/*******************/
|
||||
/* CLEAN VARIABLES */
|
||||
/*******************/
|
||||
mpfr_clears (stepsize, eps, calc, tol, NULL);
|
||||
|
||||
mpfr_set_default_prec(oldprec);
|
||||
_Polynomial::releaseCombinatorialNumbers ();
|
||||
return NULL;
|
||||
}
|
||||
|
||||
}
|
||||
#ifdef DEVELOPER
|
||||
std::cout << "x0 = " << x0.r << " step = " << stepsize << std::endl;
|
||||
#endif
|
||||
|
||||
// // /********************/
|
||||
// // /* UPDATE VARIABLES */
|
||||
// // /********************/
|
||||
x0 = MPComplex (Ix0.r.a) + MPComplex (stepsize);
|
||||
y0 = y0 + stepsize*a;
|
||||
|
||||
|
||||
|
||||
|
||||
// /*******************/
|
||||
// /* CORRECT ROOT */
|
||||
// /*******************/
|
||||
correctRoot<MPComplex> (f, x0, y0);
|
||||
|
||||
/**********************************************/
|
||||
/* CHECK WE HAVE NOT JUMPED TO ANOTHER ROOT */
|
||||
/**********************************************/
|
||||
if (!(MPIComplex(y0).subset (IY0+Ia*stepsize))) {
|
||||
mpfr_clears (stepsize, eps, calc, tol, NULL);
|
||||
printf ("error! Jumped to other thread!\n");
|
||||
mpfr_set_default_prec(oldprec);
|
||||
_Polynomial::releaseCombinatorialNumbers ();
|
||||
return NULL;
|
||||
}
|
||||
|
||||
|
||||
outputList.push_back (Cell<mpfr_t> (x0.r, y0.r, y0.i));
|
||||
|
||||
}
|
||||
|
||||
|
||||
// /***********************************************/
|
||||
/* PREPARE OUTPUT AS A BINARY ARRAY OF DOUBLES */
|
||||
/***********************************************/
|
||||
|
||||
mpfr_t *rop = new mpfr_t [((3 * outputList.size () + 1))];
|
||||
mpfr_init(rop[0]);
|
||||
mpfr_set_si(rop[0], outputList.size(), MPFR_RNDN);
|
||||
for (int i=0; i<outputList.size (); i++) {
|
||||
mpfr_init(rop[3*i+1]);
|
||||
mpfr_init(rop[3*i+2]);
|
||||
mpfr_init(rop[3*i+3]);
|
||||
mpfr_set(rop[3*i+1], outputList[i].x[0], MPFR_RNDN);
|
||||
mpfr_set(rop[3*i+2], outputList[i].x[1], MPFR_RNDN);
|
||||
mpfr_set(rop[3*i+3], outputList[i].x[2], MPFR_RNDN);
|
||||
}
|
||||
|
||||
|
||||
mpfr_clears (stepsize, eps, calc, tol, NULL);
|
||||
|
||||
mpfr_set_default_prec(oldprec);
|
||||
_Polynomial::releaseCombinatorialNumbers ();
|
||||
return rop;
|
||||
}
|
||||
|
||||
@ -3,26 +3,29 @@
|
||||
#clang C++
|
||||
#clib sirocco
|
||||
|
||||
include 'cysignals/signals.pxi'
|
||||
#include 'cysignals/signals.pxi'
|
||||
|
||||
from cysignals.signals cimport sig_on, sig_off
|
||||
from cysignals.memory cimport check_allocarray, sig_free as free
|
||||
|
||||
from sage.libs.mpfr cimport *
|
||||
from sage.rings.real_mpfr cimport RealNumber
|
||||
from sage.rings.real_mpfr import RealField
|
||||
|
||||
|
||||
cdef extern from "stdlib.h":
|
||||
void free(void* ptr)
|
||||
|
||||
|
||||
cdef extern from "sirocco.h":
|
||||
mpfr_t* homotopyPath_mp (int degree, mpfr_t *_coef, mpfr_t _y0R, mpfr_t _y0I, int prec)
|
||||
mpfr_t* homotopyPath_mp_comps (int degree, mpfr_t *_coef, mpfr_t _y0R, mpfr_t _y0I, int prec, int nothercomps, int *degreescomps, mpfr_t *_coefscomps)
|
||||
double * homotopyPath (int degree, double *_coef, double _y0R, double _y0I)
|
||||
double * homotopyPath_comps (int degree, double *_coef, double _y0R, double _y0I, int nothercomps, int *degreescomps, double *_coefscomps)
|
||||
|
||||
|
||||
|
||||
def contpath_mp(deg, values, RealNumber y0r, RealNumber y0i, int prec):
|
||||
cdef int cdeg = deg
|
||||
cdef mpfr_t* cvalues = <mpfr_t*> sage_malloc(sizeof(mpfr_t)*len(values))
|
||||
cdef mpfr_t* cvalues = <mpfr_t*> check_allocarray(len(values), sizeof(mpfr_t))
|
||||
cdef mpfr_t* rop
|
||||
cdef int j
|
||||
for i in range(len(values)):
|
||||
@ -66,16 +69,82 @@ def contpath_mp(deg, values, RealNumber y0r, RealNumber y0i, int prec):
|
||||
|
||||
|
||||
|
||||
def contpath_mp_comps(deg, values, RealNumber y0r, RealNumber y0i, int prec, otherdegrees, othercoefs):
|
||||
cdef int cdeg = deg
|
||||
cdef mpfr_t* cvalues = <mpfr_t*> check_allocarray(len(values),sizeof(mpfr_t))
|
||||
cdef mpfr_t* cothercoefs = <mpfr_t*> check_allocarray(len(othercoefs),sizeof(mpfr_t))
|
||||
cdef mpfr_t* rop
|
||||
cdef int j
|
||||
for i in range(len(values)):
|
||||
j = <int>i
|
||||
sig_on()
|
||||
mpfr_init2(cvalues[j], prec)
|
||||
mpfr_set(cvalues[j], (<RealNumber>values[j]).value, MPFR_RNDN)
|
||||
sig_off()
|
||||
|
||||
for i in range(len(othercoefs)):
|
||||
j = <int>i
|
||||
sig_on()
|
||||
mpfr_init2(cothercoefs[j], prec)
|
||||
mpfr_set(cothercoefs[j], (<RealNumber>othercoefs[j]).value, MPFR_RNDN)
|
||||
sig_off()
|
||||
|
||||
cdef int* c_otherdegrees = <int*> check_allocarray(len(otherdegrees),sizeof(int))
|
||||
for i in range(len(otherdegrees)):
|
||||
c_otherdegrees[i] = int(otherdegrees[i])
|
||||
|
||||
cdef mpfr_t y0R
|
||||
cdef mpfr_t y0I
|
||||
sig_on()
|
||||
mpfr_init2(y0R, prec)
|
||||
mpfr_set(y0R, (<RealNumber>y0r).value, MPFR_RNDN)
|
||||
mpfr_init2(y0I, prec)
|
||||
mpfr_set(y0I, (<RealNumber>y0i).value, MPFR_RNDN)
|
||||
rop = homotopyPath_mp_comps(cdeg, cvalues, y0R, y0I, prec, int(len(otherdegrees)), c_otherdegrees, cothercoefs)
|
||||
sig_off()
|
||||
for i in range(len(values)):
|
||||
j = <int>i
|
||||
sig_on()
|
||||
mpfr_clear(cvalues[j])
|
||||
sig_off()
|
||||
|
||||
for i in range(len(othercoefs)):
|
||||
j = <int>i
|
||||
sig_on()
|
||||
mpfr_clear(cothercoefs[j])
|
||||
sig_off()
|
||||
|
||||
free(cvalues)
|
||||
free(cothercoefs)
|
||||
free(c_otherdegrees)
|
||||
if rop == NULL:
|
||||
raise ValueError("libsirocco could not guarantee one step")
|
||||
cdef int n = mpfr_get_si(rop[0], MPFR_RNDN)
|
||||
|
||||
#cdef RealNumber res = RealField(53)()
|
||||
cdef double res
|
||||
l = []
|
||||
|
||||
for i in range(1, 3*n+1):
|
||||
RN = RealField(prec)()
|
||||
sig_on()
|
||||
res = mpfr_set((<RealNumber>RN).value, rop[i], MPFR_RNDN)
|
||||
mpfr_clear(rop[i])
|
||||
sig_off()
|
||||
l.append(RN)
|
||||
free(rop)
|
||||
return [(l[3*i], l[3*i+1], l[3*i+2]) for i in range(len(l)/3)]
|
||||
|
||||
|
||||
|
||||
|
||||
def contpath(deg,values,y0r,y0i):
|
||||
cdef double* rop
|
||||
cdef double* c_values = <double*> sage_malloc(sizeof(double)*len(values))
|
||||
cdef double* c_values = <double*> check_allocarray(len(values),sizeof(double))
|
||||
cdef int clen = <int> len(values)
|
||||
for i,v in enumerate(values):
|
||||
c_values[i] = values[i]
|
||||
|
||||
|
||||
|
||||
cdef double y0R = y0r
|
||||
cdef double y0I = y0i
|
||||
sig_on()
|
||||
@ -90,3 +159,37 @@ def contpath(deg,values,y0r,y0i):
|
||||
free(rop)
|
||||
free(c_values)
|
||||
return l
|
||||
|
||||
|
||||
def contpath_comps(deg, values, y0r, y0i, otherdegrees, othercomps):
|
||||
cdef double* rop
|
||||
cdef double* c_values = <double*> check_allocarray(len(values),sizeof(double))
|
||||
for i,v in enumerate(values):
|
||||
c_values[i] = values[i]
|
||||
|
||||
cdef double* c_othercomps = <double*> check_allocarray(len(othercomps),sizeof(double))
|
||||
|
||||
for i,v in enumerate(othercomps):
|
||||
c_othercomps[i] = othercomps[i]
|
||||
|
||||
cdef int* c_otherdegrees = <int *> check_allocarray(len(otherdegrees),sizeof(int))
|
||||
|
||||
for i,v in enumerate(otherdegrees):
|
||||
c_otherdegrees[i] = int(otherdegrees[i])
|
||||
|
||||
cdef double y0R = y0r
|
||||
cdef double y0I = y0i
|
||||
sig_on()
|
||||
rop = homotopyPath_comps (int(deg), c_values, y0R, y0I, int(len(otherdegrees)), c_otherdegrees, c_othercomps)
|
||||
sig_off()
|
||||
if rop == NULL:
|
||||
raise ValueError("libsirocco could not guarantee one step")
|
||||
n=int(rop[0])
|
||||
l=[0 for i in range(n)]
|
||||
for i in range(n):
|
||||
l[i]=(rop[3*i+1],rop[3*i+2],rop[3*i+3])
|
||||
free(rop)
|
||||
free(c_values)
|
||||
free(c_othercomps)
|
||||
free(c_otherdegrees)
|
||||
return l
|
||||
|
||||
@ -16,6 +16,23 @@ int main(int argc, char const *argv[]) {
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
double testCoefsotcomp[] = {
|
||||
-2.100000000000, -2.10000000000000, 0.000000000000000, 0.000000000000000,
|
||||
0.0000000000000, 0.00000000000000, 0.000000000000000, -0.000000000000000,
|
||||
0.000000000000, 0.000000000000, 0.0000000000000, 0.0000000000000,
|
||||
1.000000000000000, 1.000000000000000, 0.000000000000000, -0.000000000000000,
|
||||
0.00000000000000, 0.000000000000000, 0.000000000000000, 0.0000000000000000,
|
||||
1.0000000000000, 1.000000000000000, 0.000000000000000, 0.0000000000000000
|
||||
};
|
||||
|
||||
|
||||
|
||||
int compdegs[]={2};
|
||||
|
||||
|
||||
|
||||
double y0R = 1.414;
|
||||
double y0I = 0.0;
|
||||
double *rop = homotopyPath (2, testCoefs, y0R, y0I);
|
||||
@ -25,28 +42,46 @@ int main(int argc, char const *argv[]) {
|
||||
|
||||
delete [] rop;
|
||||
|
||||
double *rop2 = homotopyPath_comps (2, testCoefs, y0R, y0I,1, compdegs, testCoefsotcomp);
|
||||
|
||||
delete [] rop2;
|
||||
|
||||
mpfr_t mp_testCoefs[24];
|
||||
for (int i=0; i<24; ++i) {
|
||||
mpfr_init_set_d (mp_testCoefs[i], testCoefs[i],
|
||||
(i%2)? MPFR_RNDU : MPFR_RNDD);
|
||||
}
|
||||
|
||||
mpfr_t mp_othercompsCoefs[24];
|
||||
for (int i=0; i<24; ++i) {
|
||||
mpfr_init_set_d (mp_othercompsCoefs[i], testCoefsotcomp[i],
|
||||
(i%2)? MPFR_RNDU : MPFR_RNDD);
|
||||
}
|
||||
|
||||
|
||||
mpfr_t mp_y0R, mp_y0I;
|
||||
mpfr_init_set_d (mp_y0R, y0R, MPFR_RNDN);
|
||||
mpfr_init_set_d (mp_y0I, y0I, MPFR_RNDN);
|
||||
|
||||
mpfr_t *mp_rop = homotopyPath_mp (2, mp_testCoefs, mp_y0R, mp_y0I, 106);
|
||||
mpfr_t *mp_rop2 = homotopyPath_mp_comps (2, mp_testCoefs, mp_y0R, mp_y0I,102,1,compdegs, mp_othercompsCoefs);
|
||||
|
||||
int size = mpfr_get_si (mp_rop[0], MPFR_RNDN);
|
||||
|
||||
|
||||
for (int i=0; i<3*size+1; ++i)
|
||||
mpfr_clear(mp_rop[i]);
|
||||
delete [] mp_rop;
|
||||
delete [] mp_rop2;
|
||||
|
||||
|
||||
for (int i=0; i<24; ++i) {
|
||||
mpfr_clear (mp_testCoefs[i]);
|
||||
}
|
||||
|
||||
for (int i=0; i<24; ++i) {
|
||||
mpfr_clear (mp_othercompsCoefs[i]);
|
||||
}
|
||||
|
||||
mpfr_clear (mp_y0R);
|
||||
mpfr_clear (mp_y0I);
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user