blob: 35286ff95c12998f83a810e7b93755b4e426bc5b [file] [log] [blame] [raw]
/* hp2100_fp.c: HP 2100 floating point instructions
Copyright (c) 2002-2015, Robert M. Supnik
Copyright (c) 2017 J. David Bryan
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Except as contained in this notice, the names of the authors shall not be
used in advertising or otherwise to promote the sale, use or other dealings
in this Software without prior written authorization from the authors.
26-Jun-17 JDB Replaced SEXT with SEXT16
03-Jan-15 JDB Made the utility routines static
21-Jan-08 JDB Corrected fp_unpack mantissa high-word return
(from Mark Pizzolato)
01-Dec-06 JDB Reworked FFP helpers for 1000-F support, deleted f_pwr2
22-Jul-05 RMS Fixed compiler warning in Solaris (from Doug Glyn)
25-Feb-05 JDB Added FFP helpers f_pack, f_unpack, f_pwr2
11-Feb-05 JDB Fixed missing negative overflow renorm in StoreFP
26-Dec-04 RMS Separated A/B from M[0/1] for DMA IO (from Dave Bryan)
15-Jul-03 RMS Fixed signed/unsigned warning
21-Oct-02 RMS Recoded for compatibility with 21MX microcode algorithms
The HP2100 uses a unique binary floating point format:
15 14 0
+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
|S | fraction high | : A
+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
| fraction low | exponent |XS| : A + 1
+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
15 8 7 1 0
where S = 0 for plus fraction, 1 for minus fraction
fraction = s.bbbbb..., 24 binary digits
exponent = 2**+/-n
XS = 0 for plus exponent, 1 for minus exponent
Numbers can be normalized or unnormalized but are always normalized
when loaded.
Unpacked floating point numbers are stored in structure ufp
exp = exponent, 2's complement
h'l = fraction, 2's comp, left justified
This routine tries to reproduce the algorithms of the 2100/21MX
microcode in order to achieve 'bug-for-bug' compatibility. In
particular,
- The FIX code produces various results in B.
- The fraction multiply code uses 16b x 16b multiplies to produce
a 31b result. It always loses the low order bit of the product.
- The fraction divide code is an approximation that may produce
an error of 1 LSB.
- Signs are tracked implicitly as part of the fraction. Unnormalized
inputs may cause the packup code to produce the wrong sign.
- "Unclean" zeros (zero fraction, non-zero exponent) are processed
like normal operands.
Implementation notes:
1. The 2100/1000-M/E Fast FORTRAN Processor (FFP) and 1000 F-Series Floating
Point Processor (FPP) simulations require that the host compiler support
64-bit integers and the HAVE_INT64 symbol be defined during compilation.
If this symbol is defined, two-word floating-point operations are handled
in the FPP code, and this module is not used. If it is not defined, then
FFP and FPP operations are not available, and this module provides the
floating-point support.
*/
#include "hp2100_defs.h"
#include "hp2100_cpu.h"
#include "hp2100_cpu1.h"
#include "hp2100_fp.h"
#if !defined (HAVE_INT64) /* int64 support unavailable */
struct ufp { /* unpacked fp */
int32 exp; /* exp */
uint32 fr; /* frac */
};
#define FP_V_SIGN 31 /* sign */
#define FP_M_SIGN 01
#define FP_V_FR 8 /* fraction */
#define FP_M_FR 077777777
#define FP_V_EXP 1 /* exponent */
#define FP_M_EXP 0177
#define FP_V_EXPS 0 /* exp sign */
#define FP_M_EXPS 01
#define FP_SIGN (FP_M_SIGN << FP_V_SIGN)
#define FP_FR (FP_M_FR << FP_V_FR)
#define FP_EXP (FP_M_EXP << FP_V_EXP)
#define FP_EXPS (FP_M_EXPS << FP_V_EXPS)
#define FP_GETSIGN(x) (((x) >> FP_V_SIGN) & FP_M_SIGN)
#define FP_GETEXP(x) (((x) >> FP_V_EXP) & FP_M_EXP)
#define FP_GETEXPS(x) (((x) >> FP_V_EXPS) & FP_M_EXPS)
#define FP_NORM (1 << (FP_V_SIGN - 1)) /* normalized */
#define FP_LOW (1 << FP_V_FR)
#define FP_RNDP (1 << (FP_V_FR - 1)) /* round for plus */
#define FP_RNDM (FP_RNDP - 1) /* round for minus */
#define FPAB ((((uint32) AR) << 16) | ((uint32) BR))
/* Fraction shift; 0 < shift < 32 */
#define FR_ARS(v,s) (((v) >> (s)) | (((v) & FP_SIGN)? \
(((uint32) DMASK32) << (32 - (s))): 0)) & DMASK32
#define FR_NEG(v) ((~(v) + 1) & DMASK32)
/* Utility routines */
static uint32 UnpackFP (struct ufp *fop, uint32 opnd);
static void NormFP (struct ufp *fop);
static uint32 PackFP (struct ufp *fop);
static uint32 StoreFP (struct ufp *fop);
/* Floating to integer conversion */
uint32 f_fix (void)
{
struct ufp fop;
uint32 res = 0;
UnpackFP (&fop, FPAB); /* unpack op */
if (fop.exp < 0) { /* exp < 0? */
AR = 0; /* result = 0 */
return 0; /* B unchanged */
}
if (fop.exp > 15) { /* exp > 15? */
BR = AR; /* B has high bits */
AR = 077777; /* result = 77777 */
return 1; /* overflow */
}
if (fop.exp < 15) { /* if not aligned */
res = FR_ARS (fop.fr, 15 - fop.exp); /* shift right */
AR = (res >> 16) & DMASK; /* AR gets result */
}
BR = AR;
if ((AR & SIGN) && ((fop.fr | res) & DMASK)) /* any low bits lost? */
AR = (AR + 1) & DMASK; /* round up */
return 0;
}
/* Integer to floating conversion */
uint32 f_flt (void)
{
struct ufp res = { 15, 0 }; /* +, 2**15 */
res.fr = ((uint32) AR) << 16; /* left justify */
StoreFP (&res); /* store result */
return 0; /* clr overflow */
}
/* Floating point add/subtract */
uint32 f_as (uint32 opnd, t_bool sub)
{
struct ufp fop1, fop2, t;
int32 ediff;
UnpackFP (&fop1, FPAB); /* unpack A-B */
UnpackFP (&fop2, opnd); /* get op */
if (sub) { /* subtract? */
fop2.fr = FR_NEG (fop2.fr); /* negate frac */
if (fop2.fr == ((uint32) FP_SIGN)) { /* -1/2? */
fop2.fr = fop2.fr >> 1; /* special case */
fop2.exp = fop2.exp + 1;
}
}
if (fop1.fr == 0) fop1 = fop2; /* op1 = 0? res = op2 */
else if (fop2.fr != 0) { /* op2 = 0? no add */
if (fop1.exp < fop2.exp) { /* |op1| < |op2|? */
t = fop2; /* swap operands */
fop2 = fop1;
fop1 = t;
}
ediff = fop1.exp - fop2.exp; /* get exp diff */
if (ediff <= 24) {
if (ediff) fop2.fr = FR_ARS (fop2.fr, ediff); /* denorm, signed */
if ((fop1.fr ^ fop2.fr) & FP_SIGN) /* unlike signs? */
fop1.fr = fop1.fr + fop2.fr; /* eff subtract */
else { /* like signs */
fop1.fr = fop1.fr + fop2.fr; /* eff add */
if (fop2.fr & FP_SIGN) { /* both -? */
if ((fop1.fr & FP_SIGN) == 0) { /* overflow? */
fop1.fr = FP_SIGN | (fop1.fr >> 1); /* renormalize */
fop1.exp = fop1.exp + 1; /* incr exp */
}
}
else if (fop1.fr & FP_SIGN) { /* both +, cry out? */
fop1.fr = fop1.fr >> 1; /* renormalize */
fop1.exp = fop1.exp + 1; /* incr exp */
}
} /* end else like */
} /* end if ediff */
} /* end if fop2 */
return StoreFP (&fop1); /* store result */
}
/* Floating point multiply - passes diagnostic */
uint32 f_mul (uint32 opnd)
{
struct ufp fop1, fop2;
struct ufp res = { 0, 0 };
int32 shi1, shi2, t1, t2, t3, t4, t5;
UnpackFP (&fop1, FPAB); /* unpack A-B */
UnpackFP (&fop2, opnd); /* unpack op */
if (fop1.fr && fop2.fr) { /* if both != 0 */
res.exp = fop1.exp + fop2.exp + 1; /* exp = sum */
shi1 = SEXT16 (UPPER_WORD (fop1.fr)); /* mpy hi */
shi2 = SEXT16 (UPPER_WORD (fop2.fr)); /* mpc hi */
t1 = shi2 * ((int32) ((fop1.fr >> 1) & 077600)); /* mpc hi * (mpy lo/2) */
t2 = shi1 * ((int32) ((fop2.fr >> 1) & 077600)); /* mpc lo * (mpy hi/2) */
t3 = t1 + t2; /* cross product */
t4 = (shi1 * shi2) & ~1; /* mpy hi * mpc hi */
t5 = SEXT16 (UPPER_WORD (t3)) << 1; /* add in cross */
res.fr = (t4 + t5) & DMASK32; /* bit<0> is lost */
}
return StoreFP (&res); /* store */
}
/* Floating point divide - reverse engineered from diagnostic */
static uint32 divx (uint32 ba, uint32 dvr, uint32 *rem)
{
int32 sdvd = 0, sdvr = 0;
uint32 q, r;
if (ba & FP_SIGN) sdvd = 1; /* 32b/16b signed dvd */
if (dvr & SIGN) sdvr = 1; /* use old-fashioned */
if (sdvd) ba = (~ba + 1) & DMASK32; /* unsigned divides, */
if (sdvr) dvr = (~dvr + 1) & DMASK; /* as results may ovflo */
q = ba / dvr;
r = ba % dvr;
if (sdvd ^ sdvr) q = (~q + 1) & DMASK;
if (sdvd) r = (~r + 1) & DMASK;
if (rem) *rem = r;
return q;
}
uint32 f_div (uint32 opnd)
{
struct ufp fop1, fop2;
struct ufp quo = { 0, 0 };
uint32 ba, q0, q1, q2, dvrh;
UnpackFP (&fop1, FPAB); /* unpack A-B */
UnpackFP (&fop2, opnd); /* unpack op */
dvrh = (fop2.fr >> 16) & DMASK; /* high divisor */
if (dvrh == 0) { /* div by zero? */
AR = 0077777; /* return most pos */
BR = 0177776;
return 1;
}
if (fop1.fr) { /* dvd != 0? */
quo.exp = fop1.exp - fop2.exp + 1; /* exp = diff */
ba = FR_ARS (fop1.fr, 2); /* prevent ovflo */
q0 = divx (ba, dvrh, &ba); /* Q0 = dvd / dvrh */
ba = (ba & ~1) << 16; /* remainder */
ba = FR_ARS (ba, 1); /* prevent ovflo */
q1 = divx (ba, dvrh, NULL); /* Q1 = rem / dvrh */
ba = (fop2.fr & 0xFF00) << 13; /* dvrl / 8 */
q2 = divx (ba, dvrh, NULL); /* dvrl / dvrh */
ba = - SEXT16 (LOWER_WORD (q2)) * SEXT16 (LOWER_WORD (q0)); /* -Q0 * Q2 */
ba = (ba >> 16) & 0xFFFF; /* save ms half */
if (q1 & SIGN) quo.fr = quo.fr - 0x00010000; /* Q1 < 0? -1 */
if (ba & SIGN) quo.fr = quo.fr - 0x00010000; /* -Q0*Q2 < 0? */
quo.fr = quo.fr + ((ba << 2) & 0xFFFF) + q1; /* rest prod, add Q1 */
quo.fr = quo.fr << 1; /* shift result */
quo.fr = quo.fr + (q0 << 16); /* add Q0 */
} /* end if fop1.h */
return StoreFP (&quo); /* store result */
}
/* Utility routines */
/* Unpack operand */
static uint32 UnpackFP (struct ufp *fop, uint32 opnd)
{
fop->fr = opnd & FP_FR; /* get frac */
fop->exp = FP_GETEXP (opnd); /* get exp */
if (FP_GETEXPS (opnd)) fop->exp = fop->exp | ~FP_M_EXP; /* < 0? sext */
return FP_GETSIGN (opnd); /* return sign */
}
/* Normalize unpacked floating point number */
static void NormFP (struct ufp *fop)
{
if (fop->fr) { /* any fraction? */
uint32 test = (fop->fr >> 1) & FP_NORM;
while ((fop->fr & FP_NORM) == test) { /* until norm */
fop->exp = fop->exp - 1;
fop->fr = (fop->fr << 1);
}
}
else fop->exp = 0; /* clean 0 */
return;
}
/* Pack fp number */
static uint32 PackFP (struct ufp *fop)
{
return (fop->fr & FP_FR) | /* merge frac */
((fop->exp & FP_M_EXP) << FP_V_EXP) | /* and exp */
((fop->exp < 0)? (1 << FP_V_EXPS): 0); /* add exp sign */
}
/* Round fp number, store, generate overflow */
static uint32 StoreFP (struct ufp *fop)
{
uint32 sign, svfr, hi, ov = 0;
NormFP (fop); /* normalize */
svfr = fop->fr; /* save fraction */
sign = FP_GETSIGN (fop->fr); /* save sign */
fop->fr = (fop->fr + (sign? FP_RNDM: FP_RNDP)) & FP_FR; /* round */
if ((fop->fr ^ svfr) & FP_SIGN) { /* sign change? */
fop->fr = fop->fr >> 1; /* renormalize */
fop->exp = fop->exp + 1;
}
else NormFP (fop); /* check for norm */
if (fop->fr == 0) hi = 0; /* result 0? */
else if (fop->exp < -(FP_M_EXP + 1)) { /* underflow? */
hi = 0; /* store clean 0 */
ov = 1;
}
else if (fop->exp > FP_M_EXP) { /* overflow? */
hi = 0x7FFFFFFE; /* all 1's */
ov = 1;
}
else hi = PackFP (fop); /* pack mant and exp */
AR = (hi >> 16) & DMASK;
BR = hi & DMASK;
return ov;
}
/* Single-precision Fast FORTRAN Processor helpers. */
/* Pack mantissa and exponent and return fp value. */
uint32 fp_pack (OP *result, OP mantissa, int32 exponent, OPSIZE precision)
{
struct ufp fop;
uint32 val;
fop.fr = ((uint32) mantissa.fpk[0] << 16) | mantissa.fpk[1];
fop.exp = exponent;
val = PackFP (&fop);
result->fpk[0] = UPPER_WORD (val);
result->fpk[1] = LOWER_WORD (val);
return 0;
}
/* Normalize, round, and pack mantissa and exponent and return fp value. */
uint32 fp_nrpack (OP *result, OP mantissa, int32 exponent, OPSIZE precision)
{
struct ufp fop;
uint32 ovf;
fop.fr = ((uint32) mantissa.fpk[0] << 16) | mantissa.fpk[1];
fop.exp = exponent;
ovf = StoreFP (&fop);
result->fpk[0] = AR;
result->fpk[1] = BR;
return ovf;
}
/* Unpack fp number in into mantissa and exponent. */
uint32 fp_unpack (OP *mantissa, int32 *exponent, OP packed, OPSIZE precision)
{
struct ufp fop;
uint32 operand;
operand = ((uint32) packed.fpk[0] << 16) | packed.fpk[1];
UnpackFP (&fop, operand);
mantissa->fpk[0] = UPPER_WORD (fop.fr);
mantissa->fpk[1] = LOWER_WORD (fop.fr);
*exponent = fop.exp;
return 0;
}
#endif /* int64 support unavailable */