/* 
   ***************************************************************************
   <vpoisson.h> is a C++ class to manage Poisson Series
    Copyright (C) 2018. Vicente Agost Gomez & Jose Antonio Lopez Orti

    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, version 3 of the License.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <https://www.gnu.org/licenses/>. 
    **************************************************************************
*/





// This class vpoisson.h is named as class taylor_poisson.h in old versions
// (taylor_poisson.h is now vpoisson.h)



#include <vector>
#include <list>
#include <iostream>
#include <cmath>
#include <iomanip>

#define sPoisson vector<tPoisson>
#define vPoisson vector<vector<tPoisson>>

//#define PI_L acos(-1.L)
#define  PI_L 3.1415926535897932384626433L
#define _2PI_L 6.2831853071795864769252868L
#define ORDEN_TAYLOR 3  // Order of perturbation
#define NP 1            // Number of angular variables
#define TOL 1.0e-14L
#define CERO 1.0e-18L

using   namespace std;

long double  fact(const long double);
long double  binomial(const long double, int);

//__________________________________________________________________________
// CLASE Termino de Poisson

class tPoisson {
    private:
        long double m_A;
        long double m_B;
        int m_i;
        int m_j[NP];
    public:
        tPoisson(void);
        tPoisson(const tPoisson&);     //constructor de copia
        tPoisson(long double, int, int[], long double);
        tPoisson(long double);
        ~tPoisson(void);
        long double getA(void) const 	{return m_A;}
        long double getB(void) const 	{return m_B;}
        int getI(void) const		    {return m_i;}
        void setA(const long double A)	{m_A = A;}
        void setB(const long double B)	{m_B = B;}
        void setI(const int& I)		    {m_i = I;}
        int getJ(int pos) const		    {return m_j[pos];}
        void setJ(int pos, int val)	    {if(pos<NP) m_j[pos] = val;}
        friend bool operator== (const tPoisson&, const tPoisson&);
        friend bool operator!= (const tPoisson&, const tPoisson&);
        friend bool operator>  (const tPoisson&, const tPoisson&);
        friend bool operator>= (const tPoisson&, const tPoisson&);
        friend bool operator<  (const tPoisson&, const tPoisson&);
        friend bool operator<= (const tPoisson&, const tPoisson&);
};

//__________________________________________________________________________
// Cabeceras

bool menorTermino(tPoisson&, tPoisson&);
bool menorTerminoA(tPoisson&, tPoisson&);
bool mayorTerminoA(tPoisson&, tPoisson&);

//__________________________________________________________________________
// Constructores

tPoisson::tPoisson(void): m_A(0.0L), m_i(0), m_B(0.0L){
    //for(int i=0;i<NP;i++)
	for(int i=NP; i--; )
        m_j[i] = 0;
};

tPoisson::tPoisson(const tPoisson& original):
    m_A(original.m_A), m_i(original.m_i), m_B(original.m_B){
    //for(int i=0;i<NP;i++)
	for(int i=NP; i--; )
        m_j[i] = original.m_j[i];
};

tPoisson::tPoisson(long double A, int i, int j[NP], long double B):
    m_A(A), m_i(i), m_B(B){    
	//for(int i=0;i<NP;i++)
	for(int i=NP; i--; )
        m_j[i] = j[i];
};

tPoisson::tPoisson(long double A): m_A(A), m_B(0.0L), m_i(0) {
    //for(int i=0;i<NP;i++)
	for(int i=NP; i--; )
        m_j[i] = 0;
};

tPoisson::~tPoisson (void){};

//__________________________________________________________________________
// Acceso a miembros de la clase   ---> DECLARADOS INLINE
/*
long double tPoisson::getA(void) const {return m_A;};
long double tPoisson::getB(void) const {return m_B;};
int tPoisson::getI(void) const {return m_i;};
void tPoisson::setA(const long double A) {m_A = A;};
void tPoisson::setB(const long double B) {m_B = B;};
void tPoisson::setI(const int& I) {m_i = I;};
int tPoisson::getJ(int pos) const {return m_j[pos];};
void tPoisson::setJ(int pos, int valor) {if(pos<NP) m_j[pos] = valor;};
*/
//__________________________________________________________________________
// Operadores de comparacion

bool operator> (const tPoisson &tp1, const tPoisson &tp2){
    //OJO: Un term es mayor que otro si sus variables angulares son mayores
    if(tp1.m_i != tp2.m_i)
        return (tp1.m_i > tp2.m_i);
    for(int i=0;i<NP;i++)
        if(tp1.m_j[i] != tp2.m_j[i])
            return(tp1.m_j[i] > tp2.m_j[i]);
    return false; //si son iguales no es mayor
};

bool operator>= (const tPoisson& tp1, const tPoisson& tp2) {
    return !(tp1<tp2);
};

bool operator< (const tPoisson &tp1, const tPoisson &tp2){
    //OJO: Un term es menor que otro si sus variables angulares son menores
    if(tp1.m_i != tp2.m_i)
        return (tp1.m_i < tp2.m_i);
    for(int i=0;i<NP;i++)
        if(tp1.m_j[i] != tp2.m_j[i])
            return(tp1.m_j[i] < tp2.m_j[i]);
    return false; //si son iguales no es mayor
};

bool operator<= (const tPoisson& tp1, const tPoisson& tp2){
    return !(tp1>tp2);
};

bool operator== (const tPoisson &tp1, const tPoisson &tp2){
    //OJO: Igualdad de terminos indica solo mismas variables angulares
    if(tp1.m_i != tp2.m_i)
        return false;
    for(int i=0;i<NP;i++)
        if(tp1.m_j[i] != tp2.m_j[i])
            return false;
    return true;
};

bool operator!= (const tPoisson& tp1, const tPoisson& tp2){
    return !(tp1==tp2);
};

bool menorTermino(tPoisson &tp1, tPoisson &tp2){
    return (tp1 < tp2);
};

bool menorTerminoA(tPoisson &tp1, tPoisson &tp2){
    return (fabs(tp1.getA()) < fabs(tp2.getA()));
};

bool mayorTerminoA(tPoisson &tp1, tPoisson &tp2){
    return (fabs(tp1.getA()) > fabs(tp2.getA()));
};

//__________________________________________________________________________
// Declaraciones series de Poisson

long double norma_1(const sPoisson&, long double);
long double norma_1(const vPoisson&, long double);
long double norma_2(const sPoisson&, long double);
void arregla(tPoisson&);
void arregla(sPoisson&);
void show(const vPoisson&, int = 0);
void show(const sPoisson&, int = 0);
void show_sincos(const sPoisson&);
long double valor(const sPoisson&, long double =0.L);
long double valor(const sPoisson&, long double, long double[NP]);
void _multiplica(const sPoisson&, const sPoisson&, sPoisson&);
void _multiplica(const long double x, const vPoisson& v1n, vPoisson& aux);
void _multiplica(const sPoisson& s1, const vPoisson& v1n, vPoisson& aux);
void _multiplica(const vPoisson&, const vPoisson&, vPoisson&);
void _multiplica(const vPoisson&, long double, vPoisson&);
void _suma(const sPoisson&, const sPoisson&, sPoisson&);
void _suma(const vPoisson&, const vPoisson&, vPoisson&);
sPoisson& operator*(const sPoisson&,const sPoisson&);
sPoisson& operator*(const sPoisson&, const long double);
sPoisson& operator*(long double, const sPoisson&);
vPoisson& operator*(const long double, const vPoisson&);
vPoisson& operator*(const sPoisson&, const vPoisson&);
vPoisson& operator*(const vPoisson&, const vPoisson&);
sPoisson& operator/(const sPoisson&, long double);
sPoisson& operator+(const sPoisson&, const sPoisson&);
sPoisson& operator+(const sPoisson&, long double);
sPoisson& operator+(long double, const sPoisson&);
vPoisson& operator+(const vPoisson&, const vPoisson&);
sPoisson& operator-(const sPoisson&, const sPoisson&);
sPoisson& operator-(const sPoisson&);
vPoisson& operator-(const vPoisson&);
sPoisson& operator-(long double, const sPoisson&);
sPoisson& operator-(const sPoisson&, long double);
sPoisson& operator-(const sPoisson&, long double);
vPoisson& operator-(const vPoisson& , const vPoisson& );
sPoisson& compacta(sPoisson&);
vPoisson& compacta(vPoisson&);
sPoisson& sin(const sPoisson&);
vPoisson& sin(const vPoisson&);
vPoisson& sin2(const vPoisson&);
sPoisson& cos(const sPoisson&);
vPoisson& cos(const vPoisson&);
sPoisson& log(const sPoisson&);
vPoisson& log(const vPoisson&);
sPoisson& exp(const sPoisson&);
vPoisson& exp(const vPoisson&);
sPoisson& raiz(const sPoisson&, long double = 2.L);
vPoisson& raiz(const vPoisson&, long double = 2.L);
sPoisson& pow(const sPoisson&, int);
vPoisson& pow(const vPoisson&, int);

void ordenar_por_angulo(sPoisson&);
void ordenar_por_amplitud(sPoisson&);
void ordenar_por_amplitud_decreciente(sPoisson&);
sPoisson& trunca(const sPoisson&, long double = TOL);
void trunca(const sPoisson&, sPoisson&, long double = TOL);
sPoisson& chop(const sPoisson& ,long double);
void show_sincos_sec(const sPoisson& sp);
void show_sec(const sPoisson& sp);
void show_sec_char(const sPoisson &sp, int C);

//__________________________________________________________________________
// Implementaciones

long double norma_1(const sPoisson &sp1, long double T){
    long double norma = 0.0L;

    for(auto const& tp: sp1)
        if(tp.getI()==0)
            norma+=fabs( tp.getA() );
        else
            norma+=fabs( tp.getA()) * pow(fabs(T), (long double)tp.getI() );

    return norma;
};

long double norma_1(const vPoisson& vp1, long double T){
    long double norma = 0.0L;

    for(auto const& sp: vp1)
        norma += norma_1(sp, T);

    return norma;
};

long double norma_2(const sPoisson &sp1, long double T){
    long double norma = 0.0L;

    for(auto const& tp: sp1)
        if(tp.getI()==0)
            norma+=pow(fabs( tp.getA()), 2.0L);
        else
            norma+=pow(fabs( tp.getA()) * pow(fabs(T), (long double)tp.getI() ),2.0L);

    return sqrt(norma);
};

void arregla(tPoisson& tp){ 
    //Quitar primeros terminos negativos y todos nulos
    bool todosCero;
    int i,pos;

    for(i=0;i<NP;i++)
    {
       todosCero = true;
       pos = 0;
       while(pos<i && todosCero) todosCero = (tp.getJ(pos++)==0);
       if(todosCero&&tp.getJ(i)<0)
       {
          for(int j=i;j<NP;j++) tp.setJ(j, -tp.getJ(j));
          tp.setB( fmod(-tp.getB(), _2PI_L) );
          break;
       };
    };

    todosCero = true;
    pos = 0;
    while(pos<NP && todosCero) todosCero = (tp.getJ(pos++)==0);
    
	if(todosCero && fabs(tp.getB()) > CERO)
    {
        tp.setA(tp.getA()*cos(tp.getB()));
        tp.setB(0.0L);
    };
    tp.setB(fmod(tp.getB(),_2PI_L));
    while(tp.getB()<0.0L) tp.setB(tp.getB()+_2PI_L);

	if(!todosCero && tp.getA()<CERO)
    {
        tp.setA(-tp.getA());
        tp.setB(fmod(tp.getB()+PI_L,_2PI_L));
    };

    if( fabs(tp.getA())<TOL ) tp.setA(0.L);
    if( fabs(tp.getB())<TOL ) tp.setB(0.L);
};

void arregla(sPoisson& sp){

    sPoisson::iterator tp = sp.begin();

    while( tp != sp.end() )
        if( fabs(tp->getA()) < TOL )
            tp=sp.erase(tp);
        else{
            arregla(*tp);
            ++tp;
        };

    if(sp.empty())
        sp.push_back( tPoisson(0.L) );
};

void show(const sPoisson &sp, int nMostrar /* Default 0 */){

    //Muestra 鷑icamente el n鷐ero de t閞minos que se pasa como segundo par醡etro
    //y si no hay segundo par醡etro muestra todos.

    sPoisson aux(sp);

    if(nMostrar == 0)
        nMostrar = aux.size();
    else
        if(nMostrar>aux.size())
            nMostrar = aux.size();

    cout << "/-------------- Elementos a mostrar: " << nMostrar << " de " << aux.size()<< endl;

    cout.precision( 1+(int)(-log10(TOL)));

    int count=1;
    sPoisson::iterator i = aux.begin();
    while(i != aux.end() && count<=nMostrar) {
        cout << "| " << count++ <<" | "<< i->getA() <<", "<< i->getI() <<", "<< i->getB() << ", [";
        for(int j=0;j<NP;j++){
            cout << i->getJ(j);
            if(j!=NP-1) cout << ",";
        }
        cout << "]" << endl;
        ++i;
    };
    cout << "\\-------------------------------------------" << endl << endl;
};

void show_sec(const sPoisson &sp){

    sPoisson aux(sp);
    long double rad_sec=3600.0L*180.0L/PI_L;
    cout << "/------------------------- Elementos: " << aux.size() << endl;

    for(sPoisson::iterator i = aux.begin(); i != aux.end(); ++i)
        printf("%8.3Lf, %3d, %8.3Lf, [%4d, %4d]\n",i->getA()*rad_sec,i->getI(),i->getB(),i->getJ(0),i->getJ(1));

    cout << "\\-----------------------------------" << endl << endl;
};

void show_sec_char(const sPoisson &sp, int C){

    sPoisson aux(sp);
    long double rad_sec=3600.0L*180.0L/PI_L;
    cout << "/------------------------- Elementos: " << aux.size() << endl;

    for(sPoisson::iterator i = aux.begin(); i != aux.end(); ++i)
        if(fabs(i->getJ(0)+i->getJ(1))==C)
            printf("%8.3Lf, %3d, %8.3Lf, [%4d, %4d]\n",i->getA()*rad_sec,i->getI(),i->getB(),i->getJ(0),i->getJ(1));

    cout << "\\-----------------------------------" << endl << endl;
};

void show_sincos(const sPoisson &sp){

    sPoisson aux(sp);
    std::cout << "/------------------------- Elementos: " << aux.size() <<std::endl;

    for(sPoisson::iterator i = aux.begin(); i != aux.end(); ++i){

        printf("| %20.12Lf, %20.12Lf, %4i, [" ,i->getA()*cos(i->getB()),-i->getA()*sin(i->getB()), i->getI());
        for(int j=0;j<NP;j++)
        {
            if(j<NP-1)
                printf("%4i, ",i->getJ(j));
            else
                printf("%4i]\n",i->getJ(j));
        }
    };

    cout << "\\-----------------------------------" << endl << endl;
};

void show_sincos_sec(const sPoisson& sp){

    sPoisson aux(sp);
    long double rad_sec=3600.0L*180.0L/PI_L;
    std::cout << "/------------------------- Elementos: " << aux.size() <<std::endl;


    for(sPoisson::iterator i = aux.begin(); i != aux.end(); ++i){

        printf("| %20.12Lf, %20.12Lf, %4i, [" ,i->getA()*cos(i->getB())*rad_sec,-i->getA()*sin(i->getB())*rad_sec, i->getI());
        for(int j=0;j<NP;j++)
        {
            if(j<NP-1)
                printf("%4i, ",i->getJ(j));
            else
                printf("%4i]\n ",i->getJ(j));
        }
    };
    std::cout << "\\-----------------------------------" << std::endl<<std::endl;
};

long double valor(const sPoisson &sp, long double t){
    long double suma = 0.0L;

    sPoisson sp1(sp);
    ordenar_por_amplitud(sp1);

    for(sPoisson::reverse_iterator i = sp1.rbegin(); i != sp1.rend(); ++i)
        suma += i->getA() * pow(t, (long double)i->getI()) * cos(i->getB());

    return suma;
};

long double valor(const sPoisson& sp, long double t, long double M[NP]){
    long double suma = 0.0L, angulo;

    sPoisson sp1(sp);

    //Ordenar por amplitud para minimizar errores.
    ordenar_por_amplitud(sp1);

    for(sPoisson::reverse_iterator i = sp1.rbegin(); i != sp1.rend(); ++i){
        angulo = i->getB();
        for(int j=0; j<NP; j++)
            angulo += M[j]*(i->getJ(j));
        if(i->getI()==0)
            suma += i->getA() * cos(angulo);
        //else if(t!=0.0L)
		else if(fabs(t)>TOL)
            suma += i->getA() * pow(t, (long double)i->getI()) * cos(angulo);
    };
    return suma;
};

sPoisson& operator* (const sPoisson& sp1n, const sPoisson& sp2n){

    sPoisson aux;

    _multiplica(sp1n, sp2n, aux);

    sPoisson *prod = new sPoisson(aux);
    return *prod;
};

/*
sPoisson& sPoisson::operator= (const sPoisson &sp1n){

    sPoisson *copia = new sPoisson(sp1n);
    return *copia;
};
*/

void _suma(const sPoisson &sp1, const sPoisson &sp2, sPoisson &sSuma){

    sPoisson aux(sp1);
    aux.reserve(sp1.size() + sp2.size());
    aux.insert(aux.end(), sp2.begin(), sp2.end());

    if(aux.empty())
        aux.push_back( tPoisson(0.0L) );

    sSuma = compacta(aux);
};

void _multiplica(const sPoisson &sp1n, const sPoisson &sp2n, sPoisson &prod){

    sPoisson sp1, sp2;

    trunca(sp1n, sp1, CERO);
    trunca(sp2n, sp2, CERO);

    int tam_sp1 = sp1.size();
    int tam_sp2 = sp2.size();

    long int tam = 2*tam_sp1*tam_sp2;
    sPoisson aux1, aux2, aux3, aux4;

    //Vamos a trabajar con 4 procesadores.

    aux1.reserve( (int)(0.5+tam/4) );
    aux2.reserve( (int)(0.5+tam/4) );
    aux3.reserve( (int)(0.5+tam/4) );
    aux4.reserve( (int)(0.5+tam/4) );

    sPoisson sp11(sp1.begin(), sp1.begin()+(int)(tam_sp1/2));
    sPoisson sp12(sp1.begin()+(int)(tam_sp1/2), sp1.end());
    sPoisson sp21(sp2.begin(), sp2.begin()+(int)(tam_sp2/2));
    sPoisson sp22(sp2.begin()+(int)(tam_sp2/2), sp2.end());

	if(!sp11.empty()) ordenar_por_amplitud_decreciente(sp11);
	if(!sp12.empty()) ordenar_por_amplitud_decreciente(sp12);
	if(!sp21.empty()) ordenar_por_amplitud_decreciente(sp21);
	if(!sp22.empty()) ordenar_por_amplitud_decreciente(sp22);

#pragma omp parallel sections shared(sp1n,sp2n,aux1,aux2,aux3,aux4) num_threads(4)
{
    #pragma omp section
    {
        tPoisson tmp;
        for(auto const& i: sp21){
            for(auto const& j: sp11){
                if( fabs(0.5L * i.getA() * j.getA())>CERO ){
                    tmp.setA( 0.5L * i.getA() * j.getA() );
                    tmp.setI( i.getI() + j.getI() );

                    tmp.setB( i.getB() - j.getB() );
                    for(int k=0; k<NP; k++)
                        tmp.setJ(k, i.getJ(k) - j.getJ(k));
                    aux1.push_back(tmp);

                    tmp.setB( i.getB() + j.getB() );
                    for(int k=0; k<NP; k++)
                        tmp.setJ(k, i.getJ(k) + j.getJ(k));
                    aux1.push_back(tmp);
                }
                else break;
            }
        }
    }

    #pragma omp section
    {
        tPoisson tmp;
        for(auto const& i: sp21){
            for(auto const& j: sp12){
                if( fabs(0.5L * i.getA() * j.getA())>CERO ){
                    tmp.setA( 0.5L * i.getA() * j.getA() );
                    tmp.setI( i.getI() + j.getI() );

                    tmp.setB( i.getB() - j.getB() );
                    for(int k=0; k<NP; k++)
                        tmp.setJ(k, i.getJ(k) - j.getJ(k));
                    aux2.push_back(tmp);

                    tmp.setB( i.getB() + j.getB() );
                    for(int k=0; k<NP; k++)
                        tmp.setJ(k, i.getJ(k) + j.getJ(k));
                    aux2.push_back(tmp);
             	}
                else break;
            }
        }
    }

   #pragma omp section
    {
        tPoisson tmp;
        for(auto const& i: sp22){
            for(auto const& j: sp11){
                if( fabs(0.5L * i.getA() * j.getA())>CERO ){
                    tmp.setA( 0.5L * i.getA() * j.getA() );
                    tmp.setI( i.getI() + j.getI() );

                    tmp.setB( i.getB() - j.getB() );
                    for(int k=0; k<NP; k++)
                        tmp.setJ(k, i.getJ(k) - j.getJ(k));
                    aux3.push_back(tmp);

                    tmp.setB( i.getB() + j.getB() );
                    for(int k=0; k<NP; k++)
                        tmp.setJ(k, i.getJ(k) + j.getJ(k));
                    aux3.push_back(tmp);
                }
                else break;
            }
        }
    }

    #pragma omp section
    {
        tPoisson tmp;
        for(auto const& i: sp22){
            for(auto const& j: sp12){
                if( fabs(0.5L * i.getA() * j.getA())>CERO ){
                    tmp.setA( 0.5L * i.getA() * j.getA() );
                    tmp.setI( i.getI() + j.getI() );

                    tmp.setB( i.getB() - j.getB() );
                    for(int k=0; k<NP; k++)
                        tmp.setJ(k, i.getJ(k) - j.getJ(k));
                    aux4.push_back(tmp);

                    tmp.setB( i.getB() + j.getB() );
                    for(int k=0; k<NP; k++)
                        tmp.setJ(k, i.getJ(k) + j.getJ(k));
                    aux4.push_back(tmp);
                }
                else break;
            }
        }
    }
}

    prod.clear();
    long int reserva_mem = 0;

    if(!aux1.empty()){
        aux1=compacta(aux1);
        reserva_mem += aux1.size();
        prod.reserve( reserva_mem );
        prod.insert(prod.end(),
            make_move_iterator(aux1.begin()),
            make_move_iterator(aux1.end()));
        //prod.insert(prod.end(), aux1.begin(), aux1.end());
    };

    if(!aux2.empty()){
        aux2=compacta(aux2);
        reserva_mem += aux2.size();
        prod.reserve( reserva_mem );
        prod.insert(prod.end(),
            make_move_iterator(aux2.begin()),
            make_move_iterator(aux2.end()));
        //prod.insert(prod.end(), aux2.begin(), aux2.end());
    };

    if(!aux3.empty()){
        aux3=compacta(aux3);
        reserva_mem += aux3.size();
        prod.reserve( reserva_mem );
        prod.insert(prod.end(),
            make_move_iterator(aux3.begin()),
            make_move_iterator(aux3.end()));
        //prod.insert(prod.end(), aux3.begin(), aux3.end());
    };

    if(!aux4.empty()){
        aux4=compacta(aux4);
        reserva_mem += aux4.size();
        prod.reserve( reserva_mem );
        prod.insert(prod.end(),
            make_move_iterator(aux4.begin()),
            make_move_iterator(aux4.end()));
        //prod.insert(prod.end(), aux4.begin(), aux4.end());
    };

    if(prod.empty())
        prod.push_back( tPoisson(0.0L) );
    else
        prod=compacta(prod);
};

sPoisson& operator*(long double A, const sPoisson &sp){

    sPoisson aux(sp);
    sPoisson::iterator i;

    for(i = aux.begin(); i < aux.end(); ++i)
        i->setA( A*i->getA() );

    sPoisson *producto = new sPoisson( compacta(aux) );
    return *producto;
};

sPoisson& operator*(const sPoisson &sp, const long double A){
    return A*sp;
};

sPoisson& operator/(const sPoisson &sp, long double A){
    return (1.0L/A)*sp;
};

sPoisson& operator+(const sPoisson &sp1, const sPoisson &sp2){

    sPoisson aux;
    _suma(sp1,sp2,aux);

    sPoisson *suma = new sPoisson(aux);
    return *suma;
};

sPoisson& operator+ (const sPoisson& sp, long double numero){
    sPoisson *suma = new sPoisson(sp);

    suma->push_back( (tPoisson)numero );
    return compacta(*suma);
};

sPoisson& operator+ (long double numero, const sPoisson& sp){
    return sp+numero;
};

sPoisson& operator- (const sPoisson &sp1, const sPoisson &sp2){
 
    sPoisson aux(sp2);
    aux.reserve(sp1.size()+sp2.size());

    sPoisson::iterator i;
    for(i = aux.begin(); i < aux.end(); ++i)
        i->setA( -i->getA() );

    aux.insert(aux.end(),sp1.begin(),sp1.end());

    sPoisson *resta = new sPoisson( compacta(aux) );
    return *resta;
};

sPoisson& operator- (long double numero, const sPoisson &sp2){

    sPoisson aux(sp2);

    sPoisson::iterator i;
    for(i = aux.begin(); i < aux.end(); ++i)
        i->setA( -i->getA() );

    aux.push_back( tPoisson(numero) );

    sPoisson *resta = new sPoisson( compacta(aux) );
    return *resta;
};

sPoisson& operator- (const sPoisson& sp1, long double numero){
    return sp1+(-numero);
};

sPoisson& operator- (const sPoisson& sp){
    return (-1.0L)*sp;
};

sPoisson& compacta(sPoisson& spn){
    sPoisson aux;

    //arregla(spn);

    ordenar_por_angulo(spn);

    long double a, b, amp1, amp2, des1, des2;
    int terms = 0, contador = 0;

    //anterior inicializado a distinto del primero
    //sPoisson::const_iterator i = sp.begin();
    //tPoisson actual, anterior = *i;
    tPoisson actual, anterior = spn[0];
    anterior.setI(anterior.getI()+1);

    for(auto const& i: spn){
        actual = i;
        if(actual != anterior)
            aux.push_back(i);
        else{
            anterior = aux.back(); //ahora anterior es el ultimo que pusimos en la lista compactada
            amp1 = anterior.getA();
            amp2 = actual.getA();
            des1 = anterior.getB();
            des2 = actual.getB();
            a = amp1 * cos(des1) + amp2 * cos(des2);
            b = amp1 * sin(des1) + amp2 * sin(des2);
            anterior.setA(sqrt(a*a + b*b));
            anterior.setB(atan2(b,a));
            aux.pop_back();
            aux.push_back(anterior);
        };
        anterior = actual;
    };

    arregla(aux);
    return (spn = aux);
};

long double fact(const long double n){
    return (n<2.0L) ? 1.0L : n*fact(n-1);   //cuidado
};

long double binomial(const long double m, int n){
    long double valor = 1.0L;
    for(int i=0;i<n;i++)
        valor *= (m-(long double)i)/((long double)i+1.0L);
    return valor;
};

//__________________________________________________________________________
// OPERACIONES con series de Poisson
//
// La operatoria es similar en todas las funciones
// 1 - Encontrar el numero de terminos NT necesario para precision
// 2 - Aproximar por Taylor
//

sPoisson& sin(const sPoisson& sp){
    int NT = 0;
    long double errCometido = 1.0L, norma = norma_1(sp, 1.0L);

    do {
        errCometido *= norma * norma / ((3.0L+2.0L*(long double)NT)*(2.0L+ 2.0L*(long double)NT));
        NT++;
    }   while (fabs(errCometido) > TOL);

    sPoisson aux, sp_cuad(sp);
    _multiplica(sp_cuad, sp_cuad, sp_cuad);  //temp1=sp*sp  pero sin reservar memoria.

    for(int i=NT-1; i>=1; i--){
        aux.push_back( (tPoisson)(-1.0L/fact(2.0L*(long double)i+1.0L)) );
        for(auto& iter : aux)
            iter.setA( -iter.getA() );        //*sSeno = -(*sSeno);
        _multiplica(aux, sp_cuad, aux);          //aux = aux*sp*sp;
    };

    for(auto& iter : aux)
        iter.setA( -iter.getA() );   //*sSeno = -(*sSeno);

    aux.push_back( tPoisson(1.0L) );   //*sSeno = *sSeno+1
    _multiplica(aux, sp, aux);       //*sSeno = sSeno*sp  --> sin=x(1-x2/3!+x4/5!...)

    sPoisson *sSeno = new sPoisson(aux);
    return *sSeno;
};

sPoisson& cos(const sPoisson& sp){
    int NT = 1;
    long double norma = norma_1(sp, 1.0L), errCometido = norma;

    do {
        errCometido *= norma * norma / (2.0L*(long double)NT*(1.0L+ 2.0L*(long double)NT));
        NT++;
    }   while (fabs(errCometido) > TOL);

    sPoisson aux, sp_cuad(sp);
    _multiplica(sp_cuad, sp_cuad, sp_cuad);  //temp1=sp*sp  pero sin consumir memoria.

    for(int i=NT-1; i>=1; i--){
        aux.push_back( (-1.0L/fact(2.0L*(long double)i)) );
        for(auto& iter : aux)
            iter.setA( -iter.getA() );        //*sCoseno = -(*sCoseno);
        _multiplica(aux, sp_cuad, aux);          //aux = aux*sp*sp;
    };

    for(auto& iter : aux)
        iter.setA( -iter.getA() );

    aux.push_back( tPoisson(1.L) );       //*sCoseno = (-1.0L*(*sCoseno)+1.0L); //cos=1-x2/2!+x4/4!...

    sPoisson *sCoseno = new sPoisson( compacta(aux) );
    return *sCoseno;
};

sPoisson& raiz(const sPoisson& sp, long double indice /* = 2.L */){
    if(fabs(indice) < TOL) {
        cout << "Error en sqrt(sp, indice). Indice nulo. No continuar." << endl << endl;
        exit(1);
    }
    else
        indice = 1.0L/indice;

    int NT = 0;
    long double norma = norma_1(sp, 1.0L), errCometido = 1.L;

    do {
        errCometido *= norma*(indice-(long double)NT)/((long double)NT+1.0L);
        NT++;
    }   while (fabs(errCometido) > TOL);
                                                    //Revisar precision de esto.
    sPoisson aux;

    for(int i=NT; i>=1; i--){
        aux.push_back( binomial(indice, i) );   	//sRaiz = (sRaiz) + binomial(indice, i);
        _multiplica(aux, sp, aux);             		//sRaiz = (sRaiz)*sp;
    };
    aux.push_back( tPoisson(1.L) );        //sRaiz = sRaiz + 1.0L;

    sPoisson *sRaiz = new sPoisson( compacta(aux) );
    return *sRaiz;
};

sPoisson& exp(const sPoisson& sp){
    int NT = 0;
    long double errCometido = 1.0L, norma = norma_1(sp, 1.0L);

    do {
        errCometido *= norma/(1.0L+(long double)NT);
        NT++;
    }
    while (fabs(errCometido) > TOL);

	sPoisson::iterator iter;
	sPoisson::const_iterator iter1;

    sPoisson aux, sp_aux(sp);

	aux.push_back( tPoisson(1.L) );
	for(int i=NT; i>=1; i--){
		for(iter = sp_aux.begin(), iter1 = sp.begin(); iter1 != sp.end(); ++iter, ++iter1)
            iter->setA( iter1->getA()/(long double)i );

		_multiplica(aux, sp_aux, aux);
		aux.push_back( tPoisson(1.L) );

	};

//    for(int i=NT; i>=1; i--){
//        aux.push_back( 1.0L/fact((long double)i) );     //*sExp = *sExp + 1.0L/fact((long double)i);
//        _multiplica(aux, sp, aux);                       //*sExp = (*sExp) * sp;
//    };

//    aux.push_back( tPoisson(1L) );        //*sExp = *sExp + 1.0L;

    sPoisson *sExp = new sPoisson( compacta(aux) );
    return *sExp;
};

sPoisson& log(const sPoisson& sp){
    int NT = 1;
    long double errCometido = 1.0L, norma = norma_1(sp, 1.0L);

/*
    do {
        errCometido *= norma;
        if(NT>1000){
            cout << "Log no converge." << endl;
            break;
        }
		NT++;
    }
    while (fabs(errCometido)/NT > TOL);

	NT--;	//Revisar esto. Creo que sobre una iteracion
*/

	if(fabs(norma)>=1.0L){
		cout << "Log no converge" << endl;
		exit(1);
	};

    do {
        errCometido *= norma*(long double)NT/(1.0L+(long double)NT);
        NT++;
    }
    while (fabs(errCometido) > TOL);

    sPoisson aux;

    for(int i=NT; i>=1; i--){
        aux.push_back( 1.0L/((long double)i) );
        for(sPoisson::iterator iter = aux.begin(); iter != aux.end(); ++iter)
            iter->setA( -iter->getA() );        //sLog = -1.0L*(sLog);
        _multiplica(aux, sp, aux);

    };

    for(sPoisson::iterator iter = aux.begin(); iter != aux.end(); ++iter)
        iter->setA( -iter->getA() );

    sPoisson *sLog = new sPoisson(aux);
    return *sLog;
};

sPoisson& pow(const sPoisson& sp, int exp){

    sPoisson aux(sp);
    if(exp > 0)
        for(int i=1; i<exp; i++)
            _multiplica(aux,sp,aux);
    else{
        aux.clear();
        aux.insert(aux.begin(), tPoisson(0.L) );
    };

    sPoisson *sPow = new sPoisson( compacta(aux) );
    return *sPow;
};

sPoisson& trunca(const sPoisson& sp, long double err /* = TOL */){

    sPoisson aux(sp);

    sPoisson::iterator i = aux.begin();
    while(i != aux.end())
        if(fabs(i->getA()) < err)
            i = aux.erase(i);
        else
            i++;

    sPoisson *sp_truncada = new sPoisson( aux );
    return *sp_truncada;
};

void trunca(const sPoisson& sp, sPoisson& sp_truncada, long double err /* = TOL*/){

    sPoisson aux(sp);

    sPoisson::iterator i = aux.begin();
    while(i != aux.end())
        if(fabs(i->getA()) < err)
            i = aux.erase(i);
        else
            i++;

    sp_truncada = aux;
};

void ordenar_por_angulo(sPoisson &sp){

    list<tPoisson> lp(sp.begin(), sp.end());
    lp.sort(menorTermino);
    //sp.assign(lp.begin(), lp.end());
    sp.assign(make_move_iterator(lp.begin()),
              make_move_iterator(lp.end()));
};

void ordenar_por_amplitud(sPoisson &sp){

    list<tPoisson> lp(sp.begin(), sp.end());
    lp.sort(menorTerminoA);
    sp.assign(make_move_iterator(lp.begin()),
              make_move_iterator(lp.end()));
};

void ordenar_por_amplitud_decreciente(sPoisson &sp){

    list<tPoisson> lp(sp.begin(), sp.end());
    lp.sort(mayorTerminoA);
    sp.assign(make_move_iterator(lp.begin()),
              make_move_iterator(lp.end()));
};

sPoisson& chop(const sPoisson& sp1,long double err)
{
     tPoisson tp;
     sPoisson sp(sp1);
     sPoisson::iterator it;
     sPoisson *sp2=new sPoisson;
     it=sp.begin();

     while(it!=sp.end())                                        //OJO CON ESTO!!
     {
     if(fabsl(it->getA())>=err)
        (*sp2).push_back(*it);
        it++;
     }
     if((*sp2).size()==0)
     {
         tp=tPoisson(0.0L);
         (*sp2).push_back(tp);
     }
//   sp=sp2;//linea a帽adida para no tener que asignar cuando se llama
     return *sp2;// linea que se puede eleiminar haciendo la funci贸n vacia
};


/*
____________________________________________ Vectores de series
*/

void show(const vPoisson& vp, int nMostrar){

    cout << "####### Vector vPoisson de " << vp.size() << " elementos:" << endl;
    for(auto const& iter: vp)
        show(iter, nMostrar);
    cout << "#######" << endl << endl;
};

vPoisson& compacta(vPoisson& vp){

    for(auto& iter: vp)
        iter=compacta(iter);

    return vp;
};

vPoisson& operator-(const vPoisson& v1){ 
        long double i{};
	vPoisson::const_iterator it1;
        vPoisson:: iterator it2;
	vPoisson aux(v1.size());
        for(it1 = v1.begin(), it2 =aux.begin()+i;it1 != v1.end();it1++,i++)
	{
		*it2=-(*it1);
		arregla(*it2);
	}
	vPoisson *cambio_signo= new vPoisson(aux);
        return *cambio_signo;
};
	

vPoisson& operator+(const vPoisson& v1, const vPoisson& v2){

    vPoisson aux;
    aux.reserve(ORDEN_TAYLOR+1);
    _suma(v1, v2, aux);

    vPoisson *vSuma = new vPoisson(aux);
    return *vSuma;
};

vPoisson& operator-(const vPoisson& v1, const vPoisson& v2){

    vPoisson aux;
    vPoisson::const_iterator it1,it2;
    for(it1=v1.begin(),it2=v2.begin();it1 !=v1.end(),it2 !=v2.end();it1++,it2++)aux.push_back((*it1)-(*it2));	

    vPoisson *vSuma = new vPoisson(aux);
    return *vSuma;
};

void _suma(const vPoisson& v1n, const vPoisson& v2n, vPoisson& aux){

    sPoisson temp;
    temp.push_back( tPoisson(0.L) );

    aux.resize( ORDEN_TAYLOR+1, temp );

    vPoisson v1(v1n);
    v1.resize(ORDEN_TAYLOR+1, temp); //Resize borra los que sobran.

    vPoisson v2(v2n);
    v2.resize(ORDEN_TAYLOR+1, temp);

    vPoisson::const_iterator it1, it2;
    vPoisson::iterator it3;
    for(it1 = v1.begin(), it2 = v2.begin(), it3 = aux.begin();
        it1 != v1.end() && it2 != v2.end() && it3 != aux.end();
        ++it1, ++it2, ++it3 )
            _suma(*it1, *it2, *it3);

    aux=compacta(aux);
};

vPoisson& operator*(const vPoisson& v1, const vPoisson& v2){

    vPoisson aux;
    aux.reserve(ORDEN_TAYLOR+1);
    _multiplica(v1, v2, aux);

    vPoisson *vProd = new vPoisson(aux);
    return *vProd;
};

vPoisson& operator*(const sPoisson& v1, const vPoisson& v2){

    vPoisson aux;
    aux.reserve(ORDEN_TAYLOR+1);
    _multiplica(v1, v2, aux);

    vPoisson *vProd = new vPoisson(aux);
    return *vProd;
};

void _multiplica(long double x, const vPoisson& v1n, vPoisson& aux){
	long int i{};	
        vPoisson::const_iterator it1;
	vPoisson::iterator it2;
	for(it1=v1n.begin(),it2=aux.begin()+i; it1!=v1n.end(); ++i)*it2=*it1*x;
	aux=compacta(aux);
};

void _multiplica(const sPoisson& s1, const vPoisson& v1n, vPoisson& aux){
	long int i{};	
        vPoisson::const_iterator it1;
	vPoisson::iterator it2;
	for(it1=v1n.begin(),it2=aux.begin()+i; it1!=v1n.end(); ++i)*it2=*it1*s1;
	aux=compacta(aux);
};

void _multiplica(const vPoisson& v1n, const vPoisson& v2n, vPoisson& aux){

    vPoisson v1(v1n), v2(v2n);
    sPoisson temp(1, tPoisson(0.L) );

    //Inicializamos un vector de series nulas.
    aux.clear();
    aux.insert(aux.begin(), ORDEN_TAYLOR+1, temp);
    //vPoisson aux( ORDEN_TAYLOR+1, temp );

    vPoisson::const_iterator it1, it2;
    vPoisson::iterator it3;
    long int i{}, j{};

    for(it1=v1.begin(); it1!=v1.end(); ++it1, i++)
        for(it2=v2.begin(), it3=aux.begin()+i, j=0 ; it2!=v2.end(); ++it2, ++it3, j++)
                if(i+j<=ORDEN_TAYLOR){
                    _multiplica(*it1, *it2, temp);
                    _suma(*it3, temp, *it3);
                };
    aux=compacta(aux);
};

void _multiplica(const vPoisson& vp1, long double n, vPoisson& vp2){

    vp2.assign(vp1.begin(), vp1.end());

    for(auto& serie: vp2)
        for(auto& termino: serie)
            termino.setA(n*termino.getA());
};

vPoisson& sin(const vPoisson& vp){

    int NT = 0;
    long double errCometido = 1.0L, norma = norma_1(vp, 1.0L);

    do {
        errCometido *= norma * norma / ((3.0L+2.0L*(long double)NT)*(2.0L+ 2.0L*(long double)NT));
        NT++;
    }
    while (errCometido > TOL);

    sPoisson sCero(1, tPoisson(0.L));
    sPoisson sUno(1, tPoisson(1.L));
    vPoisson vUno( ORDEN_TAYLOR, sCero );
    vUno.insert(vUno.begin(), sUno);
    vPoisson vpn;
    _multiplica(vp, vp, vpn);
    _multiplica(vpn, -1.L, vpn);  //vpn=-x^2

    vPoisson suma(vUno), tmp1(vUno), tmp2;
    for(int orden=1; orden<=NT; orden++){
        _multiplica(tmp1, vpn, tmp1);
        _multiplica(tmp1, 1.L/fact(2*orden+1), tmp2);
        _suma(suma, tmp2, suma);
    };
    _multiplica(suma, vp, suma);

    vPoisson *vSin = new vPoisson(suma);
    return *vSin;
};

vPoisson& sin2(const vPoisson& vp){

    int NT = 0;
    long double errCometido = 1.0L, norma = norma_1(vp, 1.0L);

    do {
        errCometido *= norma * norma / ((1.0L+2.0L*(long double)NT)*(2.0L+ 2.0L*(long double)NT));
        NT++;
    }
    while (errCometido > TOL);

cout << "Hacen falta " << NT << " iteraciones para Taylor." << endl;

    sPoisson sCero(1, tPoisson(0.L));
    sPoisson sUno(1, tPoisson(1.L));
    vPoisson vUno( ORDEN_TAYLOR, sCero );
    vUno.insert(vUno.begin(), sUno);
    vPoisson vpn;
    _multiplica(vp, vp, vpn);

    vPoisson prod(vUno), tmp1(vUno), tmp2;
    for(int orden=1; orden<=           15000              ; orden++){
        _multiplica(vpn,-1.0L/(PI_L*PI_L*(long double)(orden*orden)), tmp2);
        _suma(vUno,tmp2,tmp2);
        _multiplica(tmp2, prod, prod);
    };
    _multiplica(prod, vp, prod);

    vPoisson *vSin = new vPoisson(prod);
    return *vSin;
};

vPoisson& cos(const vPoisson& vp){

    int NT = 1;
    long double norma = norma_1(vp, 1.0L), errCometido = norma;

    do {
        errCometido *= norma * norma / (2.0L*(long double)NT*(1.0L+ 2.0L*(long double)NT));
        NT++;
    }
    while (errCometido > TOL);

    sPoisson sCero(1, tPoisson(0.L));
    sPoisson sUno(1, tPoisson(1.L));
    vPoisson vUno( ORDEN_TAYLOR, sCero );
    vUno.insert(vUno.begin(), sUno);
    vPoisson vpn;
    _multiplica(vp, vp, vpn);
    _multiplica(vpn, -1.L, vpn);  //vpn=-x^2

    vPoisson suma(vUno), tmp1(vUno), tmp2;
    for(int orden=1; orden<=NT; orden++){
        _multiplica(tmp1, vpn, tmp1);
        _multiplica(tmp1, 1.L/fact(2*orden), tmp2);
        _suma(suma, tmp2, suma);
    };

    vPoisson *vSin = new vPoisson(suma);
    return *vSin;
};

vPoisson& exp(const vPoisson& vp){

    int NT = 0;
    long double errCometido = 1.0L, norma = norma_1(vp, 1.L);

    do {
        errCometido *= norma/(1.0L+(long double)NT);
        NT++;
    }
    while (errCometido > TOL);

    sPoisson sCero(1, tPoisson(0.L));
    sPoisson sUno(1, tPoisson(1.L));
    vPoisson vUno( ORDEN_TAYLOR, sCero );
    vUno.insert(vUno.begin(), sUno);

    vPoisson suma(vUno), tmp(vUno);

    for(int orden=1; orden<=NT; orden++){
        _multiplica(tmp, vp, tmp);
        _multiplica(tmp, 1.0L/(long double)orden, tmp);
        _suma(suma, tmp, suma);
    };

    vPoisson *vExp = new vPoisson(suma);
    return *vExp;
};

vPoisson& pow(const vPoisson& vp, int exp){

    if(exp > 0){
        vPoisson aux(vp);
        for(int i=1; i<exp; i++)
            _multiplica(aux,vp,aux);

        vPoisson *vPow = new vPoisson(aux);
        return *vPow;
    }
    else{
        sPoisson sCero {tPoisson(0.L)};
        sPoisson sUno {tPoisson(1.L)};
        vPoisson aux( ORDEN_TAYLOR, sCero );
        aux.insert(aux.begin(), sUno);

        vPoisson *vPow = new vPoisson(aux);
        return *vPow;
    }
};

vPoisson& log(const vPoisson& vp){

    int NT = 0;
    long double errCometido = 1.0L, norma = norma_1(vp, 1.L);

    do {
        errCometido *= norma;
        if(NT>1000){
            cout << "Log no converge." << endl;
            break;
        }
    }
    while (errCometido/NT++ > TOL);

    sPoisson sCero(1, tPoisson(0.L));
    sPoisson sUno(1, tPoisson(1.L));
    vPoisson vUno( ORDEN_TAYLOR, sCero );
    vUno.insert(vUno.begin(), sUno);

    vPoisson suma(vUno), tmp(vUno), aux(vUno);

    for(int orden=1, signo=-1; orden<=NT; orden++, signo=-signo){
        _multiplica(tmp, vp, tmp);
        _multiplica(tmp, (long double)signo/(long double)(orden+1), aux);
        _suma(suma, aux, suma);
    };

    _multiplica(suma, vp, suma);

    vPoisson *vLog = new vPoisson(suma);
    return *vLog;
};

vPoisson& raiz(const vPoisson& vp, long double indice /* = 2.L */){
    if(fabs(indice) < CERO) {
        cout << "Error en raiz(sp, indice). Indice nulo. No continuar." << endl;
        exit(1);
    }
    else
        indice = 1.0L/indice;

    int NT = 0;
    long double norma = norma_1(vp, 1.L), errCometido = 1.0L;

    do {
        errCometido *= norma*(indice-(long double)NT)/((long double)NT+1.0L);
        NT++;
    }   while (fabs(errCometido) > TOL);

    sPoisson sCero(1, tPoisson(0.L));
    sPoisson sUno(1, tPoisson(1.L));
    vPoisson vUno( ORDEN_TAYLOR, sCero );
    vUno.insert(vUno.begin(), sUno);

    vPoisson suma(vUno), tmp(vUno), aux(vUno);

    for(int i=1; i<=NT; i++){
        _multiplica(tmp, vp, tmp);
        _multiplica(tmp, binomial(indice, i), aux);
        _suma(suma, aux, suma);
    };

    vPoisson *vRaiz = new vPoisson(suma);
    return *vRaiz;
};