github.com/prirun.p50em
1
0
Fork 0
mirror of https://github.com/prirun/p50em.git synced 2026-02-19 05:16:58 +00:00
Code Issues Releases Wiki Activity

fp.c now fp.h, FP rewrite, more CC macros: CLEARCC, SETEQ, SETLT

Browse Source
This commit is contained in:
Jim
2007-06-25 00:00:00 -04:00
parent 31f4e2e22e
commit 29334a2945
3 changed files with 950 additions and 788 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1067
em.c
View File

File diff suppressed because it is too large Load Diff

304
fp.c
View File

@@ -1,304 +0,0 @@
/* ** general conversion doc here, or (referenced) above. ***************** */
/* doc exponent details here: (-1 & bin pt. placement by prieee) vs.
vs. (-2 by ieeepr()): note that
IEEE denormalized bin pxxxxxxx prime (relative to the 23 bits);
in ieeepr the trick is to treat normalized #'s as the special case. */
void prieee4(xlongp) /* ****** convert from Prime to IEEE. ****** */
long *xlongp;
{
long xlong;
long sign; /* sign */
long mant; /* mantissa: note this is signed: mant = -mant below */
long exp; /* exponent */
/* note: when converting form Prime to IEEE, we need to check for exponent
underflow, but not for overflow. */
xlong = *xlongp;
mant = xlong & 0xffffff00; /* 24 bits: includes sign bit */
sign = mant & 0x80000000;
mant >>= 8; /* ok to trash upper byte, which will get masked out. */
if (sign)
mant = - mant;
/* now have 24 bit # w/ leading 0. */
exp = (xlong & 0x000000ff) - 1; /* exp now in excess 127 (IEEE) form. */
if (exp < 0) {
mant >>= 1; /* tiny tiny #; will get shifted once more below. */
++exp; /* Bring exp back up to 0. */
}
else /* *** normalize the mantissa. If input Prime mantissa was already
noramalized, this loop will still execute once, unless input #
was -(2**n) in which case it will execute 0 times. This is
because Prime mantissas are 2's compliment, and powers of 2
normalize differently when negated. Note no loop for 0 exp. */
while (exp && !(mant & 0x00800000)) {
mant <<= 1;
--exp;
}
/* we now have a 24 bit unsigned # w/ a leading 1. If the exponent is > 0,
we are all set: the resulting IEEE number will be normalized. This
leading (hidden) IEEE 1 will be masked out below. */
if ( ! exp) /* result must be denormalized: shift 1 further to */
mant >>= 1; /* include IEEE leading bit, which may be 0 or 1. */
mant &= 0x007fffff;
exp <<= 23;
*xlongp = sign | exp | mant;
}
void ieeepr4(xlongp) /* ****** convert from IEEE to Prime. ****** */
long *xlongp;
{
long xlong;
long sign; /* sign */
long mant; /* mantissa: note this is signed: mant = -mant below */
long exp; /* exponent */
long templ;
/* note: when converting form Prime to IEEE, we need to check for exponent
overflow, but not for underflow. */
xlong = *xlongp;
if ( ! (xlong & 0x7fffffff))
*xlongp = 0; /* +0 or -0 */
else {
sign = xlong & 0x80000000;
mant = (xlong & 0x007fffff) << 8; /* assume denormalized, adjust below */
exp = (xlong & 0x7f800000) >> 23; /* still in excess 127 (IEEE) form. */
if (exp == 0xff) { /* NaN (not a number) or +/- infinity? */
if (mant == 0) {
/* +/- infinity. */
if (sign == 0)
*xlongp = 0x7fffffff; /* largest Prime # in for +infinity. */
else /* note: 0x800000ff cant be negated; */
*xlongp = 0x800001ff; /* use - 0x7fffffff for -infinity. */
}
else {
/* NaN (not a number) */
if (sign == 0)
*xlongp = 0x539733ff; /* use 1.11111e+38 for NaN, since */
else /* it stands out if printed. */
*xlongp = 0xac68cdff; /* use -1.11111e+38 for -NaN. */
}
}
else { /* actual number */
if (exp != 0) {
if(mant != 0x7fffff00)/* Special case of mantissa value of all 1s*/
mant = (mant >> 1) + 0x40000080; /* ieee normalized number: */
/* shift to make room for */
/* hidden leading 1 bit, 'or'*/
/* it in and round truncated */
/* LSBit up. */
mant &= 0xffffff00;
}
exp += 2; /* exp now in excess 128 (Prime) form. */
if (sign != 0) /* sign bit of mant will be 0 at this point. */
mant = - mant;
/* mant is now a 24 bit signed mantissa (shifted to Prime position) */
/* *** normalize the number. In most cases, the number will already
be normalized. Negative powers of 2 will be shifted once, since
they normalize differently. IEEE denormalized numbers can be
adjusted at most 2 bits, due to the excess 127 vs. 128 exponent and
binary point position differences (cant negate exponent). *** */
while (exp && ((~(mant^(mant<<1))) & 0x80000000)) {
mant <<= 1; /* sign and next most significant */
--exp; /* bit are equal: normalize. */
}
if (exp > 0xff) {
if ( ! sign)
*xlongp = 0x7fffffff; /* largest Prime # in for +infinity. */
else /* note: 0x800000ff cant be negated; */
*xlongp = 0x800001ff; /* use - 0x7fffffff for -infinity. */
}
else
*xlongp = mant | exp; /* mant is signed; */
}
}
}
#define l1 0
#define l0 1
/* ** general conversion doc here, or (referenced) above. ***************** */
/* doc exponent details here: (-1 & bin pt. placement by prieee) vs.
vs. (-2 by ieeepr()): note that
IEEE denormalized bin pxxxxxxx prime (relative to the 23 bits);
in ieeepr the trick is to treat normalized #'s as the special case. */
void prieee8(xlongp) /* ****** convert from Prime to IEEE. ****** */
long *xlongp;
{
long xlong1,xlong0;
long mant1,mant0; /* mantissa: these are signed: mant = -mant below */
long sign; /* sign */
long exp; /* exponent */
/* note: when converting REAL*8 form Prime to IEEE, we need to check for both
exponent underflow, and overflow. */
xlong1 = xlongp[l1]; /* high 4 bytes */
xlong0 = xlongp[l0]; /* low 4 bytes */
sign = xlong1 & 0x80000000;
mant1 = xlong1 >> 11; /* 48 bits, includes */
mant0 = (xlong1 << 21) | ((xlong0 >> 11) & 0xfffe0); /* includes sign bit */
if ( ! mant0 && ! mant1) { /* zero (dirty or otherwise)? */
xlongp[l1] = xlongp[l0] = 0;
return; /* return ****** */
}
if (sign) { /* 2's comp mant1/mant0 pair */
mant0 = - mant0;
mant1 = ~ mant1;
if ( ! mant0)
++ mant1; /* mant0==0: have a carry fr. mant0 to mant1. */
}
/* now have 48 bit # w/ leading 0. */
exp = (xlong0 & 0xffff) - 0x80; /* convert 16 bit */
if (exp & 0x8000) /* Prime exponent */
exp |= 0xffff0000; /* from "excess 128" */
else /* form to 2's */
exp &= 0x0000ffff; /* complement form. */
exp += 0x3ff; /* exp now in excess 1023 (IEEE) form. */
if (exp < -50) { /* Prime exp too small? */
xlongp[l1] = xlongp[l0] = 0; /* closest IEEE is 0. */
return; /* return ****** */
}
if (exp < 0) {
/* note: can still get 0, iff have leading (non sign bit) 0's */
mant0 = (mant1 << 32+exp) | (mant0 >> -exp); /* mant >>= -exp; */
mant1 >>= -exp; /* tiny tiny #; will get shifted once more below. */
++mant0; /* lsb will get lost below, since # is denormalized; round */
if ( ! mant0)
++mant1;
exp = 0; /* exp += -exp. */
}
else /* *** normalize the mantissa. */
while (exp && !(mant1 & 0x00100000)) {
mant1 <<= 1;
if (mant0 & 0x80000000)
mant1 |= 1;
mant0 <<= 1;
--exp;
}
if (exp > 0x7fe) {
xlongp[l1] = sign | 0x7fefffff; /* Prime exp too big; closest */
xlongp[l0] = 0xffffffff; /* closest IEEE is (+/-) max. */
}
else {
/* we now have a 48 bit unsigned # w/ a leading 1. If the exponent is
> 0, we are all set: the resulting IEEE number will be normalized.
This leading (hidden) IEEE 1 will be masked out below. */
if ( ! exp) { /* result must be denormalized: shift 1 further to */
mant0 >>= 1; /* include IEEE leading bit, which may be 0 or 1. */
if (mant1 & 1)
mant0 |= 0x80000000;
mant1 >>= 1;
}
mant1 &= 0x000fffff;
exp <<= 20;
xlongp[l1] = sign | exp | mant1;
xlongp[l0] = mant0;
}
}
void ieeepr8(xlongp) /* ****** convert from IEEE to Prime. ****** */
long *xlongp;
{
long xlong1,xlong0;
long mant1,mant0; /* mantissa: these are signed: mant = -mant below */
long sign; /* sign */
long exp; /* exponent */
/* note: when converting REAL*8 form Prime to IEEE, we dont need to check for
exponent overflow, or underflow. */
xlong1 = xlongp[l1]; /* high 4 bytes */
xlong0 = xlongp[l0]; /* low 4 bytes */
if ((xlong1 & 0x7fffffff) == 0 && xlong0 == 0) {
xlongp[l1] = xlongp[l0] = 0; /* +0 or -0 */
return; /* return ****** */
}
sign = xlong1 & 0x80000000;
exp = (xlong1 & 0x7ff00000) >> 20; /* still in excess 1023 (IEEE) form. */
mant1 = ((xlong1 & 0xfffff) << 11) | ((xlong0 >> 21) & 0x7ff);
mant0 = xlong0 << 11; /* assume denormalized, adjust below */
if (exp == 0x7ff) { /* NaN (not a number) or +/- infinity? */
if (mant1 == 0 && mant0 == 0) {
/* +/- infinity. */
if (sign == 0) {
xlongp[l1] = 0x7fffffff; /* largest Prime # in for +infinity. */
xlongp[l0] = 0xffffffff;
}
else { /* note: 0x80000000 0000ffff cant be */
xlongp[l1] = 0x80000000; /* negated, use -0x7fffffff ffffffff */
xlongp[l0] = 0x0001ffff; /* -infinity. */
}
}
else {
/* NaN (not a number) */
if (sign == 0) {
xlongp[l1] = 0x7fffffff; /* For now, use +/- infinity values */
xlongp[l0] = 0xffffffff; /* for +/- NaN. */
}
else {
xlongp[l1] = 0x80000000;
xlongp[l0] = 0x0001ffff;
}
}
}
else { /* actual number */
if (exp != 0) {
/* ieee normalized number: shift mantissa to make room for hidden */
mant0 >>= 1; /* leading 1 bit and 'or' it in. */
if (mant1 & 1)
mant0 |= 0x80000000;
else
mant0 &= 0x7fffffff;
mant1 >>= 1;
mant1 |= 0x40000000;
}
exp -= 894; /* exp now in 'excess 128' (Prime) form. */
mant0 &= 0xffff0000;
if (sign) { /* sign bit of mant will be 0 at this point. */
mant0 |= 0xffff; /* 2's comp mant1/mant0 pair */
mant0 = - mant0;
mant1 = ~ mant1;
if ( ! mant0)
++ mant1; /* mant0==0: have a carry fr. mant0 to mant1. */
}
/* mant is now a 48 bit signed mantissa (shifted to Prime position) */
/* *** normalize the number. In most cases, the number will already
be normalized. Negative powers of 2 will be shifted once, since
they normalize differently. IEEE denormalized numbers can be
adjusted at most 2 bits, due to the excess 127 vs. 128 exponent and
binary point position differences (cant negate exponent). *** */
while (exp && ((~(mant1^(mant1<<1))) & 0x80000000)) {
mant1 <<= 1; /* sign and next most significant */
if (mant0 & 0x80000000) /* bit are equal: normalize. */
mant1 |= 1;
mant0 <<= 1;
--exp;
}
mant0 &= 0xffff0000;
xlongp[l1] = mant1; /* mant is signed; no sign to 'or' in. */
xlongp[l0] = mant0 | exp;
}
}

367
fp.h Normal file
View File

@@ -0,0 +1,367 @@
/* fp.h, Jim Wilcoxson, June 2007
Floating point conversion and helper routines for the Prime emulator.
Prime DPFP format:
- 48 mantissa bits stored in 2's complement format
- 16 exponent bits stored in 2's complement with a bias of 128
- exponent follows the mantissa in both memory and registers
- some early Primes store the exponent in the 3rd halfword
- the leading "1" bit in the mantissa is only suppressed for negative
powers of 2
- Prime treats zero mantissa as 0.0, even if exponent is non-zero
- all FP operations are carried out in double precision
IEEE DPFP format:
- 1 sign bit
- 11 exponent bits with a bias of 1023
- 52 mantissa bits (sign-magnitude format)
- leading bit of mantissa has a suppressed "1" (except subnormal)
- if exponent is zero:
- and frac = 0 => positive or negative zero, depending on sign
- and frac non-zero => subnormal (aka denormal, unnormalized)
- subnormals have an implied leading "0" bit (XXX: true??)
*/
/* getdp unpacks a Prime DPFP into 48-bit sign + mantissa (left
justified in 64 bits) and a 32-bit signed exponent */
inline getdp (void *p, long long *frac64, int *exp32) {
*frac64 = *(long long *)p & 0xFFFFFFFFFFFF0000LL; /* unpack fraction */
*exp32 = *((short *)p+3); /* unpack SIGNED exponent */
}
inline putdp (void *p, long long frac64, int exp32) {
*(long long *)p = (frac64 & 0xFFFFFFFFFFFF0000LL) | (exp32 & 0xFFFF);
}
/* Conversion from Prime DPFP to IEEE DPFP
Returns true if conversion succeeded, false if it failed.
Conversions may fail because Prime exponents are 16 bits whereas
IEEE DPFP exponents are only 11 bits.
*/
int prieee8(void *dp, double *d) {
long long frac64, sign;
int exp32;
/* unpack Prime DPFP */
getdp (dp, &frac64, &exp32);
/* if negative, change to sign-magnitude */
sign = 0;
if (frac64 < 0) {
/* special case: negative power of 2 */
if (frac64 == 0x8000000000000000LL) {
exp32 += (1023-128);
if (exp32 < 0 || exp32 > 0x3ff)
return 0;
frac64 |= ((long long)exp32 << 52);
*d = *(double *)&frac64;
return 1;
} else {
sign = 0x8000000000000000LL;
frac64 = -frac64;
}
/* special case: zero */
} else if (frac64 == 0) {
*d = 0.0;
return 1;
}
/* normalize positive fraction until bit 2 is 1 */
while ((*(int *)&frac64 & 0x40000000) == 0) {
frac64 = frac64 << 1;
exp32--;
}
/* adjust exponent bias and check range */
exp32 += (1023-128) - 1;
if (exp32 < 0 || exp32 > 0x7ff)
return 0;
/* pack into an IEEE DPFP, losing the leading 1 bit in the process */
frac64 = sign | ((long long)exp32 << 52) | ((frac64 >> 10) & 0xfffffffffffffLL);
*d = *(double *)&frac64;
return 1;
}
/* Conversion from Prime SPFP to IEEE DPFP */
double prieee4(unsigned int sp) {
int frac32, sign;
long long frac64;
int exp32;
/* unpack Prime SPFP */
frac32 = sp & 0xFFFFFF00;
exp32 = sp & 0xFF;
/* if negative, change to sign-magnitude */
sign = 0;
if (frac32 < 0) {
/* special case: negative power of 2 */
if (frac32 == 0x80000000) {
exp32--;
frac64 = 0x8000000000000000LL | ((long long)exp32 << 52);
return *(double *)&frac64;
} else {
sign = 0x80000000;
frac32 = -frac32;
}
/* special case: zero */
} else if (frac32 == 0)
return 0.0;
/* normalize positive fraction until bit 2 is 1 */
while ((frac32 & 0x40000000) == 0) {
frac32 = frac32 << 1;
exp32--;
}
/* adjust exponent bias and check range */
exp32 += (1023-128) - 1;
/* pack into an IEEE DPFP, losing the leading 1 bit in the process */
frac64 = (long long)sign | ((long long)exp32 << 52) | (((long long)frac32 << 22) & 0xfffffffffffffLL);
return *(double *)&frac64;
}
/* conversion from IEEE back to Prime. Prime exponents are larger, so
this conversion cannot overflow/underflow, but precision is lost */
double ieeepr8(double d) {
long long frac64;
int exp32, neg;
/* unpack IEEE DPFP */
*(double *)&frac64 = d;
neg = frac64 < 0;
exp32 = (frac64 >> 52) & 0x7ff;
frac64 &= 0xfffffffffffffLL;
//printf("dp=%llx, neg=%d, frac64=%llx, exp32=%d, \n", *(long long *)dp, neg, frac64, exp32);
/* special case: NaN & +-infinity (these shouldn't happen!) */
if (exp32 == 0x7ff) {
if (frac64 == 0)
if (neg)
printf("em: +infinity in ieeepr8\n");
else
printf("em: -infinity in ieeepr8\n");
else
printf("em: NaN in ieeepr8\n");
return 0.0;
}
/* add back the hidden "1" bit except for the special cases +-0.0
and subnormal */
if (exp32 != 0) /* typical IEEE normalized */
frac64 |= 0x10000000000000LL;
else if (frac64 == 0) /* IEEE +-0.0 (zero exp+frac) */
return 0.0; /* IEEE and Prime zero are the same */
else
; /* subnormal: no hidden 1 bit */
/* adjust exponent, change sign-magnitude to 2's complement,
and shift fraction into Prime placement (high 48 bits) */
exp32 -= (1023-128) - 1;
if (neg)
frac64 = -frac64;
frac64 <<= 10;
/* normalize Prime DPFP */
while ((frac64 ^ (frac64 << 1)) >= 0) {
frac64 = frac64 << 1;
exp32--;
}
#if 0
if (exp32 > 32767 || exp32 < -32768) {
printf("em: exponent = %d in ieeepr8\n", exp32);
return 0.0;
}
#endif
/* round the fraction to 48 bits, ensuring no overflow */
if ((frac64 & 0x8000) && ((frac64 & 0x7fffffffffff0000LL) != 0x7fffffffffff0000LL))
/* XXX: should this be a subtract for negative numbers? */
frac64 += 0x10000;
frac64 = (frac64 & 0xffffffffffff0000LL) | (exp32 & 0xffff);
return *(double *)&frac64;
}
/* 32-bit signed integer to Prime DPFP conversion */
double fltl (int int32) {
long long frac64;
int exp32, sign32;
/* have to special case zero, or we end up with
a "dirty zero" (zero fraction, exponent of 128) */
if (int32 == 0)
return 0.0;
exp32 = 128+31;
sign32 = int32 & 0x80000000;
if ((int32 & 0xFFFF8000) == sign32>>16) {
int32 = int32<<16;
exp32 -= 16;
}
if ((int32 & 0xFF800000) == sign32>>8) {
int32 = int32<<8;
exp32 -= 8;
}
if ((int32 & 0xF8000000) == sign32>>4) {
int32 = int32<<4;
exp32 -= 4;
}
if ((int32 & 0xE0000000) == sign32>>2) {
int32 = int32<<2;
exp32 -= 2;
}
if ((int32 & 0xC0000000) == sign32>>1) {
int32 = int32<<1;
exp32 -= 1;
}
frac64 = ((long long)int32 << 32) | exp32;
return *(double *)&frac64;
}
/* Prime DPFP complement */
dfcm (void *dp) {
long long frac64;
int exp32;
CLEARC;
getdp(dp, &frac64, &exp32);
if (frac64 != 0) { /* can't normalize zero */
if (frac64 == 0x8000000000000000LL) { /* overflow case? */
frac64 = 0x4000000000000000LL; /* complement power of 2 */
exp32 += 1;
} else {
frac64 = -frac64; /* complement fraction */
while ((frac64 ^ (frac64 << 1)) >= 0) {
frac64 = frac64 << 1; /* normalize */
exp32--;
}
}
if (exp32 > 32767 || exp32 < -32768)
mathexception('f', FC_DFP_OFLOW, 0);
else
putdp(dp, frac64, exp32);
}
}
/* double precision floating point normalize
Passed a pointer to a Prime double precision variable and updates
it in place. May set the C-bit or cause a floating point
exception. */
void norm(void *dp) {
long long frac64;
int exp32;
getdp(dp, &frac64, &exp32);
while ((frac64 ^ (frac64 << 1)) >= 0) {
frac64 = frac64 << 1;
exp32--;
}
putdp(dp, frac64, exp32);
if (exp32 > 32767 || exp32 < -32768)
mathexception('f', FC_DFP_OFLOW, 0);
}
/* double->single floating point round (FRN) instruction.
Passed a pointer to a Prime double precision variable, one of the
FACC's, and updates it in place. May also take an arithmetic fault
if overflow occurs. For faults, the ea is zero because FACC's
don't have an effective address. */
void frn(void *dp) {
long long frac64;
int exp32;
int origsign, newsign;
CLEARC;
getdp(dp, &frac64, &exp32);
if (frac64 == 0)
*(long long *)dp = 0;
else {
origsign = (frac64 < 0);
if ((frac64 & 0x18000000000LL)
| ((frac64 & 0x8000000000LL) && (frac64 & 0x7FFFFF0000LL))) {
frac64 += 0x10000000000LL;
frac64 &= 0xFFFFFF0000000000LL;
norm(&frac64);
*(long long *)dp = frac64;
newsign = (frac64 < 0);
if (newsign != origsign)
/* XXX: is this fault code right? */
mathexception('f', FC_DFP_OFLOW, 0);
}
}
}
#if 0
/* Prime DPFP multiply */
dfmp(void *dp1, void *dp2, ea_t ea) {
long long frac64, frac641, frac642;
int exp32, exp321, exp322;
short fcode;
fcode = 0;
CLEARC;
getdp(dp1, &frac641, &exp321);
getdp(dp2, &frac642, &exp322);
exp32 = exp321 + exp322;
/* XXX: need to get 128-bit result to test for overflow? */
frac64 = frac641 * frac642;
if (exp32 > 32767 || exp32 < -32768)
fcode = FC_DFP_OFLOW;
/* insert (optional) rounding code here */
if (fcode == 0)
putdp(dp1, frac64, exp32);
else
mathexception('f', fcode, ea);
}
#endif