mrpower.c

大数运算库

            subdiv(_MIPP_ mr_mip->w6,2,mr_mip->w6);
        }
    }
    MR_OUT
}

void powmod2(_MIPD_ big x,big y,big a,big b,big n,big w)
{/* w=x^y.a^b mod n */
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(79)

    prepare_monty(_MIPP_ n);
    nres(_MIPP_ x,mr_mip->w2);
    nres(_MIPP_ a,mr_mip->w4);
    nres_powmod2(_MIPP_ mr_mip->w2,y,mr_mip->w4,b,w);
    redc(_MIPP_ w,w);

    MR_OUT   
}

void powmod(_MIPD_ big x,big y,big n,big w)
{ /* fast powmod, using Montgomery's method internally */

    mr_small norm;
    BOOL mty;
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(18)

    mty=TRUE;
    if (mr_mip->base!=mr_mip->base2)
    {
        if (size(n)<2 || sgcd(n->w[0],mr_mip->base)!=1) mty=FALSE;
    }
    else
        if (subdivisible(_MIPP_ n,2)) mty=FALSE;


    if (!mty)
    { /* can't use Montgomery */
        copy(y,mr_mip->w1);
        copy(x,mr_mip->w3);
        zero(w);
        if (size(mr_mip->w3)==0) 
        {
            MR_OUT
            return;
        }
        convert(_MIPP_ 1,w);
        if (size(mr_mip->w1)==0) 
        {
            MR_OUT
            return;
        }
        if (size(mr_mip->w1)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);
        if (w==n)           mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS) ;
        if (mr_mip->ERNUM)
        {
            MR_OUT
            return;
        }

        norm=normalise(_MIPP_ n,n);
        divide(_MIPP_ mr_mip->w3,n,n);
        forever
        { 
            if (mr_mip->user!=NULL) (*mr_mip->user)();

            if (subdiv(_MIPP_ mr_mip->w1,2,mr_mip->w1)!=0)
                mad(_MIPP_ w,mr_mip->w3,mr_mip->w3,n,n,w);
            if (mr_mip->ERNUM || size(mr_mip->w1)==0) break;
            mad(_MIPP_ mr_mip->w3,mr_mip->w3,mr_mip->w3,n,n,mr_mip->w3);
        }
        if (norm!=1)
        {
#ifdef MR_FP_ROUNDING
            mr_sdiv(_MIPP_ n,norm,mr_invert(norm),n);
#else
            mr_sdiv(_MIPP_ n,norm,n);
#endif
            divide(_MIPP_ w,n,n);
        }
    }
    else
    { /* optimized code for odd moduli */
        prepare_monty(_MIPP_ n); 
        nres(_MIPP_ x,mr_mip->w3);
        nres_powmod(_MIPP_ mr_mip->w3,y,w);
        redc(_MIPP_ w,w);
    }
    
    MR_OUT
}

int powltr(_MIPD_ int x, big y, big n, big w)
{
    mr_small norm;
    BOOL clean_up,mty;
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return 0;

    MR_IN(71)
    mty=TRUE;
    if (mr_mip->base!=mr_mip->base2)
    {
        if (size(n)<2 || sgcd(n->w[0],mr_mip->base)!=1) mty=FALSE;
    }
    else
    {
        if (subdivisible(_MIPP_ n,2)) mty=FALSE;
    }

/* Is Montgomery modulus in use... */

    clean_up=FALSE;
    if (mty)
    { /* still a chance to use monty... */
       if (mr_mip->modulus!=NULL)
       { /* somebody else is using it */
           if (compare(n,mr_mip->modulus)!=0) mty=FALSE;
       }
       else
       { /* its unused, so use it, but clean up after */
           clean_up=TRUE;
       }
    }
    if (!mty)
    { /* can't use Montgomery! */
        copy(y,mr_mip->w1);
        zero(w);
        if (x==0) 
        {
            MR_OUT
            return 0;
        }
        convert(_MIPP_ 1,w);
        if (size(mr_mip->w1)==0)
        {
            MR_OUT
            return 1;
        }
        
        if (size(mr_mip->w1)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);
        if (w==n)               mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS) ;
        if (mr_mip->ERNUM)
        {
            MR_OUT
            return 0;
        }
        
        norm=normalise(_MIPP_ n,n);

        expint(_MIPP_ 2,logb2(_MIPP_ mr_mip->w1)-1,mr_mip->w2);
        while (size(mr_mip->w2)!=0)
        { /* Left to Right binary method */

            if (mr_mip->user!=NULL) (*mr_mip->user)();
            mad(_MIPP_ w,w,w,n,n,w);
            if (compare(mr_mip->w1,mr_mip->w2)>=0)
            {
                premult(_MIPP_ w,x,w);
                divide(_MIPP_ w,n,n);
                subtract(_MIPP_ mr_mip->w1,mr_mip->w2,mr_mip->w1);
            }
            subdiv(_MIPP_ mr_mip->w2,2,mr_mip->w2);
        }
        if (norm!=1) 
        {
#ifdef MR_FP_ROUNDING
            mr_sdiv(_MIPP_ n,norm,mr_invert(norm),n);
#else
            mr_sdiv(_MIPP_ n,norm,n);
#endif
            divide(_MIPP_ w,n,n);
        }
    }
    else
    { /* optimized code for odd moduli */
        prepare_monty(_MIPP_ n);
        nres_powltr(_MIPP_ x,y,w);
        redc(_MIPP_ w,w);
        if (clean_up) kill_monty(_MIPPO_ );
    }
    MR_OUT 
    return (size(w));
}

BOOL nres_sqroot(_MIPD_ big x,big w)
{ /* w=sqrt(x) mod p. This depends on p being prime! */
    int i,n,e,r,cat;
    BOOL pp;
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    copy(x,mr_mip->w15);
    zero(w);
    if (size(mr_mip->w15)==0)
    {
        zero(w);
        return TRUE;
    }
    if (mr_mip->ERNUM) return FALSE;

    MR_IN(100)

    convert(_MIPP_ 1,w);
    nres(_MIPP_ w,w);
    if (compare(w,mr_mip->w15)==0)
    {
        MR_OUT
        return TRUE;
    }

    cat=remain(_MIPP_ mr_mip->modulus,8);
    switch(cat)
    {
    case 0: case 2: case 4: case 6: 
        zero(w);
        MR_OUT
        return FALSE;

    case 3: case 7:                          /* easy case */
        incr(_MIPP_ mr_mip->modulus,1,mr_mip->w14);
        subdiv(_MIPP_ mr_mip->w14,4,mr_mip->w14);
        nres_powmod(_MIPP_ x,mr_mip->w14,w);
        nres_modmult(_MIPP_ w,w,mr_mip->w14);
        MR_OUT
        if (compare(mr_mip->w14,mr_mip->w15)==0) 
            return TRUE;
        zero(w);
        return FALSE;

    case 5:
        nres_modadd(_MIPP_ mr_mip->w15,mr_mip->w15,mr_mip->w15); /* 2x */
        decr(_MIPP_ mr_mip->modulus,5,mr_mip->w14);
        subdiv(_MIPP_ mr_mip->w14,8,mr_mip->w14);
        nres_powmod(_MIPP_ mr_mip->w15,mr_mip->w14,w);
        nres_modmult(_MIPP_ w,w,mr_mip->w14);
        nres_modmult(_MIPP_ mr_mip->w15,mr_mip->w14,mr_mip->w14);
        convert(_MIPP_ 1,mr_mip->w1);
        nres(_MIPP_ mr_mip->w1,mr_mip->w1);
        nres_modsub(_MIPP_ mr_mip->w14,mr_mip->w1,mr_mip->w14);
        nres_modmult(_MIPP_ mr_mip->w14,w,w);
        if (!subdivisible(_MIPP_ mr_mip->w15,2))
            add(_MIPP_ mr_mip->w15,mr_mip->modulus,mr_mip->w15);
        subdiv(_MIPP_ mr_mip->w15,2,mr_mip->w15);                /* x */
        nres_modmult(_MIPP_ w,mr_mip->w15,w);
        nres_modmult(_MIPP_ w,w,mr_mip->w14);
        MR_OUT
        if (compare(mr_mip->w14,mr_mip->w15)==0) 
            return TRUE;
        zero(w);
        return FALSE;
        
    case 1:                      /* difficult case. Shank's method */

        decr(_MIPP_ mr_mip->modulus,1,mr_mip->w14);
        e=0;
        while (subdivisible(_MIPP_ mr_mip->w14,2))
        {
            subdiv(_MIPP_ mr_mip->w14,2,mr_mip->w14);
            e++;
        }
        for (r=2;;r++)
        {
            convert(_MIPP_ 1,mr_mip->w3);
            nres(_MIPP_ mr_mip->w3,mr_mip->w3);             /* w3=1 */
            nres_powltr(_MIPP_ r,mr_mip->w14,w);
            if (compare(w,mr_mip->w3)==0) continue;
            copy(w,mr_mip->w4);
            nres_negate(_MIPP_ mr_mip->w3,mr_mip->w1);      /* w1 = -1 */
            pp=FALSE;
            for (i=1;i<e;i++)         
            { /* check for composite modulus */
                if (mr_mip->user!=NULL) (*mr_mip->user)();

                if (compare(mr_mip->w4,mr_mip->w1)==0) pp=TRUE;
                nres_modmult(_MIPP_ mr_mip->w4,mr_mip->w4,mr_mip->w4);
                if (!pp && compare(mr_mip->w4,mr_mip->w3)==0)
                {              
                    zero(w);
                    MR_OUT
                    return FALSE;
                }
            }   
            if (compare(mr_mip->w4,mr_mip->w1)==0) break;
            if (!pp)
            {
                zero(w);
                MR_OUT
                return FALSE;
            }
        }                                                /* w= y    */
        copy(mr_mip->w15,mr_mip->w3);                    /* w3 = x  */
        nres_powmod(_MIPP_ mr_mip->w3,mr_mip->w14,mr_mip->w15);  /* w15 = w */
        incr(_MIPP_ mr_mip->w14,1,mr_mip->w14);
        subdiv(_MIPP_ mr_mip->w14,2,mr_mip->w14);
        nres_powmod(_MIPP_ mr_mip->w3,mr_mip->w14,mr_mip->w14);  /* w14 = v */
        forever
        {
            if (mr_mip->user!=NULL) (*mr_mip->user)();

            convert(_MIPP_ 1,mr_mip->w1);
            nres(_MIPP_ mr_mip->w1,mr_mip->w1);
            if (compare(w,mr_mip->w1)==0)  break;
            copy(mr_mip->w15,mr_mip->w2);
            for (n=0;compare(mr_mip->w2,mr_mip->w1)!=0;n++)
                nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w2,mr_mip->w2);
            if (n>=e)
            {
                zero(w);
                MR_OUT
                return FALSE;
            }
            r=e-n-1;
            for (i=0;i<r;i++) nres_modmult(_MIPP_ w,w,w);
            nres_modmult(_MIPP_ mr_mip->w14,w,mr_mip->w14);
            nres_modmult(_MIPP_ w,w,w);
            nres_modmult(_MIPP_ mr_mip->w15,w,mr_mip->w15);
            e=n;
        }
        copy(mr_mip->w14,w);
        nres_modmult(_MIPP_ w,w,mr_mip->w14);
        MR_OUT
        if (compare(mr_mip->w14,mr_mip->w3)==0) 
            return TRUE;
        zero(w);
        return FALSE;
        
    }
    return FALSE;
}

BOOL sqroot(_MIPD_ big x,big p,big w)
{ /* w = sqrt(x) mod p */
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return FALSE;

    MR_IN(101)

    if (subdivisible(_MIPP_ p,2))
    { /* p must be prime */
        zero(w);
        MR_OUT
        return FALSE;
    }
    prepare_monty(_MIPP_ p);
    nres(_MIPP_ x,mr_mip->w15);
    if (nres_sqroot(_MIPP_ mr_mip->w15,w))
    {
        redc(_MIPP_ w,w);
        MR_OUT
        return TRUE;
    }
    zero(w);
    MR_OUT
    return FALSE;
}