zmodn.cpp

关于大数运算的加 减 乘 除 模 等运算 运用了libtommath库的算法

/*
    Multiple-precision modular arithmetic
    BY CSK(陈士凯)
    [e-mail]
    www.csksoft.net
*/

#include "CiperLib.h"
#include "inner_support.h"

//ZModn::genPrimeArr will generate all primes <= MAX_PRIME_NUM to quickly test
//whether an Integer is a prime, implemented by ZModn::quickPrimeDec
#define MAX_PRIME_NUM  1000


//ZModn::isPrime() func will trying R_M_TRYING_TIME times to safely make a decision 
#define R_M_TRYING_TIME 100

bool ZModn::is_prime_arr_inited = false;
unsigned int ZModn::prime_array[MAX_PRIME_NUM];




void ZModn::genPrimeArr(){
    char *tmp_buff;
    int pix,i,k,prime;
    
    int max_prime;
    max_prime=(MAX_PRIME_NUM+1)/2;

    tmp_buff=new char[max_prime];

    memset(tmp_buff,1,max_prime);

    pix=1;

    for (i=0;i<max_prime;i++)
    {
        if (tmp_buff[i])
        {
            prime=i+i+3;
            for (k=i+prime;k<max_prime;k+=prime)
                tmp_buff[k]=0;
            pix++;
        }
    }

   
    prime_array[0]=2;
    pix=1;
    for (i=0;i<max_prime;i++)
        if (tmp_buff[i]) prime_array[pix++]=i+i+3;
    prime_array[pix]=0;
    delete [] tmp_buff;
    is_prime_arr_inited = true;
}

int ZModn::quickPrimeDec( Integer& src)
{
    //abs(src) <=1 ?
    static Integer tmp;
    if (src.m_acutual_len ==1 && src.m_p_data_arr[0] <=1) return 0; 
   

    if (!is_prime_arr_inited) genPrimeArr();

    for (int i=0;prime_array[i]!=0;i++)
    { /* test for divisible by first few primes */
        
        while (Integer::raw_div_int(src,prime_array[i],tmp)==0)
        { 
          //abs(tmp) == 1
           if (tmp.m_acutual_len==1 && tmp.m_p_data_arr[0]==1) return 1;
           else return 0;
        }
        //abs(tmp) <= prime_array[i]
        if (tmp.m_acutual_len==1 && tmp.m_p_data_arr[0]<=prime_array[i])
        { 
            
            return 1;
        }
    }
   
    return 2;

}

bool ZModn::setMod( Integer& src)  
{
    if (isPrime(src))
    {
        ZModn::m_mod_value = src;
        return true;
    }
    return false;
}

bool  ZModn::genPrime(Integer & start_val, Integer &  dest)
{
    static Integer tmp;
    dest = start_val;
    if (dest.m_acutual_len==1 && dest.m_p_data_arr[0]<2) 
    {
        dest = 2;
        return true;
    }
    if (Integer::raw_div_int(dest,2,tmp)==0) dest+=1;
    else                           dest+=2;
    while (!isPrime(dest)) 
    {
        dest+=2;
       
    }

    return true;
}



bool ZModn::isPrime( Integer& src)   //judge whether src is a prime
{
    int j,k,n,r,times,d;

    if (src.m_acutual_len==1 && src.m_p_data_arr[0]<=1) return false;    //abs(src) <= 1

    k=quickPrimeDec(src);
    if (k==0) 
    {
   
        return false;
    }
    if (k==1)
    {

        return true;
    }
/*
    w1 -> tmp_k
    w8 -> tmp2
    w2 -> tmp3
    w9 -> tmp4
*/
/* Miller-Rabin */
    static Integer tmp_k,tmp2,tmp3,tmp4;
    tmp_k = src - 1;
    
    k=0;
    while (Integer::raw_div_int(tmp_k,2,tmp_k)==0)
    {
        k++;

        tmp2 = tmp_k;
    }        
    times=R_M_TRYING_TIME;
    d=IntHelper::logb2(src);    /* for larger primes, reduce NTRYs */
    if (d>220) times=2+((10*times)/(d-210));
    ZModn tmpZm;
    tmpZm.m_mod_value  = src;
    tmp_k = tmp2;
    tmpZm.toMod(tmp_k);
    for (n=1;n<=times;n++)
    { 
        j=0;
        do
        {
            r=(int)IntHelper::MZ_rand()%MR_TOOBIG;
            if (src.m_acutual_len==1 && src.m_p_data_arr[0]<MR_TOOBIG) r%=src.m_p_data_arr[0];
        } while (r<2);
        
        
        tmp4.fromInt( r);
        tmpZm.toMod(tmp4);
        
        
        tmpZm.mod_pow(tmp4,tmp_k);

        
        tmp3 = src -1;
        while ((j>0 || tmp4.m_acutual_len!=1 || tmp4.m_p_data_arr[0]!=1) 
              && tmp4!=tmp3)
        {
            j++;
            if ((j>1 && (tmp4.m_acutual_len==1 && tmp4.m_p_data_arr[0]==1)) || j==k)
            { 
                return false;
            }

/*
    x = y = z = r = w9 = tmp4
    w = q = src

*/
            tmp4 *= tmp4;
            tmp4 %=src;
         
        }

    }

    return true;

}
Integer & ZModn::toMod(Integer& src)    //src = src mod m_mod_value
{
    if (!m_mod_value.isZero())
    src %= m_mod_value;
   
    return src;

}
//
Integer & ZModn::mod_add(Integer&src, Integer& dest)
{
    dest+=src;
    if (!m_mod_value.isZero())
      if (dest>=m_mod_value) dest-=m_mod_value;
    return dest;
}

Integer & ZModn::mod_mul(Integer&src, Integer& dest)
{

//TODO not the efficient way
//Using Montgomery's REDC after mul would be better
    
    dest *= src;
    if (!m_mod_value.isZero())
      dest %= m_mod_value;
    return dest;
}



Integer & ZModn::mod_pow(Integer&src, Integer& N)
{
// "Analysis of Sliding Window Techniques for Exponentiation",
// C.K. Koc, Computers. Math. & Applic. Vol. 30 pp17-24 1995. 

    int i,j,k,t,nb,nbw,nzs,n;
    Integer *int_tb[16];
    

    
    if (src.isZero() ){
        if (N.isZero()) src.zero();
        return src;
    }

    Integer tmp = N, tmpAns;
    Integer tmp1;
    if (tmp.isZero())
    {
        src = 1;
        return src;
    }

    tmpAns = 1;
    // prepare for further chgs..
    toMod(tmpAns);
    
    static Integer i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12;
    i1 = src;
    int_tb[0]=&i1;
    int_tb[1]=&i2;
    int_tb[2]=&i3;
    int_tb[3]=&i4;
    int_tb[4]=NULL;  
    int_tb[5]=&i5;
    int_tb[6]=&i6;
    int_tb[7]=&i7;
    int_tb[8]=NULL;
    int_tb[9]=NULL;
    int_tb[10]=&i8;
    int_tb[11]=&i9;
    int_tb[12]=NULL;
    int_tb[13]=&i10;
    int_tb[14]=&i11;
    int_tb[15]=&i12;
   

    tmp1 = *int_tb[0];
    mod_mul(tmp1,tmp1); //x^2

    n=15;
    j=0;
    do
    { // pre-calc
        t=1; k=j+1;
        while (int_tb[k]==NULL) {k++; t++;}
        *int_tb[k] = *int_tb[j];
        for (i=0;i<t;i++) mod_mul(tmp1,*int_tb[k]);
        j=k;
    } while (j<n);

    nb=IntHelper::logb2(tmp);
    tmpAns = *int_tb[0];
  
	if (nb>1) for (i=nb-2;i>=0;)
    { /* Left to Right method */

        n=IntHelper::sliding_window(tmp,i,&nbw,&nzs,5); 

        for (j=0;j<nbw;j++) mod_mul(tmpAns,tmpAns); //x^2
        if (n>0) mod_mul(*int_tb[n/2],tmpAns);
            
        i-=nbw;
        if (nzs)
        {
            for (j=0;j<nzs;j++) mod_mul(tmpAns,tmpAns);
            i-=nzs;
        }
     }
 
    src = tmpAns;
    return src;
}
    

Integer & ZModn::Euc(Integer&, Integer&dest)
{

    return dest;
}