mp.c

来自「开放源码的编译器open watcom 1.6.0版的源代码」· C语言 代码 · 共 1,082 行 · 第 1/2 页

/****************************************************************************
*
*                            Open Watcom Project
*
*    Portions Copyright (c) 1983-2002 Sybase, Inc. All Rights Reserved.
*
*  ========================================================================
*
*    This file contains Original Code and/or Modifications of Original
*    Code as defined in and that are subject to the Sybase Open Watcom
*    Public License version 1.0 (the 'License'). You may not use this file
*    except in compliance with the License. BY USING THIS FILE YOU AGREE TO
*    ALL TERMS AND CONDITIONS OF THE LICENSE. A copy of the License is
*    provided with the Original Code and Modifications, and is also
*    available at www.sybase.com/developer/opensource.
*
*    The Original Code and all software distributed under the License are
*    distributed on an 'AS IS' basis, WITHOUT WARRANTY OF ANY KIND, EITHER
*    EXPRESS OR IMPLIED, AND SYBASE AND ALL CONTRIBUTORS HEREBY DISCLAIM
*    ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF
*    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT OR
*    NON-INFRINGEMENT. Please see the License for the specific language
*    governing rights and limitations under the License.
*
*  ========================================================================
*
* Description:  WHEN YOU FIGURE OUT WHAT THIS FILE DOES, PLEASE
*               DESCRIBE IT HERE!
*
****************************************************************************/


/* This module contains functions for doing multiple precision, positive
 * integer arithmetic.  Allocations for these numbers increase automatically
 * as required when doing the arithmetic.  The size of a number can grow until
 * there is no more heap space.
 */
#ifndef TEST
#include header
#else
#include <malloc.h>
#define _MemoryAllocate( size )         malloc( size )
#define _MemoryFree( ptr )              free( ptr )
#endif
#include "mp.h"
#include <limits.h>
#include <string.h>
#include "fp.h"

#define two_to_the_31           0x80000000
#define two_to_the_32           (uint64)0x100000000

#define MP_GROW( num, size ) \
        if( mp_grow( num, size ) == MP_OUT_OF_MEMORY ) { \
            return MP_OUT_OF_MEMORY; \
        }
#define MP_COPY( dst, src ) \
        if( mp_copy( dst, src ) == MP_OUT_OF_MEMORY ) { \
            return MP_OUT_OF_MEMORY; \
        }
#define MP_INIT( num, val ) \
        if( mp_init( num, val ) == MP_OUT_OF_MEMORY ) { \
            return MP_OUT_OF_MEMORY; \
        }

/* initialize an mpnum */
int mp_init( mpnum *mp, uint64 value )
{
    mp->num = (uint32*)_MemoryAllocate( sizeof(uint32) * 2 );
    if( mp->num == NULL ) return MP_OUT_OF_MEMORY;
    memset( mp->num, 0, sizeof(uint32) * 2 );
    mp->allocated = 2;
    if( value >= two_to_the_32 ) {
        mp->num[0] = (uint32)value;
        mp->num[1] = value>>32;
        mp->len = 2;
    } else {
        mp->num[0] = value;
        mp->len = 1;
    }
    return MP_NO_ERR;
}

/* extend the number of uint32s allocated to src by size */
int mp_grow( mpnum *src, uint32 size )
{
    uint32 *temp;
    size = max( 3, size );
    temp = (uint32*)_MemoryAllocate( sizeof(uint32) * (src->allocated + size) );
    if( temp == NULL ) return MP_OUT_OF_MEMORY;
    memcpy( temp, src->num, src->allocated * sizeof(uint32) );
    memset( &temp[src->allocated], 0, size * sizeof(uint32) );
    _MemoryFree( src->num );
    src->num = temp;
    src->allocated += size;
    return MP_NO_ERR;
}

/* free the array of uint32s */
int mp_free( mpnum *num )
{
    if( num->allocated > 0 ) {
        _MemoryFree( num->num );
    }
    num->num = NULL;
    num->allocated = 0;
    num->len = 0;
    return MP_NO_ERR;
}

/* set num to zero */
int mp_zero( mpnum *num )
{
    memset( num->num, 0, num->allocated * sizeof(uint32) );
    num->len = 0;
    return MP_NO_ERR;
}

/* copy the src number to the destination number; dst must have already been
 * initialized */
int mp_copy( mpnum *dst, mpnum *src )
{
    if( dst->allocated < src->len ) {
        MP_GROW( dst, src->len - dst->allocated );
    }
    memcpy( dst->num, src->num, src->len * sizeof(uint32) );
    dst->len = src->len;
    return MP_NO_ERR;
}

/* compare num1 and num2
 * returns:  -1 if num1 > num2
 *            0 if num1 = num2
 *            1 if num1 < num2
 */
int mp_cmp( mpnum *num1, mpnum *num2 )
{
    uint32 i;
    uint32 minlen = min( num1->len, num2->len );

    if( num1->len > num2->len ) {
        for( i = num1->len - 1; i >= minlen; i-- ) {
            if( num1->num[i] > 0 ) return -1;
            if( i == 0 ) break;
        }
    } else if( num1->len < num2->len ) {
        for( i = num2->len - 1; i >= minlen; i-- ) {
            if( num2->num[i] > 0 ) return 1;
            if( i == 0 ) break;
        }
    }
    if( minlen == 0 ) return 0;
    for( i = minlen - 1; ; i-- ) {
        if( num1->num[i] > num2->num[i] ) {
            return -1;
        } else if( num1->num[i] < num2->num[i] ) {
            return 1;
        }
        if( i == 0 ) break;
    }
    return 0;
}

int mp_gt( mpnum *num1, mpnum *num2 )
{
    if( mp_cmp( num1, num2 ) < 0 ) return TRUE;
    return FALSE;
}

int mp_gte( mpnum *num1, mpnum *num2 )
{
    if( mp_cmp( num1, num2 ) <= 0 ) return TRUE;
    return FALSE;
}

int mp_lt( mpnum *num1, mpnum *num2 )
{
    if( mp_cmp( num1, num2 ) > 0 ) return TRUE;
    return FALSE;
}

int mp_lte( mpnum *num1, mpnum *num2 )
{
    if( mp_cmp( num1, num2 ) >= 0 ) return TRUE;
    return FALSE;
}

int mp_eq( mpnum *num1, mpnum *num2 )
{
    if( mp_cmp( num1, num2 ) == 0 ) return TRUE;
    return FALSE;
}

int mp_ne( mpnum *num1, mpnum *num2 )
{
    if( mp_cmp( num1, num2 ) == 0 ) return FALSE;
    return TRUE;
}

/* shift src 'bits' number of bits to the left */
int mp_shiftleft( mpnum *dst, mpnum *src, uint32 bits )
{
    uint32 newsize = src->len + ((bits-1) / 32) + 1;
    uint32 cursrc = src->len - 1;
    uint32 curdst = newsize - 1;
    uint32 temp1, temp2;
    uint8  shift = ((bits-1) % 32) + 1;
    uint32 i;

    if( bits == 0 || src->len == 0 ) {
        if( dst != src ) MP_COPY( dst, src );
        return MP_NO_ERR;
    }

    if( dst->allocated < newsize ) {
        MP_GROW( dst, newsize - dst->allocated );
    }
    for( i = dst->allocated-1; i >= newsize - 1; i-- ) {
        dst->num[i] = 0;
    }
    for(;;cursrc--,curdst--) {
        temp1 = src->num[cursrc] >> (32-shift);
        if( shift == 32 ) {
            temp2 = 0;
        } else {
            temp2 = src->num[cursrc] << shift;
        }
        dst->num[curdst] |= temp1;
        dst->num[curdst-1] = temp2;
        if( cursrc == 0 ) break;
    }
    for( curdst--;; curdst-- ) {
        if( curdst == 0 ) break;
        dst->num[curdst-1] = 0;
    }
    dst->len = newsize;
    if( dst->num[newsize-1] == 0 ) {
        dst->len -= 1;
    }
    return MP_NO_ERR;
}

/* shift src 'bits' number of bits to the right */
int mp_shiftright( mpnum *dst, mpnum *src, uint32 bits )
{
    uint32 newsize;
    int temp = src->len - bits / 32;
    uint32 cursrc = bits / 32;
    uint32 curdst = 0;
    uint32 temp1, temp2;
    uint8  shift = bits % 32;

    if( temp <= 0 ) {
        mp_zero( dst );
        return MP_NO_ERR;
    }
    newsize = (uint32)temp;

    if( bits == 0 || src->len == 0 ) {
        if( dst != src ) MP_COPY( dst, src );
        return MP_NO_ERR;
    }

    if( dst->allocated < newsize ) {
        MP_GROW( dst, newsize - dst->allocated );
    }
    for(;;cursrc++,curdst++) {
        if( shift == 0 ) {
            temp1 = 0;
        } else {
            temp1 = src->num[cursrc] << (32-shift);
        }
        temp2 = src->num[cursrc] >> shift;
        if( curdst > 0 ) dst->num[curdst-1] |= temp1;
        dst->num[curdst] = temp2;
        if( cursrc == src->len - 1 ) break;
    }
    for( curdst++;; curdst++ ) {
        if( curdst >= newsize ) break;
        dst->num[curdst] = 0;
    }
    dst->len = newsize;
    if( dst->num[newsize-1] == 0 ) {
        dst->len -= 1;
    }

    return MP_NO_ERR;
}

/* add a scalar value to src */
int mp_addsc( mpnum *dst, mpnum *src, uint32 scalar )
{
    uint32 maxlen = src->len;
    uint64 sum = 0;
    uint64 carry = scalar;
    uint32 i;

    if( dst->allocated < maxlen + 1 ) {
        MP_GROW( dst, maxlen + 1 - dst->allocated );
    }
    for( i = 0; ; i++ ) {
        if( i >= src->len ) {
            if( carry > 0 ) {
                dst->num[i] = carry;
            }
            break;
        } else {
            sum = (uint64)src->num[i] + carry;
        }
        if( sum >= two_to_the_32 ) {
            dst->num[i] = sum - two_to_the_32;
            carry = 1;
        } else {
            dst->num[i] = sum;
            carry = 0;
            break;
        }
    }
    if( dst->num[maxlen] == 0 ) {
        dst->len = maxlen;
    } else {
        dst->len = maxlen + 1;
    }
    return MP_NO_ERR;
}

/* add two mpnums */
int mp_add( mpnum *dst, mpnum *src1, mpnum *src2 )
{
    uint32 maxlen = max( src1->len, src2->len );
    uint64 sum = 0;
    uint64 carry = 0;
    uint32 i;

    if( dst->allocated < maxlen + 1 ) {
        MP_GROW( dst, maxlen + 1 - dst->allocated );
    }
    for( i = 0; ; i++ ) {
        if( i >= src1->len && i >= src2->len ) {
            if( carry > 0 ) {
                dst->num[i] = carry;
                dst->len = maxlen + 1;
            } else {
                dst->len = maxlen;
            }
            break;
        } else if( i >= src1->len ) {
            sum = (uint64)src2->num[i] + carry;
        } else if( i >= src2->len ) {
            sum = (uint64)src1->num[i] + carry;
        } else {
            sum = (uint64)src1->num[i] + (uint64)src2->num[i] + carry;
        }
        if( sum >= two_to_the_32 ) {
            dst->num[i] = sum - two_to_the_32;
            carry = 1;
        } else {
            dst->num[i] = sum;
            carry = 0;
        }
    }
    return MP_NO_ERR;
}

/* subtract two npnums; returns MP_NEGATIVE_RESULT if src1 < src2 */
int mp_sub( mpnum *dst, mpnum *src1, mpnum *src2 )
{
    uint32 maxlen = max( src1->len, src2->len );
    uint32 carry = 0;
    uint32 i;
    uint32 temp;

    if( mp_lt( src1, src2 ) ) {
        /* src1 < src2 => error */
        return MP_NEGATIVE_RESULT;
    }
    /* ASSERT num1 >= num2 => don't need to look for errors */
    if( dst->allocated < maxlen ) {
        MP_GROW( dst, maxlen - dst->allocated );
    }
    for( i = 0; ; i++ ) {
        if( i >= src1->len && i >= src2->len ) {
            break;
        } else if( i >= src2->len ) {
            temp = 0;
        } else {
            temp = src2->num[i];
        }
        if( src1->num[i] >= temp + carry ) {
            dst->num[i] = src1->num[i] - temp - carry;
            carry = 0;
        } else {
            dst->num[i] = (uint64)src1->num[i] + two_to_the_32 - temp - carry;
            carry = 1;
        }
    }
    for( i = maxlen - 1; ; i-- ) {
        if( dst->num[i] > 0 ) {
            dst->len = i + 1;
            break;
        }
        if( i == 0 ) {
            dst->len = 0;
            break;
        }
    }
    return MP_NO_ERR;
}

/* multiply an mpnum by a scalar value */
int mp_mulsc( mpnum *dst, mpnum *src, uint32 scalar )
{
    uint64 product = 0;
    uint64 carry = 0;
    uint32 i;

    if( dst->allocated < src->len + 1 ) {
        MP_GROW( dst, src->len + 1 - dst->allocated );
    }
    for( i = 0; ; i++ ) {
        if( i >= src->len ) {
            if( carry > 0 ) {
                dst->num[i] = carry;
                dst->len = src->len + 1;
            } else {
                dst->len = src->len;
            }
            break;
        } else {
            product = (uint64)src->num[i] * (uint64)scalar + carry;
        }
        if( product >= two_to_the_32 ) {
            dst->num[i] = (uint32)product;
            carry = product>>32;
        } else {
            dst->num[i] = product;
            carry = 0;
        }
    }
    return MP_NO_ERR;
}

/* Multiply two mpnums together using the intuitive pen-and-paper method for
 * multilpication.
 */
int mp_mul( mpnum *dst, mpnum *src1, mpnum *src2 )
{
    uint32 cur;
    uint32 sc;
    mpnum temp;
    mpnum result;
    int rc;

    if( dst->allocated < src1->len + src2->len ) {
        MP_GROW( dst, src1->len + src2->len - dst->allocated );
    }

    MP_INIT( &temp, 0 );
    MP_INIT( &result, 0 );

    for( cur = 0; cur < src2->len; cur++ ) {
        sc = src2->num[cur];
        if( sc != 1 ) {
            rc = mp_mulsc( &temp, src1, sc );
            if( rc != 0 ) goto done;
            rc = mp_shiftleft( &temp, &temp, cur * 32 );
        } else {
            rc = mp_shiftleft( &temp, src1, cur * 32 );
        }
        if( rc != 0 ) goto done;
        rc = mp_add( &result, &result, &temp );
        if( rc != 0 ) goto done;
    }

    MP_COPY( dst, &result );
    rc = MP_NO_ERR;

done:
    mp_free( &temp );
    mp_free( &result );

    return rc;
}

/* divide an mpnum by a scalar value; mode indicates ceiling or floor */
int mp_divsc( mpnum *dst, mpnum *src, uint32 scalar, int mode )
{
    uint64 dividend = 0;
    uint32 remainder = 0;
    uint32 i;

    if( scalar == 0 ) return MP_DIVIDE_BY_ZERO;
    if( scalar == 1 || src->len == 0 ) {
        if( dst != src ) MP_COPY( dst, src );
        return MP_NO_ERR;
    }
    if( dst->allocated < src->len ) {
        MP_GROW( dst, src->len - dst->allocated );
    }
    for( i = src->len - 1; ; i-- ) {
        dividend = (((uint64)remainder) << 32) | (uint64)src->num[i];
        dst->num[i] = (uint32)(dividend / scalar);
        remainder = (uint32)(dividend % scalar);
        if( i == 0 ) break;
    }
    dst->len = src->len;
    if( remainder > 0 && mode == CEILING ) {
        mp_addsc( dst, dst, 1 );
    }
    if( dst->num[src->len-1] == 0 ) dst->len--;
    return MP_NO_ERR;
}

/* This division algorithm is essentially the pen-and-paper long division
 * method.
 */
int mp_div( mpnum *qdst, mpnum *rdst, mpnum *x, mpnum *y )
{
    mpnum temp, temp2, q, r, mptmp1, mptmp2, mptmp3;
    uint32      n = x->len;
    uint32      t = y->len;
    uint32      i;
    int         rc = -1;
    uint64      tmp;
    uint32      shift = 0;

    MP_INIT( &q, 0 );
    if( q.allocated < x->len ) MP_GROW( &q, x->len - q.allocated );
    q.len = x->len;
    MP_INIT( &r, 0 );
    MP_INIT( &temp, 0 );
    MP_INIT( &temp2, 0 );
    i = 0;

    /* copy x into r */
    mp_copy( &r, x );

    /* check division by zero */
    if( mp_eq( &q, y ) ) {
        rc = MP_DIVIDE_BY_ZERO;
        goto done;