[simh.git] / PDP10 / pdp10_mdfp.c

/* pdp10_mdfp.c: PDP-10 multiply/divide and floating point simulator\r
\r
   Copyright (c) 1993-2005, Robert M Supnik\r
\r
   Permission is hereby granted, free of charge, to any person obtaining a\r
   copy of this software and associated documentation files (the "Software"),\r
   to deal in the Software without restriction, including without limitation\r
   the rights to use, copy, modify, merge, publish, distribute, sublicense,\r
   and/or sell copies of the Software, and to permit persons to whom the\r
   Software is furnished to do so, subject to the following conditions:\r
\r
   The above copyright notice and this permission notice shall be included in\r
   all copies or substantial portions of the Software.\r
\r
   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\r
   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\r
   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL\r
   ROBERT M SUPNIK BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER\r
   IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN\r
   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\r
\r
   Except as contained in this notice, the name of Robert M Supnik shall not be\r
   used in advertising or otherwise to promote the sale, use or other dealings\r
   in this Software without prior written authorization from Robert M Supnik.\r
\r
   2-Apr-04     RMS     Fixed bug in floating point unpack\r
                        Fixed bug in FIXR (found by Phil Stone, fixed by\r
                        Chris Smith)\r
   31-Aug-01    RMS     Changed int64 to t_int64 for Windoze\r
   10-Aug-01    RMS     Removed register in declarations\r
\r
   Instructions handled in this module:\r
        imul            integer multiply\r
        idiv            integer divide\r
        mul             multiply\r
        div             divide\r
        dmul            double precision multiply\r
        ddiv            double precision divide\r
        fad(r)          floating add (and round)\r
        fsb(r)          floating subtract (and round)\r
        fmp(r)          floating multiply (and round)\r
        fdv(r)          floating divide and round\r
        fsc             floating scale\r
        fix(r)          floating to fixed (and round)\r
        fltr            fixed to floating and round\r
        dfad            double precision floating add/subtract\r
        dfmp            double precision floating multiply\r
        dfdv            double precision floating divide\r
\r
   The PDP-10 stores double (quad) precision integers in sequential\r
   AC's or memory locations.  Integers are stored in 2's complement\r
   form.  Only the sign of the high order word matters; the signs\r
   in low order words are ignored on input and set to the sign of\r
   the result on output.  Quad precision integers exist only in the\r
   AC's as the result of a DMUL or the dividend of a DDIV.\r
\r
    0 00000000011111111112222222222333333\r
    0 12345678901234567890123456789012345\r
   +-+-----------------------------------+\r
   |S|      high order integer           | AC(n), A\r
   +-+-----------------------------------+\r
   |S|      low order integer            | AC(n + 1), A + 1\r
   +-+-----------------------------------+\r
   |S|      low order integer            | AC(n + 2)\r
   +-+-----------------------------------+\r
   |S|      low order integer            | AC(n + 3)\r
   +-+-----------------------------------+\r
\r
   The PDP-10 supports two floating point formats: single and double\r
   precision.  In both, the exponent is 8 bits, stored in excess\r
   128 notation.  The fraction is expected to be normalized.  A\r
   single precision floating point number has 27 bits of fraction;\r
   a double precision number has 62 bits of fraction (the sign\r
   bit of the second word is ignored and is set to zero).\r
\r
   In a negative floating point number, the exponent is stored in\r
   one's complement form, the fraction in two's complement form.\r
\r
    0 00000000 011111111112222222222333333\r
    0 12345678 901234567890123456789012345\r
   +-+--------+---------------------------+\r
   |S|exponent|      high order fraction  | AC(n), A\r
   +-+--------+---------------------------+\r
   |0|      low order fraction            | AC(n + 1), A + 1\r
   +-+------------------------------------+\r
\r
   Note that treatment of the sign is different for double precision\r
   integers and double precision floating point.  DMOVN (implemented\r
   as an inline macro) follows floating point conventions.\r
\r
   The original PDP-10 CPU (KA10) used a different format for double\r
   precision numbers and included certain instructions to make\r
   software support easier.  These instructions were phased out in\r
   the KL10 and KS10 and are treated as MUUO's.\r
\r
   The KL10 added extended precision (11-bit exponent) floating point\r
   format (so-called G floating).  These instructions were not\r
   implemented in the KS10 and are treated as MUUO's.\r
*/\r
\r
#include "pdp10_defs.h"\r
#include <setjmp.h>\r
\r
typedef struct {                                        /* unpacked fp number */\r
    int32               sign;                           /* sign */\r
    int32               exp;                            /* exponent */\r
    t_uint64    fhi;                                    /* fraction high */\r
    t_uint64    flo;                                    /* for double prec */\r
    } UFP;\r
\r
#define MSK32           0xFFFFFFFF\r
#define FIT27           (DMASK - 0x07FFFFFF)\r
#define FIT32           (DMASK - MSK32)\r
#define SFRC            TRUE                            /* frac 2's comp */\r
#define AFRC            FALSE                           /* frac abs value */\r
\r
/* In packed floating point number */\r
\r
#define FP_BIAS         0200                            /* exponent bias */\r
#define FP_N_FHI        27                              /* # of hi frac bits */\r
#define FP_V_FHI        0                               /* must be zero */\r
#define FP_M_FHI        0000777777777\r
#define FP_N_EXP        8                               /* # of exp bits */\r
#define FP_V_EXP        (FP_V_FHI + FP_N_FHI)\r
#define FP_M_EXP        0377\r
#define FP_V_SIGN       (FP_V_EXP + FP_N_EXP)           /* sign */\r
#define FP_N_FLO        35                              /* # of lo frac bits */\r
#define FP_V_FLO        0                               /* must be zero */\r
#define FP_M_FLO        0377777777777\r
#define GET_FPSIGN(x)   ((int32) (((x) >> FP_V_SIGN) & 1))\r
#define GET_FPEXP(x)    ((int32) (((x) >> FP_V_EXP) & FP_M_EXP))\r
#define GET_FPHI(x)     ((x) & FP_M_FHI)\r
#define GET_FPLO(x)     ((x) & FP_M_FLO)\r
\r
/* In unpacked floating point number */\r
\r
#define FP_N_GUARD      1                               /* # of guard bits */\r
#define FP_V_UFLO       FP_N_GUARD                      /* <35:1> */\r
#define FP_V_URNDD      (FP_V_UFLO - 1)                 /* dp round bit */\r
#define FP_V_UFHI       (FP_V_UFLO + FP_N_FLO)          /* <62:36> */\r
#define FP_V_URNDS      (FP_V_UFHI - 1)                 /* sp round bit */\r
#define FP_V_UCRY       (FP_V_UFHI + FP_N_FHI)          /* <63> */\r
#define FP_V_UNORM      (FP_V_UCRY - 1)                 /* normalized bit */\r
#define FP_UFHI         0x7FFFFFF000000000\r
#define FP_UFLO         0x0000000FFFFFFFFE\r
#define FP_UFRAC        0x7FFFFFFFFFFFFFFE\r
#define FP_URNDD        0x0000000000000001\r
#define FP_URNDS        0x0000000800000000\r
#define FP_UNORM        0x4000000000000000\r
#define FP_UCRY         0x8000000000000000\r
#define FP_ONES         0xFFFFFFFFFFFFFFFF\r
\r
#define UNEG(x)         ((~x) + 1)\r
#define DUNEG(x)        x.flo = UNEG (x.flo); x.fhi = ~x.fhi + (x.flo == 0)\r
\r
extern d10 *ac_cur;                                     /* current AC block */\r
extern int32 flags;                                     /* flags */\r
void mul (d10 a, d10 b, d10 *rs);\r
void funpack (d10 h, d10 l, UFP *r, t_bool sgn);\r
void fnorm (UFP *r, t_int64 rnd);\r
d10 fpack (UFP *r, d10 *lo, t_bool fdvneg);\r
\r
/* Integer multiply - checked against KS-10 ucode */\r
\r
d10 imul (d10 a, d10 b)\r
{\r
d10 rs[2];\r
\r
if ((a == SIGN) && (b == SIGN)) {                       /* KS10 hack */\r
    SETF (F_AOV | F_T1);                                /* -2**35 squared */\r
    return SIGN;\r
    }\r
mul (a, b, rs);                                         /* mpy, dprec result */\r
if (rs[0] && (rs[0] != ONES)) {                         /* high not all sign? */\r
    rs[1] = TSTS (a ^ b)? SETS (rs[1]): CLRS (rs[1]);   /* set sign */\r
    SETF (F_AOV | F_T1);                                /* overflow */\r
    }\r
return rs[1];\r
}\r
\r
/* Integer divide, return quotient, remainder  - checked against KS10 ucode\r
   The KS10 does not recognize -2^35/-1 as an error.  Instead, it produces\r
   2^35 (that is, -2^35) as the incorrect result.\r
*/\r
\r
t_bool idiv (d10 a, d10 b, d10 *rs)\r
{\r
d10 dvd = ABS (a);                                      /* make ops positive */\r
d10 dvr = ABS (b);\r
\r
if (dvr == 0) {                                         /* divide by 0? */\r
    SETF (F_DCK | F_AOV | F_T1);                        /* set flags, return */\r
    return FALSE;\r
    }\r
rs[0] = dvd / dvr;                                      /* get quotient */\r
rs[1] = dvd % dvr;                                      /* get remainder */\r
if (TSTS (a ^ b)) rs[0] = NEG (rs[0]);                  /* sign of result */\r
if (TSTS (a)) rs[1] = NEG (rs[1]);                      /* sign of remainder */\r
return TRUE;\r
}\r
\r
/* Multiply, return double precision result - checked against KS10 ucode */\r
\r
void mul (d10 s1, d10 s2, d10 *rs)\r
{\r
t_uint64 a = ABS (s1);\r
t_uint64 b = ABS (s2);\r
t_uint64 t, u, r;\r
\r
if ((a == 0) || (b == 0)) {                             /* operand = 0? */\r
    rs[0] = rs[1] = 0;                                  /* result 0 */\r
    return;\r
    }\r
if ((a & FIT32) || (b & FIT32)) {                       /* fit in 64b? */\r
    t = a >> 18;                                        /* no, split in half */\r
    a = a & RMASK;                                      /* "dp" multiply */\r
    u = b >> 18;\r
    b = b & RMASK;\r
    r = (a * b) + (((a * u) + (b * t)) << 18);          /* low is only 35b */\r
    rs[0] = ((t * u) << 1) + (r >> 35);                 /* so lsh hi 1 */\r
    rs[1] = r & MMASK;\r
    }\r
else {\r
    r = a * b;                                          /* fits, native mpy */\r
    rs[0] = r >> 35;                                    /* split at bit 35 */\r
    rs[1] = r & MMASK;\r
    }\r
\r
if (TSTS (s1 ^ s2)) { MKDNEG (rs); }                    /* result -? */\r
else if (TSTS (rs[0])) {                                /* result +, 2**70? */\r
    SETF (F_AOV | F_T1);                                /* overflow */\r
    rs[1] = SETS (rs[1]);                               /* consistent - */\r
    }\r
return;\r
}\r
\r
/* Divide, return quotient and remainder - checked against KS10 ucode\r
   Note that the initial divide check catches the case -2^70/-2^35;\r
   thus, the quotient can have at most 35 bits.\r
*/\r
\r
t_bool divi (int32 ac, d10 b, d10 *rs)\r
{\r
int32 p1 = ADDAC (ac, 1);\r
d10 dvr = ABS (b);                                      /* make divr positive */\r
t_int64 t;\r
int32 i;\r
d10 dvd[2];\r
\r
dvd[0] = AC(ac);                                        /* divd high */\r
dvd[1] = CLRS (AC(p1));                                 /* divd lo, clr sgn */\r
if (TSTS (AC(ac))) { DMOVN (dvd); }                     /* make divd positive */\r
if (dvd[0] >= dvr) {                                    /* divide fail? */\r
    SETF (F_AOV | F_DCK | F_T1);                        /* set flags, return */\r
    return FALSE;\r
    }\r
if (dvd[0] & FIT27) {                                   /* fit in 63b? */\r
    for (i = 0, rs[0] = 0; i < 35; i++) {               /* 35 quotient bits */\r
        dvd[0] = (dvd[0] << 1) | ((dvd[1] >> 34) & 1);\r
        dvd[1] = (dvd[1] << 1) & MMASK;                 /* shift dividend */\r
        rs[0] = rs[0] << 1;                             /* shift quotient */\r
        if (dvd[0] >= dvr) {                            /* subtract work? */\r
            dvd[0] = dvd[0] - dvr;                      /* quo bit is 1 */\r
            rs[0] = rs[0] + 1;\r
            }\r
        }\r
    rs[1] = dvd[0];                                     /* store remainder */\r
    }\r
else {\r
    t = (dvd[0] << 35) | dvd[1];                        /* concatenate */\r
    rs[0] = t / dvr;                                    /* quotient */\r
    rs[1] = t % dvr;                                    /* remainder */\r
    }\r
if (TSTS (AC(ac) ^ b)) rs[0] = NEG (rs[0]);             /* sign of result */\r
if (TSTS (AC(ac))) rs[1] = NEG (rs[1]);                 /* sign of remainder */\r
return TRUE;\r
}\r
\r
/* Double precision multiply.  This is done the old fashioned way.  Cross\r
   product multiplies would be a lot faster but would require more code.\r
*/\r
\r
void dmul (int32 ac, d10 *mpy)\r
{\r
int32 p1 = ADDAC (ac, 1);\r
int32 p2 = ADDAC (ac, 2);\r
int32 p3 = ADDAC (ac, 3);\r
int32 i;\r
d10 mpc[2], sign;\r
\r
mpc[0] = AC(ac);                                        /* mplcnd hi */\r
mpc[1] = CLRS (AC(p1));                                 /* mplcnd lo, clr sgn */\r
sign = mpc[0] ^ mpy[0];                                 /* sign of result */\r
if (TSTS (mpc[0])) { DMOVN (mpc); }                     /* get abs (mpcnd) */\r
if (TSTS (mpy[0])) { DMOVN (mpy); }                     /* get abs (mpyer) */\r
else mpy[1] = CLRS (mpy[1]);                            /* clear mpy lo sign */\r
AC(ac) = AC(p1) = AC(p2) = AC(p3) = 0;                  /* clear AC's */\r
if (((mpy[0] | mpy[1]) == 0) || ((mpc[0] | mpc[1]) == 0)) return;\r
for (i = 0; i < 71; i++) {                              /* 71 mpyer bits */\r
    if (i) {                                            /* shift res, mpy */\r
        AC(p3) = (AC(p3) >> 1) | ((AC(p2) & 1) << 34);\r
        AC(p2) = (AC(p2) >> 1) | ((AC(p1) & 1) << 34);\r
        AC(p1) = (AC(p1) >> 1) | ((AC(ac) & 1) << 34);\r
        AC(ac) = AC(ac) >> 1;\r
        mpy[1] = (mpy[1] >> 1) | ((mpy[0] & 1) << 34);\r
        mpy[0] = mpy[0] >> 1;\r
        }\r
    if (mpy[1] & 1) {                                   /* if mpy lo bit = 1 */\r
        AC(p1) = AC(p1) + mpc[1];\r
        AC(ac) = AC(ac) + mpc[0] + (TSTS (AC(p1) != 0));\r
        AC(p1) = CLRS (AC(p1));\r
        }\r
    }\r
if (TSTS (sign)) {                                      /* result minus? */\r
    AC(p3) = (-AC(p3)) & MMASK;                         /* quad negate */\r
    AC(p2) = (~AC(p2) + (AC(p3) == 0)) & MMASK;\r
    AC(p1) = (~AC(p1) + (AC(p2) == 0)) & MMASK;\r
    AC(ac) = (~AC(ac) + (AC(p1) == 0)) & DMASK;\r
    }\r
else if (TSTS (AC(ac))) SETF (F_AOV | F_T1);            /* wrong sign */\r
if (TSTS (AC(ac))) {                                    /* if result - */\r
    AC(p1) = SETS (AC(p1));                             /* make signs consistent */\r
    AC(p2) = SETS (AC(p2));\r
    AC(p3) = SETS (AC(p3));\r
    }\r
return;\r
}\r
\r
/* Double precision divide - checked against KS10 ucode */\r
\r
void ddiv (int32 ac, d10 *dvr)\r
{\r
int32 i, cryin;\r
d10 sign, qu[2], dvd[4];\r
\r
dvd[0] = AC(ac);                                        /* save dividend */\r
for (i = 1; i < 4; i++) dvd[i] = CLRS (AC(ADDAC (ac, i)));\r
sign = AC(ac) ^ dvr[0];                                 /* sign of result */\r
if (TSTS (AC(ac))) {                                    /* get abs (dividend) */\r
    for (i = 3, cryin = 1; i > 0; i--) {                /* negate quad */\r
        dvd[i] = (~dvd[i] + cryin) & MMASK;             /* comp + carry in */\r
        if (dvd[i]) cryin = 0;                          /* next carry in */\r
        }\r
    dvd[0] = (~dvd[0] + cryin) & DMASK;\r
    }\r
if (TSTS (dvr[0])) { DMOVN (dvr); }                     /* get abs (divisor) */\r
else dvr[1] = CLRS (dvr[1]);\r
if (DCMPGE (dvd, dvr)) {                                /* will divide work? */\r
    SETF (F_AOV | F_DCK | F_T1);                        /* no, set flags */\r
    return;\r
    }\r
qu[0] = qu[1] = 0;                                      /* clear quotient */\r
for (i = 0; i < 70; i++) {                              /* 70 quotient bits */\r
    dvd[0] = ((dvd[0] << 1) | ((dvd[1] >> 34) & 1)) & DMASK;;\r
    dvd[1] = ((dvd[1] << 1) | ((dvd[2] >> 34) & 1)) & MMASK;\r
    dvd[2] = ((dvd[2] << 1) | ((dvd[3] >> 34) & 1)) & MMASK;\r
    dvd[3] = (dvd[3] << 1) & MMASK;                     /* shift dividend */\r
    qu[0] = (qu[0] << 1) | ((qu[1] >> 34) & 1);         /* shift quotient */\r
    qu[1] = (qu[1] << 1) & MMASK;\r
    if (DCMPGE (dvd, dvr)) {                            /* subtract work? */\r
        dvd[0] = dvd[0] - dvr[0] - (dvd[1] < dvr[1]);\r
        dvd[1] = (dvd[1] - dvr[1]) & MMASK;             /* do subtract */\r
        qu[1] = qu[1] + 1;                              /* set quotient bit */\r
        }\r
    }\r
if (TSTS (sign) && (qu[0] | qu[1])) { MKDNEG (qu); }\r
if (TSTS (AC(ac)) && (dvd[0] | dvd[1])) { MKDNEG (dvd); }\r
AC(ac) = qu[0];                                         /* quotient */\r
AC(ADDAC(ac, 1)) = qu[1];\r
AC(ADDAC(ac, 2)) = dvd[0];                              /* remainder */\r
AC(ADDAC(ac, 3)) = dvd[1];\r
return;\r
}\r
\r
/* Single precision floating add/subtract - checked against KS10 ucode\r
   The KS10 shifts the smaller operand regardless of the exponent diff.\r
   This code will not shift more than 63 places; shifts beyond that\r
   cannot change the value of the smaller operand.\r
\r
   If the signs of the operands are the same, the result sign is the\r
   same as the source sign; the sign of the result fraction is actually\r
   part of the data.  If the signs of the operands are different, the\r
   result sign is determined by the fraction sign.\r
*/\r
\r
d10 fad (d10 op1, d10 op2, t_bool rnd, int32 inv)\r
{\r
int32 ediff;\r
UFP a, b, t;\r
\r
if (inv) op2 = NEG (op2);                               /* subtract? -b */\r
if (op1 == 0) funpack (op2, 0, &a, AFRC);               /* a = 0? result is b */\r
else if (op2 == 0) funpack (op1, 0, &a, AFRC);          /* b = 0? result is a */\r
else {\r
    funpack (op1, 0, &a, SFRC);                         /* unpack operands */\r
    funpack (op2, 0, &b, SFRC);                         /* fracs are 2's comp */\r
    ediff = a.exp - b.exp;                              /* get exp diff */\r
    if (ediff < 0) {                                    /* a < b? switch */\r
        t = a;\r
        a = b;\r
        b = t;\r
        ediff = -ediff;\r
        }\r
    if (ediff > 63) ediff = 63;                         /* cap diff at 63 */\r
    if (ediff) b.fhi = (t_int64) b.fhi >> ediff;        /* shift b (signed) */\r
    a.fhi = a.fhi + b.fhi;                              /* add fractions */\r
    if (a.sign ^ b.sign) {                              /* add or subtract? */\r
        if (a.fhi & FP_UCRY) {                          /* subtract, frac -? */\r
            a.fhi = UNEG (a.fhi);                       /* complement result */\r
            a.sign = 1;                                 /* result is - */\r
            }\r
        else a.sign = 0;                                /* result is + */\r
        }\r
    else {\r
        if (a.sign) a.fhi = UNEG (a.fhi);               /* add, src -? comp */\r
        if (a.fhi & FP_UCRY) {                          /* check for carry */\r
            a.fhi = a.fhi >> 1;                         /* flo won't be used */\r
            a.exp = a.exp + 1;\r
            }\r
        }\r
    }\r
fnorm (&a, (rnd? FP_URNDS: 0));                         /* normalize, round */\r
return fpack (&a, NULL, FALSE);\r
}\r
\r
/* Single precision floating multiply.  Because the fractions are 27b,\r
   a 64b multiply can be used for the fraction multiply.  The 27b\r
   fractions are positioned 0'frac'0000, resulting in 00'hifrac'0..0.\r
   The extra 0 is accounted for by biasing the result exponent.\r
*/\r
\r
#define FP_V_SPM        (FP_V_UFHI - (32 - FP_N_FHI - 1))\r
d10 fmp (d10 op1, d10 op2, t_bool rnd)\r
{\r
UFP a, b;\r
\r
funpack (op1, 0, &a, AFRC);                             /* unpack operands */\r
funpack (op2, 0, &b, AFRC);                             /* fracs are abs val */\r
if ((a.fhi == 0) || (b.fhi == 0)) return 0;             /* either 0?  */\r
a.sign = a.sign ^ b.sign;                               /* result sign */\r
a.exp = a.exp + b.exp - FP_BIAS + 1;                    /* result exponent */\r
a.fhi = (a.fhi >> FP_V_SPM) * (b.fhi >> FP_V_SPM);      /* high 27b of result */\r
fnorm (&a, (rnd? FP_URNDS: 0));                         /* normalize, round */\r
return fpack (&a, NULL, FALSE);\r
}\r
\r
/* Single precision floating divide.  Because the fractions are 27b, a\r
   64b divide can be used for the fraction divide.  Note that 28b-29b\r
   of fraction are developed; the code will do one special normalize to\r
   make sure that the 28th bit is not lost.  Also note the special\r
   treatment of negative quotients with non-zero remainders; this\r
   implements the note on p2-23 of the Processor Reference Manual.\r
*/\r
\r
t_bool fdv (d10 op1, d10 op2, d10 *rs, t_bool rnd)\r
{\r
UFP a, b;\r
t_uint64 savhi;\r
t_bool rem = FALSE;\r
\r
funpack (op1, 0, &a, AFRC);                             /* unpack operands */\r
funpack (op2, 0, &b, AFRC);                             /* fracs are abs val */\r
if (a.fhi >= 2 * b.fhi) {                               /* will divide work? */\r
    SETF (F_AOV | F_DCK | F_FOV | F_T1);\r
    return FALSE;\r
    }\r
if (savhi = a.fhi) {                                    /* dvd = 0? quo = 0 */\r
    a.sign = a.sign ^ b.sign;                           /* result sign */\r
    a.exp = a.exp - b.exp + FP_BIAS + 1;                /* result exponent */\r
    a.fhi = a.fhi / (b.fhi >> (FP_N_FHI + 1));          /* do divide */\r
    if (a.sign && (savhi != (a.fhi * (b.fhi >> (FP_N_FHI + 1)))))\r
        rem = TRUE;                                     /* KL/KS hack */\r
    a.fhi = a.fhi << (FP_V_UNORM - FP_N_FHI - 1);       /* put quo in place */\r
    if ((a.fhi & FP_UNORM) == 0) {                      /* normalize 1b */\r
        a.fhi = a.fhi << 1;                             /* before masking */\r
        a.exp = a.exp - 1;\r
        }\r
    a.fhi = a.fhi & (FP_UFHI | FP_URNDS);               /* mask quo to 28b */\r
    }\r
fnorm (&a, (rnd? FP_URNDS: 0));                         /* normalize, round */\r
*rs = fpack (&a, NULL, rem);                            /* pack result */\r
return TRUE;\r
}\r
\r
/* Single precision floating scale. */\r
\r
d10 fsc (d10 val, a10 ea)\r
{\r
int32 sc = LIT8 (ea);\r
UFP a;\r
\r
if (val == 0) return 0;\r
funpack (val, 0, &a, AFRC);                             /* unpack operand */\r
if (ea & RSIGN) a.exp = a.exp - sc;                     /* adjust exponent */\r
else a.exp = a.exp + sc;\r
fnorm (&a, 0);                                          /* renormalize */\r
return fpack (&a, NULL, FALSE);                         /* pack result */\r
}\r
\r
/* Float integer operand and round */\r
\r
d10 fltr (d10 mb)\r
{\r
UFP a;\r
d10 val = ABS (mb);\r
\r
a.sign = GET_FPSIGN (mb);                               /* get sign */\r
a.exp = FP_BIAS + 36;                                   /* initial exponent */\r
a.fhi = val << (FP_V_UNORM - 35);                       /* left justify op */\r
a.flo = 0;\r
fnorm (&a, FP_URNDS);                                   /* normalize, round */\r
return fpack (&a, NULL, FALSE);                         /* pack result */\r
}\r
\r
/* Fix and truncate/round floating operand */\r
\r
void fix (int32 ac, d10 mb, t_bool rnd)\r
{\r
int32 sc;\r
t_uint64 so;\r
UFP a;\r
\r
funpack (mb, 0, &a, AFRC);                              /* unpack operand */\r
if (a.exp > (FP_BIAS + FP_N_FHI + FP_N_EXP)) SETF (F_AOV | F_T1);\r
else if (a.exp < FP_BIAS) AC(ac) = 0;                   /* < 1/2? */\r
else {\r
    sc = FP_V_UNORM - (a.exp - FP_BIAS) + 1;\r
    AC(ac) = a.fhi >> sc;\r
    if (rnd) {\r
        so = a.fhi << (64 - sc);\r
        if (so >= (0x8000000000000000 + a.sign)) AC(ac) = AC(ac) + 1;\r
        }\r
    if (a.sign) AC(ac) = NEG (AC(ac));\r
    }\r
return;\r
}\r
\r
/* Double precision floating add/subtract\r
   Since a.flo is 0, adding b.flo is just a copy - this is incorporated into\r
   the denormalization step.  If there's no denormalization, bflo is zero too.\r
*/\r
\r
void dfad (int32 ac, d10 *rs, int32 inv)\r
{\r
int32 p1 = ADDAC (ac, 1);\r
int32 ediff;\r
UFP a, b, t;\r
\r
if (inv) { DMOVN (rs); }                                /* subtract? -b */\r
if ((AC(ac) | AC(p1)) == 0) funpack (rs[0], rs[1], &a, AFRC);\r
                                                        /* a == 0? sum = b */\r
else if ((rs[0] | rs[1]) == 0) funpack (AC(ac), AC(p1), &a, AFRC);\r
                                                        /* b == 0? sum = a */\r
else {\r
    funpack (AC(ac), AC(p1), &a, SFRC);                 /* unpack operands */\r
    funpack (rs[0], rs[1], &b, SFRC);\r
    ediff = a.exp - b.exp;                              /* get exp diff */\r
    if (ediff < 0) {                                    /* a < b? switch */\r
        t = a;\r
        a = b;\r
        b = t;\r
        ediff = -ediff;\r
        }\r
    if (ediff > 127) ediff = 127;                       /* cap diff at 127 */\r
    if (ediff > 63) {                                   /* diff > 63? */\r
        a.flo = (t_int64) b.fhi >> (ediff - 64);        /* b hi to a lo */\r
        b.fhi = b.sign? FP_ONES: 0;                     /* hi = all sign */\r
        }\r
    else if (ediff) {                                   /* diff <= 63 */\r
        a.flo = (b.flo >> ediff) | (b.fhi << (64 - ediff));\r
        b.fhi = (t_int64) b.fhi >> ediff;               /* shift b (signed) */\r
        }\r
    a.fhi = a.fhi + b.fhi;                              /* do add */\r
    if (a.sign ^ b.sign) {                              /* add or subtract? */\r
        if (a.fhi & FP_UCRY) {                          /* subtract, frac -? */\r
            DUNEG (a);                                  /* complement result */\r
            a.sign = 1;                                 /* result is - */\r
            }\r
        else a.sign = 0;                                /* result is + */\r
        }\r
    else {\r
        if (a.sign) { DUNEG (a); };                     /* add, src -? comp */\r
        if (a.fhi & FP_UCRY) {                          /* check for carry */\r
            a.fhi = a.fhi >> 1;                         /* flo won't be used */\r
            a.exp = a.exp + 1;\r
            }\r
        }\r
    }\r
fnorm (&a, FP_URNDD);                                   /* normalize, round */\r
AC(ac) = fpack (&a, &AC(p1), FALSE);                    /* pack result */\r
return;\r
}\r
\r
/* Double precision floating multiply\r
   The 62b fractions are multiplied, with cross products, to produce a\r
   124b fraction with two leading and two trailing 0's.  Because the\r
   product has 2 leading 0's, instead of the normal 1, an extra\r
   normalization step is needed.  Accordingly, the exponent calculation\r
   increments the result exponent, to compensate for normalization.\r
*/\r
\r
void dfmp (int32 ac, d10 *rs)\r
{\r
int32 p1 = ADDAC (ac, 1);\r
t_uint64 xh, xl, yh, yl, mid;\r
UFP a, b;\r
\r
funpack (AC(ac), AC(p1), &a, AFRC);                     /* unpack operands */\r
funpack (rs[0], rs[1], &b, AFRC);\r
if ((a.fhi == 0) || (b.fhi == 0)) {                     /* either 0? result 0 */\r
    AC(ac) = AC(p1) = 0;\r
    return;\r
    }\r
a.sign = a.sign ^ b.sign;                               /* result sign */\r
a.exp = a.exp + b.exp - FP_BIAS + 1;                    /* result exponent */\r
xh = a.fhi >> 32;                                       /* split 62b fracs */\r
xl = a.fhi & MSK32;                                     /* into 32b halves */\r
yh = b.fhi >> 32;\r
yl = b.fhi & MSK32;\r
a.fhi = xh * yh;                                        /* hi xproduct */\r
a.flo = xl * yl;                                        /* low xproduct */\r
mid = (xh * yl) + (yh * xl);                            /* fits in 64b */\r
a.flo = a.flo + (mid << 32);                            /* add mid lo to lo */\r
a.fhi = a.fhi + ((mid >> 32) & MSK32) + (a.flo < (mid << 32));\r
fnorm (&a, FP_URNDD);                                   /* normalize, round */\r
AC(ac) = fpack (&a, &AC(p1), FALSE);                    /* pack result */\r
return;\r
}\r
\r
/* Double precision floating divide\r
   This algorithm develops a full 62 bits of quotient, plus one rounding\r
   bit, in the low order 63b of a 64b number.  To do this, we must assure\r
   that the initial divide step generates a 1.  If it would fail, shift\r
   the dividend left and decrement the result exponent accordingly.\r
*/\r
\r
void dfdv (int32 ac, d10 *rs)\r
{\r
int32 p1 = ADDAC (ac, 1);\r
int32 i;\r
t_uint64 qu = 0;\r
UFP a, b;\r
\r
funpack (AC(ac), AC(p1), &a, AFRC);                     /* unpack operands */\r
funpack (rs[0], rs[1], &b, AFRC);\r
if (a.fhi >= 2 * b.fhi) {                               /* will divide work? */\r
    SETF (F_AOV | F_DCK | F_FOV | F_T1);\r
    return;\r
    }\r
if (a.fhi) {                                            /* dvd = 0? quo = 0 */\r
    a.sign = a.sign ^ b.sign;                           /* result sign */\r
    a.exp = a.exp - b.exp + FP_BIAS + 1;                /* result exponent */\r
    if (a.fhi < b.fhi) {                                /* make sure initial */\r
        a.fhi = a.fhi << 1;                             /* divide step will work */\r
        a.exp = a.exp - 1;\r
        }\r
    for (i = 0; i < 63; i++) {                          /* 63b of quotient */\r
        qu = qu << 1;                                   /* shift quotient */\r
        if (a.fhi >= b.fhi) {                           /* will div work? */\r
            a.fhi = a.fhi - b.fhi;                      /* sub, quo = 1 */\r
            qu = qu + 1;\r
            }\r
        a.fhi = a.fhi << 1;                             /* shift dividend */\r
        }\r
    a.fhi = qu;\r
    }\r
fnorm (&a, FP_URNDD);                                   /* normalize, round */\r
AC(ac) = fpack (&a, &AC(p1), FALSE);                    /* pack result */\r
return;\r
}\r
\r
/* Unpack floating point operand */\r
\r
void funpack (d10 h, d10 l, UFP *r, t_bool sgn)\r
{\r
d10 fphi, fplo;\r
\r
r->sign = GET_FPSIGN (h);\r
r->exp = GET_FPEXP (h);\r
fphi = GET_FPHI (h);\r
fplo = GET_FPLO (l);\r
r->fhi = (fphi << FP_V_UFHI) | (fplo << FP_V_UFLO);\r
r->flo = 0;\r
if (r->sign) {\r
    r->exp = r->exp ^ FP_M_EXP;                         /* 1s comp exp */\r
    if (sgn) {                                          /* signed frac? */\r
        if (r->fhi) r->fhi = r->fhi | FP_UCRY;          /* extend sign */ \r
        else {\r
            r->exp = r->exp + 1;\r
            r->fhi = FP_UCRY | FP_UNORM;\r
            }\r
        }\r
    else {                                              /* abs frac */\r
        if (r->fhi) r->fhi = UNEG (r->fhi) & FP_UFRAC;\r
        else {\r
                r->exp = r->exp + 1;\r
            r->fhi = FP_UNORM;\r
            }\r
        }\r
    }\r
return;\r
}\r
\r
/* Normalize and optionally round floating point operand */\r
 \r
void fnorm (UFP *a, t_int64 rnd)\r
{\r
int32 i;\r
static t_uint64 normmask[6] = {\r
 0x6000000000000000, 0x7800000000000000, 0x7F80000000000000,\r
 0x7FFF800000000000, 0x7FFFFFFF80000000, 0x7FFFFFFFFFFFFFFF\r
 };\r
static int32 normtab[7] = { 1, 2, 4, 8, 16, 32, 63 };\r
extern a10 pager_PC;\r
\r
if (a->fhi & FP_UCRY) {                                 /* carry set? */\r
    printf ("%%PDP-10 FP: carry bit set at normalization, PC = %o\n", pager_PC);\r
    a->flo = (a->flo >> 1) | ((a->fhi & 1) << 63);      /* try to recover */\r
    a->fhi = a->fhi >> 1;                               /* but root cause */\r
    a->exp = a->exp + 1;                                /* should be fixed! */\r
    }\r
if ((a->fhi | a->flo) == 0) {                           /* if fraction = 0 */\r
    a->sign = a->exp = 0;                               /* result is 0 */\r
    return;\r
    }\r
while ((a->fhi & FP_UNORM) == 0) {                      /* normalized? */\r
    for (i = 0; i < 6; i++) {\r
        if (a->fhi & normmask[i]) break;\r
        }\r
    a->fhi = (a->fhi << normtab[i]) | (a->flo >> (64 - normtab[i]));\r
    a->flo = a->flo << normtab[i];\r
    a->exp = a->exp - normtab[i];\r
    }\r
if (rnd) {                                              /* rounding? */\r
    a->fhi = a->fhi + rnd;                              /* add round const */\r
    if (a->fhi & FP_UCRY) {                             /* if carry out, */\r
        a->fhi = a->fhi >> 1;                           /* renormalize */\r
        a->exp = a->exp + 1;\r
        }\r
    }\r
return;\r
}\r
\r
/* Pack floating point result */\r
\r
d10 fpack (UFP *r, d10 *lo, t_bool fdvneg)\r
{\r
d10 val[2];\r
\r
if (r->exp < 0) SETF (F_AOV | F_FOV | F_FXU | F_T1);\r
else if (r->exp > FP_M_EXP) SETF (F_AOV | F_FOV | F_T1);\r
val[0] = (((((d10) r->exp) & FP_M_EXP) << FP_V_EXP) |\r
    ((r->fhi & FP_UFHI) >> FP_V_UFHI)) & DMASK;\r
if (lo) val[1] = ((r->fhi & FP_UFLO) >> FP_V_UFLO) & MMASK;\r
else val[1] = 0;\r
if (r->sign) {                                          /* negate? */\r
    if (fdvneg) {                                       /* fdvr special? */\r
        val[1] = ~val[1] & MMASK;                       /* 1's comp */\r
        val[0] = ~val[0] & DMASK;\r
		}\r
    else { DMOVN (val); }                               /* 2's comp */\r
    }\r
if (lo) *lo = val[1];\r
return val[0];\r
}\r