First Commit of my working state

[simh.git] / PDP10 / pdp10_mdfp.c

/* pdp10_mdfp.c: PDP-10 multiply/divide and floating point simulator

   Copyright (c) 1993-2005, Robert M Supnik

   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
   ROBERT M SUPNIK 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 name of Robert M Supnik shall not be
   used in advertising or otherwise to promote the sale, use or other dealings
   in this Software without prior written authorization from Robert M Supnik.

   2-Apr-04     RMS     Fixed bug in floating point unpack
                        Fixed bug in FIXR (found by Phil Stone, fixed by
                        Chris Smith)
   31-Aug-01    RMS     Changed int64 to t_int64 for Windoze
   10-Aug-01    RMS     Removed register in declarations

   Instructions handled in this module:
        imul            integer multiply
        idiv            integer divide
        mul             multiply
        div             divide
        dmul            double precision multiply
        ddiv            double precision divide
        fad(r)          floating add (and round)
        fsb(r)          floating subtract (and round)
        fmp(r)          floating multiply (and round)
        fdv(r)          floating divide and round
        fsc             floating scale
        fix(r)          floating to fixed (and round)
        fltr            fixed to floating and round
        dfad            double precision floating add/subtract
        dfmp            double precision floating multiply
        dfdv            double precision floating divide

   The PDP-10 stores double (quad) precision integers in sequential
   AC's or memory locations.  Integers are stored in 2's complement
   form.  Only the sign of the high order word matters; the signs
   in low order words are ignored on input and set to the sign of
   the result on output.  Quad precision integers exist only in the
   AC's as the result of a DMUL or the dividend of a DDIV.

    0 00000000011111111112222222222333333
    0 12345678901234567890123456789012345
   +-+-----------------------------------+
   |S|      high order integer           | AC(n), A
   +-+-----------------------------------+
   |S|      low order integer            | AC(n + 1), A + 1
   +-+-----------------------------------+
   |S|      low order integer            | AC(n + 2)
   +-+-----------------------------------+
   |S|      low order integer            | AC(n + 3)
   +-+-----------------------------------+

   The PDP-10 supports two floating point formats: single and double
   precision.  In both, the exponent is 8 bits, stored in excess
   128 notation.  The fraction is expected to be normalized.  A
   single precision floating point number has 27 bits of fraction;
   a double precision number has 62 bits of fraction (the sign
   bit of the second word is ignored and is set to zero).

   In a negative floating point number, the exponent is stored in
   one's complement form, the fraction in two's complement form.

    0 00000000 011111111112222222222333333
    0 12345678 901234567890123456789012345
   +-+--------+---------------------------+
   |S|exponent|      high order fraction  | AC(n), A
   +-+--------+---------------------------+
   |0|      low order fraction            | AC(n + 1), A + 1
   +-+------------------------------------+

   Note that treatment of the sign is different for double precision
   integers and double precision floating point.  DMOVN (implemented
   as an inline macro) follows floating point conventions.

   The original PDP-10 CPU (KA10) used a different format for double
   precision numbers and included certain instructions to make
   software support easier.  These instructions were phased out in
   the KL10 and KS10 and are treated as MUUO's.

   The KL10 added extended precision (11-bit exponent) floating point
   format (so-called G floating).  These instructions were not
   implemented in the KS10 and are treated as MUUO's.
*/

#include "pdp10_defs.h"
#include <setjmp.h>

typedef struct {                                        /* unpacked fp number */
    int32               sign;                           /* sign */
    int32               exp;                            /* exponent */
    t_uint64    fhi;                                    /* fraction high */
    t_uint64    flo;                                    /* for double prec */
    } UFP;

#define MSK32           0xFFFFFFFF
#define FIT27           (DMASK - 0x07FFFFFF)
#define FIT32           (DMASK - MSK32)
#define SFRC            TRUE                            /* frac 2's comp */
#define AFRC            FALSE                           /* frac abs value */

/* In packed floating point number */

#define FP_BIAS         0200                            /* exponent bias */
#define FP_N_FHI        27                              /* # of hi frac bits */
#define FP_V_FHI        0                               /* must be zero */
#define FP_M_FHI        0000777777777
#define FP_N_EXP        8                               /* # of exp bits */
#define FP_V_EXP        (FP_V_FHI + FP_N_FHI)
#define FP_M_EXP        0377
#define FP_V_SIGN       (FP_V_EXP + FP_N_EXP)           /* sign */
#define FP_N_FLO        35                              /* # of lo frac bits */
#define FP_V_FLO        0                               /* must be zero */
#define FP_M_FLO        0377777777777
#define GET_FPSIGN(x)   ((int32) (((x) >> FP_V_SIGN) & 1))
#define GET_FPEXP(x)    ((int32) (((x) >> FP_V_EXP) & FP_M_EXP))
#define GET_FPHI(x)     ((x) & FP_M_FHI)
#define GET_FPLO(x)     ((x) & FP_M_FLO)

/* In unpacked floating point number */

#define FP_N_GUARD      1                               /* # of guard bits */
#define FP_V_UFLO       FP_N_GUARD                      /* <35:1> */
#define FP_V_URNDD      (FP_V_UFLO - 1)                 /* dp round bit */
#define FP_V_UFHI       (FP_V_UFLO + FP_N_FLO)          /* <62:36> */
#define FP_V_URNDS      (FP_V_UFHI - 1)                 /* sp round bit */
#define FP_V_UCRY       (FP_V_UFHI + FP_N_FHI)          /* <63> */
#define FP_V_UNORM      (FP_V_UCRY - 1)                 /* normalized bit */
#define FP_UFHI         0x7FFFFFF000000000
#define FP_UFLO         0x0000000FFFFFFFFE
#define FP_UFRAC        0x7FFFFFFFFFFFFFFE
#define FP_URNDD        0x0000000000000001
#define FP_URNDS        0x0000000800000000
#define FP_UNORM        0x4000000000000000
#define FP_UCRY         0x8000000000000000
#define FP_ONES         0xFFFFFFFFFFFFFFFF

#define UNEG(x)         ((~x) + 1)
#define DUNEG(x)        x.flo = UNEG (x.flo); x.fhi = ~x.fhi + (x.flo == 0)

extern d10 *ac_cur;                                     /* current AC block */
extern int32 flags;                                     /* flags */
void mul (d10 a, d10 b, d10 *rs);
void funpack (d10 h, d10 l, UFP *r, t_bool sgn);
void fnorm (UFP *r, t_int64 rnd);
d10 fpack (UFP *r, d10 *lo, t_bool fdvneg);

/* Integer multiply - checked against KS-10 ucode */

d10 imul (d10 a, d10 b)
{
d10 rs[2];

if ((a == SIGN) && (b == SIGN)) {                       /* KS10 hack */
    SETF (F_AOV | F_T1);                                /* -2**35 squared */
    return SIGN;
    }
mul (a, b, rs);                                         /* mpy, dprec result */
if (rs[0] && (rs[0] != ONES)) {                         /* high not all sign? */
    rs[1] = TSTS (a ^ b)? SETS (rs[1]): CLRS (rs[1]);   /* set sign */
    SETF (F_AOV | F_T1);                                /* overflow */
    }
return rs[1];
}

/* Integer divide, return quotient, remainder  - checked against KS10 ucode
   The KS10 does not recognize -2^35/-1 as an error.  Instead, it produces
   2^35 (that is, -2^35) as the incorrect result.
*/

t_bool idiv (d10 a, d10 b, d10 *rs)
{
d10 dvd = ABS (a);                                      /* make ops positive */
d10 dvr = ABS (b);

if (dvr == 0) {                                         /* divide by 0? */
    SETF (F_DCK | F_AOV | F_T1);                        /* set flags, return */
    return FALSE;
    }
rs[0] = dvd / dvr;                                      /* get quotient */
rs[1] = dvd % dvr;                                      /* get remainder */
if (TSTS (a ^ b)) rs[0] = NEG (rs[0]);                  /* sign of result */
if (TSTS (a)) rs[1] = NEG (rs[1]);                      /* sign of remainder */
return TRUE;
}

/* Multiply, return double precision result - checked against KS10 ucode */

void mul (d10 s1, d10 s2, d10 *rs)
{
t_uint64 a = ABS (s1);
t_uint64 b = ABS (s2);
t_uint64 t, u, r;

if ((a == 0) || (b == 0)) {                             /* operand = 0? */
    rs[0] = rs[1] = 0;                                  /* result 0 */
    return;
    }
if ((a & FIT32) || (b & FIT32)) {                       /* fit in 64b? */
    t = a >> 18;                                        /* no, split in half */
    a = a & RMASK;                                      /* "dp" multiply */
    u = b >> 18;
    b = b & RMASK;
    r = (a * b) + (((a * u) + (b * t)) << 18);          /* low is only 35b */
    rs[0] = ((t * u) << 1) + (r >> 35);                 /* so lsh hi 1 */
    rs[1] = r & MMASK;
    }
else {
    r = a * b;                                          /* fits, native mpy */
    rs[0] = r >> 35;                                    /* split at bit 35 */
    rs[1] = r & MMASK;
    }

if (TSTS (s1 ^ s2)) { MKDNEG (rs); }                    /* result -? */
else if (TSTS (rs[0])) {                                /* result +, 2**70? */
    SETF (F_AOV | F_T1);                                /* overflow */
    rs[1] = SETS (rs[1]);                               /* consistent - */
    }
return;
}

/* Divide, return quotient and remainder - checked against KS10 ucode
   Note that the initial divide check catches the case -2^70/-2^35;
   thus, the quotient can have at most 35 bits.
*/

t_bool divi (int32 ac, d10 b, d10 *rs)
{
int32 p1 = ADDAC (ac, 1);
d10 dvr = ABS (b);                                      /* make divr positive */
t_int64 t;
int32 i;
d10 dvd[2];

dvd[0] = AC(ac);                                        /* divd high */
dvd[1] = CLRS (AC(p1));                                 /* divd lo, clr sgn */
if (TSTS (AC(ac))) { DMOVN (dvd); }                     /* make divd positive */
if (dvd[0] >= dvr) {                                    /* divide fail? */
    SETF (F_AOV | F_DCK | F_T1);                        /* set flags, return */
    return FALSE;
    }
if (dvd[0] & FIT27) {                                   /* fit in 63b? */
    for (i = 0, rs[0] = 0; i < 35; i++) {               /* 35 quotient bits */
        dvd[0] = (dvd[0] << 1) | ((dvd[1] >> 34) & 1);
        dvd[1] = (dvd[1] << 1) & MMASK;                 /* shift dividend */
        rs[0] = rs[0] << 1;                             /* shift quotient */
        if (dvd[0] >= dvr) {                            /* subtract work? */
            dvd[0] = dvd[0] - dvr;                      /* quo bit is 1 */
            rs[0] = rs[0] + 1;
            }
        }
    rs[1] = dvd[0];                                     /* store remainder */
    }
else {
    t = (dvd[0] << 35) | dvd[1];                        /* concatenate */
    rs[0] = t / dvr;                                    /* quotient */
    rs[1] = t % dvr;                                    /* remainder */
    }
if (TSTS (AC(ac) ^ b)) rs[0] = NEG (rs[0]);             /* sign of result */
if (TSTS (AC(ac))) rs[1] = NEG (rs[1]);                 /* sign of remainder */
return TRUE;
}

/* Double precision multiply.  This is done the old fashioned way.  Cross
   product multiplies would be a lot faster but would require more code.
*/

void dmul (int32 ac, d10 *mpy)
{
int32 p1 = ADDAC (ac, 1);
int32 p2 = ADDAC (ac, 2);
int32 p3 = ADDAC (ac, 3);
int32 i;
d10 mpc[2], sign;

mpc[0] = AC(ac);                                        /* mplcnd hi */
mpc[1] = CLRS (AC(p1));                                 /* mplcnd lo, clr sgn */
sign = mpc[0] ^ mpy[0];                                 /* sign of result */
if (TSTS (mpc[0])) { DMOVN (mpc); }                     /* get abs (mpcnd) */
if (TSTS (mpy[0])) { DMOVN (mpy); }                     /* get abs (mpyer) */
else mpy[1] = CLRS (mpy[1]);                            /* clear mpy lo sign */
AC(ac) = AC(p1) = AC(p2) = AC(p3) = 0;                  /* clear AC's */
if (((mpy[0] | mpy[1]) == 0) || ((mpc[0] | mpc[1]) == 0)) return;
for (i = 0; i < 71; i++) {                              /* 71 mpyer bits */
    if (i) {                                            /* shift res, mpy */
        AC(p3) = (AC(p3) >> 1) | ((AC(p2) & 1) << 34);
        AC(p2) = (AC(p2) >> 1) | ((AC(p1) & 1) << 34);
        AC(p1) = (AC(p1) >> 1) | ((AC(ac) & 1) << 34);
        AC(ac) = AC(ac) >> 1;
        mpy[1] = (mpy[1] >> 1) | ((mpy[0] & 1) << 34);
        mpy[0] = mpy[0] >> 1;
        }
    if (mpy[1] & 1) {                                   /* if mpy lo bit = 1 */
        AC(p1) = AC(p1) + mpc[1];
        AC(ac) = AC(ac) + mpc[0] + (TSTS (AC(p1) != 0));
        AC(p1) = CLRS (AC(p1));
        }
    }
if (TSTS (sign)) {                                      /* result minus? */
    AC(p3) = (-AC(p3)) & MMASK;                         /* quad negate */
    AC(p2) = (~AC(p2) + (AC(p3) == 0)) & MMASK;
    AC(p1) = (~AC(p1) + (AC(p2) == 0)) & MMASK;
    AC(ac) = (~AC(ac) + (AC(p1) == 0)) & DMASK;
    }
else if (TSTS (AC(ac))) SETF (F_AOV | F_T1);            /* wrong sign */
if (TSTS (AC(ac))) {                                    /* if result - */
    AC(p1) = SETS (AC(p1));                             /* make signs consistent */
    AC(p2) = SETS (AC(p2));
    AC(p3) = SETS (AC(p3));
    }
return;
}

/* Double precision divide - checked against KS10 ucode */

void ddiv (int32 ac, d10 *dvr)
{
int32 i, cryin;
d10 sign, qu[2], dvd[4];

dvd[0] = AC(ac);                                        /* save dividend */
for (i = 1; i < 4; i++) dvd[i] = CLRS (AC(ADDAC (ac, i)));
sign = AC(ac) ^ dvr[0];                                 /* sign of result */
if (TSTS (AC(ac))) {                                    /* get abs (dividend) */
    for (i = 3, cryin = 1; i > 0; i--) {                /* negate quad */
        dvd[i] = (~dvd[i] + cryin) & MMASK;             /* comp + carry in */
        if (dvd[i]) cryin = 0;                          /* next carry in */
        }
    dvd[0] = (~dvd[0] + cryin) & DMASK;
    }
if (TSTS (dvr[0])) { DMOVN (dvr); }                     /* get abs (divisor) */
else dvr[1] = CLRS (dvr[1]);
if (DCMPGE (dvd, dvr)) {                                /* will divide work? */
    SETF (F_AOV | F_DCK | F_T1);                        /* no, set flags */
    return;
    }
qu[0] = qu[1] = 0;                                      /* clear quotient */
for (i = 0; i < 70; i++) {                              /* 70 quotient bits */
    dvd[0] = ((dvd[0] << 1) | ((dvd[1] >> 34) & 1)) & DMASK;;
    dvd[1] = ((dvd[1] << 1) | ((dvd[2] >> 34) & 1)) & MMASK;
    dvd[2] = ((dvd[2] << 1) | ((dvd[3] >> 34) & 1)) & MMASK;
    dvd[3] = (dvd[3] << 1) & MMASK;                     /* shift dividend */
    qu[0] = (qu[0] << 1) | ((qu[1] >> 34) & 1);         /* shift quotient */
    qu[1] = (qu[1] << 1) & MMASK;
    if (DCMPGE (dvd, dvr)) {                            /* subtract work? */
        dvd[0] = dvd[0] - dvr[0] - (dvd[1] < dvr[1]);
        dvd[1] = (dvd[1] - dvr[1]) & MMASK;             /* do subtract */
        qu[1] = qu[1] + 1;                              /* set quotient bit */
        }
    }
if (TSTS (sign) && (qu[0] | qu[1])) { MKDNEG (qu); }
if (TSTS (AC(ac)) && (dvd[0] | dvd[1])) { MKDNEG (dvd); }
AC(ac) = qu[0];                                         /* quotient */
AC(ADDAC(ac, 1)) = qu[1];
AC(ADDAC(ac, 2)) = dvd[0];                              /* remainder */
AC(ADDAC(ac, 3)) = dvd[1];
return;
}

/* Single precision floating add/subtract - checked against KS10 ucode
   The KS10 shifts the smaller operand regardless of the exponent diff.
   This code will not shift more than 63 places; shifts beyond that
   cannot change the value of the smaller operand.

   If the signs of the operands are the same, the result sign is the
   same as the source sign; the sign of the result fraction is actually
   part of the data.  If the signs of the operands are different, the
   result sign is determined by the fraction sign.
*/

d10 fad (d10 op1, d10 op2, t_bool rnd, int32 inv)
{
int32 ediff;
UFP a, b, t;

if (inv) op2 = NEG (op2);                               /* subtract? -b */
if (op1 == 0) funpack (op2, 0, &a, AFRC);               /* a = 0? result is b */
else if (op2 == 0) funpack (op1, 0, &a, AFRC);          /* b = 0? result is a */
else {
    funpack (op1, 0, &a, SFRC);                         /* unpack operands */
    funpack (op2, 0, &b, SFRC);                         /* fracs are 2's comp */
    ediff = a.exp - b.exp;                              /* get exp diff */
    if (ediff < 0) {                                    /* a < b? switch */
        t = a;
        a = b;
        b = t;
        ediff = -ediff;
        }
    if (ediff > 63) ediff = 63;                         /* cap diff at 63 */
    if (ediff) b.fhi = (t_int64) b.fhi >> ediff;        /* shift b (signed) */
    a.fhi = a.fhi + b.fhi;                              /* add fractions */
    if (a.sign ^ b.sign) {                              /* add or subtract? */
        if (a.fhi & FP_UCRY) {                          /* subtract, frac -? */
            a.fhi = UNEG (a.fhi);                       /* complement result */
            a.sign = 1;                                 /* result is - */
            }
        else a.sign = 0;                                /* result is + */
        }
    else {
        if (a.sign) a.fhi = UNEG (a.fhi);               /* add, src -? comp */
        if (a.fhi & FP_UCRY) {                          /* check for carry */
            a.fhi = a.fhi >> 1;                         /* flo won't be used */
            a.exp = a.exp + 1;
            }
        }
    }
fnorm (&a, (rnd? FP_URNDS: 0));                         /* normalize, round */
return fpack (&a, NULL, FALSE);
}

/* Single precision floating multiply.  Because the fractions are 27b,
   a 64b multiply can be used for the fraction multiply.  The 27b
   fractions are positioned 0'frac'0000, resulting in 00'hifrac'0..0.
   The extra 0 is accounted for by biasing the result exponent.
*/

#define FP_V_SPM        (FP_V_UFHI - (32 - FP_N_FHI - 1))
d10 fmp (d10 op1, d10 op2, t_bool rnd)
{
UFP a, b;

funpack (op1, 0, &a, AFRC);                             /* unpack operands */
funpack (op2, 0, &b, AFRC);                             /* fracs are abs val */
if ((a.fhi == 0) || (b.fhi == 0)) return 0;             /* either 0?  */
a.sign = a.sign ^ b.sign;                               /* result sign */
a.exp = a.exp + b.exp - FP_BIAS + 1;                    /* result exponent */
a.fhi = (a.fhi >> FP_V_SPM) * (b.fhi >> FP_V_SPM);      /* high 27b of result */
fnorm (&a, (rnd? FP_URNDS: 0));                         /* normalize, round */
return fpack (&a, NULL, FALSE);
}

/* Single precision floating divide.  Because the fractions are 27b, a
   64b divide can be used for the fraction divide.  Note that 28b-29b
   of fraction are developed; the code will do one special normalize to
   make sure that the 28th bit is not lost.  Also note the special
   treatment of negative quotients with non-zero remainders; this
   implements the note on p2-23 of the Processor Reference Manual.
*/

t_bool fdv (d10 op1, d10 op2, d10 *rs, t_bool rnd)
{
UFP a, b;
t_uint64 savhi;
t_bool rem = FALSE;

funpack (op1, 0, &a, AFRC);                             /* unpack operands */
funpack (op2, 0, &b, AFRC);                             /* fracs are abs val */
if (a.fhi >= 2 * b.fhi) {                               /* will divide work? */
    SETF (F_AOV | F_DCK | F_FOV | F_T1);
    return FALSE;
    }
if (savhi = a.fhi) {                                    /* dvd = 0? quo = 0 */
    a.sign = a.sign ^ b.sign;                           /* result sign */
    a.exp = a.exp - b.exp + FP_BIAS + 1;                /* result exponent */
    a.fhi = a.fhi / (b.fhi >> (FP_N_FHI + 1));          /* do divide */
    if (a.sign && (savhi != (a.fhi * (b.fhi >> (FP_N_FHI + 1)))))
        rem = TRUE;                                     /* KL/KS hack */
    a.fhi = a.fhi << (FP_V_UNORM - FP_N_FHI - 1);       /* put quo in place */
    if ((a.fhi & FP_UNORM) == 0) {                      /* normalize 1b */
        a.fhi = a.fhi << 1;                             /* before masking */
        a.exp = a.exp - 1;
        }
    a.fhi = a.fhi & (FP_UFHI | FP_URNDS);               /* mask quo to 28b */
    }
fnorm (&a, (rnd? FP_URNDS: 0));                         /* normalize, round */
*rs = fpack (&a, NULL, rem);                            /* pack result */
return TRUE;
}

/* Single precision floating scale. */

d10 fsc (d10 val, a10 ea)
{
int32 sc = LIT8 (ea);
UFP a;

if (val == 0) return 0;
funpack (val, 0, &a, AFRC);                             /* unpack operand */
if (ea & RSIGN) a.exp = a.exp - sc;                     /* adjust exponent */
else a.exp = a.exp + sc;
fnorm (&a, 0);                                          /* renormalize */
return fpack (&a, NULL, FALSE);                         /* pack result */
}

/* Float integer operand and round */

d10 fltr (d10 mb)
{
UFP a;
d10 val = ABS (mb);

a.sign = GET_FPSIGN (mb);                               /* get sign */
a.exp = FP_BIAS + 36;                                   /* initial exponent */
a.fhi = val << (FP_V_UNORM - 35);                       /* left justify op */
a.flo = 0;
fnorm (&a, FP_URNDS);                                   /* normalize, round */
return fpack (&a, NULL, FALSE);                         /* pack result */
}

/* Fix and truncate/round floating operand */

void fix (int32 ac, d10 mb, t_bool rnd)
{
int32 sc;
t_uint64 so;
UFP a;

funpack (mb, 0, &a, AFRC);                              /* unpack operand */
if (a.exp > (FP_BIAS + FP_N_FHI + FP_N_EXP)) SETF (F_AOV | F_T1);
else if (a.exp < FP_BIAS) AC(ac) = 0;                   /* < 1/2? */
else {
    sc = FP_V_UNORM - (a.exp - FP_BIAS) + 1;
    AC(ac) = a.fhi >> sc;
    if (rnd) {
        so = a.fhi << (64 - sc);
        if (so >= (0x8000000000000000 + a.sign)) AC(ac) = AC(ac) + 1;
        }
    if (a.sign) AC(ac) = NEG (AC(ac));
    }
return;
}

/* Double precision floating add/subtract
   Since a.flo is 0, adding b.flo is just a copy - this is incorporated into
   the denormalization step.  If there's no denormalization, bflo is zero too.
*/

void dfad (int32 ac, d10 *rs, int32 inv)
{
int32 p1 = ADDAC (ac, 1);
int32 ediff;
UFP a, b, t;

if (inv) { DMOVN (rs); }                                /* subtract? -b */
if ((AC(ac) | AC(p1)) == 0) funpack (rs[0], rs[1], &a, AFRC);
                                                        /* a == 0? sum = b */
else if ((rs[0] | rs[1]) == 0) funpack (AC(ac), AC(p1), &a, AFRC);
                                                        /* b == 0? sum = a */
else {
    funpack (AC(ac), AC(p1), &a, SFRC);                 /* unpack operands */
    funpack (rs[0], rs[1], &b, SFRC);
    ediff = a.exp - b.exp;                              /* get exp diff */
    if (ediff < 0) {                                    /* a < b? switch */
        t = a;
        a = b;
        b = t;
        ediff = -ediff;
        }
    if (ediff > 127) ediff = 127;                       /* cap diff at 127 */
    if (ediff > 63) {                                   /* diff > 63? */
        a.flo = (t_int64) b.fhi >> (ediff - 64);        /* b hi to a lo */
        b.fhi = b.sign? FP_ONES: 0;                     /* hi = all sign */
        }
    else if (ediff) {                                   /* diff <= 63 */
        a.flo = (b.flo >> ediff) | (b.fhi << (64 - ediff));
        b.fhi = (t_int64) b.fhi >> ediff;               /* shift b (signed) */
        }
    a.fhi = a.fhi + b.fhi;                              /* do add */
    if (a.sign ^ b.sign) {                              /* add or subtract? */
        if (a.fhi & FP_UCRY) {                          /* subtract, frac -? */
            DUNEG (a);                                  /* complement result */
            a.sign = 1;                                 /* result is - */
            }
        else a.sign = 0;                                /* result is + */
        }
    else {
        if (a.sign) { DUNEG (a); };                     /* add, src -? comp */
        if (a.fhi & FP_UCRY) {                          /* check for carry */
            a.fhi = a.fhi >> 1;                         /* flo won't be used */
            a.exp = a.exp + 1;
            }
        }
    }
fnorm (&a, FP_URNDD);                                   /* normalize, round */
AC(ac) = fpack (&a, &AC(p1), FALSE);                    /* pack result */
return;
}

/* Double precision floating multiply
   The 62b fractions are multiplied, with cross products, to produce a
   124b fraction with two leading and two trailing 0's.  Because the
   product has 2 leading 0's, instead of the normal 1, an extra
   normalization step is needed.  Accordingly, the exponent calculation
   increments the result exponent, to compensate for normalization.
*/

void dfmp (int32 ac, d10 *rs)
{
int32 p1 = ADDAC (ac, 1);
t_uint64 xh, xl, yh, yl, mid;
UFP a, b;

funpack (AC(ac), AC(p1), &a, AFRC);                     /* unpack operands */
funpack (rs[0], rs[1], &b, AFRC);
if ((a.fhi == 0) || (b.fhi == 0)) {                     /* either 0? result 0 */
    AC(ac) = AC(p1) = 0;
    return;
    }
a.sign = a.sign ^ b.sign;                               /* result sign */
a.exp = a.exp + b.exp - FP_BIAS + 1;                    /* result exponent */
xh = a.fhi >> 32;                                       /* split 62b fracs */
xl = a.fhi & MSK32;                                     /* into 32b halves */
yh = b.fhi >> 32;
yl = b.fhi & MSK32;
a.fhi = xh * yh;                                        /* hi xproduct */
a.flo = xl * yl;                                        /* low xproduct */
mid = (xh * yl) + (yh * xl);                            /* fits in 64b */
a.flo = a.flo + (mid << 32);                            /* add mid lo to lo */
a.fhi = a.fhi + ((mid >> 32) & MSK32) + (a.flo < (mid << 32));
fnorm (&a, FP_URNDD);                                   /* normalize, round */
AC(ac) = fpack (&a, &AC(p1), FALSE);                    /* pack result */
return;
}

/* Double precision floating divide
   This algorithm develops a full 62 bits of quotient, plus one rounding
   bit, in the low order 63b of a 64b number.  To do this, we must assure
   that the initial divide step generates a 1.  If it would fail, shift
   the dividend left and decrement the result exponent accordingly.
*/

void dfdv (int32 ac, d10 *rs)
{
int32 p1 = ADDAC (ac, 1);
int32 i;
t_uint64 qu = 0;
UFP a, b;

funpack (AC(ac), AC(p1), &a, AFRC);                     /* unpack operands */
funpack (rs[0], rs[1], &b, AFRC);
if (a.fhi >= 2 * b.fhi) {                               /* will divide work? */
    SETF (F_AOV | F_DCK | F_FOV | F_T1);
    return;
    }
if (a.fhi) {                                            /* dvd = 0? quo = 0 */
    a.sign = a.sign ^ b.sign;                           /* result sign */
    a.exp = a.exp - b.exp + FP_BIAS + 1;                /* result exponent */
    if (a.fhi < b.fhi) {                                /* make sure initial */
        a.fhi = a.fhi << 1;                             /* divide step will work */
        a.exp = a.exp - 1;
        }
    for (i = 0; i < 63; i++) {                          /* 63b of quotient */
        qu = qu << 1;                                   /* shift quotient */
        if (a.fhi >= b.fhi) {                           /* will div work? */
            a.fhi = a.fhi - b.fhi;                      /* sub, quo = 1 */
            qu = qu + 1;
            }
        a.fhi = a.fhi << 1;                             /* shift dividend */
        }
    a.fhi = qu;
    }
fnorm (&a, FP_URNDD);                                   /* normalize, round */
AC(ac) = fpack (&a, &AC(p1), FALSE);                    /* pack result */
return;
}

/* Unpack floating point operand */

void funpack (d10 h, d10 l, UFP *r, t_bool sgn)
{
d10 fphi, fplo;

r->sign = GET_FPSIGN (h);
r->exp = GET_FPEXP (h);
fphi = GET_FPHI (h);
fplo = GET_FPLO (l);
r->fhi = (fphi << FP_V_UFHI) | (fplo << FP_V_UFLO);
r->flo = 0;
if (r->sign) {
    r->exp = r->exp ^ FP_M_EXP;                         /* 1s comp exp */
    if (sgn) {                                          /* signed frac? */
        if (r->fhi) r->fhi = r->fhi | FP_UCRY;          /* extend sign */ 
        else {
            r->exp = r->exp + 1;
            r->fhi = FP_UCRY | FP_UNORM;
            }
        }
    else {                                              /* abs frac */
        if (r->fhi) r->fhi = UNEG (r->fhi) & FP_UFRAC;
        else {
                r->exp = r->exp + 1;
            r->fhi = FP_UNORM;
            }
        }
    }
return;
}

/* Normalize and optionally round floating point operand */
 
void fnorm (UFP *a, t_int64 rnd)
{
int32 i;
static t_uint64 normmask[6] = {
 0x6000000000000000, 0x7800000000000000, 0x7F80000000000000,
 0x7FFF800000000000, 0x7FFFFFFF80000000, 0x7FFFFFFFFFFFFFFF
 };
static int32 normtab[7] = { 1, 2, 4, 8, 16, 32, 63 };
extern a10 pager_PC;

if (a->fhi & FP_UCRY) {                                 /* carry set? */
    printf ("%%PDP-10 FP: carry bit set at normalization, PC = %o\n", pager_PC);
    a->flo = (a->flo >> 1) | ((a->fhi & 1) << 63);      /* try to recover */
    a->fhi = a->fhi >> 1;                               /* but root cause */
    a->exp = a->exp + 1;                                /* should be fixed! */
    }
if ((a->fhi | a->flo) == 0) {                           /* if fraction = 0 */
    a->sign = a->exp = 0;                               /* result is 0 */
    return;
    }
while ((a->fhi & FP_UNORM) == 0) {                      /* normalized? */
    for (i = 0; i < 6; i++) {
        if (a->fhi & normmask[i]) break;
        }
    a->fhi = (a->fhi << normtab[i]) | (a->flo >> (64 - normtab[i]));
    a->flo = a->flo << normtab[i];
    a->exp = a->exp - normtab[i];
    }
if (rnd) {                                              /* rounding? */
    a->fhi = a->fhi + rnd;                              /* add round const */
    if (a->fhi & FP_UCRY) {                             /* if carry out, */
        a->fhi = a->fhi >> 1;                           /* renormalize */
        a->exp = a->exp + 1;
        }
    }
return;
}

/* Pack floating point result */

d10 fpack (UFP *r, d10 *lo, t_bool fdvneg)
{
d10 val[2];

if (r->exp < 0) SETF (F_AOV | F_FOV | F_FXU | F_T1);
else if (r->exp > FP_M_EXP) SETF (F_AOV | F_FOV | F_T1);
val[0] = (((((d10) r->exp) & FP_M_EXP) << FP_V_EXP) |
    ((r->fhi & FP_UFHI) >> FP_V_UFHI)) & DMASK;
if (lo) val[1] = ((r->fhi & FP_UFLO) >> FP_V_UFLO) & MMASK;
else val[1] = 0;
if (r->sign) {                                          /* negate? */
    if (fdvneg) {                                       /* fdvr special? */
        val[1] = ~val[1] & MMASK;                       /* 1's comp */
        val[0] = ~val[0] & DMASK;
                }
    else { DMOVN (val); }                               /* 2's comp */
    }
if (lo) *lo = val[1];
return val[0];
}