}

Poly pow(const Poly& f,const Big& k,const Poly& m)
{
    Poly u,t,u2,table[16];
    Big w,e;
    int i,j,nb,n,nbw,nzs;
    if (k==0) 
    {
        u.addterm((ZZn)1,0);
        return u;
    }
    u=f%m;
    if (k==1) return u;
    if (degree(m)>=FFT_BREAK_EVEN)
    { /* Uses FFT for fixed base - much faster for large m */
        setpolymod(m);

        u2=modmult(u,u,m);
        table[0]=u;
        for (i=1;i<16;i++)
            table[i]=modmult(u2,table[i-1],m);
        nb=bits(k);
        if (nb>1) for (i=nb-2;i>=0;)
        { 
            n=window(k,i,&nbw,&nzs,5);
            for (j=0;j<nbw;j++)
                u=modmult(u,u,m);
            if (n>0) u=modmult(u,table[n/2],m);
            i-=nbw;
            if (nzs)
            {
                for (j=0;j<nzs;j++) u=modmult(u,u,m);
                i-=nzs;  
            }
        }

        fft_reset();
    }
    else
    {
        e=k;
        w=pow((Big)2,bits(e)-1);
        e-=w; w/=2;
        while (w>0)
        {
            u=(u*u)%m;
            if (e>=w)
            {
               e-=w;
               u=(u*f)%m;
            }
            w/=2;
        }
    } 
    return u;
}

int degree(const Poly& p)
{
    if (p.start==NULL) return 0;
    else return p.start->n;
}


BOOL iszero(const Poly& p) 
{
    if (degree(p)==0 && p.coeff(0)==0) return TRUE;
    else return FALSE;
}

BOOL isone(const Poly& p)
{
    if (degree(p)==0 && p.coeff(0)==1) return TRUE;
    else return FALSE;
}

Poly factor(const Poly& f,int d)
{
    Poly c,u,h,g=f;
    Big r,p=get_modulus();
//    int lim=0;
    while (degree(g) > d)
    {
//       lim++;
//       if (lim>10) break;
// random monic polynomial
       u.clear();
       u.addterm((ZZn)1,2*d-1);
       for (int i=2*d-2;i>=0;i--)
       {
          r=rand(p);
          u.addterm((ZZn)r,i);
       }
       r=(pow(p,d)-1)/2;

       c=pow(u,r,g);
       c.addterm((ZZn)(-1),0);

       h=gcd(c,g);
       if (degree(h)==0 || degree(h)==degree(g)) continue;
       if (2*degree(h)>degree(g))
            g=g/h;
       else g=h;
    }
    return g;
}

void Poly::clear()
{
    term *ptr;
    while (start!=NULL)
    {   
       ptr=start->next;
       delete start;
       start=ptr;
    }
    
}

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

Poly& Poly::operator=(const ZZn& m)
{
    clear();
    if (!m.iszero()) addterm(m,0);
    return *this;
}

Poly &Poly::operator=(const Poly& p)
{
    term *ptr,*pos=NULL;
    clear();
    ptr=p.start;
    while (ptr!=NULL)
    {  
        pos=addterm(ptr->an,ptr->n,pos);
        ptr=ptr->next;
    }    
    return *this;
}

Poly operator+(const Poly& a,const ZZn& z)
{
    Poly p=a;
    p.addterm(z,0);
    return p;
}

Poly operator-(const Poly& a,const ZZn& z)
{
    Poly p=a;
    p.addterm((-z),0);
    return p;
}


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

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

Poly& Poly::operator+=(const Poly& p)
{
    term *ptr,*pos=NULL;
    ptr=p.start;
    while (ptr!=NULL)
    {  
        pos=addterm(ptr->an,ptr->n,pos);
        ptr=ptr->next;
    }    
    return *this;
}

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

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

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

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

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

Poly& Poly::operator-=(const Poly& p)
{
    term *ptr,*pos=NULL;
    ptr=p.start;
    while (ptr!=NULL)
    {  
        pos=addterm(-(ptr->an),ptr->n,pos);
        ptr=ptr->next;
    }    
    return *this;
}
 
void Poly::multerm(const ZZn& a,int power)
{
    term *ptr=start;
    while (ptr!=NULL)
    {
        ptr->an*=a;
        ptr->n+=power;
        ptr=ptr->next;
    }
}

Poly invmodxn(const Poly& a,int n)
{ // Newton's method to find 1/a mod x^n
    int i,k;
    Poly b;
    k=0; while ((1<<k)<n) k++;
    b.addterm((ZZn)1/a.coeff(0),0); // important that a0 != 0
    for (i=1;i<=k;i++) b=modxn (2*b-a*b*b,1<<i);
    b=modxn(b,n);
    return b;
}

Poly modxn(const Poly& a,int n)
{ // reduce polynomial mod x^n
    Poly b;
    term* ptr=a.start;
    term *pos=NULL;
    while (ptr!=NULL && ptr->n>=n) ptr=ptr->next;
    while (ptr!=NULL)
    {
        pos=b.addterm(ptr->an,ptr->n,pos);
        ptr=ptr->next;
    }
    return b;
}

Poly divxn(const Poly& a,int n)
{ // divide polynomial by x^n
    Poly b;
    term *ptr=a.start;
    term *pos=NULL;
    while (ptr!=NULL)
    {
        if (ptr->n>=n)
            pos=b.addterm(ptr->an,ptr->n-n,pos);
        else break;
        ptr=ptr->next;
    }
    return b;
}

Poly mulxn(const Poly& a,int n)
{ // multiply polynomial by x^n
    Poly b;
    term *ptr=a.start;
    term *pos=NULL;
    while (ptr!=NULL)
    {
        pos=b.addterm(ptr->an,ptr->n+n,pos);
        ptr=ptr->next;
    }
    return b;
}

Poly reverse(const Poly& a)
{
    term *ptr=a.start;
    int deg=degree(a);
    Poly b;
    while (ptr!=NULL)
    {
        b.addterm(ptr->an,deg-ptr->n);
        ptr=ptr->next;
    } 
    return b;
}

// add term to polynomial. The pointer pos remembers the last
// accessed element - this is faster. 
// Polynomial is stored with large powers first, down to low powers
// e.g. 9x^6 + x^4 + 3x^2 + 1

term* Poly::addterm(const ZZn& a,int power,term *pos)
{
    term* newone;  
    term* ptr;
    term *t,*iptr;
    ptr=start;
    iptr=NULL;
    if (a.iszero()) return pos;
// quick scan through to detect if term exists already
// and to find insertion point
    if (pos!=NULL) ptr=pos;      // start looking from here
    while (ptr!=NULL) 
    { 
        if (ptr->n==power)
        {
            ptr->an+=a;
            if (ptr->an.iszero()) 
            { // delete term
                if (ptr==start)
                { // delete first one
                    start=ptr->next;
                    delete ptr;
                    return start;
                }
                iptr=ptr;
                ptr=start;
                while (ptr->next!=iptr)ptr=ptr->next;
                ptr->next=iptr->next;
                delete iptr;
                return ptr;
            }
            return ptr;
        }
        if (ptr->n>power) iptr=ptr;
        else break;
        ptr=ptr->next;
    }
    newone=new term;
    newone->next=NULL;
    newone->an=a;
    newone->n=power;
    pos=newone;
    if (start==NULL)
    {
        start=newone;
        return pos;
    }

// insert at the start

    if (iptr==NULL)
    { 
        t=start;
        start=newone;
        newone->next=t;
        return pos;
    }

// insert new term

    t=iptr->next;
    iptr->next=newone;
    newone->next=t;
    return pos;    
}

// A function to differentiate a polynomial
Poly differentiate(const Poly& orig)
{
    Poly newpoly;
    term *ptr = orig.start;
    term *pos = NULL;
    int power;

    while (ptr!=NULL)
    {   
        power = ptr->n;
        if(ptr->n > 0) 
            pos = newpoly.addterm(ptr->an*power,ptr->n-1,pos);
            
        ptr = ptr->next;
    }
    
    return newpoly;
}

void makemonic(Poly& p)
{
    term *ptr = p.start;
    p.multerm((ZZn)1/ptr->an,0);
}

// The extended euclidean algorithm
// The result is returned in an array of Polys, with the gcd
// in first place, then the two coefficients
void egcd(Poly result[], const Poly& u, const Poly& v)
{       
    Poly u1, u2, u3, v1, v2, v3, t1, t2, t3, zero, q;
    u1 = 1;  
    u2 = 0;
    u3 = u;
    v1 = 0;
    v2 = 1;
    v3 = v;
    zero = 0;
    term *ptr;
            
    while(v3 != zero) {
        q = u3/v3;
        t1 = u1 - v1*q;
        t2 = u2 - v2*q;
        t3 = u3 - v3*q;

        u1 = v1;
        u2 = v2;
        u3 = v3;

        v1 = t1;
        v2 = t2;
        v3 = t3;
    }
    result[0] = u3;
    result[1] = u1;
    result[2] = u2;

}

ostream& operator<<(ostream& s,const Poly& p)
{
    BOOL first=TRUE;
    ZZn a;
    term *ptr=p.start;
    if (ptr==NULL)
    { 
        s << "0";
        return s;
    }
    while (ptr!=NULL)
    {
        a=ptr->an;
        if ((Big)a<get_modulus()/2)
        {
            if (!first) s << " + ";
        }
        else 
        {
           a=(-a);
           s << " - ";
        }
        if (ptr->n==0) 
           s << (Big)a; 
        else 
        {
            if (a!=(ZZn)1)  s << (Big)a << "*x"; 
            else            s << "x";
            if (ptr->n!=1)  s << "^" << ptr->n;
        }
        first=FALSE;
        ptr=ptr->next;
    }
    return s;
}