Script 'mail_helper' called by obssrc Hello community, here is the log from the commit of package sirocco for openSUSE:Factory checked in at 2021-07-01 07:05:46 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Comparing /work/SRC/openSUSE:Factory/sirocco (Old) and /work/SRC/openSUSE:Factory/.sirocco.new.2625 (New) ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Package is "sirocco" Thu Jul 1 07:05:46 2021 rev:2 rq:903226 version:2.1.0 Changes: -------- --- /work/SRC/openSUSE:Factory/sirocco/sirocco.changes 2020-09-15 16:27:33.322554579 +0200 +++ /work/SRC/openSUSE:Factory/.sirocco.new.2625/sirocco.changes 2021-07-01 07:06:03.375279384 +0200 @@ -1,0 +2,7 @@ +Wed Jun 30 11:20:09 UTC 2021 - Jan Engelhardt <[email protected]> + +- Update to release 2.1.0 + * This release adds support for treating the different factors + of a polynomial separately. + +------------------------------------------------------------------- Old: ---- libsirocco-2.0.2.tar.gz New: ---- libsirocco-2.1.0.tar.gz ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Other differences: ------------------ ++++++ sirocco.spec ++++++ --- /var/tmp/diff_new_pack.tRwOki/_old 2021-07-01 07:06:03.811275977 +0200 +++ /var/tmp/diff_new_pack.tRwOki/_new 2021-07-01 07:06:03.811275977 +0200 @@ -1,7 +1,7 @@ # # spec file for package sirocco # -# Copyright (c) 2020 SUSE LLC +# Copyright (c) 2021 SUSE LLC # # All modifications and additions to the file contributed by third parties # remain the property of their copyright owners, unless otherwise agreed @@ -18,7 +18,7 @@ %define lname libsirocco0 Name: sirocco -Version: 2.0.2 +Version: 2.1.0 Release: 0 Summary: Library for computing homotopy continuation of roots License: GPL-3.0-or-later ++++++ libsirocco-2.0.2.tar.gz -> libsirocco-2.1.0.tar.gz ++++++ ++++ 6569 lines of diff (skipped) ++++ retrying with extended exclude list diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/README.md new/libsirocco-2.1.0/README.md --- old/libsirocco-2.0.2/README.md 2020-02-16 20:38:32.000000000 +0100 +++ new/libsirocco-2.1.0/README.md 2021-06-30 11:56:56.000000000 +0200 @@ -10,7 +10,7 @@ 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 @@ as before. -## To do: - - - 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 +## Multiple components: +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. diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/build-aux/ar-lib new/libsirocco-2.1.0/build-aux/ar-lib --- old/libsirocco-2.0.2/build-aux/ar-lib 2020-02-08 11:56:17.000000000 +0100 +++ new/libsirocco-2.1.0/build-aux/ar-lib 2017-02-28 10:35:32.000000000 +0100 @@ -4,7 +4,7 @@ me=ar-lib scriptversion=2012-03-01.08; # UTC -# Copyright (C) 2010-2018 Free Software Foundation, Inc. +# Copyright (C) 2010-2014 Free Software Foundation, Inc. # Written by Peter Rosin <[email protected]>. # # This program is free software; you can redistribute it and/or modify @@ -18,7 +18,7 @@ # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License -# along with this program. If not, see <https://www.gnu.org/licenses/>. +# along with this program. If not, see <http://www.gnu.org/licenses/>. # As a special exception to the GNU General Public License, if you # distribute this file as part of a program that contains a diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/build-aux/compile new/libsirocco-2.1.0/build-aux/compile --- old/libsirocco-2.0.2/build-aux/compile 2020-02-08 11:56:17.000000000 +0100 +++ new/libsirocco-2.1.0/build-aux/compile 2017-02-28 10:35:32.000000000 +0100 @@ -1,9 +1,9 @@ #! /bin/sh # Wrapper for compilers which do not understand '-c -o'. -scriptversion=2018-03-07.03; # UTC +scriptversion=2012-10-14.11; # UTC -# Copyright (C) 1999-2018 Free Software Foundation, Inc. +# Copyright (C) 1999-2014 Free Software Foundation, Inc. # Written by Tom Tromey <[email protected]>. # # This program is free software; you can redistribute it and/or modify @@ -17,7 +17,7 @@ # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License -# along with this program. If not, see <https://www.gnu.org/licenses/>. +# along with this program. If not, see <http://www.gnu.org/licenses/>. # As a special exception to the GNU General Public License, if you # distribute this file as part of a program that contains a @@ -255,8 +255,7 @@ echo "compile $scriptversion" exit $? ;; - cl | *[/\\]cl | cl.exe | *[/\\]cl.exe | \ - icl | *[/\\]icl | icl.exe | *[/\\]icl.exe ) + cl | *[/\\]cl | cl.exe | *[/\\]cl.exe ) func_cl_wrapper "$@" # Doesn't return... ;; esac @@ -340,9 +339,9 @@ # Local Variables: # mode: shell-script # sh-indentation: 2 -# eval: (add-hook 'before-save-hook 'time-stamp) +# eval: (add-hook 'write-file-hooks 'time-stamp) # time-stamp-start: "scriptversion=" # time-stamp-format: "%:y-%02m-%02d.%02H" -# time-stamp-time-zone: "UTC0" +# time-stamp-time-zone: "UTC" # time-stamp-end: "; # UTC" # End: diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/build-aux/missing new/libsirocco-2.1.0/build-aux/missing --- old/libsirocco-2.0.2/build-aux/missing 2020-02-08 11:56:17.000000000 +0100 +++ new/libsirocco-2.1.0/build-aux/missing 2017-02-28 10:35:32.000000000 +0100 @@ -1,9 +1,9 @@ #! /bin/sh # Common wrapper for a few potentially missing GNU programs. -scriptversion=2018-03-07.03; # UTC +scriptversion=2013-10-28.13; # UTC -# Copyright (C) 1996-2018 Free Software Foundation, Inc. +# Copyright (C) 1996-2014 Free Software Foundation, Inc. # Originally written by Fran,cois Pinard <[email protected]>, 1996. # This program is free software; you can redistribute it and/or modify @@ -17,7 +17,7 @@ # GNU General Public License for more details. # You should have received a copy of the GNU General Public License -# along with this program. If not, see <https://www.gnu.org/licenses/>. +# along with this program. If not, see <http://www.gnu.org/licenses/>. # As a special exception to the GNU General Public License, if you # distribute this file as part of a program that contains a @@ -101,9 +101,9 @@ exit $st fi -perl_URL=https://www.perl.org/ -flex_URL=https://github.com/westes/flex -gnu_software_URL=https://www.gnu.org/software +perl_URL=http://www.perl.org/ +flex_URL=http://flex.sourceforge.net/ +gnu_software_URL=http://www.gnu.org/software program_details () { @@ -207,9 +207,9 @@ exit $st # Local variables: -# eval: (add-hook 'before-save-hook 'time-stamp) +# eval: (add-hook 'write-file-hooks 'time-stamp) # time-stamp-start: "scriptversion=" # time-stamp-format: "%:y-%02m-%02d.%02H" -# time-stamp-time-zone: "UTC0" +# time-stamp-time-zone: "UTC" # time-stamp-end: "; # UTC" # End: diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/build-aux/test-driver new/libsirocco-2.1.0/build-aux/test-driver --- old/libsirocco-2.0.2/build-aux/test-driver 2020-02-08 11:56:17.000000000 +0100 +++ new/libsirocco-2.1.0/build-aux/test-driver 2017-02-28 10:35:32.000000000 +0100 @@ -1,9 +1,9 @@ #! /bin/sh # test-driver - basic testsuite driver script. -scriptversion=2018-03-07.03; # UTC +scriptversion=2013-07-13.22; # UTC -# Copyright (C) 2011-2018 Free Software Foundation, Inc. +# Copyright (C) 2011-2014 Free Software Foundation, Inc. # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -16,7 +16,7 @@ # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License -# along with this program. If not, see <https://www.gnu.org/licenses/>. +# along with this program. If not, see <http://www.gnu.org/licenses/>. # As a special exception to the GNU General Public License, if you # distribute this file as part of a program that contains a @@ -140,9 +140,9 @@ # Local Variables: # mode: shell-script # sh-indentation: 2 -# eval: (add-hook 'before-save-hook 'time-stamp) +# eval: (add-hook 'write-file-hooks 'time-stamp) # time-stamp-start: "scriptversion=" # time-stamp-format: "%:y-%02m-%02d.%02H" -# time-stamp-time-zone: "UTC0" +# time-stamp-time-zone: "UTC" # time-stamp-end: "; # UTC" # End: diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/configure.ac new/libsirocco-2.1.0/configure.ac --- old/libsirocco-2.0.2/configure.ac 2020-02-16 20:37:54.000000000 +0100 +++ new/libsirocco-2.1.0/configure.ac 2021-06-30 11:54:31.000000000 +0200 @@ -1,4 +1,4 @@ -AC_INIT([libsirocco], [2.0.2], [[email protected]]) +AC_INIT([libsirocco], [2.1.0], [[email protected]]) AC_CONFIG_AUX_DIR([build-aux]) AC_CONFIG_MACRO_DIR([m4]) AM_INIT_AUTOMAKE([foreign -Wall]) diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/include/sirocco.h new/libsirocco-2.1.0/include/sirocco.h --- old/libsirocco-2.0.2/include/sirocco.h 2020-02-16 20:35:36.000000000 +0100 +++ new/libsirocco-2.1.0/include/sirocco.h 2021-06-30 11:32:29.000000000 +0200 @@ -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 diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/include/sirocco.hpp new/libsirocco-2.1.0/include/sirocco.hpp --- old/libsirocco-2.0.2/include/sirocco.hpp 2020-02-16 20:35:36.000000000 +0100 +++ new/libsirocco-2.1.0/include/sirocco.hpp 2020-11-20 15:42:22.000000000 +0100 @@ -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 @@ 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 @@ 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 diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/lib/sirocco.cpp new/libsirocco-2.1.0/lib/sirocco.cpp --- old/libsirocco-2.0.2/lib/sirocco.cpp 2020-02-16 20:35:36.000000000 +0100 +++ new/libsirocco-2.1.0/lib/sirocco.cpp 2020-11-20 15:37:59.000000000 +0100 @@ -4,6 +4,7 @@ #include <sirocco.hpp> #include <list.hpp> + using namespace std; // unsigned int *_Polynomial::combinatorial; @@ -69,6 +70,284 @@ 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 @@ _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; +} diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/m4/libtool.m4 new/libsirocco-2.1.0/m4/libtool.m4 --- old/libsirocco-2.0.2/m4/libtool.m4 2020-02-08 11:56:12.000000000 +0100 +++ new/libsirocco-2.1.0/m4/libtool.m4 2017-02-28 10:35:28.000000000 +0100 @@ -1708,11 +1708,6 @@ lt_cv_sys_max_cmd_len=8192; ;; - mint*) - # On MiNT this can take a long time and run out of memory. - lt_cv_sys_max_cmd_len=8192; - ;; - amigaos*) # On AmigaOS with pdksh, this test takes hours, literally. # So we just punt and use a minimum line length of 8192. @@ -2641,11 +2636,11 @@ version_type=darwin need_lib_prefix=no need_version=no - library_names_spec='$libname$release$versuffix$shared_ext $libname$release$major$shared_ext $libname$shared_ext' + library_names_spec='$libname$release$major$shared_ext $libname$shared_ext' soname_spec='$libname$release$major$shared_ext' shlibpath_overrides_runpath=yes shlibpath_var=DYLD_LIBRARY_PATH - shrext_cmds='`test .$module = .yes && echo .bundle || echo .dylib`' + shrext_cmds='`test .$module = .yes && echo .so || echo .dylib`' m4_if([$1], [],[ sys_lib_search_path_spec="$sys_lib_search_path_spec /usr/local/lib"]) sys_lib_dlsearch_path_spec='/usr/local/lib /lib /usr/lib' diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/sage-sirocco_interface.pyx new/libsirocco-2.1.0/sage-sirocco_interface.pyx --- old/libsirocco-2.0.2/sage-sirocco_interface.pyx 2020-02-16 20:38:32.000000000 +0100 +++ new/libsirocco-2.1.0/sage-sirocco_interface.pyx 2021-06-30 11:56:56.000000000 +0200 @@ -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,20 +69,118 @@ +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() + rop = homotopyPath (int(deg), c_values, y0R, y0I) + 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) + 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 (int(deg), c_values, y0R, y0I) + 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") @@ -89,4 +190,6 @@ 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 diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' --exclude Makefile.in --exclude configure --exclude config.guess --exclude '*.pot' --exclude mkinstalldirs --exclude aclocal.m4 --exclude config.sub --exclude depcomp --exclude install-sh --exclude ltmain.sh old/libsirocco-2.0.2/tests/sirocco_minimal_test.cpp new/libsirocco-2.1.0/tests/sirocco_minimal_test.cpp --- old/libsirocco-2.0.2/tests/sirocco_minimal_test.cpp 2020-02-08 11:55:41.000000000 +0100 +++ new/libsirocco-2.1.0/tests/sirocco_minimal_test.cpp 2020-11-20 15:35:29.000000000 +0100 @@ -16,6 +16,23 @@ }; + + +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 @@ 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);
