mrfast.c

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

#else
#ifdef MR_FP_ROUNDING
                imuldiv(w,data[j],(mr_small)0,prime,iprime,(mr_small *)&temp);
#else
                muldiv(w,data[j],(mr_small)0,prime,(mr_small *)&temp);
#endif
#endif
                data[j]=data[i]-temp;
                if (data[j]<0) data[j]+=prime;
                data[i] += temp;
                if (data[i]>=prime) data[i]-=prime;
            }
#endif
        }
        mmax=istep;
    }
}

static void modxn_1(_MIPD_ int n,int deg,big *x)
{  /* set X (of degree deg) =X mod x^n-1 = X%x^n + X/x^n */
    int i;
    for (i=0;n+i<=deg;i++)
    {
        nres_modadd(_MIPP_ x[i],x[n+i],x[i]);
        zero(x[n+i]);
    }
}

BOOL mr_poly_rem(_MIPD_ int dg,big *G,big *R)
{ /* G is a polynomial of degree dg - G is overwritten */
    int i,j,newn,logn,np,n;
    mr_utype p,inv,fac;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_FP_ROUNDING
    mr_large ip;
#endif

    n=mr_mip->degree;  /* degree of modulus */
    if (n==0) return FALSE; /* the preset tables have been destroyed */
    np=mr_mip->nprimes;

    newn=1; logn=0;
    while (2*n>newn) { newn<<=1; logn++; }

    for (i=0;i<np;i++)
    {
        p=mr_mip->prime[i];
#ifdef MR_FP_ROUNDING
        ip=mr_invert(p);
        for (j=n;j<=dg;j++)
            mr_mip->t[i][j-n]=mr_sdiv(_MIPP_ G[j],p,ip,mr_mip->w1); 
#else
        for (j=n;j<=dg;j++)
            mr_mip->t[i][j-n]=mr_sdiv(_MIPP_ G[j],p,mr_mip->w1); 
#endif
        for (j=dg-n+1;j<newn;j++) mr_mip->t[i][j]=0;
        mr_dif_fft(_MIPP_ logn,i,mr_mip->t[i]);
        for (j=0;j<newn;j++)
#ifdef MR_FP_ROUNDING
            imuldiv(mr_mip->t[i][j],mr_mip->s1[i][j],(mr_small)0,p,ip,(mr_small *)&mr_mip->t[i][j]);
#else
            muldiv(mr_mip->t[i][j],mr_mip->s1[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 N/2, N/4 etc */
            fac=twop(mr_mip->logN-logn);
            inv=smul(fac,inv,p);
        }
        for (j=0;j<n;j++)
#ifdef MR_FP_ROUNDING
            imuldiv(mr_mip->t[i][j+n-1],inv,(mr_small)0,p,ip,(mr_small *)&mr_mip->t[i][j+n-1]);
#else
            muldiv(mr_mip->t[i][j+n-1],inv,(mr_small)0,p,(mr_small *)&mr_mip->t[i][j+n-1]);
#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<n;j++)
    {
        for (i=0;i<np;i++)
            mr_mip->cr[i]=mr_mip->t[i][j+n-1];
        scrt(_MIPP_ &mr_mip->chin,mr_mip->cr,mr_mip->w7);
        divide(_MIPP_ mr_mip->w7,mr_mip->w6,mr_mip->w6);  /* R[j] may be too big for redc */ 
        redc(_MIPP_ mr_mip->w7,R[j]);
    }
    mr_mip->check=ON;

    for (i=0;i<np;i++)
    {
        p=mr_mip->prime[i];
#ifdef MR_FP_ROUNDING
        ip=mr_invert(p);
        for (j=0;j<n;j++)
            mr_mip->t[i][j]=mr_sdiv(_MIPP_ R[j],p,ip,mr_mip->w1);
#else
        for (j=0;j<n;j++)
            mr_mip->t[i][j]=mr_sdiv(_MIPP_ R[j],p,mr_mip->w1);
#endif
        for (j=n;j<1+newn/2;j++) mr_mip->t[i][j]=0;

        mr_dif_fft(_MIPP_ logn-1,i,mr_mip->t[i]);  /* Note: Half size */
        for (j=0;j<newn/2;j++)
#ifdef MR_FP_ROUNDING
            imuldiv(mr_mip->t[i][j],mr_mip->s2[i][j],(mr_small)0,p,ip,(mr_small *)&mr_mip->t[i][j]);
#else
            muldiv(mr_mip->t[i][j],mr_mip->s2[i][j],(mr_small)0,p,(mr_small *)&mr_mip->t[i][j]);
#endif
        mr_dit_fft(_MIPP_ logn-1,i,mr_mip->t[i]);

        inv=mr_mip->inverse[i];
        if (mr_mip->logN > logn-1)
        {
            fac=twop(mr_mip->logN-logn+1);
            inv=smul(fac,inv,p);
        }
        for (j=0;j<n;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
    }

    modxn_1(_MIPP_ newn/2,dg,G);    /* G=G mod 2^x - 1 */

    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<n;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);
        divide(_MIPP_ mr_mip->w7,mr_mip->w6,mr_mip->w6);  /* R[j] may be too big for redc */ 
        redc(_MIPP_ mr_mip->w7,R[j]);
        nres_modsub(_MIPP_ G[j],R[j],R[j]);

    }
    mr_mip->check=ON;

    return TRUE;
}

void mr_polymod_set(_MIPD_ int n, big *rf,big *f)
{ /* n is degree of f */
    int i,j,np,newn,logn,degree;
    mr_utype p;
    big *F;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_FP_ROUNDING
    mr_large ip;
#endif
    degree=2*n;
    newn=1; logn=0;
    while (degree>newn) { newn<<=1; logn++; }

    if (mr_mip->degree!=0) 
    {
        for (i=0;i<mr_mip->nprimes;i++)
        {
            mr_free(mr_mip->s1[i]);
            mr_free(mr_mip->s2[i]);
        }
        mr_free(mr_mip->s1);
        mr_free(mr_mip->s2);
    }

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

    mr_mip->degree=n;
    mr_mip->s1=(mr_utype **)mr_alloc(_MIPP_ np,sizeof(mr_utype *));
    mr_mip->s2=(mr_utype **)mr_alloc(_MIPP_ np,sizeof(mr_utype *));
    
    F=(big *)mr_alloc(_MIPP_ n+1,sizeof(big));
    for (i=0;i<=n;i++) 
    {
        F[i]=mirvar(_MIPP_ 0);
        if (f[i]!=NULL) copy(f[i],F[i]);
    }

    modxn_1(_MIPP_ newn/2,n,F);

    for (i=0;i<np;i++)
    {
      /* reserve space for precomputed tables.   */
      /* Note that s2 is half the size of s1     */
        mr_mip->s1[i]=(mr_utype *)mr_alloc(_MIPP_ newn,sizeof(mr_utype));
        mr_mip->s2[i]=(mr_utype *)mr_alloc(_MIPP_ 1+newn/2,sizeof(mr_utype));

        p=mr_mip->prime[i];
#ifdef MR_FP_ROUNDING
        ip=mr_invert(p);
#endif
        for (j=0;j<n;j++)
        {
            if (rf[j]==NULL) mr_mip->s1[i][j]=0;
#ifdef MR_FP_ROUNDING
            else mr_mip->s1[i][j]=mr_sdiv(_MIPP_ rf[j],p,ip,mr_mip->w1);
#else
            else mr_mip->s1[i][j]=mr_sdiv(_MIPP_ rf[j],p,mr_mip->w1);
#endif
        }
        mr_dif_fft(_MIPP_ logn,i,mr_mip->s1[i]);
        
       for (j=0;j<=n;j++)
#ifdef MR_FP_ROUNDING
           mr_mip->s2[i][j]=mr_sdiv(_MIPP_ F[j],p,ip,mr_mip->w1);
#else
           mr_mip->s2[i][j]=mr_sdiv(_MIPP_ F[j],p,mr_mip->w1);
#endif
       mr_dif_fft(_MIPP_ logn-1,i,mr_mip->s2[i]);
    }
    for (i=0;i<=n;i++) mr_free(F[i]);
    mr_free(F);
}

int mr_ps_zzn_mul(_MIPD_ int deg,big *x,big *y,big *z)
{
    int i,j,newn,logn,np;
    mr_utype inv,p,fac;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_FP_ROUNDING
    mr_large ip;
#endif
    newn=1; logn=0;
    while (2*deg>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;

    for (i=0;i<np;i++)
    {
        p=mr_mip->prime[i];
#ifdef MR_FP_ROUNDING
        ip=mr_invert(p);
#endif
        for (j=0;j<deg;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);
#else
            else mr_mip->wa[j]=mr_sdiv(_MIPP_ x[j],p,mr_mip->w1);
#endif
        }
        for (j=deg;j<newn;j++) mr_mip->wa[j]=0;

        mr_dif_fft(_MIPP_ logn,i,mr_mip->wa);  
        for (j=0;j<deg;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);
#else
            else mr_mip->t[i][j]=mr_sdiv(_MIPP_ y[j],p,mr_mip->w1);
#endif
        }  
        for (j=deg;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->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<deg;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);
    for (j=0;j<deg;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);
        divide(_MIPP_ mr_mip->w7,mr_mip->w6,mr_mip->w6);
        redc(_MIPP_ mr_mip->w7,z[j]);
    }
    mr_mip->check=ON;
    return np;
}

int mr_ps_big_mul(_MIPD_ int deg,big *x,big *y,big *z)
{ /* Multiply two power series with large integer parameters */
    int i,j,newn,logn,np;
    mr_utype inv,p,fac;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_FP_ROUNDING
    mr_large ip;
#endif
    newn=1; logn=0;
    while (2*deg>newn) { newn <<=1; logn++; }

    zero(mr_mip->w2);
    zero(mr_mip->w4);

/* find biggest element in each series */
    for (i=0;i<deg;i++)
    {
        if (x[i]!=NULL) 
        {
            absol(x[i],mr_mip->w3);
            if (compare(mr_mip->w3,mr_mip->w2)>0) copy(mr_mip->w3,mr_mip->w2);
        }
        if (y[i]!=NULL)
        {
            absol(y[i],mr_mip->w3);
            if (compare(mr_mip->w3,mr_mip->w4)>0) copy(mr_mip->w3,mr_mip->w4);
        }
    }
    premult(_MIPP_ mr_mip->w4,2,mr_mip->w4);     /* range is +ve and -ve */
                                          /* so extra factor of 2 included */

    np=mr_fft_init(_MIPP_ logn,mr_mip->w4,mr_mip->w2,TRUE); 
    convert(_MIPP_ 1,mr_mip->w3);
/* 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     
        mr_pmul(_MIPP_ mr_mip->w3,p,mr_mip->w3); 
        for (j=0;j<deg;j++)
        {
            if (x[j]==NULL) mr_mip->wa[j]=0;
            else 
            {
                if (size(x[j])>=0)
                {
                    copy(x[j],mr_mip->w1);
#ifdef MR_FP_ROUNDING
                    mr_mip->wa[j]=mr_sdiv(_MIPP_ mr_mip->w1,p,ip,mr_mip->w1);
#else
                    mr_mip->wa[j]=mr_sdiv(_MIPP_ mr_mip->w1,p,mr_mip->w1);
#endif
                }
                else
                {
                    negify(x[j],mr_mip->w1);
#ifdef MR_FP_ROUNDING
                    mr_mip->wa[j]=p-mr_sdiv(_MIPP_ mr_mip->w1,p,ip,mr_mip->w1);
#else
                    mr_mip->wa[j]=p-mr_sdiv(_MIPP_ mr_mip->w1,p,mr_mip->w1);
#endif
                }
            }
        }
        for (j=deg;j<newn;j++) mr_mip->wa[j]=0;

        mr_dif_fft(_MIPP_ logn,i,mr_mip->wa);  
        for (j=0;j<deg;j++)
        {
            if (y[j]==NULL) mr_mip->t[i][j]=0;
            else 
            {
                if (size(y[j])>=0)
                {
                    copy(y[j],mr_mip->w1);
#ifdef MR_FP_ROUNDING
                    mr_mip->t[i][j]=mr_sdiv(_MIPP_ mr_mip->w1,p,ip,mr_mip->w1);
#else
                    mr_mip->t[i][j]=mr_sdiv(_MIPP_ mr_mip->w1,p,mr_mip->w1);
#endif
                }
                else
                {
                    negify(y[j],mr_mip->w1);
#ifdef MR_FP_ROUNDING
                    mr_mip->t[i][j]=p-mr_sdiv(_MIPP_ mr_mip->w1,p,ip,mr_mip->w1);
#else
                    mr_mip->t[i][j]=p-mr_sdiv(_MIPP_ mr_mip->w1,p,mr_mip->w1);
#endif
                }
            }
        }  
        for (j=deg;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->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<deg;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
    }
  /* w3 is product of chinese primes */
    decr(_MIPP_ mr_mip->w3,1,mr_mip->w4);
    subdiv(_MIPP_ mr_mip->w4,2,mr_mip->w4); /* find mid-point of range */
  
    for (j=0;j<deg;j++)
    {
        for (i=0;i<np;i++)