C++ Poly::addterm方法代码示例

Poly pow(const Poly& f,int k)
{
    Poly u;
    int w,e,b;

    if (k==0)
    {
        u.addterm((ZZn)1,0);
        return u;
    }
    u=f;
    if (k==1) return u;

    e=k;
    b=0; while (k>1) {k>>=1; b++; }
    w=(1<<b);
    e-=w; w/=2;
    while (w>0)
    {
        u=(u*u);
        if (e>=w)
        {
           e-=w;
           u=(u*f);
        }
        w/=2; 
    }
    return u;
}

Poly PolyXY::F(const ZZn& y)
{
    Poly r;
    term *pos=NULL;
    int i,maxy=0;
    ZZn f;
    termXY *ptr=start;

    while (ptr!=NULL)
    {
        if (ptr->ny>maxy) maxy=ptr->ny;
        ptr=ptr->next;
    }

    // max y is max power of y present

    ZZn *pw=new ZZn[maxy+1];   // powers of y

    pw[0]=(ZZn)1;
    for (i=1; i<=maxy; i++)
        pw[i]=y*pw[i-1];

    ptr=start;
    while (ptr!=NULL)
    {
        pos=r.addterm(ptr->an*pw[ptr->ny],ptr->nx,pos);
        ptr=ptr->next;
    }

    delete [] pw;

    return r;
}

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

// Convert a polyxy object to a poly object by just taking the
// "x" value
Poly PolyXY::convert_x() const
{
    termXY *ptr=start;
    term *pos=NULL;
    Poly newpoly;
    while (ptr!=NULL)
    {
        pos = newpoly.addterm(ptr->an,ptr->nx,pos);
        ptr = ptr->next;
    }
    return newpoly;
}

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

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 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 operator*(const Poly& a,const Poly& b)
{
    int i,d,dega,degb,deg;
    BOOL squaring;
    ZZn t;
    Poly prod;
    term *iptr,*pos;
    term *ptr=b.start;

    squaring=FALSE;
    if (&a==&b) squaring=TRUE;

    dega=degree(a);
    deg=dega;
    if (!squaring)
    {
        degb=degree(b);
        if (degb<dega) deg=degb;
    }
    else degb=dega;
    if (deg>=FFT_BREAK_EVEN)      /* deg is minimum - both must be less than FFT_BREAK_EVEN */
    { // use fast methods 
        big *A,*B,*C;
        deg=dega+degb;     // degree of product
   
        A=(big *)mr_alloc(dega+1,sizeof(big));
        if (!squaring) B=(big *)mr_alloc(degb+1,sizeof(big));        
        C=(big *)mr_alloc(deg+1,sizeof(big));
        for (i=0;i<=deg;i++) C[i]=mirvar(0);
        ptr=a.start;
        while (ptr!=NULL)
        {
            A[ptr->n]=getbig(ptr->an);
            ptr=ptr->next;
        }

        if (!squaring)
        {
            ptr=b.start;
            while (ptr!=NULL)
            {
                B[ptr->n]=getbig(ptr->an);
                ptr=ptr->next;
            }
            mr_poly_mul(dega,A,degb,B,C);
        }
        else mr_poly_sqr(dega,A,C);
        pos=NULL;
        for (d=deg;d>=0;d--)
        {
            t=C[d];
            mr_free(C[d]);
            if (t.iszero()) continue;
            pos=prod.addterm(t,d,pos);
        }

        mr_free(C);
        mr_free(A);
        if (!squaring) mr_free(B);

        return prod;
    }

    if (squaring)
    { // squaring
        pos=NULL;
        while (ptr!=NULL)
        { // diagonal terms
            pos=prod.addterm(ptr->an*ptr->an,ptr->n+ptr->n,pos);
            ptr=ptr->next;
        }
        ptr=b.start;
        while (ptr!=NULL)
        { // above the diagonal
            iptr=ptr->next;
            pos=NULL;
            while (iptr!=NULL)
            {
                t=ptr->an*iptr->an;
                pos=prod.addterm(t+t,ptr->n+iptr->n,pos);
                iptr=iptr->next;
            }
            ptr=ptr->next; 
        }
    }
    else while (ptr!=NULL)
    {
        pos=NULL;
        iptr=a.start;
        while (iptr!=NULL)
        {
            pos=prod.addterm(ptr->an*iptr->an,ptr->n+iptr->n,pos);
            iptr=iptr->next;
        }
        ptr=ptr->next;
    }

    return prod;
}

//.........这里部分代码省略.........
        if (i>0)
        {
            for (j=i-1;j>0;j--)
            {
                if (!T[j].iszero())
                {
                    Acc=special(Acc,T[j]); // special karatsuba function
                    T[j].clear();          // multiply into accumulator
                }   
            }
            if (!T[0].iszero())
            {                              // check for a left-over linear poly
                Acc=Acc*T[0];
                T[0].clear();
            }
        }
        for (i=0;i<25;i++) T[i].clear();

        terms=degree(Acc);
        Float f,rem;
        Big whole;
        int nbits,maxbits=0;

        unstable=FALSE;
        for (i=terms;i>=0;i--)
        {
            f=Acc.coeff(i);
            if (f>0)
                 f+=makefloat(1,10000);
            else f-=makefloat(1,10000);
            whole=f.trunc(&rem);
            nbits=bits(whole);
            if (nbits>maxbits) maxbits=nbits;
            polly.addterm((ZZn)whole,i);
            if (fabs(rem)>makefloat(1,100)) 
            {
                unstable=TRUE; 
                break; 
            }
        }
        Acc.clear();
        if (!suppress) cout << endl;
        if (unstable) 
        {
            if (!suppress) cout << "Curve abandoned - numerical instability!" << endl;
            if (!suppress) cout << "Curve abandoned - double MIRACL precision and try again!" << endl;
            if (!suppress) cout << "finding a curve..." << endl;
            return FALSE;
        }
        if (!suppress) 
        {
            cout << polly << endl;
            cout << "Maximum precision required in bits= " << maxbits << endl;
        }
    }

// save space with smaller miracl

    mirexit();
    mip=mirsys(128,0);
    modulo(p);

    ECn pt,G;
    Big a,b,x,y;
    Big w,eps;
    int at;