#define MAX(X,Y) (X>Y)?X:Y
#define MIN(X,Y) (X<Y)?X:Y

static long int lcm(long int a, long int b);
static int getmax(token_num* token_info, int n);

static poly polycut(poly poly_old,long int time_start,long int time_end)
{
    poly result;
    int poly_num,index=0;
    poly_num = poly_old.getn();
    while(poly_old.getdata()[index].getd()<=time_start) index++;
    while(poly_old.getdata()[index].getd()<=time_end && index<poly_num)
        result = oplus(result,poly_old.getdata()[index++]);
    return result;
}

static poly serie2poly(serie serie_old, long int time_start, long int time_end)
{
    poly temp,temp_add,result,serie_p,serie_q;
    gd serie_r;
    int index,gd_p_num,gd_q_num,gd_t_num;
    if(time_start>=time_end)
        return result;
    serie_p = serie_old.getp();
    serie_q = serie_old.getq();
    serie_r = serie_old.getr();
    gd_p_num = serie_p.getn();
    gd_q_num = serie_q.getn();
    temp = oplus(serie_p,serie_q);
    temp_add = serie_q;
    index = 0;
    while(true)
    {
        gd_t_num = temp.getn();
        for(;index<gd_t_num;index++)
            if(temp.getdata()[index].getd()>time_start) break;
        if(index == gd_t_num)
        {
            if(serie_r.getd()!=0)
            {
                temp_add = otimes(temp_add,serie_r);
                std::cout<<"fff:"<<temp_add<<std::endl;
                temp = oplus(temp,temp_add);
            }
            else return result;
        }
        else break;
    }
    while(true)
    {
        gd_t_num = temp.getn();
        for(;index<gd_t_num;index++)
            if(temp.getdata()[index].getd()>time_end) break;
            else
                result = oplus(result,temp.getdata()[index]);
        if(index == gd_t_num)
        {
            if(serie_r.getd()!=0)
            {
                temp_add = otimes(temp_add,serie_r);
                temp = oplus(temp,temp_add);
            }
            else break;
        }
        else break;
    }
    temp.init(result.getdata()[result.getn()-1].getg(),infinity);
    result = oplus(result,temp);
    return result;
}

serie fracodotflat(serie serie_a,serie serie_b)
{
    float niu_a,niu_b;
    poly a_p,a_q,b_p,b_q,result_p,result_q,poly1,poly2;
    serie result,temp_serie;
    gd a_r,b_r,result_r;
    long int period_lcm,period_start_time;
    a_p = serie_a.getp();
    b_p = serie_b.getp();
    a_q = serie_a.getq();
    b_q = serie_b.getq();
    a_r = serie_a.getr();
    b_r = serie_b.getr();
    if(a_r.getd() != 0)
        niu_a = a_r.getg()/a_r.getd();
    else niu_a = 0;
    if(b_r.getd() != 0)
        niu_b = b_r.getg()/b_r.getd();
    else niu_b = 0;
    if(niu_a == 0 && niu_b == 0)            //if serie a and b are all not periodic
    {
        poly1 = oplus(a_p,a_q);
        poly2 = oplus(b_p,b_q);
        poly1 = fracodotflat(poly1,poly2);
        return poly1;
    }
    if(niu_a < niu_b)
        return result;
    else
    {
        period_lcm = lcm(a_r.getd(),b_r.getd());
        period_start_time = (a_q.getdata()[0].getd()>b_q.getdata()[0].getd())?a_q.getdata()[0].getd():b_q.getdata()[0].getd();
        poly1 = serie2poly(serie_a,-1,period_start_time+period_lcm);
        std::cout<<"poly1: "<<std::endl<<poly1<<std::endl;
        poly1 = serie2poly(serie_b,-1,period_start_time+period_lcm);
        std::cout<<"poly1: "<<std::endl<<poly1<<std::endl;
        poly1 = fracodotflat(serie2poly(serie_a,-1,period_start_time+period_lcm),serie2poly(serie_b,-1,period_start_time+period_lcm));
        std::cout<<"result1: "<<std::endl<<poly1<<std::endl;
        temp_serie.init(epsilon,poly1,e);
        result_p = serie2poly(temp_serie,-1,period_start_time);
        poly1 = fracodotflat(serie2poly(serie_a,period_start_time,period_start_time+2*period_lcm),serie2poly(serie_b,period_start_time,period_start_time+2*period_lcm));
        std::cout<<"poly1: "<<std::endl<<poly1<<std::endl;
        temp_serie.init(epsilon,poly1,e);
        result_q = serie2poly(temp_serie,period_start_time,period_start_time+period_lcm);
        result_r.init(a_r.getg()*period_lcm/a_r.getd()-b_r.getg()*period_lcm/b_r.getd(),period_lcm);
        result.init(result_p,result_q,result_r);
        return result;
    }
}

serie fracodotsharp(serie serie_a,serie serie_b)
{
    token_num* token_info;
    int token_length;
    float niu_a,niu_b;
    poly a_p,a_q,b_p,b_q,result_p,result_q,poly1,poly2,temp_result;
    serie result,temp_serie;
    gd a_r,b_r,result_r;
    long int period_lcm,period_start_time;
    int index,max_number,temp;
    poly1 = serie2poly(serie_a,-1,100);
    poly2 = serie2poly(serie_b,-1,100);
    token_info = ominus(poly1,poly2,&token_length);
    print(token_info,token_length);
    a_p = serie_a.getp();
    b_p = serie_b.getp();
    a_q = serie_a.getq();
    b_q = serie_b.getq();
    a_r = serie_a.getr();
    b_r = serie_b.getr();
    if(a_r.getd() != 0)
        niu_a = a_r.getg()/a_r.getd();
    else niu_a = 0;
    if(b_r.getd() != 0)
        niu_b = b_r.getg()/b_r.getd();
    else niu_b = 0;
    if(niu_a == 0 && niu_b == 0)            //if serie a and b are all not periodic
    {
        poly1 = oplus(a_p,a_q);
        poly2 = oplus(b_p,b_q);
        poly1 = fracodotsharp(poly1,poly2);
        return poly1;
    }
    else if(niu_a < niu_b)
    {
        period_start_time = (a_q.getdata()[0].getd()>b_q.getdata()[0].getd())?a_q.getdata()[0].getd():b_q.getdata()[0].getd();
        period_lcm = lcm(a_r.getd(),b_r.getd());
        index = 1;
        while(true)
        {
            poly1 = serie2poly(serie_a,-1,period_start_time+index*period_lcm);
            poly2 = serie2poly(serie_b,-1,period_start_time+index*period_lcm);
            temp_result = fracodotsharp(poly1,poly2);
            if(result_q == temp_result)
            {
                result.init(epsilon,result_q,e);
                return result;
            }
            else
            {
                result_q = temp_result;
                index++;
            }
        }
    }
    else
    {
        period_start_time = (a_q.getdata()[0].getd()>b_q.getdata()[0].getd())?a_q.getdata()[0].getd():b_q.getdata()[0].getd();
        period_lcm = lcm(a_r.getd(),b_r.getd());
        std::cout<<period_start_time<<"::"<<period_lcm<<std::endl;
        poly1 = serie2poly(serie_a,-1,period_start_time + period_lcm);
        poly2 = serie2poly(serie_b,-1,period_start_time + period_lcm);
        token_info = ominus(poly1,poly2,&token_length);
        max_number = getmax(token_info,token_length);
        index = 1;
        while(true)
        {
            poly1 = serie2poly(serie_a,index*period_lcm,(index+1)*period_lcm);
            poly2 = serie2poly(serie_b,index*period_lcm,(index+1)*period_lcm);
            token_info = ominus(poly1,poly2,&token_length);
            index++;
            if(max_number<getmax(token_info,token_length)) break;
        }
        poly1 = serie2poly(serie_a,-1,period_start_time+index*period_lcm);
        poly2 = serie2poly(serie_b,-1,period_start_time+index*period_lcm);
        poly1 = fracodotsharp(poly1,poly2);
        std::cout<<"poly1: "<<std::endl<<poly1<<std::endl;
        result_q = polycut(poly1,(index-1)*period_lcm,index*period_lcm);
        std::cout<<"result_q: "<<std::endl<<result_q<<std::endl;
        result_p = polycut(poly1,-1,period_start_time+(index-1)*period_lcm);
        std::cout<<"result_p: "<<std::endl<<result_p<<std::endl;
        result_r.init(a_r.getg()*period_lcm/a_r.getd()-b_r.getg()*period_lcm/b_r.getd(),period_lcm);
        result.init(result_p,result_q,result_r);
        return result;
    }
}

static long int lcm(long int a, long int b)
{
    long int c,copy_a,copy_b;
    copy_a = a;
    copy_b = b;
    while(true)
    {
        c = copy_a % copy_b;
        if(c==0)
            return (a*b/copy_b);
        else
        {
            copy_a = copy_b;
            copy_b = c;
        }
    }
}

static int getmax(token_num* token_info, int n)
{
    int result = 0,index;
    for(index=0;index<n;index++)
        if((token_info+index)->token_num>result)
            result = (token_info+index)->token_num;
    return result;
}

poly odotsharp(poly p1, poly p2)
{
    long int p_max;
    int p1_index = 0,p2_index = 0;
    long int t;
    poly result;
    gd temp_gd;
    int p1_n, p2_n;
    p1_n = p1.getn();
    p2_n = p2.getn();
    p_max = p1.getpol(p1_index).getg()-p2.getpol(p2_index).getg();
    t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
    do
    {
        if (p1_index == p1_n-1 && p2_index == p2_n-1)
        {
            t = MAX(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
            temp_gd.init(p_max,t);
            result = oplus(result,temp_gd);
            return result;
        }
        if ((p1_index == p1_n-1 || p1.getpol(p1_index).getd() > p2.getpol(p2_index).getd()) && p2_index != p2_n-1)
            p2_index++;
        else if ((p2_index == p2_n-1 || p1.getpol(p1_index).getd() < p2.getpol(p2_index).getd()) && p1_index != p1_n-1)
            p1_index++;
        else
        {
            p1_index++;
            p2_index++;
        }
        std::cout<<"p1:"<<std::endl<<p1.getpol(p1_index)<<std::endl;
        std::cout<<"p2:"<<std::endl<<p2.getpol(p2_index)<<std::endl;
        std::cout<<"t:"<<std::endl<<t<<std::endl;
        std::cout<<"p_max:"<<std::endl<<p_max<<std::endl;
        if (p_max < p1.getpol(p1_index).getg() - p2.getpol(p2_index).getg())
        {
            temp_gd.init(p_max,t);
            result = oplus(result,temp_gd);
            p_max = p1.getpol(p1_index).getg() - p2.getpol(p2_index).getg();
        }
        t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
        std::cout<<"-----------------------------------------------------------------"<<std::endl;
    }while(true);
}

poly odotflat(poly p1, poly p2)
{
    long int p_min;
    int p1_index,p2_index;
    long int t;
    gd temp_gd;
    poly result;
    int p1_n, p2_n;
    p1_n = p1.getn();
    p2_n = p2.getn();
    p1_index = p1_n-1;
    p2_index = p2_n-1;
    p_min = infinity;
    do
    {
        std::cout<<"p1.getpol(p1_index):"<<std::endl<<p1.getpol(p1_index)<<std::endl;
        std::cout<<"p2.getpol(p2_index):"<<std::endl<<p2.getpol(p2_index)<<std::endl;
        if(p_min > p1.getpol(p1_index).getg() - p2.getpol(p2_index).getg())
        {
            t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
            p_min = p1.getpol(p1_index).getg() - p2.getpol(p2_index).getg();
            temp_gd.init(p_min,t);
            result = oplus(result,temp_gd);
            std::cout<<"miao"<<std::endl;
        }
        if (p1_index == 0 && p2_index == 0)
            return result;
        if ((p1_index == 0 || p1.getpol(p1_index-1).getd() < p2.getpol(p2_index-1).getd()) && p2_index != 0)
            p2_index--;
        else if ((p2_index == 0 || p1.getpol(p1_index-1).getd() > p2.getpol(p2_index-1).getd()) && p1_index != 0)
            p1_index--;
        else
        {
            p1_index--;
            p2_index--;
        }
        std::cout<<"-----------------------------------------------------------------"<<std::endl;
    }while(true);
}

serie odotsharp(serie s1, serie s2)
{
    long int p_max;
    long int A,Pk;
    long int start_time_minus;
    long int period;                    //the period
    long int r_0;                       //the steps
    int p1_index = 0,p2_index = 0;
    long int t,t_c,i,t_p1max,t_p2max;                       //i: index of period
    serie result;
    gd temp_gd,result_r;
    poly p1,p2,result_p,result_q;
    int p1_n, p2_n;
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
    start_time_minus = MAX(s1.getq().getpol(0).getd(),s2.getq().getpol(0).getd());                  //max(T'(0), T''(0))
    period = lcm(s1.getr().getd(),s2.getr().getd());                                                //lcm(tau',tau'')
    std::cout<<"period:"<<period<<std::endl;
    r_0 = s1.getr().getg()*(period/s1.getr().getd()) - s2.getr().getg()*(period/s2.getr().getd());      //steps
    std::cout<<"r_0:"<<r_0<<std::endl;
    result_r.init(r_0,period);
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
    p1.init(0,0);
    p2.init(0,0);
    do
    {
        t_p1max = p1.getpol(p1.getn()-1).getd();
        p1 = otimes(p1,s1.getr());
        p1 = oplus(p1,s1.getq());
    }while(p1.getpol(p1.getn()-1).getd() <= start_time_minus + period);
    p1 = oplus(p1,s1.getp());
    do
    {
        t_p2max = p2.getpol(p2.getn()-1).getd();
        p2 = otimes(p2,s2.getr());
        p2 = oplus(p2,s2.getq());
    }while (p2.getpol(p2.getn()-1).getd() <= start_time_minus + period);
    p2 = oplus(p2,s2.getp());
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
    std::cout<<"p1:"<<std::endl<<p1<<std::endl;
    std::cout<<"p2:"<<std::endl<<p2<<std::endl;
    p1_n = p1.getn();
    p2_n = p2.getn();
    A = MAX(p1.getpol(p1_index).getg(),p2.getpol(p2_index).getg());
    Pk = 0;
    t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
    do
    {
        if ((p1_index == p1_n-1 || p1.getpol(p1_index).getd() > p2.getpol(p2_index).getd()) && p2_index != p2_n-1)
            p2_index++;
        else if ((p2_index == p2_n-1 || p1.getpol(p1_index).getd() < p2.getpol(p2_index).getd()) && p1_index != p1_n-1)
            p1_index++;
        else
        {
            p1_index++;
            p2_index++;
        }
        t_c = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
        if (t_c < start_time_minus || (t < start_time_minus - 1 && t_c >= start_time_minus) )                                           //if t = 0,...,t_0-1
        {
            A = MAX(A,(p1.getpol(p1_index).getg()-p2.getpol(p2_index).getg()));
        }
        else if(t_c < start_time_minus + period || (t < start_time_minus + period - 1 && t_c >= start_time_minus + period))
        {
            Pk = MAX(Pk,(p1.getpol(p1_index).getg()-p2.getpol(p2_index).getg()));
        }
        t = t_c;
        if (t >= start_time_minus + period)
            break;
//        std::cout<<"t:"<<std::endl<<t<<std::endl;
//        std::cout<<"-----------------------------------------------------------------"<<std::endl;
    }while(true);
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
    if (A > Pk)
        i = (A-Pk)/r_0+1;
    else
        i = 0;
    std::cout<<"A,Pk:"<<std::endl<<A<<","<<Pk<<std::endl;
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    p1.init(0,0);
    p2.init(0,0);
    do
    {
        t_p1max = p1.getpol(p1.getn()-1).getd();
        p1 = otimes(p1,s1.getr());
        p1 = oplus(p1,s1.getq());
    }while(p1.getpol(p1.getn()-1).getd() <= start_time_minus + (i+2)*period);
    p1 = oplus(p1,s1.getp());
    do
    {
        t_p2max = p2.getpol(p2.getn()-1).getd();
        p2 = otimes(p2,s2.getr());
        p2 = oplus(p2,s2.getq());
    }while (p2.getpol(p2.getn()-1).getd() <= start_time_minus + (i+2)*period);
    p2 = oplus(p2,s2.getp());
    std::cout<<"p1:"<<std::endl<<p1<<std::endl;
    std::cout<<"p2:"<<std::endl<<p2<<std::endl;
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
    p1_n = p1.getn();
    p2_n = p2.getn();
    p1_index = 0;p2_index = 0;
    p_max = p1.getpol(p1_index).getg()-p2.getpol(p2_index).getg();
    t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
    do
    {
        if (t > start_time_minus + (i+2)*period)
        {
            result.init(result_p,result_q,result_r);
            return result;
        }
        if ((p1_index == p1_n-1 || p1.getpol(p1_index).getd() > p2.getpol(p2_index).getd()) && p2_index != p2_n-1)
            p2_index++;
        else if ((p2_index == p2_n-1 || p1.getpol(p1_index).getd() < p2.getpol(p2_index).getd()) && p1_index != p1_n-1)
            p1_index++;
        else
        {
            p1_index++;
            p2_index++;
        }
        if (p_max < p1.getpol(p1_index).getg() - p2.getpol(p2_index).getg())
        {
            if (t < start_time_minus + (i+1)*period)
            {
                temp_gd.init(p_max,t);
                result_p = oplus(result_p,temp_gd);
                p_max = p1.getpol(p1_index).getg() - p2.getpol(p2_index).getg();
            }
            else if (t < (MAX(t_p1max,t_p2max)))
            {
                temp_gd.init(p_max,t);
                result_q = oplus(result_q,temp_gd);
                p_max = p1.getpol(p1_index).getg() - p2.getpol(p2_index).getg();
            }
        }

        t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
        std::cout<<"-----------------------------------------------------------------"<<std::endl;
    }while(true);
    return result;
}

serie odotflat(serie s1, serie s2)
{
    long int p_min;         //the minimum in the second period
    long int start_time;    //the starting time of period
    long int period;        //the period
    long int r_0;           //the steps
    int p1_index,p2_index;  //The pointers to the current monomials in both polynomials
    long int t,t_p1max,t_p2max;             //temp
    int i;                  //counter
    gd temp_gd,result_r;
    poly p1,p2,result_p,result_q;
    serie result;
    int p1_n, p2_n;
/////////////////calculating the starting time of period, period and steps///////////////////////////////
    start_time = MAX(s1.getq().getpol(0).getd(),s2.getq().getpol(0).getd());    //max(T'(0), T''(0))
    period = lcm(s1.getr().getd(),s2.getr().getd());                            //lcm(tau',tau'')
    r_0 = s1.getr().getg()*period/s1.getr().getd() - s2.getr().getg()*period/s2.getr().getd();      //steps
    result_r.init(r_0,period);
////////////////////////calculating the minimum in the second period/////////////////////////////////////
    p1.init(0,0);
    p2.init(0,0);
    do
    {
        t_p1max = p1.getpol(p1.getn()-1).getd();
        p1 = otimes(p1,s1.getr());
        p1 = oplus(p1,s1.getq());
    }while(p1.getpol(p1.getn()-1).getd() <= start_time + period);
    p1 = oplus(p1,s1.getp());
    do
    {
        t_p2max = p2.getpol(p2.getn()-1).getd();
        p2 = otimes(p2,s2.getr());
        p2 = oplus(p2,s2.getq());
    }while (p2.getpol(p2.getn()-1).getd() <= start_time + period);
    p2 = oplus(p2,s2.getp());
    std::cout<<t_p1max<<":"<<t_p2max<<std::endl;
    p_min = infinity;
    p1_n = p1.getn();
    p2_n = p2.getn();
    p1_index = p1_n-1;
    p2_index = p2_n-1;
//////////////////////////////search the minimum in the period//////////////////////////////////////////////////
    do
    {
 //       std::cout<<"p_min"<<std::endl<<p_min<<std::endl;
 //       std::cout<<"p_min2"<<std::endl<<(p1.getpol(p1_index).getg()-p2.getpol(p2_index).getg())<<std::endl;
 //       std::cout<<"p1.getpol(p1_index):"<<std::endl<<p1.getpol(p1_index)<<std::endl;
 //       std::cout<<"p2.getpol(p2_index):"<<std::endl<<p2.getpol(p2_index)<<std::endl;
        t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
        if ((p1_index == 0 && p2_index == 0) || t < start_time)    //If both polynomial is at the end
            break;
        if (t <= (MAX(t_p1max,t_p2max)))
            p_min = MIN(p_min,(p1.getpol(p1_index).getg()-p2.getpol(p2_index).getg()));
        if ((p1_index == 0 || p1.getpol(p1_index-1).getd() < p2.getpol(p2_index-1).getd()) && p2_index != 0)
            p2_index--;
        else if ((p2_index == 0 || p1.getpol(p1_index-1).getd() > p2.getpol(p2_index-1).getd()) && p1_index != 0)
            p1_index--;
        else
        {
            p1_index--;
            p2_index--;
        }
        std::cout<<"-----------------------------------------------------------------"<<std::endl;
    }while(true);
////////////////////////////////////calculating the p and q///////////////////////////////////////////////////////
    p_min += r_0;
    p1 = oplus(s1.getp(),p1);
    p2 = oplus(s2.getp(),p2);
    p1_n = p1.getn();
    p2_n = p2.getn();
//            std::cout<<"p1:"<<std::endl<<p1<<std::endl;
//        std::cout<<"p2:"<<std::endl<<p2<<std::endl;
    p1_index = p1_n-1;
    p2_index = p2_n-1;
    do
    {
//        std::cout<<"p1.getpol(p1_index):"<<std::endl<<p1.getpol(p1_index)<<std::endl;
//        std::cout<<"p2.getpol(p2_index):"<<std::endl<<p2.getpol(p2_index)<<std::endl;
        if(p_min > p1.getpol(p1_index).getg() - p2.getpol(p2_index).getg())
        {
//            std::cout<<"bingop1:"<<std::endl<<p1.getpol(p1_index)<<std::endl;
//            std::cout<<"bingop2:"<<std::endl<<p2.getpol(p2_index)<<std::endl;
            t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
            if (t <= (MAX(t_p1max,t_p2max)))
            {
                p_min = p1.getpol(p1_index).getg() - p2.getpol(p2_index).getg();
                temp_gd.init(p_min,t);
                std::cout<<"temp_gd:"<<std::endl<<temp_gd<<std::endl;
                if (t >= start_time )
                    result_q = oplus(result_q,temp_gd);
                else
                    result_p = oplus(result_p,temp_gd);
            }
        }
        if (p1_index == 0 && p2_index == 0)     //If both polynomial is at the end
            break;
        if ((p1_index == 0 || p1.getpol(p1_index-1).getd() < p2.getpol(p2_index-1).getd()) && p2_index != 0)
            p2_index--;
        else if ((p2_index == 0 || p1.getpol(p1_index-1).getd() > p2.getpol(p2_index-1).getd()) && p1_index != 0)
            p1_index--;
        else
        {
            p1_index--;
            p2_index--;
        }
        std::cout<<"-----------------------------------------------------------------"<<std::endl;
    }while(true);
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    result.init(result_p,result_q,result_r);
    return result;
}

poly odot_mark2(poly p1, poly p2)
{
    int p1_index = 0,p2_index = 0;
    long int t;
    poly result;
    gd temp_gd;
    int p1_n, p2_n;
    p1_n = p1.getn();
    p2_n = p2.getn();
    t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
    do
    {
        if (p1_index == p1_n-1 && p2_index == p2_n-1)
        {
            t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
            temp_gd.init(p1.getpol(p1_index).getg()+p2.getpol(p2_index).getg(),t);
            result = oplus(result,temp_gd);
            return result;
        }
        temp_gd.init(p1.getpol(p1_index).getg()+p2.getpol(p2_index).getg(),t);
        result = oplus(result,temp_gd);
        if ((p1_index == p1_n-1 || p1.getpol(p1_index).getd() > p2.getpol(p2_index).getd()) && p2_index != p2_n-1)
            p2_index++;
        else if ((p2_index == p2_n-1 || p1.getpol(p1_index).getd() < p2.getpol(p2_index).getd()) && p1_index != p1_n-1)
            p1_index++;
        else
        {
            p1_index++;
            p2_index++;
        }
        t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
        std::cout<<"p1:"<<std::endl<<p1.getpol(p1_index)<<std::endl;
        std::cout<<"p2:"<<std::endl<<p2.getpol(p2_index)<<std::endl;
        std::cout<<"t:"<<std::endl<<t<<std::endl;
        std::cout<<"-----------------------------------------------------------------"<<std::endl;
    }while(true);
}

serie odot_mark2(serie s1, serie s2)
{
    long int start_time;
    long int period;                    //the period
    long int r_0;                       //the steps
    int p1_index = 0,p2_index = 0;
    long int t,t_c,i,t_p1max,t_p2max;                       //i: index of period
    serie result;
    gd temp_gd,result_r;
    poly p1,p2,result_p,result_q;
    int p1_n, p2_n;
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
    start_time = MAX(s1.getq().getpol(0).getd(),s2.getq().getpol(0).getd());                  //max(T'(0), T''(0))
    period = lcm(s1.getr().getd(),s2.getr().getd());                                                //lcm(tau',tau'')
    std::cout<<"period:"<<period<<std::endl;
    r_0 = s1.getr().getg()*(period/s1.getr().getd()) + s2.getr().getg()*(period/s2.getr().getd());      //steps
    std::cout<<"r_0:"<<r_0<<std::endl;
    result_r.init(r_0,period);
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
    p1.init(0,infinity);
    p2.init(0,infinity);
    do
    {
        t_p1max = p1.getpol(p1.getn()-1).getd();
        p1 = otimes(p1,s1.getr());
        p1 = oplus(p1,s1.getq());
    }while(p1.getpol(p1.getn()-1).getd() <= start_time + period);
    p1 = oplus(p1,s1.getp());
    do
    {
        t_p2max = p2.getpol(p2.getn()-1).getd();
        p2 = otimes(p2,s2.getr());
        p2 = oplus(p2,s2.getq());
    }while (p2.getpol(p2.getn()-1).getd() <= start_time + period);
    p2 = oplus(p2,s2.getp());
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
    std::cout<<"p1:"<<std::endl<<p1<<std::endl;
    std::cout<<"p2:"<<std::endl<<p2<<std::endl;
    p1_n = p1.getn();
    p2_n = p2.getn();
    t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
    p1_index = 0;p2_index = 0;
    do
    {
        if (t >= start_time + period)
        {
            result.init(result_p,result_q,result_r);
            return result;
        }
        if (t < start_time)
        {
            temp_gd.init(p1.getpol(p1_index).getg()+p2.getpol(p2_index).getg(),t);
            result_p = oplus(result_p,temp_gd);
        }
        else
        {
            temp_gd.init(p1.getpol(p1_index).getg()+p2.getpol(p2_index).getg(),t);
            result_q = oplus(result_q,temp_gd);
        }
        if ((p1_index == p1_n-1 || p1.getpol(p1_index).getd() > p2.getpol(p2_index).getd()) && p2_index != p2_n-1)
            p2_index++;
        else if ((p2_index == p2_n-1 || p1.getpol(p1_index).getd() < p2.getpol(p2_index).getd()) && p1_index != p1_n-1)
            p1_index++;
        else
        {
            p1_index++;
            p2_index++;
        }
        t = MIN(p1.getpol(p1_index).getd(),p2.getpol(p2_index).getd());
        std::cout<<"-----------------------------------------------------------------"<<std::endl;
    }while(true);
    return result;
}
