HP2100/hp2100_fp.c

   1 /* hp2100_fp.c: HP 2100 floating point instructions
   2
   3    Copyright (c) 2002-2008, Robert M. Supnik
   4
   5    Permission is hereby granted, free of charge, to any person obtaining a
   6    copy of this software and associated documentation files (the "Software"),
   7    to deal in the Software without restriction, including without limitation
   8    the rights to use, copy, modify, merge, publish, distribute, sublicense,
   9    and/or sell copies of the Software, and to permit persons to whom the
  10    Software is furnished to do so, subject to the following conditions:
  11
  12    The above copyright notice and this permission notice shall be included in
  13    all copies or substantial portions of the Software.
  14
  15    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  16    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  17    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  18    ROBERT M SUPNIK BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
  19    IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  20    CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21
  22    Except as contained in this notice, the name of Robert M Supnik shall not be
  23    used in advertising or otherwise to promote the sale, use or other dealings
  24    in this Software without prior written authorization from Robert M Supnik.
  25
  26    21-Jan-08    JDB     Corrected fp_unpack mantissa high-word return
  27                         (from Mark Pizzolato)
  28    01-Dec-06    JDB     Reworked FFP helpers for 1000-F support, deleted f_pwr2
  29    22-Jul-05    RMS     Fixed compiler warning in Solaris (from Doug Glyn)
  30    25-Feb-05    JDB     Added FFP helpers f_pack, f_unpack, f_pwr2
  31    11-Feb-05    JDB     Fixed missing negative overflow renorm in StoreFP
  32    26-Dec-04    RMS     Separated A/B from M[0/1] for DMA IO (from Dave Bryan)
  33    15-Jul-03    RMS     Fixed signed/unsigned warning
  34    21-Oct-02    RMS     Recoded for compatibility with 21MX microcode algorithms
  35
  36
  37    Implementation note: The 2100/1000-M/E Fast FORTRAN Processor (FFP) and 1000
  38    F-Series Floating Point Processor (FPP) simulations require that the host
  39    compiler support 64-bit integers and the HAVE_INT64 symbol be defined during
  40    compilation.  If this symbol is defined, two-word floating-point operations
  41    are handled in the FPP code, and this module is not used.  If it is not
  42    defined, then FFP and FPP operations are not available, and this module
  43    provides the floating-point support.
  44
  45
  46    The HP2100 uses a unique binary floating point format:
  47
  48      15 14                                         0
  49     +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
  50     |S |               fraction high                | : A
  51     +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
  52     |     fraction low      |      exponent      |XS| : A + 1
  53     +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
  54      15                    8  7                 1  0
  55
  56         where   S = 0 for plus fraction, 1 for minus fraction
  57                 fraction = s.bbbbb..., 24 binary digits
  58                 exponent = 2**+/-n
  59                 XS = 0 for plus exponent, 1 for minus exponent
  60
  61    Numbers can be normalized or unnormalized but are always normalized
  62    when loaded.
  63
  64    Unpacked floating point numbers are stored in structure ufp
  65
  66         exp     =       exponent, 2's complement
  67         h'l     =       fraction, 2's comp, left justified
  68
  69    This routine tries to reproduce the algorithms of the 2100/21MX
  70    microcode in order to achieve 'bug-for-bug' compatibility.  In
  71    particular,
  72
  73    - The FIX code produces various results in B.
  74    - The fraction multiply code uses 16b x 16b multiplies to produce
  75      a 31b result.  It always loses the low order bit of the product.
  76    - The fraction divide code is an approximation that may produce
  77      an error of 1 LSB.
  78    - Signs are tracked implicitly as part of the fraction.  Unnormalized
  79      inputs may cause the packup code to produce the wrong sign.
  80    - "Unclean" zeros (zero fraction, non-zero exponent) are processed
  81      like normal operands.
  82 */
  83
  84 #include "hp2100_defs.h"
  85 #include "hp2100_cpu1.h"
  86 #include "hp2100_fp.h"
  87
  88 #if !defined (HAVE_INT64)                               /* int64 support unavailable */
  89
  90 struct ufp {                                            /* unpacked fp */
  91     int32       exp;                                    /* exp */
  92     uint32      fr;                                     /* frac */
  93     };
  94
  95 #define FP_V_SIGN       31                              /* sign */
  96 #define FP_M_SIGN       01
  97 #define FP_V_FR         8                               /* fraction */
  98 #define FP_M_FR         077777777
  99 #define FP_V_EXP        1                               /* exponent */
 100 #define FP_M_EXP        0177
 101 #define FP_V_EXPS       0                               /* exp sign */
 102 #define FP_M_EXPS       01
 103 #define FP_SIGN         (FP_M_SIGN << FP_V_SIGN)
 104 #define FP_FR           (FP_M_FR << FP_V_FR)
 105 #define FP_EXP          (FP_M_EXP << FP_V_EXP)
 106 #define FP_EXPS         (FP_M_EXPS << FP_V_EXPS)
 107 #define FP_GETSIGN(x)   (((x) >> FP_V_SIGN) & FP_M_SIGN)
 108 #define FP_GETEXP(x)    (((x) >> FP_V_EXP) & FP_M_EXP)
 109 #define FP_GETEXPS(x)   (((x) >> FP_V_EXPS) & FP_M_EXPS)
 110
 111 #define FP_NORM         (1 << (FP_V_SIGN - 1))          /* normalized */
 112 #define FP_LOW          (1 << FP_V_FR)
 113 #define FP_RNDP         (1 << (FP_V_FR - 1))            /* round for plus */
 114 #define FP_RNDM         (FP_RNDP - 1)                   /* round for minus */
 115
 116 #define FPAB            ((((uint32) AR) << 16) | ((uint32) BR))
 117
 118 /* Fraction shift; 0 < shift < 32 */
 119
 120 #define FR_ARS(v,s)     (((v) >> (s)) | (((v) & FP_SIGN)? \
 121                         (((uint32) DMASK32) << (32 - (s))): 0)) & DMASK32
 122
 123 #define FR_NEG(v)       ((~(v) + 1) & DMASK32)
 124
 125 extern uint16 ABREG[2];
 126
 127 uint32 UnpackFP (struct ufp *fop, uint32 opnd);
 128 void NegFP (struct ufp *fop);
 129 void NormFP (struct ufp *fop);
 130 uint32 PackFP (struct ufp *fop);
 131 uint32 StoreFP (struct ufp *fop);
 132
 133 /* Floating to integer conversion */
 134
 135 uint32 f_fix (void)
 136 {
 137 struct ufp fop;
 138 uint32 res = 0;
 139
 140 UnpackFP (&fop, FPAB);                                  /* unpack op */
 141 if (fop.exp < 0) {                                      /* exp < 0? */
 142     AR = 0;                                             /* result = 0 */
 143     return 0;                                           /* B unchanged */
 144     }
 145 if (fop.exp > 15) {                                     /* exp > 15? */
 146     BR = AR;                                            /* B has high bits */
 147     AR = 077777;                                        /* result = 77777 */
 148     return 1;                                           /* overflow */
 149     }
 150 if (fop.exp < 15) {                                     /* if not aligned */
 151     res = FR_ARS (fop.fr, 15 - fop.exp);                /* shift right */
 152     AR = (res >> 16) & DMASK;                           /* AR gets result */
 153     }
 154 BR = AR;
 155 if ((AR & SIGN) && ((fop.fr | res) & DMASK))            /* any low bits lost? */
 156     AR = (AR + 1) & DMASK;                              /* round up */
 157 return 0;
 158 }
 159
 160 /* Integer to floating conversion */
 161
 162 uint32 f_flt (void)
 163 {
 164 struct ufp res = { 15, 0 };                             /* +, 2**15 */
 165
 166 res.fr = ((uint32) AR) << 16;                           /* left justify */
 167 StoreFP (&res);                                         /* store result */
 168 return 0;                                               /* clr overflow */
 169 }
 170
 171 /* Floating point add/subtract */
 172
 173 uint32 f_as (uint32 opnd, t_bool sub)
 174 {
 175 struct ufp fop1, fop2, t;
 176 int32 ediff;
 177
 178 UnpackFP (&fop1, FPAB);                                 /* unpack A-B */
 179 UnpackFP (&fop2, opnd);                                 /* get op */
 180 if (sub) {                                              /* subtract? */
 181     fop2.fr = FR_NEG (fop2.fr);                         /* negate frac */
 182     if (fop2.fr == ((uint32) FP_SIGN)) {                /* -1/2? */
 183         fop2.fr = fop2.fr >> 1;                         /* special case */
 184         fop2.exp = fop2.exp + 1;
 185         }
 186     }
 187 if (fop1.fr == 0) fop1 = fop2;                          /* op1 = 0? res = op2 */
 188 else if (fop2.fr != 0) {                                /* op2 = 0? no add */
 189     if (fop1.exp < fop2.exp) {                          /* |op1| < |op2|? */
 190         t = fop2;                                       /* swap operands */
 191         fop2 = fop1;
 192         fop1 = t;
 193         }
 194     ediff = fop1.exp - fop2.exp;                        /* get exp diff */
 195     if (ediff <= 24) {
 196         if (ediff) fop2.fr = FR_ARS (fop2.fr, ediff);   /* denorm, signed */
 197         if ((fop1.fr ^ fop2.fr) & FP_SIGN)              /* unlike signs? */
 198             fop1.fr = fop1.fr + fop2.fr;                /* eff subtract */
 199         else {                                          /* like signs */
 200             fop1.fr = fop1.fr + fop2.fr;                /* eff add */
 201             if (fop2.fr & FP_SIGN) {                    /* both -? */
 202                 if ((fop1.fr & FP_SIGN) == 0) {         /* overflow? */
 203                     fop1.fr = FP_SIGN | (fop1.fr >> 1); /* renormalize */
 204                     fop1.exp = fop1.exp + 1;            /* incr exp */
 205                     }
 206                 }
 207             else if (fop1.fr & FP_SIGN) {               /* both +, cry out? */
 208                 fop1.fr = fop1.fr >> 1;                 /* renormalize */
 209                 fop1.exp = fop1.exp + 1;                /* incr exp */
 210                 }
 211             }                                           /* end else like */
 212         }                                               /* end if ediff */
 213     }                                                   /* end if fop2 */
 214 return StoreFP (&fop1);                                 /* store result */
 215 }
 216
 217 /* Floating point multiply - passes diagnostic */
 218
 219 uint32 f_mul (uint32 opnd)
 220 {
 221 struct ufp fop1, fop2;
 222 struct ufp res = { 0, 0 };
 223 int32 shi1, shi2, t1, t2, t3, t4, t5;
 224
 225 UnpackFP (&fop1, FPAB);                                 /* unpack A-B */
 226 UnpackFP (&fop2, opnd);                                 /* unpack op */
 227 if (fop1.fr && fop2.fr) {                               /* if both != 0 */
 228     res.exp = fop1.exp + fop2.exp + 1;                  /* exp = sum */
 229     shi1 = SEXT (fop1.fr >> 16);                        /* mpy hi */
 230     shi2 = SEXT (fop2.fr >> 16);                        /* mpc hi */
 231     t1 = shi2 * ((int32) ((fop1.fr >> 1) & 077600));    /* mpc hi * (mpy lo/2) */
 232     t2 = shi1 * ((int32) ((fop2.fr >> 1) & 077600));    /* mpc lo * (mpy hi/2) */
 233     t3 = t1 + t2;                                       /* cross product */
 234     t4 = (shi1 * shi2) & ~1;                            /* mpy hi * mpc hi */
 235     t5 = (SEXT (t3 >> 16)) << 1;                        /* add in cross */
 236     res.fr = (t4 + t5) & DMASK32;                       /* bit<0> is lost */
 237     }
 238 return StoreFP (&res);                                  /* store */
 239 }
 240
 241 /* Floating point divide - reverse engineered from diagnostic */
 242
 243 uint32 divx (uint32 ba, uint32 dvr, uint32 *rem)
 244 {
 245 int32 sdvd = 0, sdvr = 0;
 246 uint32 q, r;
 247
 248 if (ba & FP_SIGN) sdvd = 1;                             /* 32b/16b signed dvd */
 249 if (dvr & SIGN) sdvr = 1;                               /* use old-fashioned */
 250 if (sdvd) ba = (~ba + 1) & DMASK32;                     /* unsigned divides, */
 251 if (sdvr) dvr = (~dvr + 1) & DMASK;                     /* as results may ovflo */
 252 q = ba / dvr;
 253 r = ba % dvr;
 254 if (sdvd ^ sdvr) q = (~q + 1) & DMASK;
 255 if (sdvd) r = (~r + 1) & DMASK;
 256 if (rem) *rem = r;
 257 return q;
 258 }
 259
 260 uint32 f_div (uint32 opnd)
 261 {
 262 struct ufp fop1, fop2;
 263 struct ufp quo = { 0, 0 };
 264 uint32 ba, q0, q1, q2, dvrh;
 265
 266 UnpackFP (&fop1, FPAB);                                 /* unpack A-B */
 267 UnpackFP (&fop2, opnd);                                 /* unpack op */
 268 dvrh = (fop2.fr >> 16) & DMASK;                         /* high divisor */
 269 if (dvrh == 0) {                                        /* div by zero? */
 270     AR = 0077777;                                       /* return most pos */
 271     BR = 0177776;
 272     return 1;
 273     }
 274 if (fop1.fr) {                                          /* dvd != 0? */
 275     quo.exp = fop1.exp - fop2.exp + 1;                  /* exp = diff */
 276     ba = FR_ARS (fop1.fr, 2);                           /* prevent ovflo */
 277     q0 = divx (ba, dvrh, &ba);                          /* Q0 = dvd / dvrh */
 278     ba = (ba & ~1) << 16;                               /* remainder */
 279     ba = FR_ARS (ba, 1);                                /* prevent ovflo */
 280     q1 = divx (ba, dvrh, NULL);                         /* Q1 = rem / dvrh */
 281     ba = (fop2.fr & 0xFF00) << 13;                      /* dvrl / 8 */
 282     q2 = divx (ba, dvrh, NULL);                         /* dvrl / dvrh */
 283     ba = -(SEXT (q2)) * (SEXT (q0));                    /* -Q0 * Q2 */
 284     ba = (ba >> 16) & 0xFFFF;                           /* save ms half */
 285     if (q1 & SIGN) quo.fr = quo.fr - 0x00010000;        /* Q1 < 0? -1 */
 286     if (ba & SIGN) quo.fr = quo.fr - 0x00010000;        /* -Q0*Q2 < 0? */
 287     quo.fr = quo.fr + ((ba << 2) & 0xFFFF) + q1;        /* rest prod, add Q1 */
 288     quo.fr = quo.fr << 1;                               /* shift result */
 289     quo.fr = quo.fr + (q0 << 16);                       /* add Q0 */
 290     }                                                   /* end if fop1.h */
 291 return StoreFP (&quo);                                  /* store result */
 292 }
 293
 294 /* Utility routines */
 295
 296 /* Unpack operand */
 297
 298 uint32 UnpackFP (struct ufp *fop, uint32 opnd)
 299 {
 300 fop->fr = opnd & FP_FR;                                 /* get frac */
 301 fop->exp = FP_GETEXP (opnd);                            /* get exp */
 302 if (FP_GETEXPS (opnd)) fop->exp = fop->exp | ~FP_M_EXP; /* < 0? sext */
 303 return FP_GETSIGN (opnd);                               /* return sign */
 304 }
 305
 306 /* Normalize unpacked floating point number */
 307
 308 void NormFP (struct ufp *fop)
 309 {
 310 if (fop->fr) {                                          /* any fraction? */
 311     uint32 test = (fop->fr >> 1) & FP_NORM;
 312     while ((fop->fr & FP_NORM) == test) {               /* until norm */
 313         fop->exp = fop->exp - 1;
 314         fop->fr = (fop->fr << 1);
 315         }
 316     }
 317 else fop->exp = 0;                                      /* clean 0 */
 318 return;
 319 }
 320
 321 /* Pack fp number */
 322
 323 uint32 PackFP (struct ufp *fop)
 324 {
 325 return (fop->fr & FP_FR) |                              /* merge frac */
 326        ((fop->exp & FP_M_EXP) << FP_V_EXP) |            /* and exp */
 327        ((fop->exp < 0)? (1 << FP_V_EXPS): 0);           /* add exp sign */
 328 }
 329
 330 /* Round fp number, store, generate overflow */
 331
 332 uint32 StoreFP (struct ufp *fop)
 333 {
 334 uint32 sign, svfr, hi, ov = 0;
 335
 336 NormFP (fop);                                           /* normalize */
 337 svfr = fop->fr;                                         /* save fraction */
 338 sign = FP_GETSIGN (fop->fr);                            /* save sign */
 339 fop->fr = (fop->fr + (sign? FP_RNDM: FP_RNDP)) & FP_FR; /* round */
 340 if ((fop->fr ^ svfr) & FP_SIGN) {                       /* sign change? */
 341     fop->fr = fop->fr >> 1;                             /* renormalize */
 342     fop->exp = fop->exp + 1;
 343     }
 344 else NormFP (fop);                                      /* check for norm */
 345 if (fop->fr == 0) hi = 0;                               /* result 0? */
 346 else if (fop->exp < -(FP_M_EXP + 1)) {                  /* underflow? */
 347     hi = 0;                                             /* store clean 0 */
 348     ov = 1;
 349     }
 350 else if (fop->exp > FP_M_EXP) {                         /* overflow? */
 351     hi = 0x7FFFFFFE;                                    /* all 1's */
 352     ov = 1;
 353     }
 354 else hi = PackFP (fop);                                 /* pack mant and exp */
 355 AR = (hi >> 16) & DMASK;
 356 BR = hi & DMASK;
 357 return ov;
 358 }
 359
 360
 361 /* Single-precision Fast FORTRAN Processor helpers. */
 362
 363 /* Pack mantissa and exponent and return fp value. */
 364
 365 uint32 fp_pack (OP *result, OP mantissa, int32 exponent, OPSIZE precision)
 366 {
 367 struct ufp fop;
 368 uint32 val;
 369
 370 fop.fr = ((uint32) mantissa.fpk[0] << 16) | mantissa.fpk[1];
 371 fop.exp = exponent;
 372 val = PackFP (&fop);
 373 result->fpk[0] = (int16) (val >> 16);
 374 result->fpk[1] = (int16) val;
 375 return 0;
 376 }
 377
 378 /* Normalize, round, and pack mantissa and exponent and return fp value. */
 379
 380 uint32 fp_nrpack (OP *result, OP mantissa, int32 exponent, OPSIZE precision)
 381 {
 382 struct ufp fop;
 383 uint32 ovf;
 384
 385 fop.fr = ((uint32) mantissa.fpk[0] << 16) | mantissa.fpk[1];
 386 fop.exp = exponent;
 387 ovf = StoreFP (&fop);
 388 result->fpk[0] = AR;
 389 result->fpk[1] = BR;
 390 return ovf;
 391 }
 392
 393 /* Unpack fp number in into mantissa and exponent. */
 394
 395 uint32 fp_unpack (OP *mantissa, int32 *exponent, OP packed, OPSIZE precision)
 396 {
 397 struct ufp fop;
 398 uint32 operand;
 399
 400 operand = ((uint32) packed.fpk[0] << 16) | packed.fpk[1];
 401 UnpackFP (&fop, operand);
 402 mantissa->fpk[0] = (uint16) (fop.fr >> 16);
 403 mantissa->fpk[1] = (uint16) fop.fr;
 404 *exponent = fop.exp;
 405 return 0;
 406 }
 407
 408 #endif                                                  /* int64 support unavailable */