Ps_ZZn operator+(const Ps_ZZn& a,const Ps_ZZn& b)
{
    Ps_ZZn sum=a;
    sum+=b;
    return sum;
}

Ps_ZZn operator-(const Ps_ZZn& a,const Ps_ZZn& b)
{
    Ps_ZZn sum=a;
    sum-=b;
    return sum;
}

//
// In situ multiplication. Coeeficients in multiplicand
// are directly overwritten. Saves a lot of copying.
//

Ps_ZZn& Ps_ZZn::operator*=(Ps_ZZn& b)
{
    term_ps_zzn *ptr,*bptr;
    int g,d,deg,ka,kb;
    if (IsInt())
    {
        if (start!=NULL) *this = start->an*b;
        return *this;
    }
    if (b.IsInt())
    {
        if (b.start==NULL) clear();
        else *this *=(b.start->an);
        return *this;
    }
    g=igcd(pwr,b.pwr);
    deg=psN/g;


#ifdef FFT
    if (deg>=FFT_BREAK_EVEN)
    {
        big *A,*B;
        A=(big *)mr_alloc(deg,sizeof(big));
        B=(big *)mr_alloc(deg,sizeof(big));
        ka=pwr/g;
        decompress(ka);
        pad();
        ptr=start;
        while (ptr!=NULL)
        {
            d=ptr->n;
            if (d>=deg) break;
            A[d]=getbig(ptr->an);
            ptr=ptr->next;
        }
        kb=b.pwr/g;
        b.decompress(kb);
        bptr=b.start;
        while (bptr!=NULL)
        {
            d=bptr->n;
            if (d>=deg) break;
            B[d]=getbig(bptr->an);
            bptr=bptr->next;
       } 
       mr_ps_zzn_mul(deg,A,B,A);
       mr_free(B); mr_free(A);
       b.compress(kb);
    }
    else
    {
#endif
        *this=*this*b;
        return *this;
#ifdef FFT
    }
#endif
    norm();
    offset+=b.offset;
    return *this;
}

Ps_ZZn operator*(Ps_ZZn& a,Ps_ZZn& b)
{
    Ps_ZZn prod;
    ZZn t;
    term_ps_zzn *aptr,*bptr,*pos;
    int ka,kb,i,pa,pb,d,deg,g;

    if (a.IsInt())
    {
        if (a.start==NULL) return prod;
        else return (a.start->an*b);
    }

    if (b.IsInt())
    {
        if (b.start==NULL) return prod;
        else return (b.start->an*a);
    }

    g=igcd(a.pwr,b.pwr);

    deg=psN/g;

#ifdef FFT
    if (deg>=FFT_BREAK_EVEN)
    { // use fast methods
        big *A,*B,*C;
        A=(big *)mr_alloc(deg,sizeof(big));
        B=(big *)mr_alloc(deg,sizeof(big));
        C=(big *)mr_alloc(deg,sizeof(big));
        char *memc=(char *)memalloc(deg);
        for (i=0;i<deg;i++) C[i]=mirvar_mem(memc,i);
       
        ka=a.pwr/g;
        a.decompress(ka);

        aptr=a.start;
        while (aptr!=NULL)
        {
            d=aptr->n;
            if (d >= deg) break;
            A[d]=getbig(aptr->an);
            aptr=aptr->next;
        }
        kb=b.pwr/g;
        b.decompress(kb);
        bptr=b.start;
        while (bptr!=NULL)
        {
            d=bptr->n;
            if (d >= deg) break;
            B[d]=getbig(bptr->an);
            bptr=bptr->next;
        }
        mr_ps_zzn_mul(deg,A,B,C);
        pos=NULL;
        for (d=0;d<deg;d++)
        {
            t=C[d];
            if (t.iszero()) continue;
            pos=prod.addterm(t,d*g,pos);
        }
        memkill(memc,deg);
        a.compress(ka); b.compress(kb);
        mr_free(C); mr_free(B); mr_free(A); 
    }
    else
    {
#endif
        bptr=b.start;

        while (bptr!=NULL)
        {
            aptr=a.start;
            pb=bptr->n*b.pwr-b.offset;
            pos=NULL;
            while (aptr!=NULL)
            {
                pa=aptr->n*a.pwr-a.offset;
                if (pb+pa>=psN) break;
                pos=prod.addterm(aptr->an*bptr->an,pa+pb,pos);
                aptr=aptr->next;
            }
            bptr=bptr->next;
        }
#ifdef FFT
    }
#endif

    if (prod.start!=NULL) prod.offset=a.offset+b.offset;
    return prod;
}

Ps_ZZn operator/(const ZZn& num,const Ps_ZZn& b)
{
    Ps_ZZn quo;
    term_ps_zzn *pos=NULL;
    int pw=b.pwr;
    ZZn w,v0=b.cf(0);

    for (int n=0;n<psN/pw;n++)
    {
        if ((n*pw+b.offset)>=psN) break;
        if (n==0) w=num;
        else w=0;
        for (int k=0;k<n;k++)
            w-=quo.cf(k)*b.cf(n-k);
        pos=quo.addterm(w/v0,n,pos); 
    }

    quo.pwr=pw;
    quo.offset=(-b.offset);
    return quo;
}

Ps_ZZn operator/(Ps_ZZn& a,Ps_ZZn& b)
{
    Ps_ZZn quo;
    term_ps_zzn *pos=NULL;
    int ka,kb,pw;

    ka=a.pwr;
    kb=b.pwr;

    if (ka!=kb)
    {
        pw=1;
        a.decompress(a.pwr);
        b.decompress(b.pwr);
    }
    else pw=ka;

    ZZn w,v0=b.cf(0);

    for (int n=0;n<psN/pw;n++)
    {
        if (n*pw+b.offset-a.offset>=psN) break;
        w=a.cf(n);
        for (int k=0;k<n;k++)
            w-=quo.cf(k)*b.cf(n-k);
        pos=quo.addterm(w/v0,n,pos); 
    }
    if (ka!=kb)
    {
       a.compress(ka); 
       b.compress(kb);
    }
    else quo.pwr=pw;

    quo.offset=a.offset-b.offset;
    return  quo;
}

Ps_ZZn& Ps_ZZn::operator/=(const ZZn& x)
{
    term_ps_zzn *ptr=start;
    while (ptr!=NULL)
    {
        ptr->an/=x;
        ptr=ptr->next;
    }
    return *this;
}

Ps_ZZn& Ps_ZZn::operator/=(int m)
{
    term_ps_zzn *ptr=start;
    while (ptr!=NULL)
    {
        ptr->an/=m;
        ptr=ptr->next;
    }
    return *this;
}

Ps_ZZn operator/(const Ps_ZZn& a,const ZZn& b)
{
    Ps_ZZn quo;
    quo=a;
    quo/=b;
    return quo;
}

Ps_ZZn operator/(const Ps_ZZn& a,int m)
{
    Ps_ZZn quo;
    quo=a;
    quo/=m;
    return quo;
}


Ps_ZZn& Ps_ZZn::operator=(int m)
{
    clear();
    if (m!=0) addterm((ZZn)m,0);
    return *this;
}

Ps_ZZn& Ps_ZZn::operator=(const Ps_ZZn& p)
{
    term_ps_zzn *ptr;
    term_ps_zzn *pos=NULL;
    clear();
    int pw;

    pwr=p.pwr;
    offset=p.offset;
    ptr=p.start;
    while (ptr!=NULL)
    {
        pw=ptr->n*p.pwr-p.offset;
        if (pw>=psN) break;
        pos=addterm(ptr->an,pw,pos);
        ptr=ptr->next;
    }

    return *this;
}

Ps_ZZn& Ps_ZZn::operator*=(const ZZn& x)
{
    term_ps_zzn *ptr=start;
    if (x.iszero()) 
    {
        clear();
        return *this;
    }
    while (ptr!=NULL)
    {
        ptr->an*=x;
        ptr=ptr->next;
    }
    return *this;
}

Ps_ZZn operator*(const ZZn& z,const Ps_ZZn &p)
{
    Ps_ZZn r=p;
    r*=z;
    return r;
}

Ps_ZZn operator*(const Ps_ZZn &p,const ZZn& z)
{
    Ps_ZZn r=p;
    r*=z;
    return r;
}

Ps_ZZn pow(Ps_ZZn &a,int k)
{
    Ps_ZZn res;
    int w,e;

    if (k==0) 
    {
        res=1;
        return res;
    }
    res=a;
    if (k==1) return res;

    e=k; w=1;
    while (w<=e) w<<=1;
    w>>=1; 
    e-=w; w/=2;
    while (w>0)
    {
        res*=res;
        if (e>=w)
        {
            e-=w;
            res*=a;
        }
        w/=2;
    }
    res.offset*=k;
    
    return res;
}


Ps_ZZn power(const Ps_ZZn &a,int e)
{ // return f(a^e), e is +ve
    Ps_ZZn res;
    term_ps_zzn *ptr;
    term_ps_zzn *pos=NULL;

    ptr=a.start;
    while (ptr!=NULL)
    {
        if ((ptr->n*a.pwr)*e < psN)
        {
            pos=res.addterm(ptr->an,ptr->n,pos);
            ptr=ptr->next;
        }
        else break;
    }

    res.pwr=e*a.pwr;
    res.offset=e*a.offset;

    return res;
}

ostream& operator<<(ostream& s,const Ps_ZZn& p)
{
    BOOL first=TRUE;
    ZZn a;
    term_ps_zzn *ptr=p.start;
    int pw;

    if (ptr==NULL)
    {
        s << "0";
        return s;
    }

    while (ptr!=NULL)
    {
        a=ptr->an;
        if (a.iszero()) 
        {
            ptr=ptr->next;
            continue;
        }
        if ((Big)a < get_modulus()/2) 
        {
            a=(-a);
            s << " - ";
        }
        else if (!first) s << " + ";

        first=FALSE;
        pw=ptr->n*p.pwr-p.offset;
        if (pw==0)
        {
            s << a;
            ptr=ptr->next;
            continue;
        } 

        if (a==(ZZn)1) s << "x";
        else      s << a << "*x";
    
        if (pw!=1) s << "^" << pw;
        ptr=ptr->next;
    }
    return s;
}