ps_zzn.cpp

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

/*
 * C++ class to implement a power series type and to allow 
 * arithmetic on it
 *
 * WARNING: This class has been cobbled together for a specific use with
 * the MIRACL library It is not complete, and may not work in other 
 * applications
 *
 * See Knuth The Art of Computer Programming Vol.2, Chapter 4.7 
 */

#include "ps_zzn.h"

#define FFT

//
// all calulations are mod x^psN
// Power series is stored {as offset, pwr and a0+a1.x+a2.x^2+a3.x^3....}
// where the power of x is to be multiplied by pwr, and the whole power
// series is to be divided by x^offset
//

int psN;

//
// Copy Constructor
//

Ps_ZZn::Ps_ZZn(const Ps_ZZn& p)
{
    term_ps_zzn *ptr=p.start;
    term_ps_zzn *pos=NULL; 
    int pw;
    start=NULL;
    offset=p.offset;
    pwr=p.pwr;
    while (ptr!=NULL)
    {
        pw=ptr->n*p.pwr-p.offset;   // conversion needed
        if (pw>=psN) break;
        pos=addterm(ptr->an,pw,pos);
        ptr=ptr->next;
    }
}

//
// decompresses PS by reducing pwr
//

void Ps_ZZn::decompress(int m)
{ // m is divisor of current pwr
  // e.g pwr is 6 and PS = 1 + x + x^2
  // If m=2 then pwr becomes 3 and PS = 1 +x^2 +x^4.... 
    term_ps_zzn *ptr=start;
    if (start==NULL || m==1 || pwr==1) return; // it is fully decompressed
    while (ptr!=NULL)
    {
        ptr->n*=m;
        ptr=ptr->next;
    }
    pwr/=m;
}

//
//  Sets new pwr value. Parameter must exactly divide 
//  all powers in the series
//

void Ps_ZZn::compress(int p)
{
    term_ps_zzn *ptr=start;
    if (p==1) return;
    while (ptr!=NULL)
    {
        ptr->n/=p;
        ptr=ptr->next;  
    }
    pwr=p;
}

//
// insert missing terms with 0 coefficients
//

void Ps_ZZn::pad()
{ // insert any missing 0 coefficients 
    int i=0;
    term_ps_zzn *ptr=start;
    term_ps_zzn *pos=NULL;
    while (ptr!=NULL)
    {
        while (i<ptr->n)
        {
            pos=addterm((ZZn)0,i*pwr-offset,pos);
            i++;
        }
        i++;
        ptr=ptr->next;
    }
    while (i<psN/pwr) 
    {
        pos=addterm((ZZn)0,i*pwr-offset,pos);
        i++;
    }
}

//
// Find max coefficient in PS
//

int Ps_ZZn::max()
{
    int b,m=0;
    term_ps_zzn *ptr=start;
    while (ptr!=NULL)
    {
        b=bits(ptr->an);
        if (b>m) m=b;
        ptr=ptr->next;
    }
    return m;
}

//
// set new offset, so first power of x is 0
//

void Ps_ZZn::norm()
{
    int m;
    term_ps_zzn *ptr;
    if (start==NULL) return;
//
// remove any leading 0 terms
//
    while (start->an.iszero()) 
    {
        ptr=start->next;
        delete start;
        start=ptr;
        if (start==NULL) return;
    }
    ptr=start;
    m=start->n;
    if (m!=0)
    {
        offset-=m*pwr;
        while (ptr!=NULL)
        {
            ptr->n-=m;
            ptr=ptr->next;
        }
    }
}

//
// Dedekind Eta function (1-x)(1-x^2)(1-x^3)....
//

Ps_ZZn eta()
{ // simple repeating pattern
    BOOL even;
    int one,ce,co,c;
    term_ps_zzn *pos=NULL;
    Ps_ZZn n;
    n.addterm((ZZn)1,0);
    n.addterm((ZZn)-1,1);
    n.addterm((ZZn)-1,2); 
    ce=2;co=1;
    even=TRUE;
    c=2;
    one=1;
    while (c<psN)
    {
        if (even)
        {
            c+=(ce+1);
            ce+=2;
            pos=n.addterm((ZZn)one,c,pos);
            even=FALSE;
        }
        else
        {
            c+=(co+1);
            co+=1;
            pos=n.addterm((ZZn)one,c,pos);
            even=TRUE;
            one=(-one);
        }
    }
    return n;
}

//
// Checks if a power series is 0 or an integer
//

BOOL Ps_ZZn::IsInt() const
{
    if (start==NULL) return TRUE;
    if (one_term() && offset==0) return TRUE;
    return FALSE;
}

//
// return TRUE if zero or one term only in PS
//

BOOL Ps_ZZn::one_term() const
{
    int t=0;
    term_ps_zzn *ptr=start;
    if (start==NULL) return TRUE;
    while (ptr!=NULL)
    {
        if (!ptr->an.iszero()) t++;
        if (t>1) return FALSE;
        ptr=ptr->next;
    }
    return TRUE;
}

//
// add a term to a PS
//

term_ps_zzn* Ps_ZZn::addterm(const ZZn& a,int power,term_ps_zzn* pos)
{
    term_ps_zzn* newone;  
    term_ps_zzn* ptr;
    term_ps_zzn *t,*iptr;
    int dc,pw;
    ptr=start;
    iptr=NULL;

//
// intelligently determine the most compressed form to use
// for example if coefficient a=1 always, and power = -7 -5 -3 -1 1 3....
// then set pwr=2, offset=7  and PS = 1 + x + x^2
//
    pw=power+offset;
    if (one_term() && pw!=0)
    { // when PS has only one term, pwr is undefined
        if (pw<0)
            pwr=-pw;
        else pwr=pw; 
    }   

    dc=igcd(pw,pwr);

    if (dc != pwr) decompress(pwr/dc);
    power=pw/pwr;
// quick scan through to detect if term exists already
// and to find insertion point
   if (pos!=NULL) ptr=pos;
    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;
                    norm();
                    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;   // determines order
        else break;
        ptr=ptr->next;
    }
    newone=new term_ps_zzn;
    newone->next=NULL;
    newone->an=a;
    newone->n=power;
    pos=newone;
    if (start==NULL)
    {
        start=newone;
        norm();
        return pos;
    }

// insert at the start

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

// insert new term

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

//
// Destructor
//

Ps_ZZn::~Ps_ZZn()
{
    term_ps_zzn *nx;
    while (start!=NULL)
    {
        nx=start->next;
        delete start;
        start=nx;
    }
}

//
// get coefficient of actual power
//

ZZn Ps_ZZn::coeff(int power) const
{
    ZZn c=0;
    term_ps_zzn *ptr=start;
    if ((power+offset)%pwr != 0) return c;  // no such term
    power=(power+offset)/pwr;

    while (ptr!=NULL)
    {
        if (ptr->n==power)
        {
            c=ptr->an;
            return c;
        }
        ptr=ptr->next;
    }
    return c;
}

//
// get coefficient of "Internal" power
//

ZZn Ps_ZZn::cf(int power) const
{
    ZZn c=0;
    term_ps_zzn *ptr=start;

    while (ptr!=NULL)
    {
        if (ptr->n==power)
        {
            c=ptr->an;
            return c;
        }
        ptr=ptr->next;
    }
    return c;
}

//
// Zeroise PS and reclaim space
//

void Ps_ZZn::clear()
{
    term_ps_zzn *ptr;
    while (start!=NULL)
    {
        ptr=start->next;
        delete start;
        start=ptr;
    }
    offset=0;
    pwr=1;
}

// Note: real power = internal power * pwr - offset 

Ps_ZZn& Ps_ZZn::operator+=(const Ps_ZZn& p)
{
    term_ps_zzn *ptr=p.start;
    term_ps_zzn *pos=NULL;
    int pw;   
    while (ptr!=NULL)
    {
        pw=ptr->n*p.pwr-p.offset;   // convert compressed to real
        if (pw>=psN) break;
        pos=addterm(ptr->an,pw,pos);
        ptr=ptr->next;
    }
    return *this;
}

Ps_ZZn operator-(const Ps_ZZn& p)
{
    Ps_ZZn r=p;
    term_ps_zzn *ptr=r.start;
    while (ptr!=NULL)
    {
        ptr->an=(-ptr->an);
        ptr=ptr->next;
    }
    return r;
}

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

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