fpc.c

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

/****************************************************************************
*
*                            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!
*
****************************************************************************/


/* C implementation of floating point arithmetic, conversion, and comparison
 */
#include <stdio.h>
#include <limits.h>
#include <string.h>
#include "mp.h"
#include "fp.h"

/* c implementation of dadd that may be used in special cases for the non-c
 * version of dadd */
void c_dadd( void *dst, void *src1, void *src2 ) {
    uint64 f1, f2;
    int e1, e2;
    int diff;
    int emax, emin;
    mpnum num1, num2, zero;
    int isNeg1 = d_isNeg( src1 );
    int isNeg2 = d_isNeg( src2 );
    int isNeg = FALSE;
    int newexp;

    /* special cases */
    if( d_isNan( src1 ) || d_isNan( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isPosInf( src1 ) && d_isNegInf( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isNegInf( src1 ) && d_isPosInf( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isPosInf( src1 ) || d_isPosInf( src2 ) ) {
        *(uint64*)dst = DOUBLE_POS_INF;
    } else if( d_isNegInf( src1 ) || d_isNegInf( src2 ) ) {
        *(uint64*)dst = DOUBLE_NEG_INF;
    } else if( d_isPosZero( src1 ) && d_isNegZero( src2 ) ) {
        *(uint64*)dst = DOUBLE_POS_ZERO;
    } else if( d_isNegZero( src1 ) && d_isPosZero( src2 ) ) {
        *(uint64*)dst = DOUBLE_POS_ZERO;
    } else if( d_isPosZero( src1 ) && d_isPosZero( src2 ) ) {
        *(uint64*)dst = DOUBLE_POS_ZERO;
    } else if( d_isNegZero( src1 ) && d_isNegZero( src2 ) ) {
        *(uint64*)dst = DOUBLE_NEG_ZERO;
    } else if( d_isZero( src1 ) ) {
        *(uint64*)dst = *(uint64*)src2;
    } else if( d_isZero( src2 ) ) {
        *(uint64*)dst = *(uint64*)src1;
    } else {
        parseDouble( *(double*)src1, &f1, &e1 );
        parseDouble( *(double*)src2, &f2, &e2 );
        diff = abs( (e1 - e2) );
        emax = max( e1, e2 );
        emin = min( e1, e2 );
        mp_init( &num1, f1 );
        mp_init( &num2, f2 );
        mp_init( &zero, 0 );
        if( e1 > e2 ) {
            mp_shiftleft( &num1, &num1, (uint64)diff );
        } else if( e2 > e1 ) {
            mp_shiftleft( &num2, &num2, (uint64)diff );
        }
        if( (isNeg1 && isNeg2) || (!isNeg1 && !isNeg2 ) ) {
            mp_add( &num1, &num1, &num2 );
            if( isNeg1 ) isNeg = TRUE;
        } else {
            if( mp_gt( &num1, &num2 ) ) {
                mp_sub( &num1, &num1, &num2 );
                if( isNeg1 ) isNeg = TRUE;
            } else {
                mp_sub( &num1, &num2, &num1 );
                if( isNeg2 ) isNeg = TRUE;
            }
        }
        if( mp_eq( &num1, &zero ) ) {
            *(uint64*)dst = DOUBLE_POS_ZERO;
        } else {
            newexp = emax - diff - 53 + mp_bitsize( &num1 ) + double_bias;
            if( newexp <= 0 && newexp > -53 ) {
                /* denormalized result */
                uint64 temp;
                mp_binround( &num1, &num1, -(emin + double_bias) + 1 );
                mp_shiftright( &num1, &num1, -(emin + double_bias) + 1 );
                mp_touint64( &temp, &num1 );
                makeDouble( dst, temp, 0 - double_bias );
            } else {
                mp_todouble( dst, &num1 );
                changeDoubleExponent( dst, dst, emax - diff - 52 );
            }
            if( isNeg ) dneg( dst, dst );
        }
        mp_free( &num1 );
        mp_free( &num2 );
        mp_free( &zero );
    }
}

void c_dsub( void *dst, void *src1, void *src2 )
{
    double nsrc2;
    dneg( &nsrc2, src2 );
    dadd( dst, src1, &nsrc2 );
}

#ifdef C_IMPL
void dadd( void *dst, void *src1, void *src2 ) {
    c_dadd( dst, src1, src2 );
}

void dsub( void *dst, void *src1, void *src2 ) {
    c_dsub( dst, src1, src2 );
}

void dmul( void *dst, void *src1, void *src2 ) {
    uint64 f1, f2;
    int e1, e2;
    mpnum num1, num2;
    int isNeg1 = d_isNeg( src1 );
    int isNeg2 = d_isNeg( src2 );
    int isDstNeg = (isNeg1 && !isNeg2) || (!isNeg1 && isNeg2);
    int newexp;

    /* special cases */
    if( d_isNan( src1 ) || d_isNan( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isInf( src1 ) && d_isZero( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isZero( src1 ) && d_isInf( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isInf( src1 ) ) {
        *(uint64*)dst = DOUBLE_POS_INF;
        if( isDstNeg ) dneg( dst, dst );
    } else if( d_isInf( src2 ) ) {
        *(uint64*)dst = DOUBLE_POS_INF;
        if( isDstNeg ) dneg( dst, dst );
    } else if( d_isZero( src1 ) || d_isZero( src2 ) ) {
        *(uint64*)dst = DOUBLE_POS_ZERO;
        if( isDstNeg ) dneg( dst, dst );
    } else {
        parseDouble( *(double*)src1, &f1, &e1 );
        parseDouble( *(double*)src2, &f2, &e2 );
        mp_init( &num1, f1 );
        mp_init( &num2, f2 );
        mp_mul( &num1, &num1, &num2 );
        newexp = e1 + e2 - 105 + mp_bitsize( &num1 ) + double_bias;
        if( newexp <= 0 && newexp > -53 ) {
            /* denormalized result */
            uint64 temp;
            mp_binround( &num1, &num1, -newexp + 53 );
            mp_shiftright( &num1, &num1, -newexp + 53 );
            mp_touint64( &temp, &num1 );
            makeDouble( dst, temp, 0 - double_bias );
        } else {
            mp_todouble( dst, &num1 );
            changeDoubleExponent( dst, dst, e1 + e2 - 104 );
        }
        if( isDstNeg ) dneg( dst, dst );
        mp_free( &num1 );
        mp_free( &num2 );
    }
}

void ddiv( void *dst, void *src1, void *src2 ) {
    uint64 f1, f2;
    int e1, e2;
    mpnum num1, num2, rem;
    int isNeg1 = d_isNeg( src1 );
    int isNeg2 = d_isNeg( src2 );
    int isDstNeg = (isNeg1 && !isNeg2) || (!isNeg1 && isNeg2);
    int newexp, rc;
    int adjust = 0;

    /* special cases */
    if( d_isNan( src1 ) || d_isNan( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isInf( src1 ) && d_isInf( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isInf( src1 ) ) {
        *(uint64*)dst = DOUBLE_POS_INF;
        if( isDstNeg ) dneg( dst, dst );
    } else if( d_isInf( src2 ) ) {
        *(uint64*)dst = DOUBLE_POS_ZERO;
        if( isDstNeg ) dneg( dst, dst );
    } else if( d_isZero( src1 ) && d_isZero( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isZero( src1 ) ) {
        *(uint64*)dst = DOUBLE_POS_ZERO;
        if( isDstNeg ) dneg( dst, dst );
    } else if( d_isZero( src2 ) ) {
        *(uint64*)dst = DOUBLE_POS_INF;
        if( isDstNeg ) dneg( dst, dst );
    } else {
        parseDouble( *(double*)src1, &f1, &e1 );
        parseDouble( *(double*)src2, &f2, &e2 );
        mp_init( &num1, f1 );
        mp_init( &num2, f2 );
        mp_init( &rem, 0 );
        if( mp_lt( &num1, &num2 ) ) {
            adjust = 1;
            e1--;
        }
        newexp = e1 - e2 + double_bias;
        if( newexp <= 0 && newexp > -52 ) {
            mp_shiftleft( &num1, &num1, 51 + newexp + adjust );
        } else {
            mp_shiftleft( &num1, &num1, 52 + adjust );
        }
        mp_div( &num1, &rem, &num1, &num2 );
        mp_mulsc( &rem, &rem, 2 );
        rc = mp_cmp( &rem, &num2 );
        if( (rc == 0 && (rem.num[0] & 2)) || rc == -1 ) {
            /* round up */
            mp_addsc( &num1, &num1, 1 );
        }
        if( newexp <= 0 && newexp > -52 ) {
            /* denormalized result */
            uint64 temp;
            mp_touint64( &temp, &num1 );
            makeDouble( dst, temp, 0 - double_bias );
        } else {
            mp_todouble( dst, &num1 );
            changeDoubleExponent( dst, dst, e1 - e2 - 52 );
        }
        if( isDstNeg ) dneg( dst, dst );
        mp_free( &num1 );
        mp_free( &num2 );
        mp_free( &rem );
    }
}

void drem( void *dst, void *src1, void *src2 ) {
    uint64 f1, f2;
    int e1, e2;
    mpnum num1, num2, zero;
    int emin;
    int newexp;

    /* special cases */
    if( d_isNan( src1 ) || d_isNan( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isInf( src1 ) || d_isZero( src2 ) ) {
        *(uint64*)dst = DOUBLE_NAN;
    } else if( d_isZero( src1 ) || d_isInf( src2 ) ) {
        *(uint64*)dst = *(uint64*)src1;
    } else {
        parseDouble( *(double*)src1, &f1, &e1 );
        parseDouble( *(double*)src2, &f2, &e2 );
        mp_init( &num1, f1 );
        mp_init( &num2, f2 );
        mp_init( &zero, 0 );
        emin = min( e1, e2 );
        if( e2 > e1 ) {
            mp_shiftleft( &num2, &num2, e2 - e1 );
        } else {
            mp_shiftleft( &num1, &num1, e1 - e2 );
        }
        mp_div( &num1, &num2, &num1, &num2 ); /* num2 will get the remainder */
        if( mp_eq( &num2, &zero ) ) {
            *(uint64*)dst = DOUBLE_POS_ZERO;
        } else {
            newexp = emin - 53 + mp_bitsize( &num2 ) + double_bias;
            if( newexp <= 0 && newexp > -53 ) {
                /* denormalized result */
                uint64 temp;
                mp_binround( &num2, &num2, -newexp + 1 );
                mp_shiftright( &num2, &num2, -newexp + 1 );
                mp_touint64( &temp, &num2 );
                makeDouble( dst, temp, 0 - double_bias );
            } else {
                mp_todouble( dst, &num2 );
                changeDoubleExponent( dst, dst, emin - 52 );
            }
        }
        if( d_isNeg( src1 ) ) dneg( dst, dst );
        mp_free( &num1 );
        mp_free( &num2 );
        mp_free( &zero );
    }
}

void dneg( void *dst, void *src ) {
    memcpy( dst, src, 8 );
    if( !d_isNan( dst ) ) {
        ((char*)dst)[7] ^= 128;
    }
}

void fadd( void *dst, void *src1, void *src2 ) {
    uint64 f1, f2;
    int e1, e2;
    int diff;
    int emax, emin;
    mpnum num1, num2, zero;
    int isNeg1 = f_isNeg( src1 );
    int isNeg2 = f_isNeg( src2 );
    int isNeg = FALSE;
    int newexp;

    /* special cases */
    if( f_isNan( src1 ) || f_isNan( src2 ) ) {
        *(uint32*)dst = FLOAT_NAN;
    } else if( f_isPosInf( src1 ) && f_isNegInf( src2 ) ) {
        *(uint32*)dst = FLOAT_NAN;
    } else if( f_isNegInf( src1 ) && f_isPosInf( src2 ) ) {
        *(uint32*)dst = FLOAT_NAN;
    } else if( f_isPosInf( src1 ) || f_isPosInf( src2 ) ) {
        *(uint32*)dst = FLOAT_POS_INF;
    } else if( f_isNegInf( src1 ) || f_isNegInf( src2 ) ) {
        *(uint32*)dst = FLOAT_NEG_INF;
    } else if( f_isPosZero( src1 ) && f_isNegZero( src2 ) ) {
        *(uint32*)dst = FLOAT_POS_ZERO;
    } else if( f_isNegZero( src1 ) && f_isPosZero( src2 ) ) {
        *(uint32*)dst = FLOAT_POS_ZERO;
    } else if( f_isPosZero( src1 ) && f_isPosZero( src2 ) ) {
        *(uint32*)dst = FLOAT_POS_ZERO;
    } else if( f_isNegZero( src1 ) && f_isNegZero( src2 ) ) {
        *(uint32*)dst = FLOAT_NEG_ZERO;
    } else if( f_isZero( src1 ) ) {
        *(uint32*)dst = *(uint32*)src2;
    } else if( f_isZero( src2 ) ) {
        *(uint32*)dst = *(uint32*)src1;
    } else {
        parseFloat( *(float*)src1, &f1, &e1 );
        parseFloat( *(float*)src2, &f2, &e2 );
        diff = abs( (e1 - e2) );
        emax = max( e1, e2 );
        emin = min( e1, e2 );
        mp_init( &num1, f1 );
        mp_init( &num2, f2 );
        mp_init( &zero, 0 );
        if( e1 > e2 ) {
            mp_shiftleft( &num1, &num1, (uint64)diff );
        } else if( e2 > e1 ) {
            mp_shiftleft( &num2, &num2, (uint64)diff );
        }
        if( (isNeg1 && isNeg2) || (!isNeg1 && !isNeg2 ) ) {
            mp_add( &num1, &num1, &num2 );
            if( isNeg1 ) isNeg = TRUE;
        } else {
            if( mp_gt( &num1, &num2 ) ) {
                mp_sub( &num1, &num1, &num2 );
                if( isNeg1 ) isNeg = TRUE;
            } else {
                mp_sub( &num1, &num2, &num1 );
                if( isNeg2 ) isNeg = TRUE;
            }
        }
        if( mp_eq( &num1, &zero ) ) {
            *(uint32*)dst = FLOAT_POS_ZERO;
        } else {
            newexp = emax - diff - 24 + mp_bitsize( &num1 ) + float_bias;
            if( newexp <= 0 && newexp > -24 ) {
                /* denormalized result */
                uint64 temp;
                mp_binround( &num1, &num1, -(emin + float_bias) + 1 );
                mp_shiftright( &num1, &num1, -(emin + float_bias) + 1 );
                mp_touint64( &temp, &num1 );
                makeFloat( dst, temp, 0 - float_bias );
            } else {
                mp_tofloat( dst, &num1 );
                changeFloatExponent( dst, dst, emax - diff - 23 );
            }
            if( isNeg ) fneg( dst, dst );
        }
        mp_free( &num1 );
        mp_free( &num2 );
        mp_free( &zero );
    }
}

void fsub( void *dst, void *src1, void *src2 ) {
    float nsrc2;
    fneg( &nsrc2, src2 );
    fadd( dst, src1, &nsrc2 );
}

void fmul( void *dst, void *src1, void *src2 ) {
    uint64 f1, f2;
    int e1, e2;
    mpnum num1, num2;
    int isNeg1 = f_isNeg( src1 );
    int isNeg2 = f_isNeg( src2 );
    int isDstNeg = (isNeg1 && !isNeg2) || (!isNeg1 && isNeg2);
    int newexp;

    /* special cases */
    if( f_isNan( src1 ) || f_isNan( src2 ) ) {
        *(uint32*)dst = FLOAT_NAN;
    } else if( f_isInf( src1 ) && f_isZero( src2 ) ) {
        *(uint32*)dst = FLOAT_NAN;
    } else if( f_isZero( src1 ) && f_isInf( src2 ) ) {
        *(uint32*)dst = FLOAT_NAN;
    } else if( f_isInf( src1 ) ) {
        *(uint32*)dst = FLOAT_POS_INF;
        if( isDstNeg ) fneg( dst, dst );
    } else if( f_isInf( src2 ) ) {
        *(uint32*)dst = FLOAT_POS_INF;
        if( isDstNeg ) fneg( dst, dst );
    } else if( f_isZero( src1 ) || f_isZero( src2 ) ) {
        *(uint32*)dst = FLOAT_POS_ZERO;
        if( isDstNeg ) fneg( dst, dst );
    } else {
        parseFloat( *(float*)src1, &f1, &e1 );
        parseFloat( *(float*)src2, &f2, &e2 );
        mp_init( &num1, f1 );
        mp_init( &num2, f2 );
        mp_mul( &num1, &num1, &num2 );
        newexp = e1 + e2 - 47 + mp_bitsize( &num1 ) + float_bias;
        if( newexp <= 0 && newexp > -24 ) {
            /* denormalized result */
            uint64 temp;
            mp_binround( &num1, &num1, -newexp + 24 );
            mp_shiftright( &num1, &num1, -newexp + 24 );
            mp_touint64( &temp, &num1 );
            makeFloat( dst, temp, 0 - float_bias );
        } else {
            mp_tofloat( dst, &num1 );
            changeFloatExponent( dst, dst, e1 + e2 - 46 );
        }
        if( isDstNeg ) fneg( dst, dst );
        mp_free( &num1 );
        mp_free( &num2 );
    }
}

void fdiv( void *dst, void *src1, void *src2 ) {
    uint64 f1, f2;
    int e1, e2;
    mpnum num1, num2, rem;
    int isNeg1 = f_isNeg( src1 );
    int isNeg2 = f_isNeg( src2 );
    int isDstNeg = (isNeg1 && !isNeg2) || (!isNeg1 && isNeg2);
    int newexp, rc;
    int adjust = 0;

    /* special cases */
    if( f_isNan( src1 ) || f_isNan( src2 ) ) {
        *(uint32*)dst = FLOAT_NAN;
    } else if( f_isInf( src1 ) && f_isInf( src2 ) ) {
        *(uint32*)dst = FLOAT_NAN;
    } else if( f_isInf( src1 ) ) {
        *(uint32*)dst = FLOAT_POS_INF;
        if( isDstNeg ) fneg( dst, dst );
    } else if( f_isInf( src2 ) ) {
        *(uint32*)dst = FLOAT_POS_ZERO;