jwasm/tbyte.c
2014-06-18 19:16:56 +04:00

469 lines
14 KiB
C

/****************************************************************************
*
* 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: Conversion from string to tbyte format (long double)
* it doesn't depend on processor/math coprocessor type
*
****************************************************************************/
/*
It use 96bit x 96bit -> 192bit multiply for converting decimal to
80-bit extended double.
Math operation is on 96-bit wide operands (mantisa).
*/
#include <ctype.h>
#include "globals.h"
#include "tbyte.h"
#if defined(LLONG_MAX) || defined(__GNUC__) || defined(__TINYC__)
#define MAXUI64BIT 0x8000000000000000ULL
#define MAXUI64 0xffffffffffffffffULL
#else
#define MAXUI64BIT 0x8000000000000000ui64
#define MAXUI64 0xffffffffffffffffui64
#endif
union u192 {
uint_64 m64[3];
uint_32 m32[6];
};
union u96 {
uint_32 m32[3];
};
/* extended (112-bit, 96+16) long double */
struct ELD {
uint_32 m32[3];
unsigned short e;
};
#define EXPONENT_BIAS 0x3FFF
#define SET_SIGN(x,y) if( y < 0 ) x->e |= 0x8000
#define U96ISNOTZERO(x) (x.m32[0] || x.m32[1] || x.m32[2])
#define MAX_EXP_INDEX 13 /* # of items in tab_plus_exp[] and tab_minus_exp[] */
static const struct ELD tab_plus_exp[MAX_EXP_INDEX] = {
{ { 0x00000000UL, 0x00000000UL, 0xA0000000UL }, 0x4002 }, /* 1e1L */
{ { 0x00000000UL, 0x00000000UL, 0xC8000000UL }, 0x4005 }, /* 1e2L */
{ { 0x00000000UL, 0x00000000UL, 0x9C400000UL }, 0x400C }, /* 1e4L */
{ { 0x00000000UL, 0x00000000UL, 0xBEBC2000UL }, 0x4019 }, /* 1e8L */
{ { 0x00000000UL, 0x04000000UL, 0x8E1BC9BFUL }, 0x4034 }, /* 1e16L */
{ { 0xF0200000UL, 0x2B70B59DUL, 0x9DC5ADA8UL }, 0x4069 }, /* 1e32L */
{ { 0x3CBF6B71UL, 0xFFCFA6D5UL, 0xC2781F49UL }, 0x40D3 }, /* 1e64L */
{ { 0xC66F336BUL, 0x80E98CDFUL, 0x93BA47C9UL }, 0x41A8 }, /* 1e128L */
{ { 0xDDBB9018UL, 0x9DF9DE8DUL, 0xAA7EEBFBUL }, 0x4351 }, /* 1e256L */
{ { 0xCC655C4BUL, 0xA60E91C6UL, 0xE319A0AEUL }, 0x46A3 }, /* 1e512L */
{ { 0x650D3D17UL, 0x81750C17UL, 0xC9767586UL }, 0x4D48 }, /* 1e1024L*/
{ { 0xA74D28B1UL, 0xC53D5DE4UL, 0x9E8b3B5DUL }, 0x5A92 }, /* 1e2048L*/
{ { 0xC94C14F7UL, 0x8A20979AUL, 0xC4605202UL }, 0x7525 } /* 1e4096L*/
};
static const struct ELD tab_minus_exp[MAX_EXP_INDEX] = {
{ { 0xCCCCCCCDUL, 0xCCCCCCCCUL, 0xCCCCCCCCUL }, 0x3FFB }, /* 1e-1L */
{ { 0x3D70A3D7UL, 0x70A3D70AUL, 0xA3D70A3DUL }, 0x3FF8 }, /* 1e-2L */
{ { 0xD3C36113UL, 0xE219652BUL, 0xD1B71758UL }, 0x3FF1 }, /* 1e-4L */
{ { 0xFDC20D2AUL, 0x8461CEFCUL, 0xABCC7711UL }, 0x3FE4 }, /* 1e-8L */
{ { 0x4C2EBE65UL, 0xC44DE15BUL, 0xE69594BEUL }, 0x3FC9 }, /* 1e-16L */
{ { 0x67DE18E7UL, 0x453994BAUL, 0xCFB11EADUL }, 0x3F94 }, /* 1e-32L */
{ { 0x3F2398CCUL, 0xA539E9A5UL, 0xA87FEA27UL }, 0x3F2A }, /* 1e-64L */
{ { 0xAC7CB3D9UL, 0x64BCE4A0UL, 0xDDD0467CUL }, 0x3E55 }, /* 1e-128L */
{ { 0xFA911122UL, 0x637A1939UL, 0xC0314325UL }, 0x3CAC }, /* 1e-256L */
{ { 0x7132D2E4UL, 0xDB23D21CUL, 0x9049EE32UL }, 0x395A }, /* 1e-512L */
{ { 0x87A600A6UL, 0xDA57C0BDUL, 0xA2A682A5UL }, 0x32B5 }, /* 1e-1024L*/
{ { 0x4925110FUL, 0x34362DE4UL, 0xCEAE534FUL }, 0x256B }, /* 1e-2048L*/
{ { 0x2DE37E46UL, 0xD2CE9FDEUL, 0xA6DD04C8UL }, 0x0AD8 } /* 1e-4096L*/
};
static int cmp_u96_max( const union u96 *x )
/*******************************************
compare u96 with maximum value before u96
overflow after multiply by 10
*/
{
if( x->m32[2] > 0x19999999UL ) {
return( 1 );
} else if( x->m32[2] < 0x19999999UL ) {
return( -1 );
} else if( x->m32[1] > 0x99999999UL ) {
return( 1 );
} else if( x->m32[1] < 0x99999999UL ) {
return( -1 );
} else if( x->m32[0] > 0x99999998UL ) {
return( 1 );
} else if( x->m32[0] < 0x99999998UL ) {
return( -1 );
} else {
return( 0 );
}
}
static int add_check_u96_overflow( union u96 *x, unsigned int c )
/****************************************************************
test u96 overflow after multiply by 10
add one decimal digit to u96
*/
{
uint_64 cy;
int i;
if( cmp_u96_max( x ) > 0 ) {
return( 1 );
} else {
cy = c;
for( i = 0; i < 3; i++ ) {
cy += (uint_64)x->m32[i] * 10;
x->m32[i] = cy;
cy >>= 32;
}
return( 0 );
}
}
static int bitsize32( uint_32 x )
/********************************
calculate bitsize for uint_32
*/
{
int i;
for( i = 32; i > 0 ; i-- ) {
if( x & 0x80000000U ) break;
x <<= 1;
}
return( i );
}
static int bitsize64( uint_64 x )
/********************************
calculate bitsize for uint_64
*/
{
int i;
for( i = 64; i > 0 ; i-- ) {
if( x & MAXUI64BIT ) break;
x <<= 1;
}
return( i );
}
static int U96LD( const union u96 *op, struct ELD *res )
/*******************************************************
convert u96 into internal extended long double
*/
{
int bs;
int shft;
memcpy( res, op, sizeof( union u96 ) );
bs = bitsize32( res->m32[2] ) + 64;
if( bs == 64 ) {
res->m32[2] = res->m32[1];
res->m32[1] = res->m32[0];
res->m32[0] = 0;
bs = bitsize32( res->m32[2] ) + 32;
}
if( bs == 32 ) {
res->m32[2] = res->m32[1];
res->m32[1] = res->m32[0];
res->m32[0] = 0;
bs = bitsize32( res->m32[2] );
}
if( bs == 0 ) {
res->e = 0;
} else {
res->e = bs - 1 + EXPONENT_BIAS;
bs %= 32;
if( bs ) {
shft = 32 - bs;
res->m32[2] <<= shft;
res->m32[2] |= res->m32[1] >> bs;
res->m32[1] <<= shft;
res->m32[1] |= res->m32[0] >> bs;
res->m32[0] <<= shft;
}
}
return( 0 );
}
static int normalize( union u192 *res )
/**************************************
normalize internal extended long double u192
return exponent shift
*/
{
int shft;
int bs;
int bs1;
bs = bitsize64( res->m64[2] ) + 128;
if( bs == 128 ) {
res->m64[2] = res->m64[1];
res->m64[1] = res->m64[0];
res->m64[0] = 0;
bs = bitsize64( res->m64[2] ) + 64;
}
if( bs == 64 ) {
res->m64[2] = res->m64[1];
res->m64[1] = res->m64[0];
res->m64[0] = 0;
bs = bitsize64( res->m64[2] );
}
if( bs == 0 ) {
return( 0 );
}
bs1 = bs % 64;
if ( bs1 ) {
shft = 64 - bs1;
res->m64[2] <<= shft;
res->m64[2] |= res->m64[1] >> bs1;
res->m64[1] <<= shft;
res->m64[1] |= res->m64[0] >> bs1;
res->m64[0] <<= shft;
}
return( bs - 192 );
}
static int add192( union u192 *res, const uint_64 x, int pos )
/*************************************************************
add uint_64 to u192 on uint_32 position
*/
{
uint_64 cy;
int i;
cy = (uint_32)x;
for( i = pos; i < 6; i++ ) {
cy += res->m32[i];
res->m32[i] = cy;
cy >>= 32;
}
cy = x >> 32;
for( i = pos + 1; i < 6; i++ ) {
cy += res->m32[i];
res->m32[i] = cy;
cy >>= 32;
}
return( 0 );
}
static int multiply( const struct ELD *op1, const struct ELD *op2, struct ELD *res )
/***********************************************************************************
multiply u96 by u96 into u96
normalize and round result
*/
{
uint_64 x1;
union u192 r1;
int_32 exp;
exp = (int_32)(op1->e & 0x7fff) + (int_32)(op2->e & 0x7fff) - EXPONENT_BIAS + 1;
r1.m64[0] = (uint_64)(op1->m32[0]) * (uint_64)(op2->m32[0]);
r1.m64[1] = (uint_64)(op1->m32[1]) * (uint_64)(op2->m32[1]);
r1.m64[2] = (uint_64)(op1->m32[2]) * (uint_64)(op2->m32[2]);
x1 = (uint_64)(op1->m32[1]) * (uint_64)(op2->m32[0]);
add192( &r1, x1, 1 );
x1 = (uint_64)(op1->m32[0]) * (uint_64)(op2->m32[1]);
add192( &r1, x1, 1 );
x1 = (uint_64)(op1->m32[0]) * (uint_64)(op2->m32[2]);
add192( &r1, x1, 2 );
x1 = (uint_64)(op1->m32[2]) * (uint_64)(op2->m32[0]);
add192( &r1, x1, 2 );
x1 = (uint_64)(op1->m32[1]) * (uint_64)(op2->m32[2]);
add192( &r1, x1, 3 );
x1 = (uint_64)(op1->m32[2]) * (uint_64)(op2->m32[1]);
add192( &r1, x1, 3 );
exp += normalize( &r1 );
/* round result */
if( r1.m32[2] & 0x80000000U ) {
if( r1.m32[5] == 0xffffffffU && r1.m32[4] == 0xffffffffU && r1.m32[3] == 0xffffffffU ) {
r1.m32[3] = 0;
r1.m32[4] = 0;
r1.m32[5] = 0x80000000U;
exp++;
} else {
x1 = 1L;
add192( &r1, x1, 3 );
}
}
res->m32[0] = r1.m32[3];
res->m32[1] = r1.m32[4];
res->m32[2] = r1.m32[5];
res->e = exp;
return( 0 );
}
static int TB_create( union u96 *value, int_32 exponent, struct TB_LD *ld )
/**************************************************************************
create tbyte/long double from u96 value and
decimal exponent, round result
*/
{
const struct ELD *tabExp;
int i;
struct ELD res;
if( exponent < 0 ) {
exponent = -exponent;
tabExp = tab_minus_exp;
} else {
tabExp = tab_plus_exp;
}
U96LD(value, &res);
for( i = 0; i < MAX_EXP_INDEX; i++ ) {
if ( exponent & 1 ) {
multiply(&res, tabExp + i, &res);
}
exponent >>= 1;
if( exponent == 0 ) break;
}
if( exponent != 0 ) {
/* exponent overflow */
}
ld->e = res.e;
ld->m = res.m32[1] + ((uint_64)res.m32[2] << 32) ;
/* round result */
if( res.m32[0] & 0x80000000U ) {
if( ld->m == MAXUI64 ) {
ld->m = MAXUI64BIT;
ld->e++;
} else {
ld->m++;
}
}
return( 0 );
}
struct TB_LD *strtotb( const char *p, struct TB_LD *ld, char negative )
/**********************************************************************
convert string into tbyte/long double
set result sign
*/
{
int sign = +1;
int exp_sign = +1;
int_32 exp_value;
int overflow;
int_32 exp1;
int_32 exponent;
int_32 exponent_tmp;
union u96 value;
union u96 value_tmp;
while ( isspace( *p ) ) p++;
switch (*p) {
case '-':
sign = -1;
case '+':
p++;
default :
break;
}
if( negative ) {
sign = -sign;
}
memset( &value, 0, sizeof( value ) );
memset( &value_tmp, 0, sizeof( value_tmp ) );
//exponent = 0;
exp1 = 0;
exponent_tmp = 0;
overflow = 0;
while ( (unsigned int)(*p - '0') < 10u ) {
if( overflow ) {
exponent_tmp++;
exp1++;
} else {
if( add_check_u96_overflow( &value_tmp, *p - '0' ) ) {
overflow = 1;
exponent_tmp++;
exp1++;
} else if( *p != '0' ) {
memcpy( &value, &value_tmp, sizeof( value ) );
exp1 = 0;
} else if( U96ISNOTZERO(value) ) {
exp1++;
}
}
p++;
}
exponent = exp1;
if ( *p == '.' ) {
p++;
while ( (unsigned int)(*p - '0') < 10u ) {
if( overflow == 0 ) {
if( add_check_u96_overflow( &value_tmp, *p - '0' ) ) {
overflow = 1;
} else {
exponent_tmp--;
if( *p != '0' ) {
memcpy( &value, &value_tmp, sizeof( value ) );
exponent = exponent_tmp;
}
}
}
p++;
}
}
exp_value = 0;
if ( (*p | 0x20) == 'e' ) {
switch ( *++p ) {
case '-':
exp_sign = -1;
case '+': p++;
break;
case '0':
case '1':
case '2':
case '3':
case '4':
case '5':
case '6':
case '7':
case '8':
case '9':
break;
default :
ld->m = 0;
ld->e = 0;
SET_SIGN( ld, sign );
return( ld );
}
while ( (unsigned int)(*p - '0') < 10u )
exp_value = 10 * exp_value + (*p++ - '0');
if( exp_sign < 0 )
exp_value = -exp_value;
}
exp_value += exponent;
TB_create( &value, exp_value, ld );
SET_SIGN( ld, sign );
return( ld );
}