int_op.cpp

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

/*
    Big Integer Operation
    ADD/SUB/MUL/DIV/MOD operation
    SRC FILE
    BY CSK(陈士凯)
    [e-mail]
    www.csksoft.net
*/

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

#pragma warning(disable:4731)

inline void Integer::raw_add(Integer &src, const Integer &dest)
{
   unsigned int dest_length, src_length;
   
   if (src.m_acutual_len >= dest.m_acutual_len)
   {
        dest_length = src.m_acutual_len;
        src_length = dest.m_acutual_len;
   }else
   {
        dest_length = dest.m_acutual_len;
        src_length = src.m_acutual_len;
   }

   if (dest_length < Integer::ms_u_data_arr_num) dest_length++;
   src.m_dw_isOverflow = 0; //clear overflow flag

   unsigned int pos = 0;
//ADD CODE-- C-Style--slow..
#ifdef _NOTUSEASM
    unsigned int carry_flag = 0;
    
 //   unsigned int result;
    for(;pos < src_length; pos++)   //add all columns of src and dest 
    {
        unsigned int src_value;
        src_value = src.m_p_data_arr[pos];
        src.m_p_data_arr[pos] += dest.m_p_data_arr[pos] + carry_flag;
        if (src.m_p_data_arr[pos] > src_value) { //is overflow?
            carry_flag = 0;
        }
        else if (src.m_p_data_arr[pos] < dest.m_p_data_arr[pos]){
            carry_flag = 1;
        }
       
    }
    for(;pos < dest_length && carry_flag; pos++)  //continue to add, caused by the sideeffect of carry_flag
    {
        unsigned int src_value;
        src_value = src.m_p_data_arr[pos];
        src.m_p_data_arr[pos] += carry_flag;
        if (src.m_p_data_arr[pos] > src_value) { //is overflow?
            carry_flag = 0;
        }
        else if (src.m_p_data_arr[pos] < src_value){
            carry_flag = 1;
        }
    }
    //ok, the algorithm should stop here, otherwise overflow happens..

    if (carry_flag)//overflow
    {
        src.m_dw_isOverflow = 1;
    }
#endif

#ifndef _NOTUSEASM
//ADD CODE-- ASM--hard to read, but faster
    unsigned int * src_ptr,* dest_ptr, * is_overflow_ptr;
    src_ptr = src.m_p_data_arr;
    dest_ptr = const_cast<unsigned int *>(dest.m_p_data_arr);
    is_overflow_ptr = &src.m_dw_isOverflow;
    unsigned int remain_offset = dest_length - src_length;
    __asm
    {
      //  PUSHAD                      //save current register
        
        mov edi, src_ptr
        mov esi, dest_ptr           //set base addr
        xor ebx, ebx                //set counter, and clear carry flag
        mov ecx, src_length

firstl: 
        mov eax, [esi + ebx]        //load dest data
        adc [edi+ebx], eax          //Add with carry flag
        inc ebx
        inc ebx
        inc ebx
        inc ebx                     //counter ++
        dec ecx
        jnz firstl                  //loop if counter<src_length
        
        mov eax,ecx                 //set eax = 0
        mov ecx, remain_offset
        test ecx,ecx
        jz  final                   //counter < dest_length?
secondl:
        jnc final                   //carry != 0 ?
        adc [edi+ebx], eax          //src.m_p_data_arr[pos] += carry_flag;
                                    //note! eax == 0
        inc ebx
        inc ebx
        inc ebx
        inc ebx                     //counter ++
        dec ecx
        jnz secondl                 //loop if counter<dest_length

final:  
        mov eax,is_overflow_ptr
        adc [eax],0                 //is overflow?
   //     POPAD                       //recover current register
        
    }
   
#endif

    //recalc the length of src
   
    while(src.m_p_data_arr[dest_length-1] == 0 && dest_length>1) dest_length--;
    src.m_acutual_len = dest_length;
}

inline void Integer::raw_sub(Integer &src, const Integer &dest)
{
    //assume src>dest and both are positive

    unsigned int pos =0;
    unsigned int dest_length = dest.m_acutual_len;
   
#ifdef _NOTUSEASM
    unsigned int borrow_flag = 0;
    for (; pos < dest_length || borrow_flag;pos++)
    {
        unsigned int src_value;
        src_value = src.m_p_data_arr[pos];
        src.m_p_data_arr[pos] -= dest.m_p_data_arr[pos] - borrow_flag;
        if (src.m_p_data_arr[pos]<src_value){
            borrow_flag=0;
        }else if (src.m_p_data_arr[pos]>src_value){
            borrow_flag=1;
        }
    }
#endif

#ifndef _NOTUSEASM
    unsigned int * src_ptr,* dest_ptr;

    src_ptr = src.m_p_data_arr;
    dest_ptr = const_cast<unsigned int *>(dest.m_p_data_arr);
    
    __asm
    {
     //   PUSHAD                      //save current register
        
        mov edi, src_ptr
        mov esi, dest_ptr           //set base addr
        xor ebx, ebx                //set counter, and clear carry flag
        mov ecx, dest_length

firstl: 
        mov eax, [esi + ebx]        //load dest data
        sbb [edi+ebx], eax          //sub with carry flag
        inc ebx
        inc ebx
        inc ebx
        inc ebx                     //counter ++
        dec ecx
        jnz firstl                  //pos < dest_length
   
secondl:         
        jnc final                   //borrow_flag == 1?
        mov eax, [esi + ebx]        //load dest data
        sbb [edi+ebx], eax          //sub with carry flag
        inc ebx
        inc ebx
        inc ebx
        inc ebx                     //counter ++
        jmp secondl                         
final:  
     //   POPAD                       //recover current register
    }
#endif
    pos = src.m_acutual_len-1;
    while (src.m_p_data_arr[pos]==0 && pos) pos--;
    src.m_acutual_len = pos + 1;

    src.m_dw_isOverflow = 0; //never overflows
}

inline void Integer::add_sub_director(Integer &src,  Integer &dest, unsigned int intend_op)
{
/*
[0] [1] [intend_op]
-----------------------------------
 +   +     +         sig = [0]
 -   -     +
 +   -     -
 -   +     -         (old_sign_src ^ intend_op) ^ old_sign_dest = 0
-----------------------------------
 +   +     -         sig = [0] * compare
 +   -     +
 -   -     -
 -   +     +         (old_sign_src ^ intend_op) ^ old_sign_dest = 1
-----------------------------------

*/


    unsigned int old_sign_src,old_sign_dest;
    old_sign_src = src.m_dw_sign;
    old_sign_dest = dest.m_dw_sign;

    src.m_dw_sign = dest.m_dw_sign = 0; //get abs value

    

    
    unsigned int selector = (old_sign_src ^ intend_op) ^ old_sign_dest;

    
    if (selector)
    {
        int compare_val = src.compare(dest);
        if (compare_val ==0)  //must be zero
        {
            src.zero();
            dest.m_dw_sign = old_sign_dest;
            return;
        }
        if (compare_val == 1) //abs(src) > abs(dest)
        {
            raw_sub(src,dest);
            src.m_dw_sign = old_sign_src;
        }else
        {

            //tricks here, actually very bad coding style, but runs faster
            Integer tmp(dest);
            raw_sub(tmp,src);
            release_data_buf(src);
            src.m_acutual_len = tmp.m_acutual_len;
            src.m_p_data_arr = tmp.m_p_data_arr;
            tmp.m_p_data_arr = NULL;
            src.m_dw_sign = 1- old_sign_src;
        }
    }
    else    
    {
        raw_add(src,dest);
        src.m_dw_sign =old_sign_src;
    }

    dest.m_dw_sign = old_sign_dest;
    
}

void Integer::raw_mul(const Integer &left, const Integer &right, Integer &dest)
{


    
    if (left.isZero() || right.isZero()) {
        dest.zero();
        return;
    }
    unsigned int dest_sign = left.m_dw_sign ^ right.m_dw_sign;
    unsigned int ans_len = left.m_acutual_len + right.m_acutual_len;
    unsigned int buffer_length = MAX(ans_len,ms_u_data_arr_num<<1);

    static Integer tmpInt;
    static bool fisrtcall = true;
    //tricks
    if (fisrtcall){
        release_data_buf(tmpInt);
        tmpInt.m_p_data_arr = new unsigned int [buffer_length];
        fisrtcall = false;
    }
    unsigned int *ans_data = tmpInt.m_p_data_arr;
    //
    memset(ans_data,0,sizeof(unsigned int) * (buffer_length));

    unsigned int *left_buffer = const_cast<unsigned int *>(left.m_p_data_arr);
    unsigned int *right_buffer = const_cast<unsigned int *>(right.m_p_data_arr);   

#ifdef _NOTUSEASM
    __int64 tmp_ans;
    unsigned int *dw_ans = (unsigned int *)&tmp_ans;
    unsigned int i = 0 , j;
    unsigned int carry_val;
    
    for (;i<left.m_acutual_len;i++) // long multiplication
    { 
        carry_val = 0;

        for (j=0;j<right.m_acutual_len;j++)
        { 

            tmp_ans=(__int64)left_buffer[i]*right_buffer[j]+carry_val+ans_data[i+j];
            ans_data[i+j]=dw_ans[0]; // (a*b + c) mod base
            carry_val=dw_ans[1];     // (a*b + c) / base
        }
        ans_data[right.m_acutual_len+i]=carry_val;
    }
#endif

#ifndef _NOTUSEASM
//from MIRACL lib
/*
 *   MIRACL arithmetic routines 2 - multiplying and dividing BIG NUMBERS.
 *   mrarth2.c
 *
 *   Copyright (c) 1988-2003 Shamus Software Ltd.
 */
     unsigned int src_len_left = left.m_acutual_len;
     unsigned int src_len_right = right.m_acutual_len;
     for (unsigned int i=0;i<src_len_left;i++)
     { 
         __asm{
             mov ecx,src_len_right
             mov esi,i
             shl esi,2
             mov ebx,left_buffer

             add ebx,esi                //i++

             mov edi,[ebx]              //left_buffer[i] -> edi

             mov ebx,ans_data           //ans_data -> edi
             add ebx,esi                //ans_data+=i

             mov esi,right_buffer       //right_buffer -> esi
             sub ebx,esi                
             sub ebx,4                  //right_buffer[]
             push ebp
             xor ebp,ebp
       loop1:
             mov eax,[esi]              //right_buffer[j] -> eax
             add esi,4                  //j++
             mul edi                    // right_buffer[j] * left_buffer[i] -> edx:eax
             add eax,ebp                // right_buffer[j] * left_buffer[i] mod base + carry -> cf
             mov ebp,[esi+ebx]          // ans_data[j+i] -> ebp
             adc edx,0                  // right_buffer[j] * left_buffer[i].high + cf
             add ebp,eax                // right_buffer[j] * left_buffer[i] mod base + ans_data[j+i] -> cf
             adc edx,0                  // right_buffer[j] * left_buffer[i].high + cf
             mov [esi+ebx],ebp          // (a.b + c) mod base -> ans_data[j+i]
             dec ecx                    // src_len_right--
             mov ebp,edx                // carry = (a.b + c) / base
             jnz loop1

             mov [esi+ebx+4],ebp        // ans_data[right.m_acutual_len+i]=carry_val;
             pop ebp
         }
     }
#endif

    //save ans to dest
    memcpy(dest.m_p_data_arr,ans_data,sizeof(unsigned int) *ms_u_data_arr_num);
    dest.m_dw_sign = dest_sign;
    
    while(ans_data[ans_len-1] ==0 && ans_len>1) ans_len--;
    //overflow checking and calc actual len
    if (dest.m_dw_isOverflow = (ans_len>ms_u_data_arr_num)){
        dest.m_acutual_len = ms_u_data_arr_num;
      
    }
    else
    {
        dest.m_acutual_len = ans_len;
    }
    
    ///////////



}

//assume &left != &right
int Integer::raw_div(Integer &left, Integer &right, Integer &dest)
{

    if (right.isZero()) return 1; //divided by zero
    
    static Integer tmpAns;
    
    
    unsigned int d_factor = 0;
    unsigned int dw_reminder;
    unsigned int *left_buf = left.m_p_data_arr;
    unsigned int *right_buf = right.m_p_data_arr;
    unsigned int *dest_buf = dest.m_p_data_arr;
    unsigned int *tmp_buf = tmpAns.m_p_data_arr;
    unsigned int right_buff_len = right.m_acutual_len;
    unsigned int borrow_flag;
    unsigned int old_sig_left,old_sig_right;
    unsigned int attemp;
    old_sig_left = left.m_dw_sign;
    old_sig_right = right.m_dw_sign;
    unsigned int dest_sign = old_sig_left ^ old_sig_right;
    //calc with positive
    left.m_dw_sign = right.m_dw_sign = 0;
    tmpAns = left; //copy left to tmpAns

//some special conditions:
    if (left.m_acutual_len == right.m_acutual_len)
    {
        if (left.m_acutual_len == 1)    //A and B have only one "digit"
        {
            d_factor = tmp_buf[0] / right_buf[0];
            tmp_buf[0] %=  right_buf[0];
        }
        else if ((tmp_buf[tmpAns.m_acutual_len-1]/4)<right_buf[tmpAns.m_acutual_len-1])
        {
            while( tmpAns >= right )
            {
                raw_sub(tmpAns,right);
                d_factor ++;
            }


        }
    }
    if (tmpAns < right) //since tmpAns<right, tmpAns is reminder
    {
        if (left_buf != dest_buf) //calc reminder
        {
            
            left = tmpAns;
            if(!left.isZero()) left.m_dw_sign = old_sig_left;

        }
        if (right_buf != dest_buf) //calc quotient
        {
            dest.zero();
            dest_buf[0]=d_factor;
            if (d_factor) dest.m_dw_sign = dest_sign;
        }
        right.m_dw_sign = old_sig_right;
        return 0;
    }