ps_big.cpp

miracl-大数运算库,大家使用有什么问题请多多提意见

}

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

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

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

Ps_Big& Ps_Big::operator*=(Ps_Big& b)
{
    term_ps_big *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]=(ptr->an).getbig();
            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]=(bptr->an).getbig();
            bptr=bptr->next;
       } 
       mr_ps_big_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_Big operator*(Ps_Big& a,Ps_Big& b)
{
    Ps_Big prod;
    Big t;
    term_ps_big *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]=(aptr->an).getbig();
            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]=(bptr->an).getbig();
            bptr=bptr->next;
        }
        mr_ps_big_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_Big operator/(const Big& num,const Ps_Big& b)
{
    Ps_Big quo;
    term_ps_big *pos=NULL;
    int pw=b.pwr;
    Big 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_Big operator/(Ps_Big& a,Ps_Big& b)
{
    Ps_Big quo;
    term_ps_big *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;

    Big 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_Big& Ps_Big::operator/=(const Big& x)
{
    term_ps_big *ptr=start;
    while (ptr!=NULL)
    {
        ptr->an/=x;
        ptr=ptr->next;
    }
    return *this;
}

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

Ps_Big operator/(const Ps_Big& a,const Big& b)
{
    Ps_Big quo;
    quo=a;
    quo/=b;
    return quo;
}

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

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

Ps_Big& Ps_Big::operator=(const Ps_Big& p)
{
    term_ps_big *ptr,*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_Big& Ps_Big::operator*=(const Big& x)
{
    term_ps_big *ptr=start;
    if (x.iszero()) 
    {
        clear();
        return *this;
    }
    while (ptr!=NULL)
    {
        ptr->an*=x;
        ptr=ptr->next;
    }
    return *this;
}

Ps_Big operator*(const Big& z,const Ps_Big &p)
{
    Ps_Big r=p;
    r*=z;
    return r;
}

Ps_Big operator*(const Ps_Big &p,const Big& z)
{
    Ps_Big r=p;
    r*=z;
    return r;
}

Ps_Big pow(Ps_Big &a,int k)
{
    Ps_Big 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_Big power(const Ps_Big &a,int e)
{ // return f(a^e), e is +ve
    Ps_Big res;
    term_ps_big *ptr,*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_Big& p)
{
    BOOL first=TRUE;
    Big a;
    term_ps_big *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 (a < (Big)0) 
        {
            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==1) s << "x";
        else      s << a << "*x";
    
        if (pw!=1) s << "^" << pw;
        ptr=ptr->next;
    }
    return s;
}