/* 
   ***************************************************************************
   <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/>. 
    **************************************************************************
*/

#ifndef HEADER_4ABD3F9DC4D3F4B6
#define HEADER_4ABD3F9DC4D3F4B6

#include <stack>
#include <vector>
#include <list>
#include <stdio.h>
#include <cstdlib>
#include <iostream>
#include <cmath>

#define sPoisson vector<tPoisson>
#define Pila stack< pair<operadores, sPoisson> >

using   namespace std;

const   int NP = 2;
const   double TOL = 1e-8;

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

enum    operadores {
        ninguno=0,      suma=1,         resta=2,        producto=4,
        division=8,     potencia=16,    raiz=32,        seno=64,
        coseno=128,     logaritmo=256,  exponencial=512
};

//__________________________________________________________________________
// CLASE Termino de Poisson

class tPoisson {
    private:
        double m_A;
        double m_B;
        int m_i;
        int m_j[NP];
    public:
        tPoisson (void);
        tPoisson (const tPoisson&);     //constructor de copia
        tPoisson (double, int, int[], double);
        tPoisson (double);
        ~tPoisson (void);
        double getA(void);
        double getB(void);
        int getI(void);
        void setA(const double&);
        void setB(const double&);
        void setI(const int&);
        int getJ(int);
        void setJ(int, int);
        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&);

//__________________________________________________________________________
// Constructores

tPoisson::tPoisson(void): m_A(0.0), m_i(0), m_B(0.0){
    for(int i=0;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++)
        m_j[i] = original.m_j[i];
};

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

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

tPoisson::~tPoisson (void){};

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

//__________________________________________________________________________
// Operadores de comparación

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());
};

//__________________________________________________________________________
// Clase seriePoisson es clase DERIVADA de sPoisson

class seriePoisson : public sPoisson{
    public:
        seriePoisson& trunca(double);
};

//__________________________________________________________________________
// Declaraciones

double norma_1(const sPoisson&, double);
double norma_2(const sPoisson&, double);
void arregla(sPoisson&);
void show(const sPoisson&);
void show_sincos(const sPoisson&);
double valor(const sPoisson&, double);
double valor(const sPoisson&, double, double[NP]);
sPoisson& operator*(const sPoisson&,const sPoisson&);
sPoisson& operator*(const sPoisson&, double);
sPoisson& operator*(double, const sPoisson&);
sPoisson& operator/(const sPoisson &, double);
sPoisson& operator+(const sPoisson&, const sPoisson&);
sPoisson& operator+(const sPoisson&, double);
sPoisson& operator+(double, const sPoisson&);
sPoisson& operator-(const sPoisson&, const sPoisson&);
sPoisson& operator-(const sPoisson&);
sPoisson& operator- (const double , const sPoisson&);
sPoisson& compacta(const sPoisson&);
sPoisson& sin(sPoisson&);
sPoisson& cos(sPoisson&);
sPoisson& log(sPoisson&);
sPoisson& exp(sPoisson&);
sPoisson& sqrt(sPoisson&, double);
sPoisson& pow(sPoisson&, int);

void operar1(Pila&);
sPoisson& operarTodo(Pila&);
sPoisson& operar(sPoisson&, operadores);
sPoisson& operar(sPoisson&, sPoisson&, operadores);
void apila(Pila&, sPoisson&, operadores);
void apila(Pila&, operadores);

//void quicksort( sPoisson &, int, int);
//int particion( sPoisson &, int, int);
//void swap(tPoisson * const, tPoisson * const);
void ordenar_por_angulo(sPoisson&);
void ordenar_por_amplitud(sPoisson&);
void ordenar_por_amplitud_decreciente(sPoisson&);
sPoisson& trunca(const sPoisson&, double);
sPoisson& chop(const sPoisson& ,double);
void show_sincos_sec(const sPoisson& sp);
void show_sec(const sPoisson& sp);
void show_sec_char(const sPoisson &sp, int C);

//__________________________________________________________________________
// Implementaciones

double norma_1(const sPoisson &sp1, double T){
    double norma = 0.0;
    sPoisson sp(sp1);
    sPoisson::iterator termino;
    arregla(sp);
    sp=compacta(sp);
    termino = sp.begin();
    while(termino != sp.end()){
        if(termino->getI()==0)  norma+=fabs( termino->getA());
        else norma += fabs( termino->getA()) * pow(fabs(T), termino->getI() );
        termino++;
    };
    return norma;
};
double norma_2(const sPoisson &sp1, double T){
    double norma = 0.0;
    sPoisson sp(sp1);
    sPoisson::iterator termino = sp.begin();
    arregla(sp);
    sp=compacta(sp);
    while(termino != sp.end()){
        if(termino->getI()==0)  norma+=pow(fabs( termino->getA()),2);
        else norma += fabs( termino->getA()) * pow(fabs(T), termino->getI() );
        termino++;
    };
    return sqrt(norma);
};
/*
void arregla(sPoisson& sp){
    //Quitar primeros terminos negativos y todos nulos
    bool todosCero;
    int pos;
    sPoisson::iterator termino = sp.begin();

    while(termino != sp.end()){
        if(termino->getJ(0)<0){
            for(int j=0;j<NP;j++)
                termino->setJ(j, -termino->getJ(j));
            termino->setB( fmod(-termino->getB(), 2.*M_PI) );
        };
        todosCero = true;
        pos = 0;
        while(pos<NP && todosCero)
            todosCero = (termino->getJ(pos++)==0);
        if(todosCero && termino->getB()!=0.0) {
            termino->setA(termino->getA()*cos(termino->getB()));
            termino->setB(0.0);
        };
        termino++;
    };
};
*/
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(), 2.*M_PI) );
              break;
           };
        }
        todosCero = true;
        pos = 0;
        while(pos<NP && todosCero) todosCero = (tp.getJ(pos++)==0);
        if(todosCero && tp.getB()!=0.0)
        {
            tp.setA(tp.getA()*cos(tp.getB()));
            tp.setB(0.0);
        };
        tp.setB(fmod(tp.getB(),6.2831853071796));
        while(tp.getB()<0.)tp.setB(tp.getB()+6.2831853071796);
        if(!todosCero&&tp.getA()<0)
        {
        	tp.setA(-tp.getA());

tp.setB(fmod(tp.getB()+M_PI,2.*M_PI));
        }

};

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


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

    while(termino != sp.end())
    {
        for(i=0;i<NP;i++)
        {
           todosCero = true;
           pos = 0;
           while(pos<i && todosCero) todosCero = (termino->getJ(pos++)==0);
           if(todosCero&&termino->getJ(i)<0)
           {
              for(int j=i;j<NP;j++)termino->setJ(j, -termino->getJ(j));
              termino->setB( fmod(-termino->getB(), 2.*M_PI) );
              break;
           };
        }
        todosCero = true;
        pos = 0;
        while(pos<NP && todosCero) todosCero = (termino->getJ(pos++)==0);
        if(todosCero && termino->getB()!=0.0)
        {
            termino->setA(termino->getA()*cos(termino->getB()));
            termino->setB(0.0);
        };
        termino->setB(fmod(termino->getB(),6.2831853071796));
        while(termino->getB()<0.)termino->setB(termino->getB()+6.2831853071796);
        if(!todosCero&&termino->getA()<0)
        {
        	termino->setA(-termino->getA());
                termino->setB(fmod(termino->getB()+M_PI,2.*M_PI));
        }
        termino++;

    };
};

void show(const sPoisson &sp){
    //sPoisson *aux = new sPoisson(sp);
    sPoisson aux(sp);
    cout << "/------------------------- Elementos: " << aux.size() << endl;

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

void show_sec(const sPoisson &sp){
    //sPoisson *aux = new sPoisson(sp);
    sPoisson aux(sp);
    double rad_sec=3600.*180./M_PI;
    cout << "/------------------------- Elementos: " << aux.size() << endl;

    sPoisson::iterator i = aux.begin();
    while(i != aux.end()) {
        printf("%8.3lf, %3d, %8.3lf, [%4d, %4d]\n",i->getA()*rad_sec,i->getI(),i->getB(),i->getJ(0),i->getJ(1));
        i++;
    };
    cout << "\\-----------------------------------" << endl << endl;
    //delete aux;
};

void show_sec_char(const sPoisson &sp, int C){
    //sPoisson *aux = new sPoisson(sp);
    sPoisson aux(sp);
    double rad_sec=3600.*180./M_PI;
    cout << "/------------------------- Elementos: " << aux.size() << endl;

    sPoisson::iterator i = aux.begin();
    while(i != aux.end()) {
        if(abs(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));
        i++;
    };
    cout << "\\-----------------------------------" << endl << endl;
    //delete aux;
};

void show_sincos(const sPoisson &sp){
    //sPoisson *aux = new sPoisson(sp);
    sPoisson aux(sp);
    std::cout << "/------------------------- Elementos: " << aux.size() <<std::endl;

    sPoisson::iterator i = aux.begin();
    while(i != aux.end()) {
  /*      std::cout << "| " << i->getA()*cos(i->getB()) << ", " <<
            -i->getA()*sin(i->getB()) << ", " << i->getI() << ", [";
        for(int j=0;j<NP;j++)
            std::cout << i->getJ(j) << ",";
        std::cout << "]" << std::endl; */

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

    };
    std::cout << "\\-----------------------------------" << std::endl<<std::endl;
    //delete aux;
};

void show_sincos_sec(const sPoisson& sp){
    //sPoisson *aux = new sPoisson(sp);
    sPoisson aux(sp);
    double rad_sec=3600.*180./M_PI;
    std::cout << "/------------------------- Elementos: " << aux.size() <<std::endl;

    sPoisson::iterator i = aux.begin();

    while(i != aux.end()) {
     /*   std::cout << "| " << i->getA()*cos(i->getB())*rad_sec << ", " <<
            -i->getA()*sin(i->getB())*rad_sec << ", " << i->getI() << ", [";
        for(int j=0;j<NP;j++)
            std::cout << i->getJ(j) << ",";
        std::cout << "]" << std::endl;*/
        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++)
            std::cout << i->getJ(j) << ",";
        std::cout << "]" << std::endl;*/
        for(int j=0;j<NP;j++)
        {
            if(j<NP-1)
                printf("%4i, ",i->getJ(j));
            else
                printf("%4i]\n ",i->getJ(j));
        }
        i++;
    };
    std::cout << "\\-----------------------------------" << std::endl<<std::endl;
    //delete aux;
};

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

    sPoisson sp1(sp);
    ordenar_por_amplitud(sp1);

    sPoisson::iterator i = sp1.end();
    while(i != sp1.begin()){
        i--;
        suma += i->getA() * pow(t, i->getI()) * cos(i->getB());

    };
    return suma;
};

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

    sPoisson sp1(sp);

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

    sPoisson::iterator i = sp1.end();
    do{
        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)
        suma += i->getA() * pow(t, i->getI()) * cos(angulo);

    }while(i != sp1.begin());
    return suma;
};

sPoisson& operator* (const sPoisson& sp1n,const sPoisson& sp2n){
    //Producto de 2 series de Poisson
    tPoisson tmp;
    sPoisson sp1(sp1n);
    sPoisson sp2(sp2n);
    //sPoisson *prod = new sPoisson;
    sPoisson *prod=new sPoisson;
    double error = 2*sp1.size()*sp2.size();
    sPoisson::iterator j, i;
 //   error = (error>0)? TOL/error : TOL;
    error = 1.e-14;
  //  printf("estoy aqui\n");
 //   ordenar_por_amplitud_decreciente(sp1);
//    ordenar_por_amplitud_decreciente(sp2);
    i= sp2.begin();
 //   printf("estoy aqui 2\n");
    while(i != sp2.end()){
        j = sp1.begin();
    //    if(fabs(i->getA() * j->getA()) < 1.e-14){ break;}
        while(j != sp1.end()){
    //            if(fabs(i->getA() * j->getA()) < 1.e-14)break;
                tmp.setA( .5 * 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));
                (*prod).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));
                (*prod).push_back(tmp);
            //};
            j++;
       };
       *prod=compacta(*prod);
       i++;
    };
//    printf("estoy aqui 4\n");
    return compacta(*prod);
};

sPoisson& operator*(double A, const sPoisson &sp){
    //sPoisson *producto = new sPoisson(sp);
    sPoisson producto(sp);
    sPoisson::iterator i = producto.begin();

    while(i != producto.end()) {
        i->setA( A*i->getA() );
        i++;
    };
    return compacta(producto);
};

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

sPoisson& operator/(const sPoisson &sp, double A){
    sPoisson spaux(sp);
    double B;
    B=1.0/A;
    return B*spaux;
};

sPoisson& operator+(const sPoisson &sp1, const sPoisson &sp2){
    //sPoisson *suma = new sPoisson(sp1);
    sPoisson suma(sp1);
    sPoisson::const_iterator i = sp2.begin();

    while(i != sp2.end()) {
        suma.push_back(*i);
        i++;
    };
    return compacta(suma);
};

sPoisson& operator+ (const sPoisson& sp, double numero){
    //sPoisson *aux = new sPoisson(sp);
    sPoisson aux(sp);
    aux.push_back((tPoisson)numero);
    return compacta(aux);
};

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

sPoisson& operator- (const sPoisson &sp1, const sPoisson &sp2){
    sPoisson aux;
    aux=sp1+(-1.0)*sp2;
    return compacta(aux);
};

sPoisson& operator- (const double numero, const sPoisson &sp2){
    sPoisson aux(sp2);
    aux.push_back(tPoisson(-numero));
    aux=-aux;
    return compacta(aux);
};

sPoisson& operator- (const sPoisson& sp){
    //return (-1.0)*sp; //??
    //sPoisson *spaux = new sPoisson(sp);
    sPoisson spaux(sp);
    sPoisson::iterator i = spaux.begin();
    while(i!=spaux.end()) {
       i->setA(-i->getA());
       i++;
    };
    return compacta(spaux);
};
/*
sPoisson& compacta(const sPoisson& spn){

    //sPoisson *sp = new sPoisson(spn);
    sPoisson sp(spn);
    sPoisson *aux = new sPoisson;
    int control=0,n1[NP],n2[NP],j;
    arregla(sp);

    if(sp->size() < 2) return spn; //OJO falta delete

    //sp.sort(menorTermino);
    //quicksort(sp, 0, sp.size()-1);
    ordenar_por_angulo(sp);

    //tPoisson anterior, actual;
    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;
    anterior.setI(anterior.getI()+1);

    while(i != sp.end()){
        actual = *i;
        control=0;
        for(j=0;j<NP;j++)
        {
            n1[j]=anterior.getJ(j);
            n2[j]=actual.getJ(j);
            if(n1[j]!=n2[j]) control=1;
        }
        if(control!=0)
            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;
        i++;
    };

*/
sPoisson& compacta(const sPoisson& spn){
    //sPoisson *sp = new sPoisson(spn);
    sPoisson sp(spn);
    sPoisson *aux = new sPoisson;
    arregla(sp);

    //if(sp->size() < 2) return spn; //OJO falta delete

    //sp.sort(menorTermino);
    //quicksort(sp, 0, sp.size()-1);
    ordenar_por_angulo(sp);

    //tPoisson anterior, actual;
    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;
    anterior.setI(anterior.getI()+1);

    while(i != sp.end()){
        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;
        i++;
    };

//    delete sp;
    return *aux;
};

double fact(double n){
    return (n<2) ? 1. : n*fact(n-1);
};

double binomial(double m, int n){
    double valor = 1.;
    for(int i=0;i<n;i++)
        valor *= (m-i)/(i+1);
    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(sPoisson& sp){
    int NT = 0;
    double errCometido = 1., norma = norma_1(sp, .1);

    do {
        errCometido *= norma * norma / ((1.+2.*NT)*(2.+ 2.*NT));
        NT++;
        if(NT>1.E+4) cout << "Seno: no convergencia!" << endl;
    }   while (errCometido > TOL);

    //sPoisson *sSeno = new sPoisson;
    sPoisson sSeno;

    for(int i=NT-1; i>=1; i--){
        sSeno = -1.0*(sSeno) + (1.0/fact(2.*i+1.));
        sSeno = (sSeno)*sp*sp;
    };
    sSeno = (-1.0*(sSeno) + 1.0) * sp; //sin=x(1-x2/3!+x4/5!...)

    return compacta(sSeno);
};

sPoisson& cos(sPoisson& sp){
    int NT = 1;
    double norma = norma_1(sp, .1), errCometido = norma;

    do {
        errCometido *= norma * norma / (2.0*NT*(1.+ 2.*NT));
        NT++;
    }   while (errCometido > TOL);

    //sPoisson *sCoseno = new sPoisson;
    sPoisson sCoseno;

    for(int i=NT-1; i>=1; i--){
        sCoseno = -1.0*(sCoseno) + (1./fact(2.*i));
        sCoseno = (sCoseno)*sp*sp;
    };
    sCoseno = (-1.0*(sCoseno)+1.); //cos=1-x2/2!+x4/4!...

    return compacta(sCoseno);
};

sPoisson& sqrt(sPoisson& sp, double indice=2.0){
    if(indice==0.0) {
        cout << "Error en sqrt(sp, indice). Indice nulo. No continuar." << endl;
        exit(1);
        //sPoisson *unidad = new sPoisson;
        //sPoisson unidad;
        //unidad.push_back((tPoisson)1.0);
        //return( unidad );
    }
    else
        indice = 1./indice;

    int NT = 0;
    double norma = norma_1(sp, 0.1), errCometido = 1.0;

    do {
        errCometido *= norma*(indice-NT)/(NT+1);
        NT++;
    }   while (fabs(errCometido) > TOL);

    //sPoisson *sRaiz = new sPoisson;
    sPoisson sRaiz;

    for(int i=NT; i>=1; i--){
        sRaiz = (sRaiz) + binomial(indice, i);
        sRaiz = (sRaiz)*sp;
    };
    sRaiz = sRaiz + 1.0;

    return compacta(sRaiz);
};

sPoisson& exp(sPoisson& sp){
    int NT = 0;
    double errCometido = 1., norma = norma_1(sp, .1);

    do errCometido *= norma/(1.+NT++);
    while (errCometido > TOL);

    //sPoisson *sExp = new sPoisson;
    sPoisson sExp;

    for(int i=NT; i>=1; i--){
        sExp = sExp + 1./fact(i);
        sExp = (sExp) * sp;
    };
    sExp = sExp + 1.0;

    return compacta(sExp);
};

sPoisson& log(sPoisson& sp){
    int NT = 0;
    double errCometido = 1., norma = norma_1(sp, .1);

    do errCometido *= norma;
    while (errCometido/NT++ > TOL);

    //sPoisson *sLog = new sPoisson;
    sPoisson sLog;

    for(int i=NT; i>=1; i--){
        sLog = -1.0*(sLog) + 1.0/i;
        sLog = (sLog)*sp;
    };

    return compacta(sLog);
};

sPoisson& pow(sPoisson& sp, int exp){
    //sPoisson *sTmp = new sPoisson(sp);
    sPoisson sTmp(sp);
    for(int i=1;i< exp; i++)
        sTmp = (sTmp)*sp;
    return compacta(sTmp);
};

void operar1(Pila& pila){
    pair<operadores, sPoisson > tupla1 = pila.top();
    operadores oper1 = tupla1.first;
    sPoisson sTmp1 = tupla1.second;
    pila.pop();

    pair<operadores, sPoisson > tupla2 = pila.top();
    operadores oper2 = tupla2.first;
    sPoisson sTmp2 = tupla2.second;
    pila.pop();

    if(oper1 == ninguno && (oper2 & 0x1e0) ){
        //tupla1.first = ninguno; //ya vale 0
        tupla1.second = operar(sTmp1, oper2);
        pila.push(tupla1); //memoria new??
    }
    else{
        pair<operadores, sPoisson > tupla3 = pila.top();
        operadores oper3 = tupla3.first;
        sPoisson sTmp3 = tupla3.second;
        pila.pop();

        if(oper1 == ninguno && oper2 == ninguno && (oper3 & 0x1f)){
            //tupla1.first = ninguno; //ya vale 0
            tupla1.second = operar(sTmp1, sTmp2, oper3);
            pila.push(tupla1); //memoria new??
        }
        else
            std::cout << "Pila no esta en Notacion Polaca Inversa" << std::endl;
    };
};

sPoisson& operarTodo(Pila& pila){
    while( pila.size() >= 2 ) //Mientras haya al menos 2
        operar1(pila);

    pair<operadores,sPoisson> *tupla=new pair<operadores, sPoisson>(pila.top());
    pila.pop();

    return tupla->second;
};

sPoisson& operar(sPoisson &sp, operadores oper){
    switch (oper){
        case raiz:          return sqrt(sp);
        case seno:          return sin(sp);
        case coseno:        return cos(sp);
        case logaritmo:     return log(sp);
        case exponencial:   return exp(sp);
    };
};

sPoisson& operar(sPoisson& sp1, sPoisson& sp2, operadores oper){
    switch (oper){
        case suma:      return sp1+sp2;
        case resta:     return sp1-sp2;
        case producto:  return sp1*sp2;
        case division:  return sp1/valor(sp2);
        case potencia:  return pow(sp1, (int)valor(sp2));
    };
};

void apila(Pila &pila, sPoisson& sp, operadores oper = ninguno){
    pila.push(make_pair(oper, sp));
};

void apila(Pila &pila, operadores oper){
    sPoisson *sp = new sPoisson;
    pila.push(make_pair(oper, *sp));
};

/*
void quicksort( sPoisson & vector, int first, int last ){
  //int particion( sPoisson &, int, int );
  int posicionActual;
  if ( first >= last )
    return;
  posicionActual = particion( vector, first, last );
  quicksort( vector, first, posicionActual - 1 );
  quicksort( vector, posicionActual + 1, last );
};

int particion( sPoisson& vector, int izda, int dcha ){
  int posicion = izda;
  while ( true ) {
    while ( vector[ posicion ] <= vector[ dcha ] && posicion != dcha )
      --dcha;
    if ( posicion == dcha )
      return posicion;
    if ( vector[ posicion ] > vector[ dcha ]) {
      swap( &vector[ posicion ], &vector[ dcha ] );
      posicion = dcha;
    }
    while ( vector[ izda ] <= vector[ posicion ] && izda != posicion )
      ++izda;
    if ( posicion == izda )
      return posicion;
    if ( vector[ izda ] > vector[ posicion ] ) {
      swap( &vector[ posicion ], &vector[ izda ] );
      posicion = izda;
    }
  }
};

void swap(tPoisson * const ptr1, tPoisson * const ptr2 ){
    tPoisson temp = *ptr1;
    *ptr1 = *ptr2;
    *ptr2 = temp;
};
*/
seriePoisson& seriePoisson::trunca(double err = TOL){

    seriePoisson *aux = new seriePoisson;

    err /= aux->size() ? aux->size() : 1.;

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

    return *aux;
};

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

    sPoisson *aux = new sPoisson(sp);

    err /= aux->size() ? aux->size() : 1.;

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

    return *aux;
};

void ordenar_por_angulo(sPoisson &sp){
    //Pasa de un vector a una lista para ordenar por amplitud
    list<tPoisson> lp;

    sPoisson::iterator i = sp.begin();
    while(i != sp.end()){
        lp.push_back(*i);
        i++;
    };

    lp.sort(menorTermino);

    sp.clear();
    list<tPoisson>::iterator j = lp.begin();
    while(j != lp.end()){
        sp.push_back(*j);
        j++;
    };
};

void ordenar_por_amplitud(sPoisson &sp){
    //Pasa de un vector a una lista para ordenar por amplitud
    list<tPoisson> lp;
/*
    lp <- sp
    ordenar lp
    sp <- lp
*/
    sPoisson::iterator i = sp.begin();
    while(i != sp.end()){
        lp.push_back(*i);
        i++;
    };

    lp.sort(menorTerminoA);

    sp.clear();
    list<tPoisson>::iterator j = lp.begin();
    while(j != lp.end()){
        sp.push_back(*j);
        j++;
    };
};

sPoisson& chop(const sPoisson& sp1,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(fabs(it->getA())>=err)(*sp2).push_back(*it);
     it++;
     }
     if((*sp2).size()==0)
     {
         tp=tPoisson(0.);
         (*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
};
void ordenar_por_amplitud_decreciente(sPoisson &sp){
    //Pasa de un vector a una lista para ordenar por amplitud
    list<tPoisson> lp;
/*
    lp <- sp
    ordenar lp
    sp <- lp
*/
    sPoisson::iterator i = sp.begin();
    while(i != sp.end()){
        lp.push_back(*i);
        i++;
    };

    lp.sort(menorTerminoA);

    sp.clear();
    list<tPoisson>::iterator j = lp.end();
    j--;
    sp.push_back(*j);
    while(j != lp.begin()){
            j--;
        sp.push_back(*j);


    };
};

#endif // header guard