From 3cf6b0d1dc0a7212ef761171f400d21e5df1c315 Mon Sep 17 00:00:00 2001 From: Miguel Marco Date: Sat, 25 Feb 2017 13:37:27 +0100 Subject: [PATCH] Initial commit --- Makefile.am | 6 + README.md | 211 +++++++++++++ ZVK.py | 426 +++++++++++++++++++++++++ configure.ac | 10 + include/icomplex.hpp | 104 ++++++ include/interval.hpp | 71 +++++ include/list.hpp | 80 +++++ include/mp_complex.hpp | 86 +++++ include/mp_icomplex.hpp | 134 ++++++++ include/mp_interval.hpp | 108 +++++++ include/polynomial.hpp | 475 ++++++++++++++++++++++++++++ include/sirocco.h | 14 + include/sirocco.hpp | 73 +++++ lib/Makefile.am | 5 + lib/icomplex.cpp | 389 +++++++++++++++++++++++ lib/interval.cpp | 368 ++++++++++++++++++++++ lib/mp_complex.cpp | 389 +++++++++++++++++++++++ lib/mp_icomplex.cpp | 483 ++++++++++++++++++++++++++++ lib/mp_interval.cpp | 559 +++++++++++++++++++++++++++++++++ lib/polynomial.cpp | 297 ++++++++++++++++++ lib/sirocco.cpp | 556 ++++++++++++++++++++++++++++++++ sage-sirocco_interface.pyx | 99 ++++++ tests/Makefile.am | 11 + tests/sirocco_minimal_test.cpp | 67 ++++ tests/sirocco_test.cpp | 95 ++++++ 25 files changed, 5116 insertions(+) create mode 100644 Makefile.am create mode 100644 README.md create mode 100644 ZVK.py create mode 100644 configure.ac create mode 100644 include/icomplex.hpp create mode 100644 include/interval.hpp create mode 100644 include/list.hpp create mode 100644 include/mp_complex.hpp create mode 100644 include/mp_icomplex.hpp create mode 100644 include/mp_interval.hpp create mode 100644 include/polynomial.hpp create mode 100644 include/sirocco.h create mode 100644 include/sirocco.hpp create mode 100644 lib/Makefile.am create mode 100644 lib/icomplex.cpp create mode 100644 lib/interval.cpp create mode 100644 lib/mp_complex.cpp create mode 100644 lib/mp_icomplex.cpp create mode 100644 lib/mp_interval.cpp create mode 100644 lib/polynomial.cpp create mode 100644 lib/sirocco.cpp create mode 100644 sage-sirocco_interface.pyx create mode 100644 tests/Makefile.am create mode 100644 tests/sirocco_minimal_test.cpp create mode 100644 tests/sirocco_test.cpp diff --git a/Makefile.am b/Makefile.am new file mode 100644 index 0000000..b45e9fa --- /dev/null +++ b/Makefile.am @@ -0,0 +1,6 @@ +SUBDIRS = lib tests +ACLOCAL_AMFLAGS = -I m4 +dist-hook: + cp README.md $(distdir) + cp ZVK.py $(distdir) + cp sage-sirocco_interface.pyx $(distdir) diff --git a/README.md b/README.md new file mode 100644 index 0000000..9c832d3 --- /dev/null +++ b/README.md @@ -0,0 +1,211 @@ +# Sirocco Is a ROot Certified COntinuator + +## General description. + +This is a library for computing homotopy continuation of a given root of one +dimensional sections of bivariate complex polynomials. The output is a +piecewise linear aproximation of the path followed by the root, with the +property that there is a tubular neighborhood, with square transversal section, +that contains the actual path, and there is a three times thicker tubular +neigborhood guaranted to contain no other root of the polynomial. This second +property ensures that the piecewise linear approximation computed from all roots +of a polynomial form a topologically correct deformation of the actual braid, +since the inner tubular neighborhoods cannot intersect. + +The first version of the library was written in C, and provided only double +precission. This is the second version, rewritten in C++, and provides the +possibility of using arbitrary precission computations. Hence, it deppends on +the [MPFR](http://www.mpfr.org/) and [GMP](https://gmplib.org/) libraries. + +## Authors: + + - Miguel Marco + - Marcos Rodriguez + +## Usage: + +After installing the library, you can call it from any .cpp program by calling +the homotopyPath or homotopyPath_mp functions. In the file test/test_example you can see a simple +example + +- Double precission: + +The input of the homotopyPath function is the following: + +- An integer representing the total degree of the polynomial +- A malloc'ed array of double precission floats. Each four consecutive entries +of this list represent a coefficient, in the form of a rectangular complex +interval. These four numbers represent the following: + - the lower bound of the real part + - the upper bound of the real part + - the lower bound of the imaginary part + - the upper bound of the imaginary part + + These coefficients are written in increasing deglex order. That is, the first +coefficients are the corresponding to the monomials $1, x, y, x^2, x y, y^2, +x^3, x^2 y, x y^2, y^3...$ + +- A double precission float corresponding to the real part of an approximate +solution y0 +- A double precission float corresponding to the imaginary part of an +approximate solution y0 + +The use of complex intervals allows us to represent polynomials with non +rational coefficients (even transcendental): just give an interval that +contains it. + +y0 must be an approximation of a simple root of the univariate polynomial +f(0,y). + +For example, the input 1, [2.0, 2.0, 0.0, 0.0, 0.0, 0.0, 3.2, 3.2, -4.5, -4.5, +0.7, 0.7], 6.1, -3.3 would represent the polynomial of degree 1 f(x,y)=(2) + +(3.2*i)*x + (-4.5 + 0.7)*y , and the approximate root y0= 6.1-3.3*i + +The output is a malloc'ed array of double floats. The first entry represents +the number of points in the piecewise linear path. Then, each three entries +correspond to one of these points. The first entry of these triplets is the x +coordinate of the point (will be a real number between 0 and 1), the next two +are the real and imaginary parts of the y coordinate. + +This path represents the movement of the root y0 of f(x,y) when x varies from 0 +to 1. So for each t between 0 and 1, there is a value y_t which is the linear +interpolation of the returned values. The computation certifies that for each +t, there is a square neighborhood centered around y_t that contains one and +only one root of f(t,y). That is, there is a tube around the piecewise linear +path that contains the actual path followed by y0. Moreover, the corresponding +tubes of two different roots computed this way will not intersect. + +Note that the method might fail if there is a singularity (i.e. a value xt for +which f(xt, y) has multiple roots). In that case, you will either get a NULL +value, or an infinite running time. + +- Arbitrary precission + +The input of the homotopyPath_mp function is the following: + +- An integer representing the total degree of the polynomial +- A malloc'ed array of mpfr_t numbers. As before, each four entries +represent an (intervalar) coefficient of a polynomial, following +increasing deglex order. +- A mpfr_t number corresponding to the real part of an approximate +solution y0 +- A mpfr_t number corresponding to the imaginary part of an +approximate solution y0 +- An integer representing the precission to be used in the internal +computations (and the output). + +The output is a malloc'ed array of mpfr_t numbers as before, being the first +equal to the length of the list of points. + + +### From [SageMath](https://www.sagemath.org) + +If you have the library and headers installed correctly, you can call this +library from a [sagemath](https://www.sagemath.org) session. The relevant files +for that are `sage-sirocco_interface.pyx` and `ZVK.py`. The first one contains +the interface code, and the second one several functions to compute the +fundamental group of the complement of a curve in the complex plane, using this +library. Reffer to the comments in `ZVK.py` for a more detailed explanation of +the functions available there. + +An example of session using this library can be the following + +``` +sage: load("sage-sirocco_interface.pyx") +Compiling ./sage-sirocco_interface.pyx... +sage: load("ZVK.py") +sage: R.=QQ[] +sage: f = x^4*y - x^2*y - y^2 + y +sage: followstrand(f, 0, 2, 1) +[(0.0, 1.0, 0.0), + (0.08838834764831845, 0.9697265625, 0.0), + (0.18495073095668707, 0.8818945614497176, 0.0), + (0.2830354627446311, 0.7822431249849692, 0.0), + (0.36442245301616083, 0.7509743697977961, 0.0), + (0.4162134713480532, 0.7872237656367932, 0.0), + (0.460346631682259, 0.8708786602233393, 0.0), + (0.5011524089991131, 1.0046362459291838, 0.0), + (0.5406631526749315, 1.1979177532866128, 0.0), + (0.5801060727487398, 1.465872049076292, 0.0), + (0.6203344167290801, 1.830059589826348, 0.0), + (0.6620003579483648, 2.319953607100401, 0.0), + (0.7056383582405797, 2.975174768739414, 0.0), + (0.7517106691371234, 3.848570611753438, 0.0), + (0.8006345979664276, 5.010356224191542, 0.0), + (0.8528001455202311, 6.553629601012362, 0.0), + (0.908582197461905, 8.601688407025646, 0.0), + (0.9683494644216709, 11.317718460782247, 0.0), + (1.0, 13.0, 0.0)] +sage: fundamental_group(f) +Finitely presented group < x0, x1 | x1*x0*x1^-1*x0^-1 > +``` + +The command followstrand shows the path followed y the root y0=1 of f(x,y) when +x varies from 0 o 2. It accepts an optional argument called `prec`, for the precission +to use. By default it is set to 53 (and uses the double precission implementation). If +the computation fails, it falls back to a bigger precission until it succeeds). +The command fundamental_group uses this approach to +compute the braid monodromy associated to this polynomial, and then uses it to +compute the fundamental group of its complement, using Zariski-VanKampen method. + +At some point we expect to have this included in SageMath as an optional (or +maybe even standard) package. When that happens, we will probably provide a +more natural interface for the fundamental group and/or braid monodromy +computation. + +## Dependencies: + +The deppendencies are: + + - autotools/libtools/autoconf if you want to generate the configuration +scripts (only needed if you get the source code from a git repository) + - /make/gcc for the build process (you need them to build from a source +tarball) + - standard c libraries (math, io and std) + - mpfr + - gmp + + +## Installation: + +### From the git repository + +You can get the git repository by running + +``` +git clone https://github.com/miguelmarco/SIROCCO2.git +``` + +And then cd in the directory and run + +``` +autoreconf --install +``` + +to create the configure script then from there you can continue as usual: + +``` +./configure +make +make install +``` + +This last line might require superuser privileges for writing in system folders + +### From a tar.gz package + +If you got a tarball produced by a `make dist`command, just unpack it and run + +``` +./configure +make +make install +``` + +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 + diff --git a/ZVK.py b/ZVK.py new file mode 100644 index 0000000..09a655c --- /dev/null +++ b/ZVK.py @@ -0,0 +1,426 @@ +""" +Zariski-Van Kampen method implementation + +This file contains functions to compute the fundamental group of +the complement of a curve in the complex affine plane, using +Zariski-Van Kampen approach. It deppends on the package... + +AUTHORS: + +- Miguel Marco (2015-09-30): Initial version + + +EXAMPLES:: + + sage: R. = QQ[] + sage: f = y^3 + x^3 -1 + sage: fundamental_group(f) + Finitely presented group < x0 | > +""" + +#***************************************************************************** +# Copyright (C) 2015 Miguel Marco +# +# 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 +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# http://www.gnu.org/licenses/ +#***************************************************************************** + + + +from sage.groups.braid import BraidGroup +from sage.groups.perm_gps.permgroup_named import SymmetricGroup +from sage.rings.rational_field import QQ +from sage.rings.qqbar import QQbar +from scipy.spatial import Voronoi +from numpy import array, vstack +from sage.rings.all import CC, CIF +from sage.rings.complex_field import ComplexField +from sage.rings.complex_interval_field import ComplexIntervalField + + + +def braid_from_piecewise(strands): + """ + Compute the braid corresponding to the piecewise linear curves strands. + strands is a list whose elements are a list of pairs (x,y), where x is + a real number between 0 and 1, and y is a complex number. + + INPUT: + + - A list of lists of tuples (t, c), where t is a number bewteen 0 and 1, + and c is a complex number. + + OUTPUT: + + The braid formed by the piecewise linear strands. + + EXAMPLES:: + + sage: paths = [[(0, I), (0.2, -1 - 0.5*I), (0.8, -1), (1, -I)], [(0, -1), (0.5, -I), (1, 1)], [(0, 1), (0.5, 1 + I), (1, I)]] + sage: braid_from_piecewise(paths) + s0*s1 + + """ + L = strands + i = min([L[k][1][0] for k in range(len(L))]) + totalpoints = [[[a[0][1].real(), a[0][1].imag()]] for a in L] + indices = [1 for a in range(len(L))] + while i < 1: + for j in range(len(L)): + if L[j][indices[j]][0] > i: + xaux = L[j][indices[j] - 1][1] + yaux = L[j][indices[j]][1] + aaux = L[j][indices[j] - 1][0] + baux = L[j][indices[j]][0] + interpola = xaux + (yaux - xaux)*(i - aaux)/(baux - aaux) + totalpoints[j].append([interpola.real(), interpola.imag()]) + else: + totalpoints[j].append([L[j][indices[j]][1].real(), L[j][indices[j]][1].imag()]) + indices[j] = indices[j] + 1 + i = min([L[k][indices[k]][0] for k in range(len(L))]) + for j in range(len(L)): + totalpoints[j].append([L[j][-1][1].real(), L[j][-1][1].imag()]) + braid = [] + G = SymmetricGroup(len(totalpoints)) + for i in range(len(totalpoints[0]) - 1): + l1 = [totalpoints[j][i] for j in range(len(L))] + l2 = [totalpoints[j][i+1] for j in range(len(L))] + M = [[l1[s], l2[s]] for s in range(len(l1))] + M.sort() + l1 = [a[0] for a in M] + l2 = [a[1] for a in M] + cruces = [] + for j in range(len(l2)): + for k in range(j): + if l2[j] < l2[k]: + t = (l1[j][0] - l1[k][0])/(l2[k][0] - l1[k][0] + l1[j][0] - l2[j][0]) + s = cmp(l1[k][1]*(1 - t) + t*l2[k][1], l1[j][1]*(1 - t) + t*l2[j][1]) + cruces.append([t, k, j, -s]) + if len(cruces) > 0: + cruces.sort() + P = G(Permutation([])) + while cruces: + # we select the crosses in the same t + crucesl = [c for c in cruces if c[0]==cruces[0][0]] + crossesl = [(P(c[2]+1) - P(c[1]+1),c[1],c[2],c[3]) for c in crucesl] + cruces = cruces[len(crucesl):] + while crossesl: + crossesl.sort() + c = crossesl.pop(0) + braid.append(c[3]*min(map(P, [c[1] + 1, c[2] + 1]))) + P = G(Permutation([(c[1] + 1, c[2] + 1)]))*P + crossesl = [(P(c[2]+1) - P(c[1]+1),c[1],c[2],c[3]) for c in crossesl] + B = BraidGroup(len(L)) + return B(braid) + + +def discrim(f): + """ + Return the points in the discriminant of f. f must be a polynomial in two variables with coefficients in a + number field with a fixed embedding in QQbar. The result is the set of values of the first variable for which + two roots in the second variable coincide. + + INPUT: + + - ``f`` -- a polynomial in two variables, with coefficients in a number field with a fixed embedding in QQbar. + + OUTPUT: + + A list with the values of the discriminant in QQbar + + EXAMPLES:: + + sage: R. = QQ[] + sage: f = (y^3 + x^3 -1) * (x+y) + sage: discrim(f) + [1, + -0.500000000000000? - 0.866025403784439?*I, + -0.500000000000000? + 0.866025403784439?*I] + + """ + x, y = f.variables() + F = f.base_ring() + disc = F[x](f.discriminant(y).resultant(f, y)).roots(QQbar, multiplicities = False) + return disc + + +def segments(points): + """ + Return the bounded segments of the Voronoi diagram of the given points. + + INPUT: + + - A list of complex points. + + OUTPUT: + + A list of pairs (p1, p2) where p1 and p2 are the endpoints of the segments in the Voronoi diagram + + EXAMPLES:: + + sage: R. = QQ[] + sage: f = y^3 + x^3 -1 + sage: disc = discrim(f) + sage: segments(disc) + [(-2.84740787203333 - 2.84740787203333*I, + -2.14285714285714 + 1.11022302462516e-16*I), + (-2.84740787203333 + 2.84740787203333*I, + -2.14285714285714 + 1.11022302462516e-16*I), + (2.50000000000000 + 2.50000000000000*I, + 1.26513881334184 + 2.19128470333546*I), + (2.50000000000000 + 2.50000000000000*I, + 2.50000000000000 - 2.50000000000000*I), + (1.26513881334184 + 2.19128470333546*I, 0.000000000000000), + (0.000000000000000, 1.26513881334184 - 2.19128470333546*I), + (2.50000000000000 - 2.50000000000000*I, + 1.26513881334184 - 2.19128470333546*I), + (-2.84740787203333 + 2.84740787203333*I, + 1.26513881334184 + 2.19128470333546*I), + (-2.14285714285714 + 1.11022302462516e-16*I, 0.000000000000000), + (-2.84740787203333 - 2.84740787203333*I, + 1.26513881334184 - 2.19128470333546*I)] + + """ + discpoints = array([(CC(a).real(), CC(a).imag()) for a in points]) + added_points = 3*abs(discpoints).max() + 1.0 + configuration = vstack([discpoints, array([[added_points, 0], [-added_points, 0], [0, added_points],[0, -added_points]])]) + V = Voronoi(configuration) + res = [] + for rv in V.ridge_vertices: + if not -1 in rv: + p1 = CC(list(V.vertices[rv[0]])) + p2 = CC(list(V.vertices[rv[1]])) + res.append((p1, p2)) + return res + +def followstrand(f, x0, x1, y0a, prec=53): + """ + Return a piecewise linear aproximation of the homotopy continuation of the root y0a + from x0 to x1 + + INPUT: + + - ``f`` -- a polynomial in two variables + - ``x0`` -- a complex value, where the homotopy starts + - ``x1`` -- a complex value, where the homotopy ends + - ``y0a`` -- an approximate solution of the polynomial $F(y)=f(x_0,y)$ + - ``prec`` -- the precission to use + + OUTPUT: + + A list of values (t, ytr, yti) such that: + - ``t`` is a real number between zero and one + - $f(t\cdot x_1+(1-t)\cdot x_0, y_{tr} + I\cdot y_{ti})$ is zero (or a good enough aproximation) + - the piecewise linear path determined by the points has a tubular + neighborhood where the actual homotopy continuation path lies, and + no other root intersects it. + + EXAMPLES:: + + sage: R. = QQ[] + sage: f = x^2 + y^3 + sage: x0 = CC(1, 0) + sage: x1 = CC(1, 0.5) + sage: followstrand(f, x0, x1, -1.0) # optional + [(0.0, -1.0, 0.0), + (0.063125, -1.0001106593601545, -0.02104011456212618), + (0.12230468750000001, -1.0004151100003031, -0.04075695242737829), + (0.17778564453125, -1.0008762007617709, -0.059227299979315154), + (0.28181243896484376, -1.0021948141328698, -0.09380038023464156), + (0.3728358840942383, -1.0038270754728402, -0.123962953123039), + (0.4524813985824585, -1.005613540368227, -0.15026634875747985), + (0.5221712237596512, -1.0074443351354077, -0.17320066690539515), + (0.5831498207896948, -1.009246118007067, -0.1931978258913501), + (0.636506093190983, -1.0109719597443307, -0.21063630045261386), + (0.6831928315421101, -1.0125937110987449, -0.2258465242053158), + (0.7648946236565826, -1.0156754074572174, -0.2523480191006915), + (0.8261709677424369, -1.0181837235093538, -0.2721208327884435), + (0.8721282258068277, -1.0201720738092597, -0.2868892148703381), + (0.9410641129034139, -1.0233210275568283, -0.3089391941950436), + (1.0, -1.026166099551513, -0.3276894025360433)] + + """ + CIF = ComplexIntervalField(prec) + CC = ComplexField(prec) + G = f.change_ring(QQbar).change_ring(CIF) + (x, y) = G.variables() + g = G.subs({x: (1-x)*CIF(x0) + x*CIF(x1)}) + coefs = [] + deg = g.total_degree() + for d in range(deg + 1): + for i in range(d + 1): + c = CIF(g.coefficient({x: d-i, y: i})) + cr = c.real() + ci = c.imag() + coefs += list(cr.endpoints()) + coefs += list(ci.endpoints()) + yr = CC(y0a).real() + yi = CC(y0a).imag() + try: + if prec == 53: + points = contpath(deg, coefs, yr, yi) + else: + points = contpath_mp(deg, coefs, yr, yi, prec) + return points + except: + return followstrand(f, x0, x1, y0a, 2*prec) + +@parallel +def braid_in_segment(f, x0, x1): + """ + Return the braid formed by the y roots of f when x moves from + x0 to x1 + + INPUT: + + - ``f`` -- a polynomial in two variables + - ``x0`` -- a complex number + - ``x1`` -- a complex number + + - OUTPUT: + + A braid + + EXAMPLES:: + + sage: R. = QQ[] + sage: f = x^2 + y^3 + sage: x0 = CC(1,0) + sage: x1 = CC(1, 0.5) + sage: braid_in_segment(f, x0, x1) # optional - + s1 + + """ + CC = ComplexField(64) + (x, y) = f.variables() + I = QQbar.gen() + X0 = QQ(x0.real()) + I*QQ(x0.imag()) + X1 = QQ(x1.real()) + I*QQ(x1.imag()) + F0 = QQbar[y](f(x=X0)) + y0s = F0.roots(multiplicities=False) + strands = [followstrand(f, x0, x1, CC(a)) for a in y0s] + complexstrands = [[(a[0], CC(a[1], a[2])) for a in b] for b in strands] + centralbraid = braid_from_piecewise(complexstrands) + initialstrands = [] + y0aps = [c[0][1] for c in complexstrands] + used = [] + for y0ap in y0aps: + distances = [((y0ap - y0).norm(), y0) for y0 in y0s] + y0 = sorted(distances)[0][1] + if y0 in used: + raise ValueError("different roots are too close") + used.append(y0) + initialstrands.append([(0, CC(y0)), (1, y0ap)]) + initialbraid = braid_from_piecewise(initialstrands) + F1 = QQbar[y](f(x=X1)) + y1s = F1.roots(multiplicities=False) + finalstrands = [] + y1aps = [c[-1][1] for c in complexstrands] + used = [] + for y1ap in y1aps: + distances = [((y1ap - y1).norm(), y1) for y1 in y1s] + y1 = sorted(distances)[0][1] + if y1 in used: + raise ValueError("different roots are too close") + used.append(y1) + finalstrands.append([(0, y1ap), (1, CC(y1))]) + finallbraid = braid_from_piecewise(finalstrands) + return initialbraid*centralbraid*finallbraid + +def fundamental_group(f, simplified=True, projective = False): + """ + Return a presentation of the fundamental group of the complement of the algebraic + set defined by the polynomial f. + + INPUT: + + - ``f`` -- a polynomial in two variables, with coefficients in either the rationals + or a number field with a fixed embedding in QQbar. + + - ``simplified`` -- boolean (default: True). If set to True the presentation will be + simplified. Otherwise, the returned presentation has as many generators + as degree of the polynomial times the points in the base used to create + the segments that surround the discriminant. In this case, the generators + are granted to be meridians of the curve. + + -- ``projective`` -- boolean (default: False). If set to True, the fundamental group of + the complement of the projective completion of the curve will be computed. + Otherwise, the fundamental group of the complement in the affine plane will + be computed. + + OUTPUT: + + A presentation of the fundamental group of the complement of the curve defined by ``f``. + + EXAMPLES:: + + sage: R. = QQ[] + sage: f = x^2 + y^3 + sage: fundamental_group(f) # optional - + Finitely presented group < x0, x1 | x0^-1*x1*x0*x1*x0^-1*x1^-1 > + sage: fundamental_group(f, simplified=False) # optional + Finitely presented group < x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11 | x3*x9^-1, x4*x5*x4^-1*x10^-1, x4*x11^-1, x0*x1*x0^-1*x3^-1, x0*x4^-1, x2*x5^-1, x0*x6^-1, x2*x7^-1, x2^-1*x1*x2*x8^-1, x7*x9^-1, x7^-1*x6*x7*x10^-1, x8*x11^-1 > + + + :: + + sage: R. = QQ[] + sage: f = y^3 + x^3 + sage: fundamental_group(f) # optional + Finitely presented group < x0, x1, x2 | x1*x2*x1^-1*x0^-1*x2^-1*x0, x2*x0*x1*x0^-1*x2^-1*x1^-1 > + + It is also possible to have coefficients in a number field with a fixed embedding in QQbar:: + + sage: zeta = QQbar['x']('x^2+x+1').roots(multiplicities = False)[0] + sage: zeta + -0.50000000000000000? - 0.866025403784439?*I + sage: F = NumberField(zeta.minpoly(), 'zeta', embedding = zeta) + sage: F.inj + F.inject_variables F.injvar + sage: F.inject_variables() + Defining zeta + sage: R. = F[] + sage: f = y^3 + x^3 +zeta *x + 1 + sage: fundamental_group(f) + Finitely presented group < x0 | > + + + + + """ + (x, y) = f.variables() + F = f.base_ring() + g = f.factor().radical().prod() + d = g.degree(y) + while not g.coefficient(y**d) in F or (projective and g.total_degree() > d): + g = g.subs({x: x + y}) + d = g.degree(y) + disc = discrim(g) + segs = segments(disc) + vertices = list(set(flatten(segs))) + Faux = FreeGroup(d) + F = FreeGroup(d * len(vertices)) + rels = [] + if projective: + rels.append(prod(F.gen(i) for i in range(d))) + braidscomputed = braid_in_segment([(g, seg[0], seg[1]) for seg in segs]) + #braidscomputed = [(((g, seg[0], seg[1]), ), braid_in_segment(g,seg[0], seg[1])) for seg in segs] + for braidcomputed in braidscomputed: + seg = (braidcomputed[0][0][1], braidcomputed[0][0][2]) + b = braidcomputed[1] + i = vertices.index(seg[0]) + j = vertices.index(seg[1]) + for k in range(d): + el1 = Faux([k + 1]) * b.inverse() + el2 = k + 1 + w1 = F([sign(a)*d*i + a for a in el1.Tietze()]) + w2 = F([d*j + el2]) + rels.append(w1/w2) + G = F/rels + if simplified: + return G.simplified() + else: + return G diff --git a/configure.ac b/configure.ac new file mode 100644 index 0000000..64f87bb --- /dev/null +++ b/configure.ac @@ -0,0 +1,10 @@ +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 +LT_INIT +AC_CONFIG_FILES([Makefile lib/Makefile tests/Makefile]) +AC_OUTPUT diff --git a/include/icomplex.hpp b/include/icomplex.hpp new file mode 100644 index 0000000..9847afa --- /dev/null +++ b/include/icomplex.hpp @@ -0,0 +1,104 @@ +#ifndef _ICOMPLEX_HPP +#define _ICOMPLEX_HPP + + +#include +#include +#include + +class IComplex { +public: + + Interval r; + Interval i; + + // CONSTRUCTORS + IComplex (double = 0); + IComplex (const Interval &, const Interval &); + IComplex (double, double, double, double); + IComplex (const std::complex &); + IComplex (const IComplex &); + ~IComplex (); + + // ASSIGNMENT + IComplex & operator= (const IComplex &); + IComplex & operator= (const Interval &); + IComplex & operator= (const std::complex &); + IComplex & operator= (double); + + // COMPARISON + int operator== (const IComplex &) const; + int operator== (const std::complex &) const; + int operator== (double) const; + int operator!= (const IComplex &) const; + int operator!= (const std::complex &) const; + int operator!= (double) const; + + // ARITHMETIC + IComplex operator+ (const IComplex &) const; + IComplex operator+ (const Interval &) const; + IComplex operator+ (const std::complex &) const; + IComplex operator+ (double) const; + + IComplex operator- () const; + IComplex operator- (const IComplex &) const; + IComplex operator- (const std::complex &) const; + IComplex operator- (const Interval &) const; + IComplex operator- (double) const; + + IComplex operator* (const IComplex &) const; + IComplex operator* (const std::complex &) const; + IComplex operator* (const Interval &) const; + IComplex operator* (double) const; + + IComplex operator/ (const IComplex &) const; + IComplex operator/ (const std::complex &) const; + IComplex operator/ (const Interval &) const; + IComplex operator/ (double) const; + + + + // MISC + IComplex inverse () const; + + int containsZero () const; + int contains (double) const; + int contains (const std::complex &); + int subset (const IComplex &) const; + int subsetEq (const IComplex &) const; + + + +}; + + +// NON MEMBER FUNCTIONS + +std::ostream & operator<< (std::ostream &, const IComplex &); + + +IComplex operator+ (const Interval &, const IComplex &); +IComplex operator+ (const std::complex &, const IComplex &); +IComplex operator+ (double, const IComplex &); + +IComplex operator- (const Interval &, const IComplex &); +IComplex operator- (const std::complex &, const IComplex &); +IComplex operator- (double, const IComplex &); + +IComplex operator* (const Interval &, const IComplex &); +IComplex operator* (const std::complex &, const IComplex &); +IComplex operator* (double, const IComplex &); + +IComplex operator/ (const Interval &, const IComplex &); +IComplex operator/ (const std::complex &, const IComplex &); +IComplex operator/ (double, const IComplex &); + + + +// MISCELLANEOUS FUNCTIONS FOR INVERT COMPLEX INTERVALS +IComplex inverse_horizontal (double, double, double); +IComplex inverse_vertical (double, double, double); +double f (double, double); +double g (double, double); + +#endif diff --git a/include/interval.hpp b/include/interval.hpp new file mode 100644 index 0000000..d3b7624 --- /dev/null +++ b/include/interval.hpp @@ -0,0 +1,71 @@ +#ifndef _INTERVAL_HPP +#define _INTERVAL_HPP + +#include + + +// IMPLEMENTS [a,b], +// a,b double + +class Interval { + // FRIEND CLASS + +public: + + double a; + double b; + + + // CONSTRUCTORS + // DEFAULT CONSTRUCTOR + Interval (double = 0); + // CONSTRUCTOR BY ELEMENTS + Interval (double, double); + // COPY CONSTRUCTOR + Interval (const Interval &); + // DESTRUCTOR + ~Interval (); + + // ASSIGNMENT + Interval & operator= (const Interval &); + Interval & operator= (double); + + // COMPARISON + int operator== (const Interval &) const; + int operator!= (const Interval &) const; + + // ARITHMETIC + Interval operator+ (const Interval &) const; + Interval operator+ (double) const; + + Interval operator- () const; + Interval operator- (const Interval &) const; + Interval operator- (double) const; + + Interval operator* (const Interval &) const; + Interval operator* (double) const; + + Interval operator/ (const Interval &) const; + Interval operator/ (double) const; + + + + // MISC + Interval inverse () const; + + int containsZero () const; + int contains (double) const; + int subset (const Interval &) const; + int subsetEq (const Interval &) const; + + +}; + +// EXTERNAL FUNCTIONS +std::ostream & operator<< (std::ostream &, const Interval &); + +Interval operator+ (double, const Interval &); +Interval operator- (double, const Interval &); +Interval operator* (double, const Interval &); +Interval operator/ (double, const Interval &); +#endif \ No newline at end of file diff --git a/include/list.hpp b/include/list.hpp new file mode 100644 index 0000000..4fbbe98 --- /dev/null +++ b/include/list.hpp @@ -0,0 +1,80 @@ +#ifndef _LIST_HPP +#define _LIST_HPP + +#include +#include +#include + + +template +class Cell { + T x[3]; + Cell () {} + Cell (T data[3]) {} + Cell (T x0, T y0R, T y0I) {} + ~Cell () {} +}; + +template <> +class Cell { +public: + double x[3]; + Cell () {} + Cell (double data[3]) {x[0] = data[0]; x[1] = data[1], x[2] = data[2];} + Cell (double x0, double y0R, double y0I) {x[0] = x0; x[1] = y0R; x[2] = y0I;} + ~Cell () {} +}; +template <> +class Cell { +public: + mpfr_t x[3]; + Cell () {} + Cell (mpfr_t x0, mpfr_t y0R, mpfr_t y0I) { + mpfr_inits (x[0], x[1], x[2], NULL); + mpfr_set (x[0], x0, MPFR_RNDN); + mpfr_set (x[1], y0R, MPFR_RNDN); + mpfr_set (x[2], y0I, MPFR_RNDN); + } + + Cell (const Cell &op) { + mpfr_inits (x[0], x[1], x[2], NULL); + mpfr_set (x[0], op.x[0], MPFR_RNDN); + mpfr_set (x[1], op.x[1], MPFR_RNDN); + mpfr_set (x[2], op.x[2], MPFR_RNDN); + } + + ~Cell () { + mpfr_clears (x[0], x[1], x[2], NULL); + } +}; + +std::ostream & operator<< (std::ostream &output, const Cell &op) { + output << op.x[0] << " "; + output << op.x[1] << " "; + output << op.x[2]; + + return output; +} + +std::ostream & operator<< (std::ostream &output, const Cell &op) { + output << op.x[0] << " "; + output << op.x[1] << " "; + output << op.x[2]; + + return output; +} + +std::ostream & operator<< (std::ostream &output, const std::vector > &op) { + for (int i=0; i > &op) { + for (int i=0; i +#include +#include + + + +class MPComplex { +public: + + mpfr_t r; + mpfr_t i; + + MPComplex (double = 0); + MPComplex (double, double); + MPComplex (const std::complex &); + MPComplex (const mpfr_t); + MPComplex (const mpfr_t, const mpfr_t); + ~MPComplex (); + + // ASSIGNMENT + MPComplex & operator= (double); + MPComplex & operator= (const std::complex &); + MPComplex & operator= (const mpfr_t ); + MPComplex & operator= (const MPComplex &); + + // COMPARISON + int operator == (double) const; + int operator == (const std::complex &) const; + int operator == (const MPComplex &) const; + int operator != (double) const; + int operator != (const std::complex &) const; + int operator != (const MPComplex &) const; + + // ARITHMETIC + MPComplex operator+ (double) const; + MPComplex operator+ (const std::complex &) const; + MPComplex operator+ (const mpfr_t ) const; + MPComplex operator+ (const MPComplex &) const; + + MPComplex operator- () const; + MPComplex operator- (double) const; + MPComplex operator- (const std::complex &) const; + MPComplex operator- (const mpfr_t ) const; + MPComplex operator- (const MPComplex &) const; + + MPComplex operator* (double) const; + MPComplex operator* (const std::complex &) const; + MPComplex operator* (const mpfr_t ) const; + MPComplex operator* (const MPComplex &) const; + + MPComplex operator/ (double) const; + MPComplex operator/ (const std::complex &) const; + MPComplex operator/ (const mpfr_t ) const; + MPComplex operator/ (const MPComplex &) const; + + +}; + +// NON MEMBER FUNCTIONS + +std::ostream & operator<< (std::ostream &, const MPComplex &); + +MPComplex operator+ (double, const MPComplex &); +MPComplex operator+ (const std::complex &, const MPComplex &); +MPComplex operator+ (const mpfr_t , const MPComplex &); + +MPComplex operator- (double, const MPComplex &); +MPComplex operator- (const std::complex &, const MPComplex &); +MPComplex operator- (const mpfr_t , const MPComplex &); + +MPComplex operator* (double, const MPComplex &); +MPComplex operator* (const std::complex &, const MPComplex &); +MPComplex operator* (const mpfr_t , const MPComplex &); + +MPComplex operator/ (double, const MPComplex &); +MPComplex operator/ (const std::complex &, const MPComplex &); +MPComplex operator/ (const mpfr_t , const MPComplex &); + + +void abs (mpfr_t rop, const MPComplex &op); + + +#endif diff --git a/include/mp_icomplex.hpp b/include/mp_icomplex.hpp new file mode 100644 index 0000000..35cccdf --- /dev/null +++ b/include/mp_icomplex.hpp @@ -0,0 +1,134 @@ +#ifndef _MP_ICOMPLEX_HPP +#define _MP_ICOMPLEX_HPP + + +#include +#include +#include +#include +#include +#include + +class MPIComplex { +public: + + + // CONSTRUCTORS + MPIComplex (double = 0); + MPIComplex (double, double, double, double); + MPIComplex (const mpfr_t); + MPIComplex (const MPComplex &); + MPIComplex (const mpfr_t, const mpfr_t, const mpfr_t, const mpfr_t); + MPIComplex (const Interval &, const Interval &); + MPIComplex (const MPInterval &, const MPInterval &); + MPIComplex (const IComplex &); + MPIComplex (const MPIComplex &); + ~MPIComplex (); + + // ASSIGNMENT + MPIComplex & operator= (const MPIComplex &); + MPIComplex & operator= (const IComplex &); + MPIComplex & operator= (const MPInterval &); + MPIComplex & operator= (const Interval &); + MPIComplex & operator= (const mpfr_t); + MPIComplex & operator= (double); + + // COMPARISON + int operator== (double) const; + int operator== (const std::complex &) const; + int operator== (const Interval &) const; + int operator== (const MPInterval &) const; + int operator== (const MPIComplex &) const; + int operator!= (double) const; + int operator!= (const std::complex &) const; + int operator!= (const Interval &) const; + int operator!= (const MPInterval &) const; + int operator!= (const MPIComplex &) const; + + // ARITHMETIC + MPIComplex operator+ (const MPIComplex &) const; + MPIComplex operator+ (const IComplex &) const; + MPIComplex operator+ (const MPInterval &) const; + MPIComplex operator+ (const Interval &) const; + MPIComplex operator+ (const mpfr_t) const; + MPIComplex operator+ (double) const; + + MPIComplex operator- () const; + MPIComplex operator- (const MPIComplex &) const; + MPIComplex operator- (const IComplex &) const; + MPIComplex operator- (const MPInterval &) const; + MPIComplex operator- (const Interval &) const; + MPIComplex operator- (const mpfr_t) const; + MPIComplex operator- (double) const; + + MPIComplex operator* (const MPIComplex &) const; + MPIComplex operator* (const IComplex &) const; + MPIComplex operator* (const MPInterval &) const; + MPIComplex operator* (const Interval &) const; + MPIComplex operator* (const mpfr_t) const; + MPIComplex operator* (double) const; + + + MPIComplex operator/ (const MPIComplex &) const; + MPIComplex operator/ (const IComplex &) const; + MPIComplex operator/ (const MPInterval &) const; + MPIComplex operator/ (const Interval &) const; + MPIComplex operator/ (const mpfr_t) const; + MPIComplex operator/ (double) const; + + + + // MISC + MPIComplex inverse () const; + + int containsZero () const; + int contains (double) const; + int contains (const mpfr_t) const; + int subset (const MPIComplex &) const; + // int subsetEq (const MPIComplex &) const; + + + MPInterval r; + MPInterval i; +}; + + +// NON MEMBER FUNCTIONS + +std::ostream & operator<< (std::ostream &, const MPIComplex &); + + +MPIComplex operator+ (const IComplex &, const MPIComplex &); +MPIComplex operator+ (const MPInterval &, const MPIComplex &); +MPIComplex operator+ (const Interval &, const MPIComplex &); +MPIComplex operator+ (const mpfr_t, const MPIComplex &); +MPIComplex operator+ (double, const MPIComplex &); + +MPIComplex operator- (const IComplex &, const MPIComplex &); +MPIComplex operator- (const MPInterval &, const MPIComplex &); +MPIComplex operator- (const Interval &, const MPIComplex &); +MPIComplex operator- (const mpfr_t, const MPIComplex &); +MPIComplex operator- (double, const MPIComplex &); + +MPIComplex operator* (const IComplex &, const MPIComplex &); +MPIComplex operator* (const MPInterval &, const MPIComplex &); +MPIComplex operator* (const Interval &, const MPIComplex &); +MPIComplex operator* (const mpfr_t, const MPIComplex &); +MPIComplex operator* (double, const MPIComplex &); + +MPIComplex operator/ (const IComplex &, const MPIComplex &); +MPIComplex operator/ (const MPInterval &, const MPIComplex &); +MPIComplex operator/ (const Interval &, const MPIComplex &); +MPIComplex operator/ (const mpfr_t, const MPIComplex &); +MPIComplex operator/ (double, const MPIComplex &); + + + + +// MISCELLANEOUS FUNCTIONS FOR INVERT COMPLEX INTERVALS +MPIComplex inverse_horizontal (const mpfr_t, const mpfr_t, const mpfr_t); +MPIComplex inverse_vertical (const mpfr_t, const mpfr_t, const mpfr_t); +void f (mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t); +void g (mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t); + +#endif diff --git a/include/mp_interval.hpp b/include/mp_interval.hpp new file mode 100644 index 0000000..ff3ba27 --- /dev/null +++ b/include/mp_interval.hpp @@ -0,0 +1,108 @@ +#ifndef _MP_INTERVAL_HPP +#define _MP_INTERVAL_HPP + +#include +#include +#include + + +// IMPLEMENTS [a,b], +// a,b mpfr + +class MPInterval { + // FRIEND CLASS +public: + + mpfr_t a; + mpfr_t b; + + + // CONSTRUCTORS + // DEFAULT CONSTRUCTOR + MPInterval (double = 0); + // CONSTRUCTOR BY double ELEMENTS + MPInterval (double, double); + // CONSTRUCTOR BY string VALUE + MPInterval (const char *); + // CONSTRUCTOR BY string PAIR OF VALUES + MPInterval (const char *, const char *); + // CONSTRUCTOR BY mpfr ELEMENT + MPInterval (const mpfr_t); + // CONSTRUCTOR BY mpfr PAIR OF ELEMENTS + MPInterval (const mpfr_t, const mpfr_t); + // CONSTRUCTOR BY double Invterval + MPInterval (const Interval &); + // COPY CONSTRUCTOR + MPInterval (const MPInterval &); + // DESTRUCTOR + ~MPInterval (); + + // COMPARISON + int operator== (double) const; + int operator== (const Interval &) const; + int operator== (const MPInterval &) const; + int operator!= (double) const; + int operator!= (const Interval &) const; + int operator!= (const MPInterval &) const; + + // ASSIGNMENT + MPInterval & operator= (const MPInterval &); + MPInterval & operator= (const Interval &); + MPInterval & operator= (const mpfr_t); + MPInterval & operator= (double); + + // ARITHMETIC + MPInterval operator+ (const MPInterval &) const; + MPInterval operator+ (const Interval &) const; + MPInterval operator+ (const mpfr_t) const; + MPInterval operator+ (double) const; + + MPInterval operator- () const; + MPInterval operator- (const MPInterval &) const; + MPInterval operator- (const Interval &) const; + MPInterval operator- (const mpfr_t) const; + MPInterval operator- (double) const; + + MPInterval operator* (const MPInterval &) const; + MPInterval operator* (const Interval &) const; + MPInterval operator* (const mpfr_t) const; + MPInterval operator* (double) const; + + MPInterval operator/ (const MPInterval &) const; + MPInterval operator/ (const Interval &) const; + MPInterval operator/ (const mpfr_t) const; + MPInterval operator/ (double) const; + + + + // MISC + MPInterval inverse () const; + + int containsZero () const; + int contains (double) const; + int contains (const mpfr_t) const; + int subset (const MPInterval &) const; + int subset (const Interval &) const; + + +}; + +// EXTERNAL FUNCTIONS +std::ostream & operator<< (std::ostream &, const MPInterval &); + +MPInterval operator+ (double, const MPInterval &); +MPInterval operator+ (const mpfr_t, const MPInterval &); +MPInterval operator+ (const Interval &, const MPInterval &); + +MPInterval operator- (double, const MPInterval &); +MPInterval operator- (const mpfr_t, const MPInterval &); +MPInterval operator- (const Interval &, const MPInterval &); + +MPInterval operator* (double, const MPInterval &); +MPInterval operator* (const mpfr_t, const MPInterval &); +MPInterval operator* (const Interval &, const MPInterval &); + +MPInterval operator/ (double, const MPInterval &); +MPInterval operator/ (const mpfr_t, const MPInterval &); +MPInterval operator/ (const Interval &, const MPInterval &); +#endif diff --git a/include/polynomial.hpp b/include/polynomial.hpp new file mode 100644 index 0000000..11202fc --- /dev/null +++ b/include/polynomial.hpp @@ -0,0 +1,475 @@ +#ifndef _POLYNOMIAL_HPP +#define _POLYNOMIAL_HPP + + +#include +#include // std::complex +#include +#include +#include + + +#define MIN(a,b) (((a) < (b))? (a) : (b)) +#define MAX(a,b) (((a) > (b))? (a) : (b)) + +class _Polynomial { +public: + _Polynomial() {} + ~_Polynomial() {} + + static unsigned int *combinatorial; + static void startCombinatorialNumbers (int degree); + static void releaseCombinatorialNumbers () { + delete [] _Polynomial::combinatorial; + } +}; + + +template +class Polynomial : _Polynomial { +public: + + // CONSTRUCTOR + Polynomial (int = 0); + Polynomial (int, const T *); + Polynomial (const Polynomial &); + ~Polynomial () {delete [] coef;} + + // ASSIGNMENT + Polynomial & operator= (const Polynomial &); + Polynomial & operator= (const T &); + Polynomial & operator= (double); + + T operator() (const T &, const T &) const; + T diffX (const T &, const T &) const; + T diffY (const T &, const T &) const; + T diffXX (const T &, const T &) const; + T diffXY (const T &, const T &) const; + T diffYY (const T &, const T &) const; + + + // MISC FUNCTIONS FOR HORNER + T evalPolClassic (const T &, const T &) const; + T evalIPolHornerXY (const T &, const T &) const {return T();} + T evalIPolHornerYX (const T &, const T &) const {return T();} + + T evalPolYClassic (const T &, const T &) const; + T evalIPolYHornerXY (const T &, const T &) const {return T();} + T evalIPolYHornerYX (const T &, const T &) const {return T();} + + + // CHANGE OF VARIABLE + Polynomial IlinearVarChange (const T &a) const; + Polynomial IlinearVarChange2 (const T &x0) const; + + // MEMBERS + T *coef; + int degree; + + +}; + +template +std::ostream & operator<< (std::ostream &output, const Polynomial &op) { + int allZero = 1; + for (int i=0; i<=op.degree; ++i) + 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) + output << "-"; + else if (op.coef[i*(i+1)/2 + j] != 1.0) + output << "+(" << op.coef[i*(i+1)/2 + j] << ")*"; + else + output << "+"; + if (i == 0) + if (op.coef[0] == T(1)) + output << "1"; + if (i != 0) + if (i-j == 0) + if (j == 1) + output << "y"; + else + output << "y^" << j; + else if (i-j == 1) + if (j == 0) + output << "x"; + else if (j == 1) + output << "x*y"; + else + output << "x*y^" << j; + 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; +} + +// CONSTRUCTOR + +// DEFAULT CONSTRUCTOR FOR ZERO POLYNOMIAL +template +Polynomial::Polynomial (int op) { + degree = op; + coef = new T [((op+1)*(op+2))/2]; + for (int i=0; i<((op+1)*(op+2))/2; ++i) + coef[i] = 0.0; +} + +// CONSTRUCTOR FROM ARRAY OF COMPLEX or ICOMPLEX +template +Polynomial::Polynomial (int op1, const T *op2) { + degree = op1; + int nCoef = ((degree+1)*(degree+2)) / 2; + coef = new T [nCoef]; + for (int i=0; i +Polynomial::Polynomial (const Polynomial &op) { + degree = op.degree; + int nCoef = ((degree+1)*(degree+2)) / 2; + coef = new T [nCoef]; + for (int i=0; i +Polynomial & Polynomial::operator= (const Polynomial &op) { + // CHECKS WHETHER ASSIGNED POL IS TARGET POL + if (this == &op) return (*this); + + // GETS FINAL NUMBER OF COEFICIENTS + int nCoef = ((op.degree+1)*(op.degree+2)) / 2; + // IF NUMBER OF COEFFICIENTS ARE THE SAME, NO NEED TO REALLOCATE + if (this->degree != op.degree) { + delete [] this->coef; + this->coef = new T[nCoef]; + this->degree = op.degree; + } + + // UPDATE VALUES OF COEFICIENTS + for (int i=0; icoef[i] = op.coef[i]; + + + return (*this); + +} + +template +Polynomial & Polynomial::operator= (const T &op) { + // CHECKS IF TARGET POL IS ALSO A CONSTANT POL + // IN THAT CASE THERE IS NO NEED TO REALLOCATE MEMORY + if (this->degree != 0) { + delete [] this->coef; + this->coef = new T[1]; + this->degree = 0; + } + + // UPDATE INDEPENDENT TERM + this->coef[0] = op; + + return (*this); +} + +template +Polynomial & Polynomial::operator= (double op) { + // USE PREVIOUS DEFINITION + return ((*this) = T(op)); +} + + + + +// EVALUATION + + +template +T Polynomial::evalPolClassic (const T &x, const T &y) const { + // coef[(i*(i+1))/2 + j] is coeficient of monomial of degree 'i', + // with degree of second variable 'j' -> coef of x^(i-j) y^j + // i; // counter along degree of polynomial + // j; // counter for monomials of equal degree + + // sets output to 0+0i + T rop; + + // compute all powers of x and y (from 0 to degree) + T xPower[this->degree+1], yPower[this->degree+1]; + xPower[0] = yPower[0] = 1.0; + xPower[1] = x; yPower[1] = y; + for (int j=1; jdegree; ++j) { + xPower[j+1] = xPower[j] * x; + yPower[j+1] = yPower[j] * y; + } + + // loop for each monomial of each degree + for (int i=0; i<=this->degree; ++i) + for (int j=0; j<=i; ++j) + rop = rop + + ((this->coef)[(i*(i+1))/2 + j]) * xPower[i-j] * yPower[j]; + + + return rop; +} + + +template +T Polynomial::operator() (const T &x, const T &y) const { + return this->evalPolClassic (x,y); +} + +template +T Polynomial::diffX (const T &x, const T &y) const { + // coef[(i*(i+1))/2 + j] is coeficient of monomial of degree 'i', + // with degree of second variable 'j' -> coef of x^(i-j) y^j + // i; // counter along degree of polynomial + // j; // counter for monomials of equal degree + + // sets output to 0+0i + T rop; + + // compute all powers of x and y + // (from 0 to degree-1) (from 0 to degree) resp + T xPower[this->degree], yPower[this->degree+1]; + xPower[0] = yPower[0] = 1.0; + xPower[1] = x; yPower[1] = y; + for (int j=1; jdegree-1; ++j) { + xPower[j+1] = xPower[j] * x; + yPower[j+1] = yPower[j] * y; + } + yPower[this->degree] = yPower[this->degree-1] * y; + + // loop for each monomial of each degree + for (int i=1; i<=this->degree; ++i) + for (int j=0; jcoef)[(i*(i+1))/2 + j]) * + xPower[i-j-1] * yPower[j]; + + return rop; + +} + +template +T Polynomial::evalPolYClassic (const T &x, const T &y) const { + // coef[(i*(i+1))/2 + j] is coeficient of monomial of degree 'i', + // with degree of second variable 'j' -> coef of x^(i-j) y^j + // i; // counter along degree of polynomial + // j; // counter for monomials of equal degree + + // sets output to 0+0i + T rop; + + // compute all powers of x and y + // (from 0 to degree) (from 0 to degree-1) resp + T xPower[this->degree], yPower[this->degree]; + xPower[0] = yPower[0] = 1.0; + xPower[1] = x; yPower[1] = y; + for (int j=1; jdegree-1; ++j) { + xPower[j+1] = xPower[j] * x; + yPower[j+1] = yPower[j] * y; + } + + // loop for each monomial of each degree + for (int i=1; i<=this->degree; ++i) + for (int j=1; j<=i; ++j) + rop = rop + + T(j)*((this->coef)[(i*(i+1))/2 + j]) * + xPower[i-j] * yPower[j-1]; + + return rop; + +} + + +template +T Polynomial::diffY (const T &x, const T &y) const { + return this->evalPolYClassic (x,y); +} + + +template +T Polynomial::diffXX (const T &x, const T &y) const { + // coef[(i*(i+1))/2 + j] is coeficient of monomial of degree 'i', + // with degree of second variable 'j' -> coef of x^(i-j) y^j + // i; // counter along degree of polynomial + // j; // counter for monomials of equal degree + + // sets output to 0+0i + T rop; + + // compute all powers of x and y + // (from 0 to degree-2) (from 0 to degree-2) resp + T xPower[this->degree-1], yPower[this->degree-1]; + xPower[0] = yPower[0] = 1.0; + xPower[1] = x; yPower[1] = y; + for (int j=1; jdegree-2; ++j) { + xPower[j+1] = xPower[j] * x; + yPower[j+1] = yPower[j] * y; + } + + // loop for each monomial of each degree + for (int i=2; i<=this->degree; ++i) + for (int j=0; jcoef)[(i*(i+1))/2 + j]) * + xPower[i-j-2] * yPower[j]; + + return rop; + +} + +template +T Polynomial::diffXY (const T &x, const T &y) const { + // coef[(i*(i+1))/2 + j] is coeficient of monomial of degree 'i', + // with degree of second variable 'j' -> coef of x^(i-j) y^j + // i; // counter along degree of polynomial + // j; // counter for monomials of equal degree + + // sets output to 0+0i + T rop; + + // compute all powers of x and y + // (from 0 to degree-2) (from 0 to degree-2) resp + T xPower[this->degree-1], yPower[this->degree-1]; + xPower[0] = yPower[0] = 1.0; + xPower[1] = x; yPower[1] = y; + for (int j=1; jdegree-2; ++j) { + xPower[j+1] = xPower[j] * x; + yPower[j+1] = yPower[j] * y; + } + + // loop for each monomial of each degree + for (int i=2; i<=this->degree; ++i) + for (int j=1; jcoef)[(i*(i+1))/2 + j]) * + xPower[i-j-1] * yPower[j-1]; + + return rop; + + +} + +template +T Polynomial::diffYY (const T &x, const T &y) const { + // coef[(i*(i+1))/2 + j] is coeficient of monomial of degree 'i', + // with degree of second variable 'j' -> coef of x^(i-j) y^j + // i; // counter along degree of polynomial + // j; // counter for monomials of equal degree + + // sets output to 0+0i + T rop; + + // compute all powers of x and y + // (from 0 to degree-2) (from 0 to degree-2) resp + T xPower[this->degree-1]; + T yPower[this->degree-1]; + xPower[0] = yPower[0] = 1.0; + xPower[1] = x; yPower[1] = y; + for (int j=1; jdegree-2; ++j) { + xPower[j+1] = xPower[j] * x; + yPower[j+1] = yPower[j] * y; + } + + // loop for each monomial of each degree + for (int i=2; i<=this->degree; ++i) + for (int j=2; j<=i; ++j) + rop = rop + + T(j)*T(j-1)*((this->coef)[(i*(i+1))/2 + j]) * + xPower[i-j] * yPower[j-2]; + + return rop; +} + + + + +// CHANGE OF VARIABLE + +// RETURNS q = p(x, y + ax); +// q = sum D(i,j)x^(i-j)y^j +template +Polynomial Polynomial::IlinearVarChange (const T &a) const { + + + int i; // counter for degree + int j; // counter for exponent of y + int k; // counter for elements modifying x^(i-j) y^j + + // for each degree i, only coefficients of monomials with j in {0, ..., i-1} are changed + // for each j, only monomials with degree i and k in {j+1, ..., i} modify it + // for each k, coef[(i*(i+1))/2+k] x^(i-k) (y+ax)^k + // from binomial expansion, the term with x^(i-j) y^j has multiplier: + // coef[(i*(i+1))/2+k] a^(k-j) (k over j) + T aPower[this->degree+1]; + aPower[0] = 1; + aPower[1] = a; + for (i=1; idegree; ++i) + aPower[i+1] = aPower[i] * a; + + unsigned int *combinatorial = _Polynomial::combinatorial; + + Polynomial rop(this->degree); + + // rop = sum_{i=0}^{degree} sum_{j=0}^{i} D(i,j)x^{i-j}y^j + // D(i,j) = sum_{k=j}^{i} C(i,k)(k over j) a^{k-j} + for (i=0; i<=this->degree; ++i) for (j=0; j<=i; ++j) { + rop.coef[(i*(i+1))/2 + j] = this->coef[(i*(i+1))/2 + j]; + for (k=j+1; k<=i; ++k) + rop.coef[(i*(i+1))/2 + j] = rop.coef[(i*(i+1))/2 + j] + + combinatorial[(k*(k+1))/2+j] * aPower[k-j] * this->coef[(i*(i+1))/2+k]; + } + + return rop; + +} +// RETURNS q = p(x+x0, y); +template +Polynomial Polynomial::IlinearVarChange2 (const T &x0) const { + int i; // counter for degree + int j; // counter for exponent of y + int k; // counter for elements modifying x^(i-j) y^j + + // for each degree i, only coefficients of monomials with j in {0, ..., i-1} are changed + // for each j, only monomials with degree i and k in {j+1, ..., i} modify it + // for each k, coef[(i*(i+1))/2+k] x^(i-k) (y+ax)^k + // from binomial expansion, the term with x^(i-j) y^j has multiplier: + // coef[(i*(i+1))/2+k] a^(k-j) (k over j) + T x0Power[this->degree+1]; + x0Power[0] = 1; + x0Power[1] = x0; + for (i=1; idegree; ++i) + x0Power[i+1] = x0Power[i] * x0; + + unsigned int *combinatorial = _Polynomial::combinatorial; + + Polynomial rop(this->degree); + + // rop = sum_{i=0}^{degree} sum_{j=0}^{i} D(i,j)x^{i-j}y^j + // D(i,j) = sum_{k=i}^{degree} (k-j over i-j) C(k,j) x0^(k-i) + for (i=0; i<=this->degree; ++i) for (j=0; j<=i; ++j) { + rop.coef[(i*(i+1))/2+j] = this->coef[(i*(i+1))/2+j]; + for (k=i+1; k<=this->degree; ++k) + rop.coef[(i*(i+1))/2+j] = rop.coef[(i*(i+1))/2+j] + + combinatorial[((k-j)*(k-j+1))/2 + i-j] * this->coef[(k*(k+1))/2+j] * x0Power[k-i]; + } + + return rop; + +} + + +#endif diff --git a/include/sirocco.h b/include/sirocco.h new file mode 100644 index 0000000..53f74bd --- /dev/null +++ b/include/sirocco.h @@ -0,0 +1,14 @@ +#ifndef __SIROCCO_H__ +#define __SIROCCO_H__ + +#include +#include + +// ACTIVATE ROUND MODE CONTROL +#pragma STDC FENV_ACCESS on + + +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); + +#endif diff --git a/include/sirocco.hpp b/include/sirocco.hpp new file mode 100644 index 0000000..536a054 --- /dev/null +++ b/include/sirocco.hpp @@ -0,0 +1,73 @@ +#ifndef _SIROCCO_HPP_ +#define _SIROCCO_HPP_ + +#include +#include +#include +#include + + +// ACTIVATE ROUND MODE CONTROL +#pragma STDC FENV_ACCESS on + +typedef std::complex 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); + + +template // T -> Complex, IComplex, MPComplex or MPIComplex +void correctRoot (const Polynomial, const T &x0, T &y0) {} + + +template // T -> Complex or MPComplex +T getA (const Polynomial f, const T &x0, const T &y0) { + return -(f.diffX (x0, y0) / f.diffY (x0, y0)); +} + + +template // T -> IComplex or MPIComplex +int newtonTest (const Polynomial &If, const T &Ix0, const T &Iy0, const T &IY0) { + + T N, aux; + + N = If(Ix0, Iy0); + + + aux = If.diffY (Ix0, IY0); + + // check we do not perform division by an interval containing 0. + // in this case, return fail the test + if (aux.containsZero ()) return 0; + + N = Iy0 - (N / aux); + + + + if (N.subset (IY0)) + return 1; + return 0; +} + + +// // RETURNS WHETHER NEWTON EVALUATION N IS CONTAINED IN THE INTERVAL ENCLOSURE IY0 OF Iy0. +// // BESIDES REPEATS WITH A ENCLOSURE THREE TIMES BIGER. RETURNS SUCCESS IF BOTH ARE SATISFIED +template // T -> IComplex or MPIComplex +int validate (const Polynomial &If, const T &Ix0, const T &Iy0, const T &IY0) { + + + T IY1; + if (!newtonTest (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 (If, Ix0, Iy0, IY1)) return 0; + + + return 1; +} +#endif diff --git a/lib/Makefile.am b/lib/Makefile.am new file mode 100644 index 0000000..ca6f06e --- /dev/null +++ b/lib/Makefile.am @@ -0,0 +1,5 @@ +lib_LTLIBRARIES = libsirocco.la +libsirocco_la_SOURCES = icomplex.cpp interval.cpp $(top_srcdir)/include/list.hpp $(top_srcdir)/include/mp_complex.hpp $(top_srcdir)/include/mp_icomplex.hpp $(top_srcdir)/include/mp_interval.hpp $(top_srcdir)/include/polynomial.hpp $(top_srcdir)/include/sirocco.hpp $(top_srcdir)/include/icomplex.hpp $(top_srcdir)/include/interval.hpp mp_complex.cpp mp_icomplex.cpp mp_interval.cpp polynomial.cpp sirocco.cpp +libsirocco_la_LDFLAGS = -lm -lmpfr -lgmp +AM_CPPFLAGS = -I$(top_srcdir)/include +include_HEADERS = $(top_srcdir)/include/sirocco.h diff --git a/lib/icomplex.cpp b/lib/icomplex.cpp new file mode 100644 index 0000000..9ca65ef --- /dev/null +++ b/lib/icomplex.cpp @@ -0,0 +1,389 @@ +#include +#include +#include "icomplex.hpp" + +#define MIN(a,b) (((a) < (b))? (a) : (b)) +#define MAX(a,b) (((a) > (b))? (a) : (b)) + +#pragma STDC FENV_ACCESS on + + +// CONSTRUCTORS +IComplex::IComplex (double op) + : r(op, op), i() {} +IComplex::IComplex (double op1, double op2, double op3, double op4) + : r(op1, op2), i(op3, op4) {} +IComplex::IComplex (const Interval &op1, const Interval &op2) + : r(op1), i(op2) {} + +IComplex::IComplex (const std::complex &op) + : r(op.real()), i(op.imag()) {} + // COPY CONSTRUCTOR +IComplex::IComplex (const IComplex &op) + : r(op.r), i(op.i) {} +IComplex::~IComplex () {} + + + +// // ASSIGNMENT +IComplex & IComplex::operator= (const IComplex &op) { + // CHECK WHETER WE ARE ASIGNING SOMETHING OVER ITSELF + if (this == &op) return (*this); + // UPDATE ELEMENTS + this->r = op.r; this->i = op.i; + return (*this); +} +IComplex & IComplex::operator= (const std::complex &op) { + // UPDATE ELEMENTS + this->r = op.real (); this->i = op.imag (); + return (*this); +} + +IComplex & IComplex::operator= (const Interval &op) { + // UPDATE ELEMENTS + this->r = op; this->i = 0.0; + return (*this); +} +IComplex & IComplex::operator= (double op) { + // UPDATE ELEMENTS + this->r = op; this->i = 0.0; + return (*this); +} + +// // COMPARISON +int IComplex::operator== (const IComplex &op) const { + return (this->r == op.r) && (this->i == op.i); +} +int IComplex::operator== (const std::complex &op) const { + return (this->r == op.real ()) && (this->i == op.imag ()); +} +int IComplex::operator== (double op) const { + return (this->r == op) && (this->i == 0); +} +int IComplex::operator!= (const IComplex &op) const { + return !((*this) == op); +} +int IComplex::operator!= (const std::complex &op) const { + return !((*this) == op); +} +int IComplex::operator!= (double op) const { + return !((*this) == op); +} + + +// // ARITHMETIC +// A INTERVAL REAL PART +// B INTERNAL IMAG PART +// (A + iB) + (C + iD) = (A+C) + i(B+D) +IComplex IComplex::operator+ (const IComplex &op) const { + return IComplex (this->r + op.r, this->i + op.i); +} + +// (A + iB) + (C + iD) = (A+C) + i(B+D) +IComplex IComplex::operator+ (const std::complex &op) const { + return IComplex (this->r + op.real (), this->i + op.imag ()); +} + +// (A + iB) + C = (A+C) + iB +IComplex IComplex::operator+ (const Interval &op) const { + return IComplex (this->r + op, this->i); +} + +// (A + iB) + C = (A+C) + iB +IComplex IComplex::operator+ (double op) const { + return IComplex (this->r + op, this->i); +} + + +// -(A + iB) = (-A) + i(-B) +IComplex IComplex::operator- () const { + return IComplex (-(this->r), -(this->i)); +} + +// (A + iB) - (C + iD) = (A-C) + i(B-D) +IComplex IComplex::operator- (const IComplex &op) const { + return IComplex (this->r - op.r, this->i - op.i); +} + +// (A + iB) - (C + iD) = (A-C) + i(B-D) +IComplex IComplex::operator- (const std::complex &op) const { + return IComplex (this->r - op.real (), this->i - op.imag ()); +} + +// (A + iB) - (C + iD) = (A-C) + i(B-D) +IComplex IComplex::operator- (const Interval &op) const { + return IComplex (this->r - op, this->i); +} + +// (A + iB) - C = (A-C) + iB +IComplex IComplex::operator- (double op) const { + return IComplex (this->r - op, this->i); +} + + +// (A + iB) * (C + iD) = (A*C - B*D) + i(A*D + B*C) +IComplex IComplex::operator* (const IComplex &op) const { + return IComplex (this->r * op.r - this->i * op.i, + this->r * op.i + this->i * op.r); +} + +// (A + iB) * (C + iD) = (A*C - B*D) + i(A*D + B*C) +IComplex IComplex::operator* (const std::complex &op) const { + return IComplex (this->r * op.real () - this->i * op.imag (), + this->r * op.imag () + this->i * op.real ()); +} + +// (A + iB) * C = (A*C) + i(B*C) +IComplex IComplex::operator* (const Interval &op) const { + return IComplex (this->r * op, this->i * op); +} + +// (A + iB) * C = (A*C) + i(B*C) +IComplex IComplex::operator* (double op) const { + return IComplex (this->r * op, this->i * op); +} + + + +IComplex IComplex::operator/ (const IComplex &op) const { + return (*this) * (op.inverse ()); +} +IComplex IComplex::operator/ (const std::complex &op) const { + return (*this) * (IComplex(op).inverse ()); +} +IComplex IComplex::operator/ (const Interval &op) const { + return IComplex (this->r / op, this->i / op); +} +IComplex IComplex::operator/ (double op) const { + return IComplex (this->r / op, this->i / op); +} + + + +// // MISC +IComplex IComplex::inverse () const { + if (this->containsZero()) + return IComplex (NAN, NAN, NAN, NAN); + IComplex inv_edge[4], rop; + inv_edge[0] = inverse_horizontal (this->r.a, this->r.b, this->i.a); + inv_edge[1] = inverse_horizontal (this->r.a, this->r.b, this->i.b); + inv_edge[2] = inverse_vertical (this->i.a, this->i.b, this->r.a); + inv_edge[3] = inverse_vertical (this->i.a, this->i.b, this->r.b); + + rop.r.a = MIN (MIN (MIN (inv_edge[0].r.a, + inv_edge[1].r.a), inv_edge[2].r.a), inv_edge[3].r.a); + rop.r.b = MAX (MAX (MAX (inv_edge[0].r.b, + inv_edge[1].r.b), inv_edge[2].r.b), inv_edge[3].r.b); + rop.i.a = MIN (MIN (MIN (inv_edge[0].i.a, + inv_edge[1].i.a), inv_edge[2].i.a), inv_edge[3].i.a); + rop.i.b = MAX (MAX (MAX (inv_edge[0].i.b, + inv_edge[1].i.b), inv_edge[2].i.b), inv_edge[3].i.b); + + return rop; +} + +// CHECKS WHETHER 0 IS IN A + iB +int IComplex::containsZero () const { + return ((this->r).containsZero () && (this->i).containsZero()); +} + +// CHECKS WHETHER A VALUE IS in A + iB +int IComplex::contains (double op) const { + return ((this->r).contains (op) && (this->i).contains(op)); +} +int IComplex::contains (const std::complex &op) { + return ((this->r).contains(op.real()) && (this->i).contains(op.imag ())); +} + + +// CHECK WHETHER A COMPLEX INTERVAL IS INSIDE ANOTHER WITHOUT TOUCHING THE EDGES +int IComplex::subset (const IComplex &op) const { + return (this->r).subset (op.r) && (this->i).subset (op.i); +} + +// CHECK WHETHER A COMPLEX INTERVAL IS INSIDE ANOTHER +int IComplex::subsetEq (const IComplex &op) const { + return (this->r).subsetEq (op.r) && (this->i).subsetEq (op.i); +} + + + + +// // NON MEMBER FUNCTIONS + +// OVERLOAD OUTPUT STREAM +std::ostream & operator<< (std::ostream &output, const IComplex &op) { + output << op.r << " + i*" << op.i; + return output; +} + + +// C + (A + iB) = (A+C) + iB +IComplex operator+ (const std::complex &op1, const IComplex &op2) { + return IComplex (op1.real () + op2.r, op1.imag () + op2.i); +} + +// C + (A + iB) = (A+C) + iB +IComplex operator+ (const Interval &op1, const IComplex &op2) { + return IComplex (op1 + op2.r, op2.i); +} + +// C + (A + iB) = (A+C) + iB +IComplex operator+ (double op1, const IComplex &op2) { + return IComplex (op1 + op2.r, op2.i); +} + +// C - (A + iB) = (C-A) + i(-B) +IComplex operator- (const std::complex &op1, const IComplex &op2) { + return IComplex (op1.real () - op2.r, op1.imag () -op2.i); +} +// C - (A + iB) = (C-A) + i(-B) +IComplex operator- (const Interval &op1, const IComplex &op2) { + return IComplex (op1 - op2.r, -op2.i); +} + +// C - (A + iB) = (C-A) + i(-B) +IComplex operator- (double op1, const IComplex &op2) { + return IComplex (op1 - op2.r, -op2.i); +} + +// C * (A + iB) = (C*A) + i(C*B) +IComplex operator* (const std::complex &op1, const IComplex &op2) { + return (IComplex (op1.real () * op2.r - op1.imag () * op2.i, + op1.real () * op2.i + op1.imag () * op2.r)); +} + +// C * (A + iB) = (C*A) + i(C*B) +IComplex operator* (const Interval &op1, const IComplex &op2) { + return (IComplex (op1 * op2.r, op1 * op2.i)); +} + +// C * (A + iB) = (C*A) + i(C*B) +IComplex operator* (double op1, const IComplex &op2) { + return (IComplex (op1 * op2.r, op1 * op2.i)); +} + +// C / (A + iB) = C * (A + iB)^-1 +IComplex operator/ (const std::complex &op1, const IComplex &op2) { + return op1 * (op2.inverse ()); +} + +// C / (A + iB) = C * (A + iB)^-1 +IComplex operator/ (const Interval &op1, const IComplex &op2) { + return op1 * (op2.inverse ()); +} + +// C / (A + iB) = C * (A + iB)^-1 +IComplex operator/ (double op1, const IComplex &op2) { + return op1 * (op2.inverse ()); +} + + +// MISC METHODS FOR INVERSION + +// returns inverse of horizontal segment [a,b] + ci +IComplex inverse_horizontal (double a, double b, double c) { + IComplex rop; + int oldRoundingMode = fegetround (); + // case c = 0 + if (c == 0) { + if (a < 0 && b > 0) { + rop.r.a = rop.r.b = rop.i.a = rop.i.b = NAN; + } else { + rop.i.a = rop.i.b = 0.0; + fesetround (FE_DOWNWARD); + rop.r.a = 1.0 / b; + fesetround (FE_UPWARD); + rop.r.b = 1.0 / a; + } + // case c > 0 + } else if (c > 0) { + // case a<0 (-a)) + rop.i.b = g (b, c); + else + rop.i.b = g (a, c); + } + + // case 0<=a +#include + + + +// ACTIVATE ROUND MODE CONTROL +#pragma STDC FENV_ACCESS on + +// MACROS FOR MIN AND MAX +#define MIN(a,b) (((a) < (b))? (a) : (b)) +#define MAX(a,b) (((a) > (b))? (a) : (b)) + +// DEFAULT CONSTRUCTOR +Interval::Interval (double op) { + a = b = op; +} + +// CONSTRUCTOR BY ELEMENT +Interval::Interval (double op1, double op2) { + a = op1; b = op2; +} + +// COPY CONSTRUCTOR +Interval::Interval (const Interval &op) { + a = op.a; b = op.b; +} + +// DESTRUCTOR. NOTHING TO DO... +Interval::~Interval () {} + + +// ASSIGNMENT +Interval & Interval::operator= (const Interval &op) { + // CHECK WHETER WE ARE ASIGNING SOMETHING OVER ITSELF + if (this == &op) return (*this); + // UPDATE VALUES + this->a = op.a; this->b = op.b; + return (*this); +} + +Interval & Interval::operator= (double op) { + // UPDATE VALUES + this->a = op; this->b = op; + return (*this); +} + + +// COMPARISON +int Interval::operator== (const Interval &op) const { + // BOTH END POINTS MUST BE THE SAME TO BE EQUAL + if (this->a == op.a && this->b == op.b) return 1; + return 0; +} +int Interval::operator!= (const Interval &op) const { + // BOTH END POINTS MUST BE THE SAME TO BE NOT DIFFERENT + if (this->a == op.a && this->b == op.b) return 0; + return 1; +} + + +// ARITHMETIC + +// [a,b] + [c,d] = [a+c, b+d] +Interval Interval::operator+ (const Interval &op) const { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + fesetround (FE_DOWNWARD); + rop.a = this->a + op.a; + + fesetround (FE_UPWARD); + rop.b = this->b + op.b; + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + + return rop; +} +// [a,b] + c = [a+c, b+c] +Interval Interval::operator+ (double op) const { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + fesetround (FE_DOWNWARD); + rop.a = this->a + op; + + fesetround (FE_UPWARD); + rop.b = this->b + op; + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + + return rop; +} + +// CHANGE OF SIGN OF AN INTERVAL +// -[a,b] = [-b,a] +Interval Interval::operator- () const { + Interval rop; + rop.a = -this->b; + rop.b = -this->a; + return rop; +} + +// [a,b] - [c,d] = [a-d, b-c] +Interval Interval::operator- (const Interval &op) const { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + fesetround (FE_DOWNWARD); + rop.a = this->a - op.b; + + fesetround (FE_UPWARD); + rop.b = this->b - op.a; + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + + return rop; +} + +// [a,b] - c = [a-c, b-c] +Interval Interval::operator- (double op) const { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + fesetround (FE_DOWNWARD); + rop.a = this->a - op; + + fesetround (FE_UPWARD); + rop.b = this->b - op; + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + + return rop; +} + +// [a,b]*[c,d] = [min(a*c,a*d,b*c,b*d), max(a*c,a*d,b*c,b*d)] +Interval Interval::operator* (const Interval &op) const { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + + fesetround (FE_DOWNWARD); + rop.a = MIN(MIN(MIN(this->a*op.a, this->a*op.b), this->b*op.a), this->b*op.b); + + fesetround (FE_UPWARD); + rop.b = MAX(MAX(MAX(this->a*op.a, this->a*op.b), this->b*op.a), this->b*op.b); + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + return rop; + + +} + +// [a,b]*c = [min(a*c,b*c), max(a*c,b*c)] +Interval Interval::operator* (double op) const { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + + fesetround (FE_DOWNWARD); + rop.a = MIN(this->a*op, this->b*op); + + fesetround (FE_UPWARD); + rop.b = MAX(this->a*op, this->b*op); + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + return rop; +} + + +// [a,b]/[c,d] = [min(a/c,a/d,b/c,b/d), max(a/c,a/d,b/c,b/d)] +Interval Interval::operator/ (const Interval &op) const { + if (op.containsZero()) return Interval (NAN); + + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + fesetround (FE_DOWNWARD); + rop.a = MIN(MIN(MIN(this->a/op.a, this->a/op.b), this->b/op.a), this->b/op.b); + + fesetround (FE_UPWARD); + rop.b = MAX(MAX(MAX(this->a/op.a, this->a/op.b), this->b/op.a), this->b/op.b); + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + return rop; +} + +// [a,b]/c = [min(a/c,b/c), max(a/c,b/c)] +Interval Interval::operator/ (double op) const { + if (op == 0.0) return Interval (NAN); + + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + + fesetround (FE_DOWNWARD); + rop.a = MIN(this->a/op, this->b/op); + + fesetround (FE_UPWARD); + rop.b = MAX(this->a/op, this->b/op); + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + return rop; +} + + + + +// MISC +Interval Interval::inverse () const { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + if (0a || this->b<0) { + fesetround (FE_DOWNWARD); + rop.a = 1.0 / (this->b); + fesetround (FE_UPWARD); + rop.b = 1.0 / (this->a); + } else { + rop.a = rop.b = NAN; + } + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + return rop; +} + +// CHECKS WHETHER 0 in [a,b] +int Interval::containsZero () const { + if (this->a <= 0 && this->b >= 0) return 1; + return 0; +} +// CHECKS WHETHER a value is in [a,b] +int Interval::contains (double op) const { + if (this->a <= op && this->b >= op) return 1; + return 0; + +} +// CHECKS WHETHER [a,b] subset of [c,d] with a!=c and b!=d +int Interval::subset (const Interval &op) const { + if (this->a > op.a && this->b < op.b) return 1; + return 0; +} + +// CHECKS WHETHER [a,b] subseteq of [c,d] +int Interval::subsetEq (const Interval &op) const { + if (this->a >= op.a && this->b <= op.b) return 1; + return 0; +} + + + + + +//EXTERNAL FUNCTIONS +// OVERLOAD OUTPUT STREAM +std::ostream & operator<< (std::ostream &output, const Interval &op) { + output << "[" << op.a << ", " << op.b << "]"; + return output; +} + + +// c + [a,b] = [a+c, b+c] +Interval operator+ (double op1, const Interval &op2) { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + fesetround (FE_DOWNWARD); + rop.a = op1 + op2.a; + + fesetround (FE_UPWARD); + rop.b = op1 + op2.b; + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + + return rop; + +} + +// c - [a,b] = [c-b, c-a] +Interval operator- (double op1, const Interval &op2) { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + fesetround (FE_DOWNWARD); + rop.a = op1 - op2.b; + + fesetround (FE_UPWARD); + rop.b = op1 - op2.a; + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + + return rop; +} + + +// c*[a,b] = [min(a*c,b*c), max(a*c,b*c)] +Interval operator* (double op1, const Interval &op2) { + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + + fesetround (FE_DOWNWARD); + rop.a = MIN(op2.a*op1, op2.b*op1); + + fesetround (FE_UPWARD); + rop.b = MAX(op2.a*op1, op2.b*op1); + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + return rop; +} + + +// c/[a,b] = [min(c/a,c/b), max(c/a,c/b)] +Interval operator/ (double op1, const Interval &op2) { + if (op2.containsZero ()) return Interval (NAN); + + Interval rop; + + // STORE OLD ROUNDING MODE + int oldRound = fegetround (); + + + fesetround (FE_DOWNWARD); + rop.a = MIN(op1 / op2.a, op1 / op2.b); + + fesetround (FE_UPWARD); + rop.b = MAX(op1 / op2.a, op1 / op2.b); + + // RESTORE OLD ROUNDING MODE + fesetround (oldRound); + return rop; +} \ No newline at end of file diff --git a/lib/mp_complex.cpp b/lib/mp_complex.cpp new file mode 100644 index 0000000..9847402 --- /dev/null +++ b/lib/mp_complex.cpp @@ -0,0 +1,389 @@ +#include "mp_complex.hpp" + +MPComplex::MPComplex (double op) { + mpfr_init (r); + mpfr_set_d (r, op, MPFR_RNDN); + mpfr_init_set_si (i, 0, MPFR_RNDN); +} +MPComplex::MPComplex (double op1, double op2) { + mpfr_init_set_d (r, op1, MPFR_RNDN); + mpfr_init_set_d (i, op2, MPFR_RNDN); +} +MPComplex::MPComplex (const std::complex &op) { + mpfr_init_set_d (r, op.real (), MPFR_RNDN); + mpfr_init_set_d (i, op.imag (), MPFR_RNDN); +} +MPComplex::MPComplex (const mpfr_t op) { + mpfr_init_set (r, op, MPFR_RNDN); + mpfr_init_set_si (i, 0, MPFR_RNDN); +} +MPComplex::MPComplex (const mpfr_t op1, const mpfr_t op2) { + mpfr_init_set (r, op1, MPFR_RNDN); + mpfr_init_set (i, op2, MPFR_RNDN); +} +MPComplex::~MPComplex () { + mpfr_clears (r, i, NULL); +} + + // ASSIGNMENT +MPComplex & MPComplex::operator= (double op) { + mpfr_set_d (this->r, op, MPFR_RNDN); + mpfr_set_si (this->i, 0, MPFR_RNDN); + return (*this); +} +MPComplex & MPComplex::operator= (const std::complex &op) { + mpfr_set_d (this->r, op.real (), MPFR_RNDN); + mpfr_set_d (this->i, op.imag (), MPFR_RNDN); + return (*this); +} +MPComplex & MPComplex::operator= (const mpfr_t op) { + mpfr_set (this->r, op, MPFR_RNDN); + mpfr_set_si (this->i, 0, MPFR_RNDN); + return (*this); +} +MPComplex & MPComplex::operator= (const MPComplex &op) { + mpfr_set (this->r, op.r, MPFR_RNDN); + mpfr_set (this->i, op.i, MPFR_RNDN); + return (*this); +} + +// COMPARISON + +int MPComplex::operator == (double op) const { + return (mpfr_cmp_d (this->r, op) == 0) && (mpfr_cmp_si (this->i, 0) == 0); +} +int MPComplex::operator == (const std::complex &op) const { + return (mpfr_cmp_d (this->r, op.real ()) == 0) && (mpfr_cmp_d (this->i, op.imag ()) == 0); +} +int MPComplex::operator == (const MPComplex &op) const { + return (mpfr_cmp (this->r, op.r) == 0) && (mpfr_cmp (this->i, op.i) == 0); +} +int MPComplex::operator != (double op) const { + return (mpfr_cmp_d (this->r, op) != 0) || (mpfr_cmp_si (this->i, 0) != 0); +} +int MPComplex::operator != (const std::complex &op) const { + return (mpfr_cmp_d (this->r, op.real ()) != 0) || (mpfr_cmp_d (this->i, op.imag ()) != 0); +} +int MPComplex::operator != (const MPComplex &op) const { + return (mpfr_cmp (this->r, op.r) != 0) || (mpfr_cmp (this->i, op.i) != 0); +} + +// ARITHMETIC +MPComplex MPComplex::operator+ (double op) const { + MPComplex rop; + mpfr_add_d (rop.r, this->r, op, MPFR_RNDN); + mpfr_set (rop.i, this->i, MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator+ (const std::complex &op) const { + MPComplex rop; + mpfr_add_d (rop.r, this->r, op.real (), MPFR_RNDN); + mpfr_add_d (rop.i, this->i, op.imag (), MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator+ (const mpfr_t op) const { + MPComplex rop; + mpfr_add (rop.r, this->r, op, MPFR_RNDN); + mpfr_set (rop.i, this->i, MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator+ (const MPComplex &op) const { + MPComplex rop; + mpfr_add (rop.r, this->r, op.r, MPFR_RNDN); + mpfr_add (rop.i, this->i, op.i, MPFR_RNDN); + return rop; +} + +MPComplex MPComplex::operator- () const { + MPComplex rop; + mpfr_neg (rop.r, this->r, MPFR_RNDN); + mpfr_neg (rop.i, this->i, MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator- (double op) const { + MPComplex rop; + mpfr_sub_d (rop.r, this->r, op, MPFR_RNDN); + mpfr_set (rop.i, this->i, MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator- (const std::complex &op) const { + MPComplex rop; + mpfr_sub_d (rop.r, this->r, op.real (), MPFR_RNDN); + mpfr_sub_d (rop.i, this->i, op.imag (), MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator- (const mpfr_t op) const { + MPComplex rop; + mpfr_sub (rop.r, this->r, op, MPFR_RNDN); + mpfr_set (rop.i, this->i, MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator- (const MPComplex &op) const { + MPComplex rop; + mpfr_sub (rop.r, this->r, op.r, MPFR_RNDN); + mpfr_sub (rop.i, this->i, op.i, MPFR_RNDN); + return rop; +} + +MPComplex MPComplex::operator* (double op) const { + MPComplex rop; + mpfr_mul_d (rop.r, this->r, op, MPFR_RNDN); + mpfr_mul_d (rop.i, this->i, op, MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator* (const std::complex &op) const { + MPComplex rop; + mpfr_t aux; mpfr_init (aux); + mpfr_mul_d (rop.r, this->r, op.real (), MPFR_RNDN); + mpfr_mul_d (aux, this->i, op.imag (), MPFR_RNDN); + mpfr_sub (rop.r, rop.r, aux, MPFR_RNDN); + mpfr_mul_d (rop.i, this->r, op.imag (), MPFR_RNDN); + mpfr_mul_d (aux, this->i, op.real (), MPFR_RNDN); + mpfr_add (rop.i, rop.i, aux, MPFR_RNDN); + + mpfr_clear (aux); + return rop; +} +MPComplex MPComplex::operator* (const mpfr_t op) const { + MPComplex rop; + mpfr_mul (rop.r, this->r, op, MPFR_RNDN); + mpfr_mul (rop.i, this->i, op, MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator* (const MPComplex &op) const { + MPComplex rop; + mpfr_t aux; mpfr_init (aux); + mpfr_mul (rop.r, this->r, op.r, MPFR_RNDN); + mpfr_mul (aux, this->i, op.i, MPFR_RNDN); + mpfr_sub (rop.r, rop.r, aux, MPFR_RNDN); + mpfr_mul (rop.i, this->r, op.i, MPFR_RNDN); + mpfr_mul (aux, this->i, op.r, MPFR_RNDN); + mpfr_add (rop.i, rop.i, aux, MPFR_RNDN); + + mpfr_clear (aux); + return rop; +} + + +MPComplex MPComplex::operator/ (double op) const { + MPComplex rop; + mpfr_div_d (rop.r, this->r, op, MPFR_RNDN); + mpfr_div_d (rop.i, this->i, op, MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator/ (const std::complex &op) const { + MPComplex rop; + mpfr_t aux; mpfr_init (aux); + double factor = 1.0 / ((op.real() * op.real()) + (op.imag() * op.imag())); + mpfr_mul_d (rop.r, this->r, op.real (), MPFR_RNDN); + mpfr_mul_d (aux, this->i, op.imag (), MPFR_RNDN); + mpfr_add (rop.r, rop.r, aux, MPFR_RNDN); + mpfr_mul_d (rop.r, rop.r, factor, MPFR_RNDN); + mpfr_mul_d (rop.i, this->i, op.real (), MPFR_RNDN); + mpfr_mul_d (aux, this->r, op.imag (), MPFR_RNDN); + mpfr_sub (rop.i, rop.i, aux, MPFR_RNDN); + mpfr_mul_d (rop.i, rop.i, factor, MPFR_RNDN); + mpfr_clear (aux); + return rop; +} +MPComplex MPComplex::operator/ (const mpfr_t op) const { + MPComplex rop; + mpfr_div (rop.r, this->r, op, MPFR_RNDN); + mpfr_div (rop.i, this->i, op, MPFR_RNDN); + return rop; +} +MPComplex MPComplex::operator/ (const MPComplex &op) const { + MPComplex rop; + mpfr_t aux, factor; mpfr_inits (aux, factor, NULL); + + mpfr_mul (factor, op.r, op.r, MPFR_RNDN); + mpfr_mul (aux, op.i, op.i, MPFR_RNDN); + mpfr_add (factor, factor, aux, MPFR_RNDN); + mpfr_si_div (factor, 1, factor, MPFR_RNDN); + + mpfr_mul (rop.r, this->r, op.r, MPFR_RNDN); + mpfr_mul (aux, this->i, op.i, MPFR_RNDN); + mpfr_add (rop.r, rop.r, aux, MPFR_RNDN); + mpfr_mul (rop.r, rop.r, factor, MPFR_RNDN); + mpfr_mul (rop.i, this->i, op.r, MPFR_RNDN); + mpfr_mul (aux, this->r, op.i, MPFR_RNDN); + mpfr_sub (rop.i, rop.i, aux, MPFR_RNDN); + mpfr_mul (rop.i, rop.i, factor, MPFR_RNDN); + + mpfr_clears (aux, factor, NULL); + return rop; +} + +// NON MEMBER FUNCTIONS + +std::ostream & operator<< (std::ostream &output, const MPComplex &op) { + + mpfr_exp_t exponent; + char *string = mpfr_get_str (NULL, &exponent, 10, 0, op.r, MPFR_RNDD); + if (mpfr_signbit (op.r) != 0) + output << "-" << string[1] << "." << string+2 << "e" << exponent-1; + else + output << string[0] << "." << string+1 << "e" << exponent-1; + mpfr_free_str (string); + + string = mpfr_get_str (NULL, &exponent, 10, 0, op.i, MPFR_RNDU); + if (mpfr_signbit (op.i) != 0) + output << "-" << string[1] << "." << string+2 << "e" << exponent-1 << "i"; + else + output << "+" << string[0] << "." << string+1 << "e" << exponent-1 << "i"; + + mpfr_free_str (string); + return output; +} + + + +MPComplex operator+ (double op1, const MPComplex &op2) { + MPComplex rop; + mpfr_add_d (rop.r, op2.r, op1, MPFR_RNDN); + mpfr_set (rop.i, op2.i, MPFR_RNDN); + return rop; + +} +MPComplex operator+ (const std::complex &op1, const MPComplex &op2) { + MPComplex rop; + mpfr_add_d (rop.r, op2.r, op1.real (), MPFR_RNDN); + mpfr_add_d (rop.i, op2.i, op1.imag (), MPFR_RNDN); + return rop; +} +MPComplex operator+ (const mpfr_t op1, const MPComplex &op2) { + MPComplex rop; + mpfr_add (rop.r, op2.r, op1, MPFR_RNDN); + mpfr_set (rop.i, op2.i, MPFR_RNDN); + return rop; +} + + +MPComplex operator- (double op1, const MPComplex &op2) { + MPComplex rop; + mpfr_d_sub (rop.r, op1, op2.r, MPFR_RNDN); + mpfr_neg (rop.i, op2.i, MPFR_RNDN); + return rop; + +} +MPComplex operator- (const std::complex &op1, const MPComplex &op2) { + MPComplex rop; + mpfr_d_sub (rop.r, op1.real (), op2.r, MPFR_RNDN); + mpfr_d_sub (rop.i, op1.imag (), op2.i, MPFR_RNDN); + return rop; +} +MPComplex operator- (const mpfr_t op1, const MPComplex &op2) { + MPComplex rop; + mpfr_sub (rop.r, op1, op2.r, MPFR_RNDN); + mpfr_neg (rop.i, op2.i, MPFR_RNDN); + return rop; +} + + + + + +MPComplex operator* (double op1, const MPComplex &op2) { + MPComplex rop; + mpfr_mul_d (rop.r, op2.r, op1, MPFR_RNDN); + mpfr_mul_d (rop.i, op2.i, op1, MPFR_RNDN); + return rop; + +} +MPComplex operator* (const std::complex &op1, const MPComplex &op2) { + MPComplex rop; + mpfr_t aux; mpfr_init (aux); + mpfr_mul_d (rop.r, op2.r, op1.real (), MPFR_RNDN); + mpfr_mul_d (aux, op2.i, op1.imag (), MPFR_RNDN); + mpfr_sub (rop.r, rop.r, aux, MPFR_RNDN); + mpfr_mul_d (rop.i, op2.r, op1.imag (), MPFR_RNDN); + mpfr_mul_d (aux, op2.r, op1.real (), MPFR_RNDN); + mpfr_add (rop.i, rop.i, aux, MPFR_RNDN); + + mpfr_clear (aux); + return rop; + +} +MPComplex operator* (const mpfr_t op1, const MPComplex &op2) { + MPComplex rop; + mpfr_mul (rop.r, op2.r, op1, MPFR_RNDN); + mpfr_mul (rop.i, op2.i, op1, MPFR_RNDN); + return rop; + +} + + + +MPComplex operator/ (double op1, const MPComplex &op2) { + MPComplex rop; + mpfr_t factor; + mpfr_init (factor); + mpfr_mul (factor, op2.r, op2.r, MPFR_RNDN); + mpfr_mul (rop.r, op2.i, op2.i, MPFR_RNDN); + mpfr_add (factor, factor, rop.r, MPFR_RNDN); + mpfr_si_div (factor, 1, factor, MPFR_RNDN); + + mpfr_mul_d (rop.r, op2.r, op1, MPFR_RNDN); + mpfr_mul (rop.r, rop.r, factor, MPFR_RNDN); + mpfr_mul_d (rop.i, op2.i, op1, MPFR_RNDN); + mpfr_neg (rop.i, rop.i, MPFR_RNDN); + mpfr_mul (rop.i, rop.i, factor, MPFR_RNDN); + + mpfr_clear (factor); + return rop; + +} +MPComplex operator/ (const std::complex &op1, const MPComplex &op2) { + MPComplex rop; + mpfr_t factor, aux; + mpfr_inits (factor, aux, NULL); + mpfr_mul (factor, op2.r, op2.r, MPFR_RNDN); + mpfr_mul (aux, op2.i, op2.i, MPFR_RNDN); + mpfr_add (factor, factor, aux, MPFR_RNDN); + mpfr_si_div (factor, 1, factor, MPFR_RNDN); + + mpfr_mul_d (rop.r, op2.r, op1.real (), MPFR_RNDN); + mpfr_mul_d (aux, op2.i, op1.imag (), MPFR_RNDN); + mpfr_add (rop.r, rop.r, aux, MPFR_RNDN); + mpfr_mul (rop.r, rop.r, factor, MPFR_RNDN); + mpfr_mul_d (rop.i, op2.r, op1.imag (), MPFR_RNDN); + mpfr_mul_d (aux, op2.i, op1.real (), MPFR_RNDN); + mpfr_sub (rop.i, rop.i, aux, MPFR_RNDN); + mpfr_neg (rop.i, rop.i, MPFR_RNDN); + mpfr_mul (rop.i, rop.i, factor, MPFR_RNDN); + + mpfr_clears (factor, aux, NULL); + return rop; + +} +MPComplex operator/ (const mpfr_t op1, const MPComplex &op2) { + MPComplex rop; + mpfr_t factor; + mpfr_init (factor); + mpfr_mul (factor, op2.r, op2.r, MPFR_RNDN); + mpfr_mul (rop.r, op2.i, op2.i, MPFR_RNDN); + mpfr_add (factor, factor, rop.r, MPFR_RNDN); + mpfr_si_div (factor, 1, factor, MPFR_RNDN); + + mpfr_mul (rop.r, op2.r, op1, MPFR_RNDN); + mpfr_mul (rop.r, rop.r, factor, MPFR_RNDN); + mpfr_mul (rop.i, op2.i, op1, MPFR_RNDN); + mpfr_neg (rop.i, rop.i, MPFR_RNDN); + mpfr_mul (rop.i, rop.i, factor, MPFR_RNDN); + + mpfr_clear (factor); + return rop; +} + +void abs (mpfr_t rop, const MPComplex &op) { + mpfr_t aux; + mpfr_init (aux); + mpfr_mul (rop, op.r, op.r, MPFR_RNDN); + mpfr_mul (aux, op.i, op.i, MPFR_RNDN); + mpfr_add (rop, rop, aux, MPFR_RNDN); + mpfr_sqrt (rop, rop, MPFR_RNDN); + + + mpfr_clear (aux); + return; +} diff --git a/lib/mp_icomplex.cpp b/lib/mp_icomplex.cpp new file mode 100644 index 0000000..658c059 --- /dev/null +++ b/lib/mp_icomplex.cpp @@ -0,0 +1,483 @@ +#include "mp_icomplex.hpp" + + // CONSTRUCTORS +MPIComplex::MPIComplex (double op) + : r (op), i (0.0) {} +MPIComplex::MPIComplex (double op1, double op2, double op3, double op4) + : r (op1, op2), i (op3, op4) {} +MPIComplex::MPIComplex (const mpfr_t op) + : r (op), i (0.0) {} +MPIComplex::MPIComplex (const MPComplex &op) + : r (op.r), i (op.i) {} +MPIComplex::MPIComplex (const mpfr_t op1, const mpfr_t op2, const mpfr_t op3, const mpfr_t op4) + : r (op1, op2), i (op3, op4) {} +MPIComplex::MPIComplex (const Interval &op1, const Interval &op2) + : r (op1), i (op2) {} +MPIComplex::MPIComplex (const MPInterval &op1, const MPInterval &op2) + : r (op1), i (op2) {} +MPIComplex::MPIComplex (const IComplex &op) + : r (op.r), i (op.i) {} + // COPY CONSTRUCTOR +MPIComplex::MPIComplex (const MPIComplex &op) + : r (op.r), i (op.i) {} +MPIComplex::~MPIComplex () {} + + // ASSIGNMENT +MPIComplex & MPIComplex::operator= (const MPIComplex &op) { + // CHECK WHETER WE ARE ASIGNING SOMETHING OVER ITSELF + if (this == &op) return (*this); + // UPDATE ELEMENTS + this->r = op.r; this->i = op.i; + return (*this); +} +MPIComplex & MPIComplex::operator= (const IComplex &op) { + // UPDATE ELEMENTS + this->r = op.r; this->i = op.i; + return (*this); +} +MPIComplex & MPIComplex::operator= (const MPInterval &op) { + // UPDATE ELEMENTS + this->r = op; this->i = 0.0; + return (*this); +} +MPIComplex & MPIComplex::operator= (const Interval &op) { + // UPDATE ELEMENTS + this->r = op; this->i = 0.0; + return (*this); +} +MPIComplex & MPIComplex::operator= (const mpfr_t op) { + // UPDATE ELEMENTS + this->r = op; this->i = 0.0; + return (*this); +} +MPIComplex & MPIComplex::operator= (double op) { + // UPDATE ELEMENTS + this->r = op; this->i = 0.0; + return (*this); +} + +// // COMPARISON + +int MPIComplex::operator== (double op) const { + return (this->r == op) && (this->i == 0); +} +int MPIComplex::operator== (const std::complex &op) const { + return (this->r == op.real() && this->i == op.imag ()); +} +int MPIComplex::operator== (const Interval &op) const { + return (this->r == op) && (this->i == 0); +} +int MPIComplex::operator== (const MPInterval &op) const { + return (this->r == op) && (this->i == 0); +} +int MPIComplex::operator== (const MPIComplex &op) const { + return (this->r == op.r) && (this->i == op.i); +} +int MPIComplex::operator!= (double op) const { + return (this->r != op) || (this->i != 0); +} +int MPIComplex::operator!= (const std::complex &op) const { + return (this->r != op.real() || this->i != op.imag ()); +} +int MPIComplex::operator!= (const Interval &op) const { + return (this->r != op) || (this->i != 0); +} +int MPIComplex::operator!= (const MPInterval &op) const { + return (this->r != op) || (this->i != 0); +} +int MPIComplex::operator!= (const MPIComplex &op) const { + return (this->r != op.r) || (this->i != op.i); +} + + + +// // ARITHMETIC +// A INTERVAL REAL PART +// B INTERNAL IMAG PART +// (A + iB) + (C + iD) = (A+C) + i(B+D) +MPIComplex MPIComplex::operator+ (const MPIComplex &op) const { + return MPIComplex (this->r + op.r, this->i + op.i); +} + +// (A + iB) + (C + iD) = (A+C) + i(B+D) +MPIComplex MPIComplex::operator+ (const IComplex &op) const { + return MPIComplex (this->r + op.r, this->i + op.i); +} + +// (A + iB) + C = (A+C) + iB +MPIComplex MPIComplex::operator+ (const MPInterval &op) const { + return MPIComplex (this->r + op, this->i); +} + +// (A + iB) + C = (A+C) + iB +MPIComplex MPIComplex::operator+ (const Interval &op) const { + return MPIComplex (this->r + op, this->i); +} + +// (A + iB) + C = (A+C) + iB +MPIComplex MPIComplex::operator+ (const mpfr_t op) const { + return MPIComplex (this->r + op, this->i); +} + +// (A + iB) + C = (A+C) + iB +MPIComplex MPIComplex::operator+ (double op) const { + return MPIComplex (this->r + op, this->i); +} + + + +// -(A + iB) = (-A) + i(-B) +MPIComplex MPIComplex::operator- () const { + return MPIComplex (-(this->r), -(this->i)); +} + +// (A + iB) - (C + iD) = (A-C) + i(B-D) +MPIComplex MPIComplex::operator- (const MPIComplex &op) const { + return MPIComplex (this->r - op.r, this->i - op.i); +} + +// (A + iB) - (C + iD) = (A-C) + i(B-D) +MPIComplex MPIComplex::operator- (const IComplex &op) const { + return MPIComplex (this->r - op.r, this->i - op.i); +} + +// (A + iB) - C = (A-C) + iB +MPIComplex MPIComplex::operator- (const MPInterval &op) const { + return MPIComplex (this->r - op, this->i); +} + +// (A + iB) - C = (A-C) + iB +MPIComplex MPIComplex::operator- (const Interval &op) const { + return MPIComplex (this->r - op, this->i); +} + +// (A + iB) - C = (A-C) + iB +MPIComplex MPIComplex::operator- (const mpfr_t op) const { + return MPIComplex (this->r - op, this->i); +} + +// (A + iB) - C = (A-C) + iB +MPIComplex MPIComplex::operator- (double op) const { + return MPIComplex (this->r - op, this->i); +} + + + +// (A + iB) * (C + iD) = (A*C - B*D) + i(A*D + B*C) +MPIComplex MPIComplex::operator* (const MPIComplex &op) const { + return MPIComplex (this->r * op.r - this->i * op.i, + this->r * op.i + this->i * op.r); +} + +// (A + iB) * (C + iD) = (A*C - B*D) + i(A*D + B*C) +MPIComplex MPIComplex::operator* (const IComplex &op) const { + return MPIComplex (this->r * op.r - this->i * op.i, + this->r * op.i + this->i * op.r); +} + +// (A + iB) * C = (A*C) + i(B*C) +MPIComplex MPIComplex::operator* (const MPInterval &op) const { + return MPIComplex (this->r * op, this->i * op); +} + +// (A + iB) * C = (A*C) + i(B*C) +MPIComplex MPIComplex::operator* (const Interval &op) const { + return MPIComplex (this->r * op, this->i * op); +} + +// (A + iB) * C = (A*C) + i(B*C) +MPIComplex MPIComplex::operator* (const mpfr_t op) const { + return MPIComplex (this->r * op, this->i * op); +} + +// (A + iB) * C = (A*C) + i(B*C) +MPIComplex MPIComplex::operator* (double op) const { + return MPIComplex (this->r * op, this->i * op); +} + + +MPIComplex MPIComplex::operator/ (const MPIComplex &op) const { + return (*this) * (op.inverse ()); +} +MPIComplex MPIComplex::operator/ (const IComplex &op) const { + return (*this) * (op.inverse ()); +} +MPIComplex MPIComplex::operator/ (const MPInterval &op) const { + return (*this) * (op.inverse ()); +} +MPIComplex MPIComplex::operator/ (const Interval &op) const { + return (*this) * (op.inverse ()); +} +MPIComplex MPIComplex::operator/ (const mpfr_t op) const { + return (*this) * (MPIComplex(op).inverse ()); +} +MPIComplex MPIComplex::operator/ (double op) const { + return (*this) * (IComplex(op).inverse ()); +} + + + +// // MISC +MPIComplex MPIComplex::inverse () const { + if (this->containsZero()) { + mpfr_t aux; + mpfr_init (aux); + mpfr_set_nan (aux); + MPIComplex rop (aux); + mpfr_clear (aux); + return rop; + } + + MPIComplex inv_edge[4], rop; + inv_edge[0] = inverse_horizontal (this->r.a, this->r.b, this->i.a); + inv_edge[1] = inverse_horizontal (this->r.a, this->r.b, this->i.b); + inv_edge[2] = inverse_vertical (this->i.a, this->i.b, this->r.a); + inv_edge[3] = inverse_vertical (this->i.a, this->i.b, this->r.b); + + mpfr_min (rop.r.a, inv_edge[0].r.a, inv_edge[1].r.a, MPFR_RNDD); + mpfr_min (rop.r.a, rop.r.a, inv_edge[2].r.a, MPFR_RNDD); + mpfr_min (rop.r.a, rop.r.a, inv_edge[3].r.a, MPFR_RNDD); + + + mpfr_max (rop.r.b, inv_edge[0].r.b, inv_edge[1].r.b, MPFR_RNDU); + mpfr_max (rop.r.b, rop.r.b, inv_edge[2].r.b, MPFR_RNDU); + mpfr_max (rop.r.b, rop.r.b, inv_edge[3].r.b, MPFR_RNDU); + + mpfr_min (rop.i.a, inv_edge[0].i.a, inv_edge[1].i.a, MPFR_RNDD); + mpfr_min (rop.i.a, rop.i.a, inv_edge[2].i.a, MPFR_RNDD); + mpfr_min (rop.i.a, rop.i.a, inv_edge[3].i.a, MPFR_RNDD); + + + mpfr_max (rop.i.b, inv_edge[0].i.b, inv_edge[1].i.b, MPFR_RNDU); + mpfr_max (rop.i.b, rop.i.b, inv_edge[2].i.b, MPFR_RNDU); + mpfr_max (rop.i.b, rop.i.b, inv_edge[3].i.b, MPFR_RNDU); + + + return rop; + +} + +int MPIComplex::containsZero () const { + return ((this->r).containsZero () && (this->i).containsZero()); +} +int MPIComplex::contains (double op) const { + return ((this->r).contains (op) && (this->i).contains(op)); +} +int MPIComplex::contains (const mpfr_t op) const { + return ((this->r).contains (op) && (this->i).contains(op)); +} +int MPIComplex::subset (const MPIComplex &op) const { + return (this->r).subset (op.r) && (this->i).subset (op.i); +} + + + + +// // NON MEMBER FUNCTIONS + +std::ostream & operator<< (std::ostream &output, const MPIComplex &op) { + output << op.r << " + i" << op.i; + return output; +} + + +MPIComplex operator+ (const IComplex &op1, const MPIComplex &op2) { + return MPIComplex (op1.r + op2.r, op1.i + op2.i); +} +MPIComplex operator+ (const MPInterval &op1, const MPIComplex &op2) { + return MPIComplex (op2.r, op1 + op2.i); +} +MPIComplex operator+ (const Interval &op1, const MPIComplex &op2) { + return MPIComplex (op1 + op2.r, op2.i); +} +MPIComplex operator+ (const mpfr_t op1, const MPIComplex &op2) { + return MPIComplex (op1 + op2.r, op2.i); +} +MPIComplex operator+ (double op1, const MPIComplex &op2) { + return MPIComplex (op1 + op2.r, op2.i); +} + +MPIComplex operator- (const IComplex &op1, const MPIComplex &op2) { + return MPIComplex (op1.r - op2.r, op1.i - op2.i); +} +MPIComplex operator- (const MPInterval &op1, const MPIComplex &op2) { + return MPIComplex (op1 - op2.r, -op2.i); +} +MPIComplex operator- (const Interval &op1, const MPIComplex &op2) { + return MPIComplex (op1 - op2.r, -op2.i); +} +MPIComplex operator- (const mpfr_t op1, const MPIComplex &op2) { + return MPIComplex (op1 - op2.r, -op2.i); +} +MPIComplex operator- (double op1, const MPIComplex &op2) { + return MPIComplex (op1 - op2.r, -op2.i); +} + +MPIComplex operator* (const IComplex &op1, const MPIComplex &op2) { + return MPIComplex (op1.r * op2.r - op1.i * op2.i, + op1.r * op2.i + op1.i * op2.r); +} +MPIComplex operator* (const MPInterval &op1, const MPIComplex &op2) { + return MPIComplex (op1 * op2.r, op1 * op2.i); +} +MPIComplex operator* (const Interval &op1, const MPIComplex &op2) { + return MPIComplex (op1 * op2.r, op1 * op2.i); +} +MPIComplex operator* (const mpfr_t op1, const MPIComplex &op2) { + return MPIComplex (op1 * op2.r, op1 * op2.i); +} +MPIComplex operator* (double op1, const MPIComplex &op2) { + return MPIComplex (op1 * op2.r, op1 * op2.i); +} + +MPIComplex operator/ (const IComplex &op1, const MPIComplex &op2) { + return op1 * (op2.inverse ()); +} +MPIComplex operator/ (const MPInterval &op1, const MPIComplex &op2) { + return op1 * (op2.inverse ()); +} +MPIComplex operator/ (const Interval &op1, const MPIComplex &op2) { + return op1 * (op2.inverse ()); +} +MPIComplex operator/ (const mpfr_t op1, const MPIComplex &op2) { + return op1 * (op2.inverse ()); +} +MPIComplex operator/ (double op1, const MPIComplex &op2) { + return op1 * (op2.inverse ()); +} + + + + +// // MISCELLANEOUS FUNCTIONS FOR INVERT COMPLEX INTERVALS + +// returns inverse of horizontal segment [a,b] + ci +MPIComplex inverse_horizontal (const mpfr_t a, const mpfr_t b, const mpfr_t c) { + MPIComplex rop; + mpfr_t aux; + mpfr_init (aux); + // case c = 0 + if (mpfr_cmp_si (c, 0) == 0) { + if (mpfr_cmp_si (a, 0) < 0 && mpfr_cmp_si (b, 0) > 0) { + mpfr_set_nan (rop.r.a); + mpfr_set_nan (rop.r.b); + mpfr_set_nan (rop.i.a); + mpfr_set_nan (rop.i.b); + } else { + mpfr_si_div (rop.r.a, 1, b, MPFR_RNDD); + mpfr_si_div (rop.r.b, 1, a, MPFR_RNDU); + mpfr_set_si (rop.i.a, 0, MPFR_RNDD); + mpfr_set_si (rop.i.b, 0, MPFR_RNDU); + } + // case c > 0 + } else if (mpfr_cmp_si (c, 0) > 0) { + // case a<0 0) { + mpfr_si_div (rop.i.a, -1, c, MPFR_RNDD); + + if (mpfr_cmpabs (b, a) > 0) + g (rop.i.b, b, c, MPFR_RNDU); + else + g (rop.i.b, a, c, MPFR_RNDU); + + + if (mpfr_cmpabs (a, c) > 0) { + mpfr_mul_si (aux, c, 2, MPFR_RNDD); + mpfr_si_div (rop.r.a, -1, aux, MPFR_RNDD); + } else + f (rop.r.a, a, c, MPFR_RNDD); + + + if (mpfr_cmp (c, b) < 0) { + mpfr_mul_si (aux, c, 2, MPFR_RNDU); + mpfr_si_div (rop.r.b, 1, aux, MPFR_RNDU); + } else { + f (rop.r.b, b, c, MPFR_RNDU); + } + } + + // case 0<=a= 0) { + g (rop.i.a, a, c, MPFR_RNDD); + + g (rop.i.b, b, c, MPFR_RNDU); + + f (aux, a, c, MPFR_RNDD); + f (rop.r.a, b, c, MPFR_RNDD); + mpfr_min (rop.r.a, rop.r.a, aux, MPFR_RNDD); + + if (mpfr_cmp (b, c) < 0) { + f (rop.r.b, b, c, MPFR_RNDU); + } else if (mpfr_cmp (c, a) < 0) { + f (rop.r.b, a, c, MPFR_RNDU); + } else { + mpfr_mul_si (rop.r.b, c, 2, MPFR_RNDU); + mpfr_si_div (rop.r.b, 1, rop.r.b, MPFR_RNDU); + } + } + + // case a 0) { + f (rop.r.a, a, c, MPFR_RNDD); + } else if (mpfr_cmpabs (b, c) > 0) { + f (rop.r.a, b, c, MPFR_RNDD); + } else { + mpfr_mul_si (rop.r.a, c, 2, MPFR_RNDD); + mpfr_si_div (rop.r.a, -1, rop.r.a, MPFR_RNDD); + } + } + // case c < 0 + } else { + mpfr_neg (aux, c, MPFR_RNDN); // ??? + rop = inverse_horizontal (a, b, aux); + rop.i = -rop.i; + } + mpfr_clear (aux); + return rop; + +} + +// returns inverse of verticla segment c + [a,b]i +MPIComplex inverse_vertical (const mpfr_t a, const mpfr_t b, const mpfr_t c) { + mpfr_t aux1, aux2; + mpfr_inits (aux1, aux2, NULL); + mpfr_neg (aux1, b, MPFR_RNDD); + mpfr_neg (aux2, a, MPFR_RNDU); + MPIComplex rop; + rop = inverse_horizontal (aux1, aux2, c); + + MPInterval aux (rop.i); + rop.i = rop.r; + rop.r = -aux; + + mpfr_clears (aux1, aux2, NULL); + return rop; + + +} +void f (mpfr_t rop, const mpfr_t x, const mpfr_t c, mpfr_rnd_t roundMode) { + mpfr_t aux; + mpfr_init(aux); + mpfr_mul (aux, c, c, roundMode); + mpfr_mul (rop, x, x, roundMode); + mpfr_add (rop, rop, aux, roundMode); + mpfr_div (rop, x, rop, roundMode); + mpfr_clear (aux); +} +void g (mpfr_t rop, const mpfr_t x, const mpfr_t c, mpfr_rnd_t roundMode) { + mpfr_t aux; + mpfr_init (aux); + mpfr_mul (aux, c, c, roundMode); + mpfr_mul (rop, x, x, roundMode); + mpfr_add (rop, rop, aux, roundMode); + mpfr_neg (aux, c, roundMode); + mpfr_div (rop, aux, rop, roundMode); + mpfr_clear (aux); +} diff --git a/lib/mp_interval.cpp b/lib/mp_interval.cpp new file mode 100644 index 0000000..22f5062 --- /dev/null +++ b/lib/mp_interval.cpp @@ -0,0 +1,559 @@ +#include "mp_interval.hpp" + +// MACROS FOR MIN AND MAX + +#define MIN(a,b) (((a) < (b))? (a) : (b)) +#define MAX(a,b) (((a) > (b))? (a) : (b)) + +// DEFAULT CONSTRUCTOR +MPInterval::MPInterval (double op) { + mpfr_init_set_d (a, op, MPFR_RNDD); + mpfr_init_set_d (b, op, MPFR_RNDU); +} +// CONSTRUCTOR BY double ELEMENTS +MPInterval::MPInterval (double op1, double op2) { + mpfr_init_set_d (a, op1, MPFR_RNDD); + mpfr_init_set_d (b, op2, MPFR_RNDU); +} +// CONSTRUCTOR BY string VALUE +MPInterval::MPInterval (const char *op) { + mpfr_init_set_str (a, op, 10, MPFR_RNDD); + mpfr_init_set_str (b, op, 10, MPFR_RNDU); +} +// CONSTRUCTOR BY string PAIR OF VALUES +MPInterval::MPInterval (const char *op1, const char *op2) { + mpfr_init_set_str (a, op1, 10, MPFR_RNDD); + mpfr_init_set_str (b, op2, 10, MPFR_RNDU); + +} +// CONSTRUCTOR BY mpfr ELEMENT +MPInterval::MPInterval (const mpfr_t op) { + mpfr_init_set (a, op, MPFR_RNDD); + mpfr_init_set (b, op, MPFR_RNDU); +} +// CONSTRUCTOR BY mpfr PAIR OF ELEMENTS +MPInterval::MPInterval (const mpfr_t op1, const mpfr_t op2) { + mpfr_init_set (a, op1, MPFR_RNDD); + mpfr_init_set (b, op2, MPFR_RNDU); +} +// CONSTRUCTOR BY double Invterval +MPInterval::MPInterval (const Interval &op) { + mpfr_init_set_d (a, op.a, MPFR_RNDD); + mpfr_init_set_d (b, op.b, MPFR_RNDU); +} +// COPY CONSTRUCTOR +MPInterval::MPInterval (const MPInterval &op) { + mpfr_init_set (a, op.a, MPFR_RNDD); + mpfr_init_set (b, op.b, MPFR_RNDU); +} +// DESTRUCTOR +MPInterval::~MPInterval () { + mpfr_clears (a, b, NULL); +} + + + +int MPInterval::operator== (double op) const { + return (mpfr_cmp_d(this->a, op) == 0) && (mpfr_cmp_d(this->b, op) == 0); +} +int MPInterval::operator== (const Interval &op) const { + return (mpfr_cmp_d (this->a, op.a) == 0) && (mpfr_cmp_d (this->b, op.b) == 0); +} +int MPInterval::operator== (const MPInterval &op) const { + return (mpfr_cmp(this->a, op.a) == 0) && (mpfr_cmp(this->b, op.b) == 0); +} +int MPInterval::operator!= (double op) const { + return (mpfr_cmp_d(this->a, op) != 0) || (mpfr_cmp_d(this->b, op) != 0); +} +int MPInterval::operator!= (const Interval &op) const { + return (mpfr_cmp_d (this->a, op.a) != 0) || (mpfr_cmp_d (this->b, op.b) != 0); +} +int MPInterval::operator!= (const MPInterval &op) const { + return (mpfr_cmp(this->a, op.a) != 0) || (mpfr_cmp(this->a, op.b) != 0); +} + + +// ASSIGNMENT +MPInterval & MPInterval::operator= (const MPInterval &op) { + // CHECK WHETER WE ARE ASIGNING SOMETHING OVER ITSELF + if (this == &op) return (*this); + // UPDATE VALUES + mpfr_set (this->a, op.a, MPFR_RNDD); + mpfr_set (this->b, op.b, MPFR_RNDU); +} +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); +} +MPInterval & MPInterval::operator= (const mpfr_t op) { + // UPDATE VALUES + mpfr_set (this->a, op, MPFR_RNDD); + mpfr_set (this->b, op, MPFR_RNDU); +} +MPInterval & MPInterval::operator= (double op) { + // UPDATE VALUES + mpfr_set_d (this->a, op, MPFR_RNDD); + mpfr_set_d (this->b, op, MPFR_RNDU); +} + + +// ARITHMETIC + +// [a,b] + [c,d] = [a+c, b+d] +MPInterval MPInterval::operator+ (const MPInterval &op) const { + MPInterval rop (*this); + mpfr_add (rop.a, rop.a, op.a, MPFR_RNDD); + mpfr_add (rop.b, rop.b, op.b, MPFR_RNDU); + return rop; +} + +// [a,b] + [c,d] = [a+c, b+d] +MPInterval MPInterval::operator+ (const Interval &op) const { + MPInterval rop (*this); + mpfr_add_d (rop.a, rop.a, op.a, MPFR_RNDD); + mpfr_add_d (rop.b, rop.b, op.b, MPFR_RNDU); + return rop; +} + +// [a,b] + c = [a+c, b+c] +MPInterval MPInterval::operator+ (const mpfr_t op) const { + MPInterval rop (*this); + mpfr_add (rop.a, rop.a, op, MPFR_RNDD); + mpfr_add (rop.b, rop.b, op, MPFR_RNDU); + return rop; +} + +// [a,b] + c = [a+c, b+c] +MPInterval MPInterval::operator+ (double op) const { + MPInterval rop (*this); + mpfr_add_d (rop.a, rop.a, op, MPFR_RNDD); + mpfr_add_d (rop.b, rop.b, op, MPFR_RNDU); + return rop; +} + +// CHANGE OF SIGN OF AN INTERVAL +// -[a,b] = [-b,a] +MPInterval MPInterval::operator- () const { + MPInterval rop; + mpfr_neg (rop.a, this->b, MPFR_RNDD); + mpfr_neg (rop.b, this->a, MPFR_RNDU); + return rop; +} + +// [a,b] - [c,d] = [a-d, b-c] +MPInterval MPInterval::operator- (const MPInterval &op) const { + MPInterval rop(*this); + mpfr_sub (rop.a, rop.a, op.b, MPFR_RNDD); + mpfr_sub (rop.b, rop.b, op.a, MPFR_RNDU); + return rop; +} +// [a,b] - [c,d] = [a-d, b-c] +MPInterval MPInterval::operator- (const Interval &op) const { + MPInterval rop(*this); + mpfr_sub_d (rop.a, rop.a, op.b, MPFR_RNDD); + mpfr_sub_d (rop.b, rop.b, op.a, MPFR_RNDU); + return rop; +} + +// [a,b] - c = [a-c, b-c] +MPInterval MPInterval::operator- (const mpfr_t op) const { + MPInterval rop(*this); + mpfr_sub (rop.a, rop.a, op, MPFR_RNDD); + mpfr_sub (rop.b, rop.b, op, MPFR_RNDU); + return rop; +} + +// [a,b] - c = [a-c, b-c] +MPInterval MPInterval::operator- (double op) const { + MPInterval rop(*this); + mpfr_sub_d (rop.a, rop.a, op, MPFR_RNDD); + mpfr_sub_d (rop.b, rop.b, op, MPFR_RNDU); + return rop; +} + +// [a,b]*[c,d] = [min(a*c,a*d,b*c,b*d), max(a*c,a*d,b*c,b*d)] +MPInterval MPInterval::operator* (const MPInterval &op) const { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_mul (rop.a, this->a, op.a, MPFR_RNDD); + mpfr_mul (aux, this->a, op.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_mul (aux, this->b, op.a, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_mul (aux, this->b, op.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_mul (rop.b, this->a, op.a, MPFR_RNDU); + mpfr_mul (aux, this->a, op.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_mul (aux, this->b, op.a, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_mul (aux, this->b, op.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; + +} + +// [a,b]*[c,d] = [min(a*c,a*d,b*c,b*d), max(a*c,a*d,b*c,b*d)] +MPInterval MPInterval::operator* (const Interval &op) const { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_mul_d (rop.a, this->a, op.a, MPFR_RNDD); + mpfr_mul_d (aux, this->a, op.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_mul_d (aux, this->b, op.a, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_mul_d (aux, this->b, op.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_mul_d (rop.b, this->a, op.a, MPFR_RNDU); + mpfr_mul_d (aux, this->a, op.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_mul_d (aux, this->b, op.a, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_mul_d (aux, this->b, op.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; +} + + +// [a,b]*c = [min(a*c,b*c), max(a*c,b*c)] +MPInterval MPInterval::operator* (const mpfr_t op) const { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_mul (rop.a, this->a, op, MPFR_RNDD); + mpfr_mul (aux, this->b, op, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_mul (rop.b, this->a, op, MPFR_RNDU); + mpfr_mul (aux, this->b, op, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; +} + +// [a,b]*c = [min(a*c,b*c), max(a*c,b*c)] +MPInterval MPInterval::operator* (double op) const { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_mul_d (rop.a, this->a, op, MPFR_RNDD); + mpfr_mul_d (aux, this->b, op, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_mul_d (rop.b, this->a, op, MPFR_RNDU); + mpfr_mul_d (aux, this->b, op, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; +} + + + +// [a,b]/[c,d] = [min(a/c,a/d,b/c,b/d), max(a/c,a/d,b/c,b/d)] +MPInterval MPInterval::operator/ (const MPInterval &op) const { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_div (rop.a, this->a, op.a, MPFR_RNDD); + mpfr_div (aux, this->a, op.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_div (aux, this->b, op.a, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_div (aux, this->b, op.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_div (rop.b, this->a, op.a, MPFR_RNDU); + mpfr_div (aux, this->a, op.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_div (aux, this->b, op.a, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_div (aux, this->b, op.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; +} + +// [a,b]/[c,d] = [min(a/c,a/d,b/c,b/d), max(a/c,a/d,b/c,b/d)] +MPInterval MPInterval::operator/ (const Interval &op) const { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_div_d (rop.a, this->a, op.a, MPFR_RNDD); + mpfr_div_d (aux, this->a, op.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_div_d (aux, this->b, op.a, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_div_d (aux, this->b, op.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_div_d (rop.b, this->a, op.a, MPFR_RNDU); + mpfr_div_d (aux, this->a, op.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_div_d (aux, this->b, op.a, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_div_d (aux, this->b, op.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; +} + +// [a,b]/c = [min(a/c,b/c), max(a/c,b/c)] +MPInterval MPInterval::operator/ (const mpfr_t op) const { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_div (rop.a, this->a, op, MPFR_RNDD); + mpfr_div (aux, this->b, op, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_div (rop.b, this->a, op, MPFR_RNDU); + mpfr_div (aux, this->b, op, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; +} + +// [a,b]/c = [min(a/c,b/c), max(a/c,b/c)] +MPInterval MPInterval::operator/ (double op) const { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_div_d (rop.a, this->a, op, MPFR_RNDD); + mpfr_div_d (aux, this->b, op, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_div_d (rop.b, this->a, op, MPFR_RNDU); + mpfr_div_d (aux, this->b, op, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; +} + + + +// MISC +// [a,b] ^ -1 = [1/b, 1/a] +MPInterval MPInterval::inverse () const { + MPInterval rop; + + // CHECK WHETER 0 IS IN THE INTERVAL + if (this->containsZero ()) { + mpfr_si_div (rop.a, 1, this->b, MPFR_RNDD); + mpfr_si_div (rop.b, 1, this->a, MPFR_RNDU); + } else { + mpfr_set_nan (rop.a); + mpfr_set_nan (rop.b); + } + + return rop; + +} + +// CHECK WHETER 0 IS IN THE INTERVAL +int MPInterval::containsZero () const { + if (mpfr_cmp_si (this->a, 0) <= 0 && mpfr_cmp_si (this->b, 0) >= 0) + return 1; + return 0; +} + +// CHECK WHETER A VALUE IS IN THE INTERVAL +int MPInterval::contains (double op) const { + if (mpfr_cmp_d (this->a, op) <= 0 && mpfr_cmp_d (this->b, op) >= 0) + return 1; + return 0; +} + + +// CHECK WHETER A VALUE IS IN THE INTERVAL +int MPInterval::contains (const mpfr_t op) const { + if (mpfr_cmp (this->a, op) <= 0 && mpfr_cmp (this->b, op) >= 0) + return 1; + return 0; +} + +// CHECKS WHETER THE INTERVAL [a,b] IS SUBSET OF ANOTHER [c,d]. +// ca, op.a) > 0 && mpfr_cmp (this->b, op.b) < 0) + return 1; + return 0; +} + +// CHECKS WHETER THE INTERVAL [a,b] IS SUBSET OF ANOTHER [c,d]. +// ca, op.a) > 0 && mpfr_cmp_d (this->b, op.b) < 0) + return 1; + return 0; +} + + + +// EXTERNAL FUNCTIONS +// OVERLOAD OUTPUT STREAM +std::ostream & operator<< (std::ostream &output, const MPInterval &op) { + mpfr_exp_t exponent; + char *string = mpfr_get_str (NULL, &exponent, 10, 0, op.a, MPFR_RNDD); + if (mpfr_signbit (op.a) != 0) + output << "[-" << string[1] << "." << string+2 << "e" << exponent-1 << ", "; + else + output << "[" << string[0] << "." << string+1 << "e" << exponent-1 << ", "; + mpfr_free_str (string); + + string = mpfr_get_str (NULL, &exponent, 10, 0, op.b, MPFR_RNDU); + if (mpfr_signbit (op.b) != 0) + output << "-" << string[1] << "." << string+2 << "e" << exponent-1 << "]"; + else + output << string[0] << "." << string+1 << "e" << exponent-1 << "]"; + mpfr_free_str (string); + return output; +} + +// c + [a,b] = [a+c, b+c] = [a,b] + c +MPInterval operator+ (double op1, const MPInterval &op2) { + return op2 + op1; +} +// c + [a,b] = [a+c, b+c] = [a,b] + c +MPInterval operator+ (const mpfr_t op1, const MPInterval &op2) { + return op2 + op1; +} +// [c,d] + [a,b] = [a+c, b+c] = [a,b] + [c,d] +MPInterval operator+ (const Interval &op1, const MPInterval &op2) { + return op2 + op1; +} + +// c - [a,b] = [c-b, c-a] +MPInterval operator- (double op1, const MPInterval &op2) { + MPInterval rop; + mpfr_d_sub (rop.a, op1, op2.b, MPFR_RNDD); + mpfr_d_sub (rop.b, op1, op2.a, MPFR_RNDU); + return rop; +} + +// c - [a,b] = [c-b, c-a] +MPInterval operator- (const mpfr_t op1, const MPInterval &op2) { + MPInterval rop; + mpfr_sub (rop.a, op1, op2.b, MPFR_RNDD); + mpfr_sub (rop.b, op1, op2.a, MPFR_RNDU); + return rop; +} + +// [c,d]- [a,b] = [c-b, d-a] +MPInterval operator- (const Interval &op1, const MPInterval &op2) { + MPInterval rop; + mpfr_d_sub (rop.a, op1.a, op2.b, MPFR_RNDD); + mpfr_d_sub (rop.b, op1.b, op2.a, MPFR_RNDU); + return rop; +} + +// c*[a,b] = [a,b]*c +MPInterval operator* (double op1, const MPInterval &op2) { + return op2 * op1; +} + +// c*[a,b] = [a,b]*c +MPInterval operator* (const mpfr_t op1, const MPInterval &op2) { + return op2 * op1; +} +// [c,d] * [a,b] = [a,b] * [c,d] +MPInterval operator* (const Interval &op1, const MPInterval &op2) { + return op2 * op1; +} + + +// c / [a,b] = [min (c/a, c/b), max (c/a, c/b)] +MPInterval operator/ (double op1, const MPInterval &op2) { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_d_div (rop.a, op1, op2.a, MPFR_RNDD); + mpfr_d_div (aux, op1, op2.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_d_div (rop.b, op1, op2.a, MPFR_RNDU); + mpfr_d_div (aux, op1, op2.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; +} + +// c / [a,b] = [min (c/a, c/b), max (c/a, c/b)] +MPInterval operator/ (const mpfr_t op1, const MPInterval &op2) { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_div (rop.a, op1, op2.a, MPFR_RNDD); + mpfr_div (aux, op1, op2.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_div (rop.b, op1, op2.a, MPFR_RNDU); + mpfr_div (aux, op1, op2.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; +} + +// [a,b]/[c,d] = [min(a/c,a/d,b/c,b/d), max(a/c,a/d,b/c,b/d)] +MPInterval operator/ (const Interval &op1, const MPInterval &op2) { + mpfr_t aux; + mpfr_init (aux); + MPInterval rop; + + mpfr_d_div (rop.a, op1.a, op2.a, MPFR_RNDD); + mpfr_d_div (aux, op1.a, op2.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_d_div (aux, op1.b, op2.a, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + mpfr_d_div (aux, op1.b, op2.b, MPFR_RNDD); + mpfr_min (rop.a, rop.a, aux, MPFR_RNDD); + + + mpfr_d_div (rop.b, op1.a, op2.a, MPFR_RNDU); + mpfr_d_div (aux, op1.a, op2.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_d_div (aux, op1.b, op2.a, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + mpfr_d_div (aux, op1.b, op2.b, MPFR_RNDU); + mpfr_max (rop.b, rop.b, aux, MPFR_RNDU); + + mpfr_clear (aux); + return rop; + +} diff --git a/lib/polynomial.cpp b/lib/polynomial.cpp new file mode 100644 index 0000000..19f66b9 --- /dev/null +++ b/lib/polynomial.cpp @@ -0,0 +1,297 @@ +#include "polynomial.hpp" + + +unsigned int *_Polynomial::combinatorial; + + +void _Polynomial::startCombinatorialNumbers (int degree) { + _Polynomial::combinatorial = new unsigned int [((degree+2)*(degree+1))/2]; + unsigned int *combinatorial = _Polynomial::combinatorial; + unsigned int i; // counter for m + unsigned int j; // counter for n + + for (i=0; i<=degree; i++) { + combinatorial[(i*(i+1))/2] = combinatorial[(i*(i+3))/2] = 1; // FIRST AND LAST = 1 + + // THE REST IS COMPUTED FROM THE RECURRENCE OF COMBINATORIAL NUMBERS. + for (j=1; j +IComplex Polynomial::evalIPolHornerXY (const IComplex &x, const IComplex &y) const { + + //evaluate polynomial as an element of K[x][y] in Horner form + //that is, each coefficient of y is a polynomial on x + int i; //i will represent the on x + int j; //j will represent the degree on y + int d = this->degree; + IComplex coef; + IComplex rop; + rop = (this->coef)[(d+2)*(d+1)/2-1]; + for (j=d-1; j>=0; --j){ + coef = (this->coef)[d*(d+1)/2+j]; + for (i=d-1-j; i>=0; --i){ + coef = x * coef + (this->coef)[(j+i)*(j+i+1)/2+j]; + } + rop = rop * y + coef; + } + return rop; +} + + + +template <> +IComplex Polynomial::evalIPolHornerYX (const IComplex &x, const IComplex &y) const { + //evaluate polynomial as an element of K[y][x] in Horner form + //that is, each coefficient of x is a polynomial on y + int i; //i will represent the degree on x + int j; //j will represent the on y + int d = this->degree; + IComplex coef; + IComplex rop; + rop = (this->coef)[(d)*(d+1)/2]; + for (i=d-1; i>=0; --i){ + coef = (this->coef)[d*(d+1)/2+d-i]; + for (j=d-i-1; j>=0; --j) { + coef = coef * y + (this->coef)[(i+j)*(i+j+1)/2+j]; + } + rop = rop * x + coef; + } + + return rop; +} + +template <> +IComplex Polynomial::operator () (const IComplex &x, const IComplex &y) const{ + IComplex rop, ropx, ropy; + + ropx = this->evalIPolHornerXY (x, y); + ropy = this->evalIPolHornerYX (x, y); + rop = this->evalPolClassic (x, y); + + + rop.r.a = MAX(MAX(ropx.r.a, ropy.r.a), rop.r.a); + rop.r.b = MIN(MIN(ropx.r.b, ropy.r.b), rop.r.b); + rop.i.a = MAX(MAX(ropx.i.a, ropy.i.a), rop.i.a); + rop.i.b = MIN(MIN(ropx.i.b, ropy.i.b), rop.i.b); + + return rop; +} + + + +template <> +IComplex Polynomial::evalIPolYHornerXY (const IComplex &x, const IComplex &y) const { + //evaluate polynomial as an element of K[x][y] in Horner form + //that is, each coefficient of y is a polynomial on x + int i; //i will represent the on x + int j; //j will represent the degree on y + int d = this->degree; + IComplex coef; + IComplex rop; + rop = (d) * (this->coef)[(d+2)*(d+1)/2-1]; + + for (j=d-2; j>=0; --j){ + coef = (this->coef)[d*(d+1)/2+j+1]; + for (i=d-1; i>=j+1; --i){ + coef = x * coef + (this->coef)[(i*(i+1))/2+j+1]; + } + rop = rop * y + (j+1.0) * coef; + } + return rop; + +} + +template <> +IComplex Polynomial::evalIPolYHornerYX (const IComplex &x, const IComplex &y) const { + + //evaluate polynomial as an element of K[x][y] in Horner form + //that is, each coefficient of y is a polynomial on x + int i; //i will represent the on x + int j; //j will represent the degree on y + int d = this->degree; + IComplex coef; + IComplex rop; + rop = (this->coef)[(d)*(d+1)/2+1]; + + for (i=d-2; i>=0; --i){ + coef = ((double) d-i)*(this->coef)[(d*(d+1))/2+d-i]; + for (j=d-1; j>=i+1; --j){ + coef = y * coef + (j-i)*(this->coef)[(j*(j+1))/2+j-i]; + } + rop = rop * x + coef; + } + return rop; + +} + + + +template <> +IComplex Polynomial::diffY (const IComplex &x, const IComplex &y) const{ + IComplex rop, ropx, ropy; + + ropx = this->evalIPolYHornerXY (x, y); + ropy = this->evalIPolYHornerYX (x, y); + rop = this->evalPolYClassic (x, y); + + + rop.r.a = MAX(MAX(ropx.r.a, ropy.r.a), rop.r.a); + rop.r.b = MIN(MIN(ropx.r.b, ropy.r.b), rop.r.b); + rop.i.a = MAX(MAX(ropx.i.a, ropy.i.a), rop.i.a); + rop.i.b = MIN(MIN(ropx.i.b, ropy.i.b), rop.i.b); + + return rop; +} + + + +template <> +MPIComplex Polynomial::evalIPolHornerXY (const MPIComplex &x, const MPIComplex &y) const { + + //evaluate polynomial as an element of K[x][y] in Horner form + //that is, each coefficient of y is a polynomial on x + int i; //i will represent the on x + int j; //j will represent the degree on y + int d = this->degree; + MPIComplex coef; + MPIComplex rop; + rop = (this->coef)[(d+2)*(d+1)/2-1]; + for (j=d-1; j>=0; --j){ + coef = (this->coef)[d*(d+1)/2+j]; + for (i=d-1-j; i>=0; --i){ + coef = x * coef + (this->coef)[(j+i)*(j+i+1)/2+j]; + } + rop = rop * y + coef; + } + return rop; +} + + + +template <> +MPIComplex Polynomial::evalIPolHornerYX (const MPIComplex &x, const MPIComplex &y) const { + //evaluate polynomial as an element of K[y][x] in Horner form + //that is, each coefficient of x is a polynomial on y + int i; //i will represent the degree on x + int j; //j will represent the on y + int d = this->degree; + MPIComplex coef; + MPIComplex rop; + rop = (this->coef)[(d)*(d+1)/2]; + for (i=d-1; i>=0; --i){ + coef = (this->coef)[d*(d+1)/2+d-i]; + for (j=d-i-1; j>=0; --j) { + coef = coef * y + (this->coef)[(i+j)*(i+j+1)/2+j]; + } + rop = rop * x + coef; + } + + return rop; +} + +template <> +MPIComplex Polynomial::operator () (const MPIComplex &x, const MPIComplex &y) const{ + MPIComplex rop, ropx, ropy; + + ropx = this->evalIPolHornerXY (x, y); + ropy = this->evalIPolHornerYX (x, y); + rop = this->evalPolClassic (x, y); + + + mpfr_max (rop.r.a, rop.r.a, ropx.r.a, MPFR_RNDD); + mpfr_max (rop.r.a, rop.r.a, ropy.r.a, MPFR_RNDD); + mpfr_max (rop.i.a, rop.i.a, ropx.i.a, MPFR_RNDD); + mpfr_max (rop.i.a, rop.i.a, ropy.i.a, MPFR_RNDD); + + mpfr_min (rop.r.b, rop.r.b, ropx.r.b, MPFR_RNDU); + mpfr_min (rop.r.b, rop.r.b, ropy.r.b, MPFR_RNDU); + mpfr_min (rop.i.b, rop.i.b, ropx.i.b, MPFR_RNDU); + mpfr_min (rop.i.b, rop.i.b, ropy.i.b, MPFR_RNDU); + + // rop.r.a = MAX(MAX(ropx.r.a, ropy.r.a), rop.r.a); + // rop.r.b = MIN(MIN(ropx.r.b, ropy.r.b), rop.r.b); + // rop.i.a = MAX(MAX(ropx.i.a, ropy.i.a), rop.i.a); + // rop.i.b = MIN(MIN(ropx.i.b, ropy.i.b), rop.i.b); + + return rop; +} + + + +template <> +MPIComplex Polynomial::evalIPolYHornerXY (const MPIComplex &x, const MPIComplex &y) const { + //evaluate polynomial as an element of K[x][y] in Horner form + //that is, each coefficient of y is a polynomial on x + int i; //i will represent the on x + int j; //j will represent the degree on y + int d = this->degree; + MPIComplex coef; + MPIComplex rop; + rop = (d) * (this->coef)[(d+2)*(d+1)/2-1]; + + for (j=d-2; j>=0; --j){ + coef = (this->coef)[d*(d+1)/2+j+1]; + for (i=d-1; i>=j+1; --i){ + coef = x * coef + (this->coef)[(i*(i+1))/2+j+1]; + } + rop = rop * y + (j+1.0) * coef; + } + return rop; + +} + +template <> +MPIComplex Polynomial::evalIPolYHornerYX (const MPIComplex &x, const MPIComplex &y) const { + + //evaluate polynomial as an element of K[x][y] in Horner form + //that is, each coefficient of y is a polynomial on x + int i; //i will represent the on x + int j; //j will represent the degree on y + int d = this->degree; + MPIComplex coef; + MPIComplex rop; + rop = (this->coef)[(d)*(d+1)/2+1]; + + for (i=d-2; i>=0; --i){ + coef = ((double) d-i)*(this->coef)[(d*(d+1))/2+d-i]; + for (j=d-1; j>=i+1; --j){ + coef = y * coef + (j-i)*(this->coef)[(j*(j+1))/2+j-i]; + } + rop = rop * x + coef; + } + return rop; + +} + + + +template <> +MPIComplex Polynomial::diffY (const MPIComplex &x, const MPIComplex &y) const{ + MPIComplex rop, ropx, ropy; + + ropx = this->evalIPolYHornerXY (x, y); + ropy = this->evalIPolYHornerYX (x, y); + rop = this->evalPolYClassic (x, y); + + + /*rop.r.a = MAX(MAX(ropx.r.a, ropy.r.a), rop.r.a); + rop.r.b = MIN(MIN(ropx.r.b, ropy.r.b), rop.r.b); + rop.i.a = MAX(MAX(ropx.i.a, ropy.i.a), rop.i.a); + rop.i.b = MIN(MIN(ropx.i.b, ropy.i.b), rop.i.b);*/ + mpfr_max (rop.r.a, rop.r.a, ropx.r.a, MPFR_RNDD); + mpfr_max (rop.r.a, rop.r.a, ropy.r.a, MPFR_RNDD); + mpfr_max (rop.i.a, rop.i.a, ropx.i.a, MPFR_RNDD); + mpfr_max (rop.i.a, rop.i.a, ropy.i.a, MPFR_RNDD); + + mpfr_min (rop.r.b, rop.r.b, ropx.r.b, MPFR_RNDU); + mpfr_min (rop.r.b, rop.r.b, ropy.r.b, MPFR_RNDU); + mpfr_min (rop.i.b, rop.i.b, ropx.i.b, MPFR_RNDU); + mpfr_min (rop.i.b, rop.i.b, ropy.i.b, MPFR_RNDU); + + return rop; +} diff --git a/lib/sirocco.cpp b/lib/sirocco.cpp new file mode 100644 index 0000000..1657a1d --- /dev/null +++ b/lib/sirocco.cpp @@ -0,0 +1,556 @@ + +#include +#include +#include +#include + +using namespace std; + +// unsigned int *_Polynomial::combinatorial; + + + +/****************************************** +// f(z1,z2) = z1^3+2*z1*z2-3*z2-z1; +// COEFS = {{1,0,0,0}, {0,2,0}, {-1,3}, {0}} +// z1 = x + iy +// z2 = t(1+0.01i) +******************************************/ + +#pragma STDC FENV_ACCESS ON + +// PERFORM A MINIMUM OF 5 ITERATIONS OF NEWTON METHOD TO CORRECT THE APPROXIMATION OF THE ROOT (SUPPOSED TO BE GOOD ENOUGH) +template <> +void correctRoot (const Polynomial f, const Complex &x0, Complex &y0) { + Complex increment; + int count = 0; + double last_increment, new_increment; // sizes of corrections + new_increment = 9999.0; + do { + last_increment = new_increment; + increment = f(x0,y0) / (f.diffY(x0,y0)); + new_increment = abs (increment); + y0 = y0 - increment; + ++count; + } while ((fabs(last_increment - new_increment) > 1.0e-14 || count < 6 ) && count < 50); + return; +} +template <> +void correctRoot (const Polynomial f, const MPComplex &x0, MPComplex &y0) { + MPComplex increment; + int count = 0; + mpfr_t last_increment, new_increment, tol; // sizes of corrections and tolerance + mpfr_inits (last_increment, new_increment, tol, NULL); + + mpfr_set_d (tol, mpfr_get_default_prec(), MPFR_RNDN); + mpfr_sub_si (tol, tol, 5, MPFR_RNDN); + mpfr_neg (tol, tol, MPFR_RNDN); + mpfr_exp2 (tol, tol, MPFR_RNDN); + // std::cout << "tol = "; + // mpfr_out_str (stdout, 10, 20, tol, MPFR_RNDN); + // std::cout << std::endl; + + mpfr_set_d (new_increment, 9999.0, MPFR_RNDN); + do { + mpfr_set (last_increment, new_increment, MPFR_RNDN); + increment = f(x0,y0) / (f.diffY(x0,y0)); + abs (new_increment, increment); + y0 = y0 - increment; + ++count; + + // compute difference of increments and store in "last_increment" + mpfr_sub (last_increment, last_increment, new_increment, MPFR_RNDN); + mpfr_abs (last_increment, last_increment, MPFR_RNDN); + } while ((mpfr_cmp (last_increment, tol) > 0 || count < 6 ) && count < 50); + + mpfr_clears (last_increment, new_increment, tol, NULL); + return; +} + +double * homotopyPath (int degree, double *_coef, double _y0R, double _y0I) { + + _Polynomial::startCombinatorialNumbers (degree); + +// /************************************************/ +// /* 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 > outputList; + + + + // INITIALIZE COEFICIENT DATA FROM ARGUMENT INTO APPROPIATE DATA STRUCTURES + for (int i=0; i f (degree, coef); // complex polynomials (with fp arithmetic) + // to be used in Newton method correction + Polynomial If (degree, Icoef); // complex polynomials (with Interval arithmetic) + Polynomial Ig (degree); // to be used in validated homotopy + Polynomial Ih (degree); + + + + // 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 (0.0, y0.real (), y0.imag ())); + + /*******************/ + /* CORRECT ROOT */ + /*******************/ + + + correctRoot (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 (If, Ix0, Iy0, IY0)) { +#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 (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); + + + while (!validate (Ih, Ix1, Iy0, IY0)) { +#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 (Ih, Ix1, Iy0, IY0)); + } + + +// /********************/ +// /* 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 (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 (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; + + + + // INITIALIZE COEFICIENT DATA FROM ARGUMENT INTO APPROPRIATE DATA STRUCTURES + for (int i=0; i f (degree, coef); // complex polynomials (with fp arithmetic) + // to be used in Newton method correction + Polynomial If (degree, Icoef); // complex polynomials (with Interval arithmetic) + Polynomial Ig (degree); // to be used in validated homotopy + Polynomial Ih (degree); + + + +// // 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 (x0.r, y0.r, y0.i)); + +// /*******************/ +// /* CORRECT ROOT */ +// /*******************/ + + + correctRoot (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 (If, Ix0, Iy0, IY0)) { +#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 (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); + + + while (!validate (Ih, Ix1, Iy0, IY0)) { +#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 (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 (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 sage_malloc(sizeof(mpfr_t)*len(values)) + cdef mpfr_t* rop + cdef int j + for i in range(len(values)): + j = i + sig_on() + mpfr_init2(cvalues[j], prec) + mpfr_set(cvalues[j], (values[j]).value, MPFR_RNDN) + sig_off() + cdef mpfr_t y0R + cdef mpfr_t y0I + sig_on() + mpfr_init2(y0R, prec) + mpfr_set(y0R, (y0r).value, MPFR_RNDN) + mpfr_init2(y0I, prec) + mpfr_set(y0I, (y0i).value, MPFR_RNDN) + rop = homotopyPath_mp(cdeg, cvalues, y0R, y0I, prec) + sig_off() + for i in range(len(values)): + j = i + sig_on() + mpfr_clear(cvalues[j]) + sig_off() + free(cvalues) + 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((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 = sage_malloc(sizeof(double)*len(values)) + cdef int clen = 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 diff --git a/tests/Makefile.am b/tests/Makefile.am new file mode 100644 index 0000000..099dd5c --- /dev/null +++ b/tests/Makefile.am @@ -0,0 +1,11 @@ +check_PROGRAMS = sirocco_test sirocco_minimal_test +sirocco_test_SOURCES = sirocco_test.cpp +sirocco_test_LDADD = @top_builddir@/lib/libsirocco.la +sirocco_test_DEPENDENCIES = @top_builddir@/lib/libsirocco.la +sirocco_test_LDFLAGS = -lm -lmpfr -lgmp +sirocco_minimal_test_SOURCES = sirocco_minimal_test.cpp +sirocco_minimal_test_LDADD = @top_builddir@/lib/libsirocco.la +sirocco_minimal_test_DEPENDENCIES = @top_builddir@/lib/libsirocco.la +sirocco_minimal_test_LDFLAGS = -lm -lmpfr -lgmp +TESTS = sirocco_test sirocco_minimal_test +AM_CPPFLAGS = -I$(top_srcdir)/include diff --git a/tests/sirocco_minimal_test.cpp b/tests/sirocco_minimal_test.cpp new file mode 100644 index 0000000..bd45ac2 --- /dev/null +++ b/tests/sirocco_minimal_test.cpp @@ -0,0 +1,67 @@ +#include +#include + + +int main(int argc, char const *argv[]) { + mpfr_set_default_prec (162); + + + double testCoefs[] = {0.000000000000000, -0.000000000000000, 0.000000000000000, +-0.000000000000000, 0.000000000000000, -0.000000000000000, +0.000000000000000, -0.000000000000000, -128.000000000000, +-128.000000000000, -96.0000000000000, -96.0000000000000, +0.000000000000000, -0.000000000000000, 0.000000000000000, +-0.000000000000000, 72.0000000000000, 72.0000000000000, +304.000000000000, 304.000000000000, 12.0000000000000, +12.0000000000000, 180.000000000000, 180.000000000000, +0.000000000000000, -0.000000000000000, 0.000000000000000, +-0.000000000000000, 90.0000000000000, 90.0000000000000, +-182.500000000000, -182.500000000000, 133.500000000000, +133.500000000000, -190.500000000000, -190.500000000000, +42.0000000000000, 42.0000000000000, -44.0000000000000, +-44.0000000000000, 0.000000000000000, -0.000000000000000, +0.000000000000000, -0.000000000000000, -40.2500000000000, +-40.2500000000000, 18.2500000000000, 18.2500000000000, +-72.0000000000000, -72.0000000000000, 21.0000000000000, +21.0000000000000, -38.5000000000000, -38.5000000000000, +5.50000000000000, 5.50000000000000, -6.00000000000000, +-6.00000000000000, -0.000000000000000, 0.000000000000000}; + + double y0R = 4.0; + double y0I = -4.0; + double *rop = homotopyPath (4, testCoefs, y0R, y0I); + + int len = (int) rop[0]; + int i; + + delete [] rop; + + mpfr_t mp_testCoefs[60]; + for (int i=0; i<60; ++i) { + mpfr_init_set_d (mp_testCoefs[i], testCoefs[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 (4, mp_testCoefs, mp_y0R, mp_y0I, 5300); + + 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; + + + for (int i=0; i<60; ++i) { + mpfr_clear (mp_testCoefs[i]); + } + + mpfr_clear (mp_y0R); + mpfr_clear (mp_y0I); + + + + return 0; +} diff --git a/tests/sirocco_test.cpp b/tests/sirocco_test.cpp new file mode 100644 index 0000000..6ee227a --- /dev/null +++ b/tests/sirocco_test.cpp @@ -0,0 +1,95 @@ +#include +#include + + + +using namespace std; + + +std::ostream & operator<< (std::ostream &output, const std::complex op) { + if (op == 0.0) + output << "0"; + else if (op.imag() > 0) + if (op.real () != 0.0) + output << op.real () << "+" << op.imag() << "i"; + else + output << op.imag() << "i"; + else if (op.imag() == 0) + output << op.real (); + else + if (op.real () != 0.0) + output << op.real () << op.imag() << "i"; + else + output << op.imag() << "i"; + return output; +} + + + +int main(int argc, char const *argv[]) { + mpfr_set_default_prec (162); + + + double testCoefs[] = {0.000000000000000, -0.000000000000000, 0.000000000000000, +-0.000000000000000, 0.000000000000000, -0.000000000000000, +0.000000000000000, -0.000000000000000, -128.000000000000, +-128.000000000000, -96.0000000000000, -96.0000000000000, +0.000000000000000, -0.000000000000000, 0.000000000000000, +-0.000000000000000, 72.0000000000000, 72.0000000000000, +304.000000000000, 304.000000000000, 12.0000000000000, +12.0000000000000, 180.000000000000, 180.000000000000, +0.000000000000000, -0.000000000000000, 0.000000000000000, +-0.000000000000000, 90.0000000000000, 90.0000000000000, +-182.500000000000, -182.500000000000, 133.500000000000, +133.500000000000, -190.500000000000, -190.500000000000, +42.0000000000000, 42.0000000000000, -44.0000000000000, +-44.0000000000000, 0.000000000000000, -0.000000000000000, +0.000000000000000, -0.000000000000000, -40.2500000000000, +-40.2500000000000, 18.2500000000000, 18.2500000000000, +-72.0000000000000, -72.0000000000000, 21.0000000000000, +21.0000000000000, -38.5000000000000, -38.5000000000000, +5.50000000000000, 5.50000000000000, -6.00000000000000, +-6.00000000000000, -0.000000000000000, 0.000000000000000}; + + double y0R = 4.0; + double y0I = -4.0; + double *rop = homotopyPath (4, testCoefs, y0R, y0I); + + int len = (int) rop[0]; + int i; + for (i=0; i