mrfast.c

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

            mr_mip->cr[i]=mr_mip->t[i][j];
        scrt(_MIPP_ &mr_mip->chin,mr_mip->cr,z[j]); /* N*3*np*np/2  */
        if (compare(z[j],mr_mip->w4)>=0)
        { /* In higher half of range, so number is negative */
            subtract(_MIPP_ mr_mip->w3,z[j],z[j]);
            negify(z[j],z[j]);
        }
    }            /* np*np*N/4 */
    return np;
}

int mr_poly_mul(_MIPD_ int degx,big *x,int degy,big *y,big *z)
{ /*  Multiply two polynomials. The big arrays are of size degree */
    int i,j,newn,logn,np,degree;
    mr_utype inv,p,fac;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_FP_ROUNDING
    mr_large ip;
#endif
    degree=degx+degy;
    if (x==y) 
    {
        mr_poly_sqr(_MIPP_ degx,x,z);
        return degree;
    }
    newn=1; logn=0;
    while (degree+1>newn) { newn<<=1; logn++; }

    if (mr_mip->logN<logn)
        np=mr_fft_init(_MIPP_ logn,mr_mip->modulus,mr_mip->modulus,TRUE);
    else np=mr_mip->nprimes;

/* compute coefficients modulo fft primes */    
    for (i=0;i<np;i++)
    {                                              
        p=mr_mip->prime[i];
#ifdef MR_FP_ROUNDING
        ip=mr_invert(p);
#endif
        for (j=0;j<=degx;j++)
        {
            if (x[j]==NULL) mr_mip->wa[j]=0;
#ifdef MR_FP_ROUNDING
            else mr_mip->wa[j]=mr_sdiv(_MIPP_ x[j],p,ip,mr_mip->w1);      /* np*np*N/2 muldivs */
#else
            else mr_mip->wa[j]=mr_sdiv(_MIPP_ x[j],p,mr_mip->w1);      /* np*np*N/2 muldivs */
#endif
        }
        for (j=degx+1;j<newn;j++) mr_mip->wa[j]=0;
        mr_dif_fft(_MIPP_ logn,i,mr_mip->wa);              /* np*N*lgN     */

        for (j=0;j<=degy;j++)
        {
            if (y[j]==NULL) mr_mip->t[i][j]=0;
#ifdef MR_FP_ROUNDING
            else mr_mip->t[i][j]=mr_sdiv(_MIPP_ y[j],p,ip,mr_mip->w1);   /* np*np*N/2   */
#else
            else mr_mip->t[i][j]=mr_sdiv(_MIPP_ y[j],p,mr_mip->w1);   /* np*np*N/2   */
#endif
        }
        for (j=degy+1;j<newn;j++) mr_mip->t[i][j]=0;
        mr_dif_fft(_MIPP_ logn,i,mr_mip->t[i]);           /* np*N*lgN */

   /* multiply FFTs */
        for (j=0;j<newn;j++)                       /* np*N */
#ifdef MR_FP_ROUNDING
            imuldiv(mr_mip->wa[j],mr_mip->t[i][j],(mr_small)0,p,ip,(mr_small *)&mr_mip->t[i][j]);
#else
            muldiv(mr_mip->wa[j],mr_mip->t[i][j],(mr_small)0,p,(mr_small *)&mr_mip->t[i][j]);
#endif
        mr_dit_fft(_MIPP_ logn,i,mr_mip->t[i]);           /* np*N*lgN */

        inv=mr_mip->inverse[i];
        if (mr_mip->logN > logn)
        {
            fac=twop(mr_mip->logN-logn);
            inv=smul(fac,inv,p);
        }
        for (j=0;j<=degree;j++)                   /* np*N */
#ifdef MR_FP_ROUNDING
            imuldiv(mr_mip->t[i][j],inv,(mr_small)0,p,ip,(mr_small *)&mr_mip->t[i][j]);
#else
            muldiv(mr_mip->t[i][j],inv,(mr_small)0,p,(mr_small *)&mr_mip->t[i][j]);
#endif
    }
    mr_mip->check=OFF;
    mr_shift(_MIPP_ mr_mip->modulus,(int)mr_mip->modulus->len,mr_mip->w6);
       /* w6 = N.R */
    for (j=0;j<=degree;j++)
    {
        for (i=0;i<np;i++)
            mr_mip->cr[i]=mr_mip->t[i][j];
        scrt(_MIPP_ &mr_mip->chin,mr_mip->cr,mr_mip->w7); /* N*3*np*np/2 */
        divide(_MIPP_ mr_mip->w7,mr_mip->w6,mr_mip->w6);  /* z[j] may be too big for redc */ 
        redc(_MIPP_ mr_mip->w7,z[j]);
    }                                        /* np*np*N/4 */
    mr_mip->check=ON;
    return degree;
}

int mr_poly_sqr(_MIPD_ int degx,big *x,big *z)
{ /*  Multiply two polynomials. The big arrays are of size degree */
    int i,j,newn,logn,np,degree;
    mr_utype inv,p,fac;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_FP_ROUNDING
    mr_large ip;
#endif
    degree=2*degx;
    newn=1; logn=0;
    while (degree+1>newn) { newn<<=1; logn++; }

    if (mr_mip->logN<logn)
        np=mr_fft_init(_MIPP_ logn,mr_mip->modulus,mr_mip->modulus,TRUE);
    else np=mr_mip->nprimes;

/* compute coefficients modulo fft primes */    
    for (i=0;i<np;i++)
    {
        p=mr_mip->prime[i];
#ifdef MR_FP_ROUNDING
        ip=mr_invert(p);
#endif
        for (j=0;j<=degx;j++)
        {
            if (x[j]==NULL) mr_mip->t[i][j]=0;
#ifdef MR_FP_ROUNDING
            else mr_mip->t[i][j]=mr_sdiv(_MIPP_ x[j],p,ip,mr_mip->w1);
#else
            else mr_mip->t[i][j]=mr_sdiv(_MIPP_ x[j],p,mr_mip->w1);
#endif
        }
        for (j=degx+1;j<newn;j++) mr_mip->t[i][j]=0;
        mr_dif_fft(_MIPP_ logn,i,mr_mip->t[i]);

   /* multiply FFTs */
        for (j=0;j<newn;j++)
#ifdef MR_FP_ROUNDING
            imuldiv(mr_mip->t[i][j],mr_mip->t[i][j],(mr_small)0,p,ip,(mr_small *)&mr_mip->t[i][j]);
#else
            muldiv(mr_mip->t[i][j],mr_mip->t[i][j],(mr_small)0,p,(mr_small *)&mr_mip->t[i][j]);
#endif
        mr_dit_fft(_MIPP_ logn,i,mr_mip->t[i]);

        inv=mr_mip->inverse[i];
        if (mr_mip->logN > logn)
        { /* adjust 1/N log p for smaller N  */
            fac=twop(mr_mip->logN-logn);
            inv=smul(fac,inv,p);
        }
        for (j=0;j<=degree;j++)
#ifdef MR_FP_ROUNDING
            imuldiv(mr_mip->t[i][j],inv,(mr_small)0,p,ip,(mr_small *)&mr_mip->t[i][j]);
#else
            muldiv(mr_mip->t[i][j],inv,(mr_small)0,p,(mr_small *)&mr_mip->t[i][j]);
#endif
    }
    mr_mip->check=OFF;
    mr_shift(_MIPP_ mr_mip->modulus,(int)mr_mip->modulus->len,mr_mip->w6);
       /* w6 = N.R */
    for (j=0;j<=degree;j++)
    { /* apply CRT to each column */
        for (i=0;i<np;i++)
            mr_mip->cr[i]=mr_mip->t[i][j];

        scrt(_MIPP_ &mr_mip->chin,mr_mip->cr,mr_mip->w7);
        divide(_MIPP_ mr_mip->w7,mr_mip->w6,mr_mip->w6);  /* z[j] may be too big for redc */ 
        redc(_MIPP_ mr_mip->w7,z[j]);
    }
    mr_mip->check=ON;

    return degree;
}

static BOOL init_it(_MIPD_ int logn)
{ /* find primes, table of roots, inverses etc for new N */
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif

#ifdef MR_ITANIUM
    mr_small tm;
#endif
   
    zero(mr_mip->w16);
    mr_mip->w16->len=2; mr_mip->w16->w[0]=0; mr_mip->w16->w[1]=1;

    if (mr_fft_init(_MIPP_ logn,mr_mip->w16,mr_mip->w16,FALSE)!=3) return FALSE;

    mr_mip->const1=invers(mr_mip->prime[0],mr_mip->prime[1]);
    mr_mip->const2=invers(mr_mip->prime[0],mr_mip->prime[2]);
    mr_mip->const3=invers(mr_mip->prime[1],mr_mip->prime[2]);
    if (mr_mip->base==0) 
    {
#ifndef MR_NOFULLWIDTH
        mr_mip->msw=muldvd(mr_mip->prime[0],mr_mip->prime[1],(mr_small)0,&mr_mip->lsw);
#endif
    }
    else mr_mip->msw=muldiv(mr_mip->prime[0],mr_mip->prime[1],(mr_small)0,mr_mip->base,&mr_mip->lsw);
    mr_mip->logN=logn;
    return TRUE;
}

BOOL fastmultop(_MIPD_ int n,big x,big y,big z)
{ /* only return top n words... assumes x and y are n in length */
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    int len;
    mr_mip->check=OFF;
    fft_mult(_MIPP_ x,y,mr_mip->w0);
    mr_mip->check=ON;
    len=mr_lent(mr_mip->w0);
    mr_shift(_MIPP_ mr_mip->w0,n-len,mr_mip->w0);
    copy(mr_mip->w0,z);
    if (len<2*n) return TRUE;
    return FALSE;
}

void fft_mult(_MIPD_ big x,big y,big z)
{ /* "fast" O(n.log n) multiplication */
    int i,pr,xl,yl,zl,newn,logn;
    mr_small v1,v2,v3,ic,c1,c2,p,fac,inv;

#ifdef MR_ITANIUM
    mr_small tm;
#endif

#ifdef MR_FP
    mr_small dres;
#endif
    mr_unsign32 sz;
    mr_utype *w[3],*wptr,*dptr,*d0,*d1,*d2,t;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_FP_ROUNDING
    mr_large ip;
#endif
    if (mr_mip->ERNUM) return;
    if (y->len==0 || x->len==0) 
    {
        zero(z);
        return;
    }

    MR_IN(72)

    if (mr_notint(x) || mr_notint(y))
    {
        mr_berror(_MIPP_ MR_ERR_INT_OP);
        MR_OUT
        return;
    }  
    sz=((x->len&MR_MSBIT)^(y->len&MR_MSBIT));
    xl=(int)(x->len&MR_OBITS);
    yl=(int)(y->len&MR_OBITS);
    zl=xl+yl;
    if (xl<512 || yl<512)    /* should be 512 */
    { /* not worth it! */
        multiply(_MIPP_ x,y,z);
        MR_OUT
        return;
    }
    if (zl>mr_mip->nib && mr_mip->check)
    {
        mr_berror(_MIPP_ MR_ERR_OVERFLOW);
        MR_OUT
        return;
    }
    newn=1; logn=0;
    while (zl>newn) { newn<<=1; logn++;}
    if (logn>mr_mip->logN)     /*  2^(N+1) settings can be used for 2^N */
    { /* numbers too big for current settings */
        if (!init_it(_MIPP_ logn)) 
        {
            mr_berror(_MIPP_ MR_ERR_OUT_OF_MEMORY);
            MR_OUT
            return;
        }
    }
    if (newn>2*mr_mip->nib)
    {
        mr_berror(_MIPP_ MR_ERR_OVERFLOW);
        MR_OUT
        return;
    }

    d0=mr_mip->t[0]; d1=mr_mip->t[1]; d2=mr_mip->t[2];
    w[0]=mr_mip->wa; w[1]=mr_mip->wb; w[2]=mr_mip->wc;

    fac=twop(mr_mip->logN-logn);
    
    for (pr=0;pr<3;pr++)
    { /* multiply mod each prime */
        p=mr_mip->prime[pr];
        inv=mr_mip->inverse[pr];  
#ifdef MR_FP_ROUNDING
        ip=mr_invert(p);
#endif         
        if (fac!=1) inv=smul(fac,inv,p);  /* adjust 1/N mod p */
        
        dptr=mr_mip->t[pr];
        wptr=w[pr];

        for (i=0;i<xl;i++) dptr[i]=MR_REMAIN(x->w[i],p);
        for (i=xl;i<newn;i++) dptr[i]=0;
                                   
        mr_dif_fft(_MIPP_ logn,pr,dptr);


        if (x!=y)
        {
            if (!mr_mip->same || !mr_mip->first_one)
            {   
                for (i=0;i<yl;i++) wptr[i]=MR_REMAIN(y->w[i],p);
                for (i=yl;i<newn;i++) wptr[i]=0;
                mr_dif_fft(_MIPP_ logn,pr,wptr);
            } 
        }
        else wptr=dptr;

        for (i=0;i<newn;i++)
        {  /* "multiply" Fourier transforms */
#ifdef MR_FP_ROUNDING
            imuldiv(dptr[i],wptr[i],(mr_small)0,p,ip,(mr_small *)&dptr[i]);
#else
            muldiv(dptr[i],wptr[i],(mr_small)0,p,(mr_small *)&dptr[i]);
#endif
        }

        mr_dit_fft(_MIPP_ logn,pr,dptr);

        for (i=0;i<newn;i++)
        {  /* convert to mixed radix form     *
            * using chinese remainder thereom *
            * but first multiply by 1/N mod p */
#ifdef MR_FP_ROUNDING
            imuldiv(dptr[i],inv,(mr_small)0,p,ip,(mr_small *)&dptr[i]); 
#else
            muldiv(dptr[i],inv,(mr_small)0,p,(mr_small *)&dptr[i]); 
#endif
            if (pr==1)
            {
                t=d1[i]-d0[i];
                while (t<0) t+=mr_mip->prime[1];
                muldiv(t,mr_mip->const1,(mr_small)0,mr_mip->prime[1],(mr_small *)&d1[i]);
            } 
            if (pr==2)
            {
                t=d2[i]-d0[i];
                while (t<0) t+=mr_mip->prime[2];
                muldiv(t,mr_mip->const2,(mr_small)0,mr_mip->prime[2],(mr_small *)&t);
                t-=d1[i];
                while (t<0) t+=mr_mip->prime[2];
                muldiv(t,mr_mip->const3,(mr_small)0,mr_mip->prime[2],(mr_small *)&d2[i]);
            }
        }    
    }

    mr_mip->first_one=TRUE;

    zero(z);
    c1=c2=0;
    if (mr_mip->base==0) for (i=0;i<zl;i++)
    { /* propagate carries */
#ifndef MR_NOFULLWIDTH
        v1=d0[i];
        v2=d1[i];
        v3=d2[i];
        v2=muldvd(v2,mr_mip->prime[0],v1,&v1);
        c1+=v1;
        if (c1<v1) v2++;
        ic=c2+muldvd(mr_mip->lsw,v3,c1,&z->w[i]);
        c2=muldvd(mr_mip->msw,v3,ic,&c1);
        c1+=v2;
        if (c1<v2) c2++;
#endif
    }
    else for (i=0;i<zl;i++)
    { /* propagate carries */
        v1=d0[i];
        v2=d1[i];
        v3=d2[i];
#ifdef MR_FP_ROUNDING
        v2=imuldiv(v2,mr_mip->prime[0],v1+c1,mr_mip->base,mr_mip->inverse_base,&v1);
        ic=c2+imuldiv(mr_mip->lsw,v3,v1,mr_mip->base,mr_mip->inverse_base,&z->w[i]);
        c2=imuldiv(mr_mip->msw,v3,v2+ic,mr_mip->base,mr_mip->inverse_base,&c1);
#else
        v2=muldiv(v2,mr_mip->prime[0],(mr_small)(v1+c1),mr_mip->base,&v1);
        ic=c2+muldiv(mr_mip->lsw,v3,v1,mr_mip->base,&z->w[i]);
        c2=muldiv(mr_mip->msw,v3,(mr_small)(v2+ic),mr_mip->base,&c1);
#endif
    }
    z->len=(sz|zl); /* set length and sign of result */
    mr_lzero(z);
    MR_OUT
}

#endif

/*
main()
{
    big x,y,z,w;
    int i,j,k;
    miracl *mip=mirsys(1024,0);
    x=mirvar(0);
    y=mirvar(0);
    z=mirvar(0);
    w=mirvar(0);

    mip->IOBASE=16;
    bigbits(512*MIRACL,x);
    bigbits(512*MIRACL,y);


    multiply(x,x,z);

    cotnum(z,stdout);

    fft_mult(x,x,w);

    cotnum(w,stdout);
    if (compare(z,w)!=0) printf("Problems\n");

}
*/