/* KN10OPS.C - PDP-10 arithmetic operations
*/
/* $Id: kn10ops.c,v 2.3 2001/11/10 21:28:59 klh Exp $
*/
/*  Copyright © 1991, 1992, 1993, 2001 Kenneth L. Harrenstien
**  All Rights Reserved
**
**  This file is part of the KLH10 Distribution.  Use, modification, and
**  re-distribution is permitted subject to the terms in the file
**  named "LICENSE", which contains the full text of the legal notices
**  and should always accompany this Distribution.
**
**  This software is provided "AS IS" with NO WARRANTY OF ANY KIND.
**
**  This notice (including the copyright and warranty disclaimer)
**  must be included in all copies or derivations of this software.
*/
/*
 * $Log: kn10ops.c,v $
 * Revision 2.3  2001/11/10 21:28:59  klh
 * Final 2.0 distribution checkin
 *
 */

/*	Contains functions to provide 36-bit PDP-10 arithmetic operations
**	on a machine which does not support them.
**
** KN10OPS is primarily intended for use by the KLH10 PDP-10 emulator, but
** can also be used by KCC (PDP-10 C compiler) or other programs which want
** to emulate PDP-10 arithmetic operations.  Note there are a number of
** artifacts defined here for KCC which aren't needed by the KLH10.
**
** Supported PDP-10 machine ops are "op10*()" where * is:
**	AND/IOR/XOR/SETCM		Single word logical operations
**	MOVN/DMOVN			Negation
**	LSH/LSHC/ASH/ASHC/ROT/ROTC	Shifting
**	ADD/SUB/MUL/DIV/IDIV		Single precision fixed point
**	DADD/DSUB/DMUL/DDIV		Double precision fixed point
**	FAD/FSB/FMP/FDV			Single precision floating point
**	FADR/FSBR/FMPR/FDVR		Single precision floating point, Round
**	DFAD/DFSB/DFMP/DFDV		Double precision floating point
** plus:
**	NEG/DNEG	Negation, single & double fixed point
**	INC/DINC	Increment, single & double fixed point
**	CMP/DCMP	Comparison (CAM*) single & double
**	TST/DTST	Testing (SKIP*)
*/

/* By default this file assumes it is being compiled in its own right
** as a module for the KLH10 emulator program.
**
** However, it can also be included in another file to provide other
** programs with their own interface to PDP-10 operation functions.
** Examples of this are KCC (PDP-10 C Compiler) and OPTEST (a program
** to test these functions).
**	To do this, the macro OP10_INCLUDEFLAG must be defined prior to
**	inclusion.
**
** Note: KCC should only include this if being built on a non-10 platform
** as a cross-compiler to generate PDP-10 code.
*/

#ifndef OP10_INCLUDEFLAG
  /* Assume being compiled for KLH10 as its own module */
# include "klh10.h"	/* KLH10 Config params */
# include "kn10def.h"	/* Emulated machine defs */

#ifdef RCSID
 RCSID(kn10ops_c,"$Id: kn10ops.c,v 2.3 2001/11/10 21:28:59 klh Exp $")
#endif

# if 0			/* FYI: OP10_PCFSET definition from klh10.h */
#  define OP10_PCFSET(f) (cpu.mr_pcflags |= (f), \
	(((f)&(PCF_TR1|PCF_TR2)) ? INSBRKSET() : 0))
# endif

# if 0 /* KLH10_DEBUG */ /* OPTEST.C much better for debugging now */
    int op10debug = 0;
#  define OP10_DBGPRINT(list) if (op10debug) printf list
#  include <stdio.h>
# endif

    /* Determine what values are OK for IDIVI. */
# if KLH10_CPU_KLX		/* KLX accepts +1 */
#  define OP10_IDIVI_OK(w) op10m_skipg(w)
# elif KLH10_CPU_KS		/* KS accepts both +1 and -1 */
#  define OP10_IDIVI_OK(w) op10m_skipn(w)
# endif

    /* Determine whether to include G format instructions */
# ifndef OP10_GFMT
#  define OP10_GFMT KLH10_CPU_KL
# endif

    /* Determine behavior of FDV */
# ifndef OP10_FDV_KL
#  define OP10_FDV_KL (KLH10_CPU_KL || KLH10_CPU_KS)
# endif

#endif

/* Includer-defined macros.  The following should be defined appropriately
** by either kn10def.h or the includer, if desired.  They will default
** to something reasonable if not defined.
**
** OP10_PCFSET(flgs) - Must be defined if the facilities here are to return
**	PC flag results; the kn10ops.h file also checks for this in order to
**	define various things.
**	If not defined, no flag operations are compiled.
** OP10_DBGPRINTF(list) - Defined if debug output capability is desired.
**	If not defined, no debug code is compiled.
** OP10_GFMT - Defined 1 to say that G format ops should be included.
**	If 0 or undefined, none are included.
** OP10_KCCOPS - Defined 1 to include emulation of KCC C ops.
**	If 0 or undefined, not included.
** OP10_IDIVI_OK(w)  - Defined to specify action on a failing IDIVI; test
**	is made on divisor and returns TRUE if IDIVI is OK (ie should not
**	set no-divide and trap flags).
**	If not defined, acts like KA/KI and always fails for cases of 1 and +1.
** OP10_FDV_KL - Defined 1 to specify KL+KS behavior on FDV.
**	If 0 or undefined, uses KA+KI behavior.
*/

#include "kn10ops.h"	/* Get defs and macro facilities */

/* Resolve any defaults for OP10_ feature request macros */

#ifndef OP10_PCFSET		/* If not bothering with flags, */
# define IFFLAGS 0
# define OP10_PCFSET(a)		/* make this be no-op */
#else				/* Ah, wants flags! */
# define IFFLAGS 1
#endif

/* Debug output.  Invoke with nested parens, eg DEBUGPRF(("format", args)); */
#ifndef OP10_DBGPRINT
# define DEBUGPRF(list)		/* No debugging, so make no-op */
#else
# define DEBUGPRF(list) OP10_DBGPRINT(list)
# define DBGV_W(w) (long)LHGET(w),(long)RHGET(w)
# define DBGV_D(d) DBGV_W((d).HI),DBGV_W((d).LO)
# define DBGV_Q(q) DBGV_D((q).D0),DBGV_D((q).D1)
#endif

#ifndef OP10_GFMT		/* Default is not to include G-fmt stuff */
# define OP10_GFMT 0
#endif
#ifndef OP10_KCCOPS		/* Default is not to include KCC stuff */
# define OP10_KCCOPS 0
#endif


/* Determine whether DIV() or LDIV() is available from C */
#if defined(__STDC__) && __STDC__
# include <stdlib.h>		/* Get div, ldiv */

# ifdef WORD10_ITX32	/* Try to define DIV32 operation and result */
#  if WORD10_ITX32 == WORD10_TX_INT
#   define OP10_DIV32(n,d) div(n,d)
#   define OP10_DIV32T	div_t
#  elif WORD10_ITX32 == WORD10_TX_LONG
#   define OP10_DIV32(n,d) ldiv(n,d)
#   define OP10_DIV32T	ldiv_t
#  endif
# endif

# ifdef WORD10_ITX64	/* Try to define DIV64 operation and result */
#  if WORD10_ITX64 == WORD10_TX_INT
#   define OP10_DIV64(n,d) div(n,d)
#   define OP10_DIV64T	div_t
#  elif WORD10_ITX64 == WORD10_TX_LONG
#   define OP10_DIV64(n,d) ldiv(n,d)
#   define OP10_DIV64T	ldiv_t
#  endif
# endif

#endif

/* Defs for compatibility with older versions of this code */
#define HBITS H10BITS		/* Number of bits in halfword */
#define HSIGN H10SIGN		/* Sign bit for halfword */
#define HMASK H10MASK		/* Halfword Mask */
#define OP10XWD(l,r) XWDINIT(l,r)
#define OP10DXWD(h0,h1,h2,h3) DXWDINIT(h0,h1,h2,h3)
#define OP10QXWD(h0,h1,h2,h3,l0,l1,l2,l3) QXWDINIT(h0,h1,h2,h3,l0,l1,l2,l3)
#define HI w[0]
#define LO w[1]
#define D0 d[0]
#define D1 d[1]

/* Exported data */
w10_t w10one = OP10XWD(0,1);
w10_t w10mask = OP10XWD(HMASK,HMASK);	/* Word of all ones */
w10_t w10zero = OP10XWD(0,0);		/* Word of all zeros */
dw10_t dw10zero = OP10DXWD(0,0,0,0);

/* Internal data */
static w10_t w10maxneg = OP10XWD(HSIGN,0);
static dw10_t dw10one = OP10DXWD(0,0,0,1);
static dw10_t dw10maxneg = OP10DXWD(HSIGN,0,HSIGN,0); /* FixPt */
static qw10_t qw10maxneg = OP10QXWD(HSIGN,0,HSIGN,0,	   /* FixPt */
				    HSIGN,0,HSIGN,0);


static w10_t w10loobits[37] = {
# define lmask(n) (((uint18)1<<(n))-1)
# define WBMLO(n) W10XWD(0,lmask(n))
# define WBMHI(n) W10XWD(lmask((n)-18),H10MASK)

	WBMLO(0),  WBMLO(1),  WBMLO(2),  WBMLO(3),  WBMLO(4),  WBMLO(5),
	WBMLO(6),  WBMLO(7),  WBMLO(8),  WBMLO(9),  WBMLO(10), WBMLO(11),
	WBMLO(12), WBMLO(13), WBMLO(14), WBMLO(15), WBMLO(16), WBMLO(17),

	WBMHI(18), WBMHI(19), WBMHI(20), WBMHI(21), WBMHI(22), WBMHI(23),
	WBMHI(24), WBMHI(25), WBMHI(26), WBMHI(27), WBMHI(28), WBMHI(29),
	WBMHI(30), WBMHI(31), WBMHI(32), WBMHI(33), WBMHI(34), WBMHI(35),

	W10XWD(H10MASK, H10MASK)	/* Entry 36 is all ones */
# undef WBMLO
# undef WBMHI
# undef lmask
};

#if 0	/* Not needed at moment */

static w10_t w10hiobits[37] = {
# define amask(b) ((((uint18)1<<(HBITS-(b)))-1) ^ HMASK)
# define WBMLO(n) W10XWD(amask(n), 0)
# define WBMHI(n) W10XWD(H10MASK,amask((n)-18))

	WBMLO(0),  WBMLO(1),  WBMLO(2),  WBMLO(3),  WBMLO(4),  WBMLO(5),
	WBMLO(6),  WBMLO(7),  WBMLO(8),  WBMLO(9),  WBMLO(10), WBMLO(11),
	WBMLO(12), WBMLO(13), WBMLO(14), WBMLO(15), WBMLO(16), WBMLO(17),

	WBMHI(18), WBMHI(19), WBMHI(20), WBMHI(21), WBMHI(22), WBMHI(23),
	WBMHI(24), WBMHI(25), WBMHI(26), WBMHI(27), WBMHI(28), WBMHI(29),
	WBMHI(30), WBMHI(31), WBMHI(32), WBMHI(33), WBMHI(34), WBMHI(35),

	W10XWD(H10MASK, H10MASK)	/* Entry 36 is all ones */
# undef WBMLO
# undef WBMHI
};
static w10_t w10hizbits[37] = {
# define smask(b) (((uint18)H10MASK)>>(b))
# define WBMLO(n) W10XWD(smask(n), H10MASK)
# define WBMHI(n) W10XWD(0,smask((n)-18))

	WBMLO(0),  WBMLO(1),  WBMLO(2),  WBMLO(3),  WBMLO(4),  WBMLO(5),
	WBMLO(6),  WBMLO(7),  WBMLO(8),  WBMLO(9),  WBMLO(10), WBMLO(11),
	WBMLO(12), WBMLO(13), WBMLO(14), WBMLO(15), WBMLO(16), WBMLO(17),

	WBMHI(18), WBMHI(19), WBMHI(20), WBMHI(21), WBMHI(22), WBMHI(23),
	WBMHI(24), WBMHI(25), WBMHI(26), WBMHI(27), WBMHI(28), WBMHI(29),
	WBMHI(30), WBMHI(31), WBMHI(32), WBMHI(33), WBMHI(34), WBMHI(35),

	W10XWD(0, 0)	/* Entry 36 is all zeros */
# undef WBMLO
# undef WBMHI
};

/* Halfword bit masks */

static uint18 h10hiobits[H10BITS+1] = {
	amask(0), amask(1), amask(2), amask(3), amask(4), amask(5),
	amask(6), amask(7), amask(8), amask(9), amask(10), amask(11),
	amask(12), amask(13), amask(14), amask(15), amask(16), amask(17),
	H10MASK		/* Entry 18 is all ones */
};

#endif /* 0 */

/* Internal routines
**	Routines called x_name() are provided in order to carry out
**	arithmetic operations without the side-effects (e.g. flag setting)
**	of an actual instruction.
*/
static int adffo(dw10_t);	/* Arith Double FFO */
static qw10_t x_qneg(qw10_t);	/* Negate quadword */
static qw10_t x_qash1(qw10_t);	/* pseudo-ASH quadword left by 1 */
static int x_div(dw10_t *, w10_t);	/* Aux for single division */
static int ddivstep(dw10_t *, w10_t, int, int);	/* " " "     "     */
static dw10_t op10dfinc(dw10_t);	/* Increment DF int */
static qw10_t x_dmul(dw10_t, dw10_t);	/* For unsigned-double-int mult */
static w10_t ffadr(w10_t, w10_t, int);	/* Aux for FAD/FSB/FADR/FSBR */
static w10_t ffmpr(w10_t, w10_t, int);	/* Aux for FMP/FMPR */
static w10_t sfnorm(int, w10_t, int);	/* Float & normalize single-prec int */
static dw10_t dfad(dw10_t, dw10_t, int);	/* Aux for DFAD/DFSB */
static int qdivstep(qw10_t *, dw10_t, int);	/* For double division */

#if IFFLAGS
static void x_ashflg(w10_t, int);
static w10_t x_ash(w10_t, int);
static dw10_t x_ashc(dw10_t, int);
static dw10_t x_dadd(dw10_t, dw10_t);
static dw10_t x_dsub(dw10_t, dw10_t);
#else
# define x_ash(w,h)  op10ash(w, (h10_t)(long)imm8op(h))
# define x_ashc(d,h) op10ashc(d, (h10_t)(long)imm8op(h))
# define x_dadd(d,e) op10dadd(d,e)
# define x_dsub(d,e) op10dsub(d,e)
#endif


/* Certain instructions use a signed 8-bit immediate operand,
**	specifically ASH/C, LSH/C, ROT/C, and FSC.
** The sign bit is bit 18 (high bit of RH), applied to the low 8 bits.
** Note that a negative-signed 0 becomes -256, not 0!
** Also, these macros explicitly return an "int" result.  This is
** important as they are used to convert h10_t values into ints.
*/
#define imm8op(h)  (((h)&HSIGN) ? ((int)(h)|~0377) : ((int)(h)&0377))
#define imm8neg(h) (-((int)(h)|~0377))	/* Negate when known negative */


/* Simple macros, mainly for backward compatibility.
**	wskipl(w) - TRUE if word is negative (sign bit set).
**	wskipn(w) - TRUE if any bits in word are set.
**	wmagskipn(w) - Similar but ignores sign bit, checks magnitude only.
**	dskipn(d) - TRUE if any bits in doubleword are set
**	dlomagskipn(d) - Similar but ignores LOW word sign bit; HI sign
**			is still checked.
*/
#define wskipl(w)	op10m_skipl(w)
#define wskipn(w)	op10m_skipn(w)
#define wmagskipn(w)	op10m_magstst(w)
#define dskipn(d)	(wskipn((d).HI) || wskipn((d).LO))
#define dlomagskipn(d)	(wskipn((d).HI) || wmagskipn((d).LO))

/* Simple auxiliaries, primarily for KCC */

/* Convert unsigned base integer to w10_t */
w10_t op10utow(uint32 u)
{
    register w10_t w;

    LRHSET(w, ((u >> H10BITS) & H10MASK), (u & H10MASK));
    return w;
}

/* Convert w10_t to signed base integer */
int32 op10wtos(register w10_t w)
{
    return ((((int32)LHGET(w))|(~H10MASK)) << HBITS) | RHGET(w);
}


#if 0 /* Not used, and incorrect anyway */

w10_t op10stow(s)	/* Convert signed base integer to w10_t */
int32 s;		/* Defaults to signed! */
{
    register w10_t w;
    RHSET(w, s & HMASK);
    LHSET(w, (s >> HBITS) & HMASK);
    return w;
}
#endif


/* JFFO-type routines.
**	A return value of 0 means bit 0 (sign).
**	A value of 36 (or 72) means no bit found.
*/
/* Find First One.  Note returns int, not word! */
int op10ffo(register w10_t w)
{
    register int i;
    register uint18 reg;		/* Need unsigned for shifting */

    if (reg = LHGET(w)) i = 17;		/* Find right halfword and set up */
    else if (reg = RHGET(w)) i = 17+18;
    else return 36;

    while (reg >>= 1) --i;
    return i;

#if 0	/* Old version of this code */
    register uint32 reg;		/* Need unsigned for shifting */
    if (reg = LHGET(w)) i = 0;		/* Find right halfword and set up */
    else if (reg = RHGET(w)) i = 18;
    else return 36;
    while (!((reg <<= 1) & (H10SIGN<<1))) ++i;	/* Shift until found bit */
    return i;
#endif
}

/* ADFFO - Arithmetic Double Find First One.
**	IGNORES the low word sign bit, regardless of its setting.
**	Returns # of bits between high end and first one-bit.
** A return value of:
**	0    - sign bit was set
**	>=36 - high word clear, bit is in low word magnitude
**	>=71 - No bits were set in high word or non-sign bits of low word.
*/
static int adffo(register dw10_t d)
{
    register int i;
    register uint18 reg;		/* Need unsigned for shifting */

    if (reg = LHGET(d.HI)) i = 17;	/* Find right halfword and set up */
    else if (reg = RHGET(d.HI)) i = 17+18;
    else if (reg = (LHGET(d.LO)&H10MAGS)) i = 17+36-1;
    else if (reg = RHGET(d.LO)) i = 17+36+18-1;
    else return 36+36-1;

    while (reg >>= 1) --i;
    return i;
}

/* Logical operations and ones-complement negation */

w10_t op10and(register w10_t a, register w10_t b)
{
    op10m_and(a,b);
    return a;
}

w10_t op10ior(register w10_t a, register w10_t b)
{
    op10m_ior(a,b);
    return a;
}

w10_t op10xor(register w10_t a, register w10_t b)
{
    op10m_xor(a,b);
    return a;
}

w10_t op10setcm(register w10_t w)
{
    op10m_setcm(w);
    return w;
}

/* Negation - Single and Double */

/* MOVM - for completeness.  Returns magnitude, sets flags if
**	trying to make the max negative # positive.
*/
w10_t op10movm(register w10_t w)
{
    return op10m_skipl(w) ? op10movn(w) : w;
}

/* MOVN - If arg 0, set CRY0+CRY1.  If SETZ, set TRP1+OV+CRY1.
*/
w10_t op10movn(register w10_t w)
{
#if IFFLAGS
    op10mf_movn(w);
#else
    op10m_movn(w);
#endif
    return w;
}

/* DMOVN - Note that bit 0 of low word is always cleared, to follow
**	double floating-point convention.  Use DNEG to follow
**	the double fixed-point convention instead.
**	As for MOVN, if arg 0, set CRY0+CRY1; if SETZ, set TRP1+OV+CRY1.
*/
dw10_t
op10dmovn(register dw10_t d)
{
#if IFFLAGS
    op10mf_dmovn(d);
#else
    op10m_dmovn(d);
#endif
    return d;
}

/* QNEG - artificial "instr" to negate a quad-length fixed-point
**	value.  Follows fixed-point convention of copying sign bit
**	to low-order words.
*/
static qw10_t x_qneg(register qw10_t q)
{
    op10m_movn(q.D1.LO);	/* Negate lowest word */
    if (wmagskipn(q.D1.LO)) {
	op10m_setcm(q.D1.HI);	/* Low word has stuff, no carry into high */
	op10m_setcm(q.D0.LO);
	op10m_setcm(q.D0.HI);
    } else {
	op10m_movn(q.D1.HI);	/* Low word clear, so carry into high word */
	if (wmagskipn(q.D1.HI)) {	/* If has stuff, */
	    op10m_setcm(q.D0.LO);	/* no carry into higher words */
	    op10m_setcm(q.D0.HI);
	} else
	    op10m_dmovn(q.D0);	/* Lower word clear, so carry up */
    }
    /* Now make sign bit uniform */
    if (wskipl(q.D0.HI)) {
	op10m_signset(q.D0.LO);
	op10m_signset(q.D1.HI);
	op10m_signset(q.D1.LO);
    } else {
	op10m_signclr(q.D0.LO);
	op10m_signclr(q.D1.HI);
	op10m_signclr(q.D1.LO);
    }
    return q;
}

/* Shift instructions (Logical and Arithmetic) */

#ifndef OP10_USEMAC
#  define OP10_USEMAC 1
#endif
#if OP10_USEMAC
#  define LSHIFTM(w,s) op10m_lshift(w,s)
#  define RSHIFTM(w,s) op10m_rshift(w,s)
#else
#  define LSHIFTM(w,s) ((w) = op10lsh(w,s))
#  define RSHIFTM(w,s) ((w) = op10lsh(w,-(s)))
#endif


/* Logical shift */

w10_t op10lsh(register w10_t w,
	      register h10_t h)
{
    if (h & HSIGN) {		/* Right shift */
	h = imm8neg(h);		/* Get magnitude (was neg, so > 0) */
#if OP10_USEMAC
	op10m_rshift(w, h);
#else
	if (h > HBITS) {
	    if (h > HBITS*2) w.rh = 0;
	    else w.rh = (w.lh >> (h-HBITS));
	    w.lh = 0;
	} else {
	    w.rh = (w.rh >> h) | ((w.lh << (HBITS-h)) & HMASK);
	    w.lh >>= h;
	}
#endif
    } else {		/* Left shift */
	if (!(h &= 0377)) ;		/* Do nothing if 0 */
	else
#if OP10_USEMAC
	  op10m_lshift(w, h);
#else
	  if (h > HBITS) {
	    if (h > HBITS*2) w.lh = 0;
	    else w.lh = (w.rh << (h-HBITS)) & HMASK;
	    w.rh = 0;
	} else {
	    w.lh = ((w.lh << h)&HMASK) | (w.rh >> (HBITS-h));
	    w.rh = (w.rh << h) & HMASK;
	}
#endif
    }
    return w;
}

/* Double Logical shift */

dw10_t op10lshc(register dw10_t d,
		register h10_t h)
{
    if (h & HSIGN) {	/* Right shift */
	h = imm8neg(h);			/* Get magnitude (was neg, so > 0) */
	if (h >= HBITS*2) {
	    if (h > HBITS*4) op10m_setz(d.LO);
	    else {
		d.LO = d.HI;		/* Set up for macro */
		h -= HBITS*2;		/* See how much farther to go */
		RSHIFTM(d.LO, h);	/* Right-shift remaining distance */
	    }
	    op10m_setz(d.HI);
	} else {
	    register w10_t tmp;
	    tmp = d.HI;			/* Save old value of high wd */
	    RSHIFTM(d.HI, h);		/* Right-shift high by amount */
	    RSHIFTM(d.LO, h);		/* ditto for low wd */
	    h = HBITS*2 - h;		/* Find left-shift for lost hi bits */
	    LSHIFTM(tmp, h);		/* Move into position for IOR */
	    op10m_ior(d.LO, tmp);	/* IOR the bits back in */
	}
    } else {
	if (!(h &= 0377)) ;		/* Pos, do nothing if 0 */
	else if (h >= HBITS*2) {
	    if (h >= HBITS*4) op10m_setz(d.HI);
	    else {
		d.HI = d.LO;		/* Set up for macro */
		h -= HBITS*2;
		LSHIFTM(d.HI, h);	/* Do positive left-shift of word */
	    }
	    op10m_setz(d.LO);
	} else {
	    register w10_t tmp;
	    tmp = d.LO;			/* Save old val of low wd */
	    LSHIFTM(d.HI, h);		/* Shift high word left in place */
	    LSHIFTM(d.LO, h);		/* Shift low word left in place */
	    h = HBITS*2 - h;		/* Find right-shift for lost hi bits */
	    RSHIFTM(tmp, h);		/* Move into position for IOR */
	    op10m_ior(d.HI, tmp);	/* IOR the bits back */
	}
    }

    return d;
}


/* ASH, ASHC - If any bits of significance (a 1 if positive, 0 if negative)
**	are shifted out, set overflow and Trap 1.
*/

/* Full-word null-bit masks, indexed by # bits for a shift
**	Entry 1 has high bit set, 2 has high 2 bits, etc.
*/

static w10_t ashwmask[37] = {
# define amask(b) ((((uint18)1<<(HBITS-(b)))-1) ^ HMASK)
# define WBMLO(n) W10XWD(amask(n), 0)
# define WBMHI(n) W10XWD(H10MASK,amask((n)-18))

	WBMLO(0),  WBMLO(1),  WBMLO(2),  WBMLO(3),  WBMLO(4),  WBMLO(5),
	WBMLO(6),  WBMLO(7),  WBMLO(8),  WBMLO(9),  WBMLO(10), WBMLO(11),
	WBMLO(12), WBMLO(13), WBMLO(14), WBMLO(15), WBMLO(16), WBMLO(17),

	WBMHI(18), WBMHI(19), WBMHI(20), WBMHI(21), WBMHI(22), WBMHI(23),
	WBMHI(24), WBMHI(25), WBMHI(26), WBMHI(27), WBMHI(28), WBMHI(29),
	WBMHI(30), WBMHI(31), WBMHI(32), WBMHI(33), WBMHI(34), WBMHI(35),

	W10XWD(H10MASK, H10MASK)	/* Entry 36 is all ones */
};

/* Macro to do right-shift ASH of negative number, when shift is
**	known to be 0 <= s <= 35
*/
#define NASH_RSHIFTM(w, s) \
	(RSHIFTM(w, s), op10m_ior(w, (ashwmask+1)[s]))

#if IFFLAGS
	/* Similar table, but for halfwords only */
static int32 ashmask[HBITS] = {
	amask(0), amask(1), amask(2), amask(3), amask(4), amask(5),
	amask(6), amask(7), amask(8), amask(9), amask(10), amask(11),
	amask(12), amask(13), amask(14), amask(15), amask(16), amask(17)
};
#endif /* IFFLAGS */
#undef amask


#if IFFLAGS
/* Auxiliary to set flags if shift will lose a bit */

static void
x_ashflg(register w10_t w,
	 register int shift)
{
    if (wskipl(w)) {		/* If negative, */
	if (shift > 35) {	/* see if shifting entire word */
	    /* Always fail here, cuz 0-bits coming in from the right will
	    ** trigger overflow when they leave bit 1 of a negative number.
	    */
	    OP10_PCFSET(PCF_ARO|PCF_TR1);
	    return;
	}
	op10m_setcm(w);		/* invert so looking for 1-bits */
    }
    if (shift >= 35) {
	if (!wskipn(w))
	    return;
    } else if (shift >= 17) {
	if (!LHGET(w) && !(RHGET(w) & ashmask[shift-17]))
	    return;
    } else if (!(LHGET(w) & ashmask[shift+1]))
	return;
    OP10_PCFSET(PCF_ARO|PCF_TR1);
}
#endif /* IFFLAGS */

/* Arithmetic shift */

w10_t op10ash(register w10_t w, register h10_t h)
{
    register int i = imm8op(h);

#if IFFLAGS
    if (i > 0)
	x_ashflg(w, i);	/* Shifting left, so check for overflow */
    return x_ash(w, i);
}

/* Arithmetic shift */

static w10_t x_ash(register w10_t w,
		   register int i)
{
#endif /* IFFLAGS */
    register int32 r;	/* Must be signed! */

    if (i >= 0) {		/* Left shift -- similar to logical */
	i &= 0377;
	r = wskipl(w);		/* Remember sign bit */
	LSHIFTM(w, i);		/* Do the shift */
	if (r) op10m_signset(w);	/* Set sign to same as before */
	else op10m_signclr(w);
	return w;
    }
    /* Do right shift */
    i = -i;			/* Get magnitude (was neg, so > 0) */
    if ((r = LHGET(w)) & HSIGN)
	r |= ~HMASK;		/* Extend sign bit of LH */
    if (i >= HBITS) {
	if (i > (HBITS*2-1)) i = HBITS*2-1;
	RHSET(w, ((r >> (i-HBITS)) & HMASK));
	LHSET(w, (wskipl(w) ? HMASK : 0));
    } else {
	RHSET(w, ((RHGET(w) >> i) | ((LHGET(w) << (HBITS-i)) & HMASK)));
	LHSET(w, ((r >> i) & HMASK));
    }
    return w;
}


/* Double Arithmetic shift */

dw10_t op10ashc(register dw10_t d,
		register h10_t h)
{
    register int i = imm8op(h);

#if IFFLAGS
    if (i > 0) {
	if (i > 35) {
	    /* If high word will be completely lost, special-case it */
	    i -= 35;			/* Adjust for remaining code */
	    if (wskipl(d.HI)) {
		op10m_signset(d.LO);	/* Must copy hi sign */
		op10m_setcm(d.HI);
	    } else
		op10m_signclr(d.LO);	/* Must copy hi sign */

	    /* Now test high and low words for losing non-sign bits */
	    if (wskipn(d.HI))		/* Losing from high? */
		OP10_PCFSET(PCF_ARO|PCF_TR1);
	    else
		x_ashflg(d.LO, i);	/* If not, see if losing from low */

	    /* Flag testing done, do the shift */
	    d.HI = x_ash(d.LO, i);	/* Shift, plop result into high */
	    op10m_tlz(d.LO, H10MAGS);	/* Clear low except for sign */
	    RHSET(d.LO, 0);
	    return d;			/* and that's it */
	}
	x_ashflg(d.HI, i);	/* Losing less than a word, just check high */
    }
    return x_ashc(d, i);	/* Do negative or within-high positive shift */
}

/* Internal ASHC.
**	NOTE!!!  Shift argument has type signed int (not halfword!)
**	and no modulo operation is applied to it.  If one is needed it
**	must be done to the argument prior to the call!
*/
static dw10_t
x_ashc(/*register*/ dw10_t d,
       register int i)
{
#endif /* IFFLAGS */

    register int r;

    if (i >= 0) {		/*  Left shift -- similar to logical */
	if (i == 0) return d;	/* If no shift, no change to operand */
	r = wskipl(d.HI);		/* Remember sign bit */
	if (i >= (HBITS*2-1)) {
	    if (i >= (HBITS*4)-2) op10m_setz(d.HI);
	    else {
		i -= (HBITS*2-1);
		d.HI = d.LO;		/* Set up for macro */
		LSHIFTM(d.HI, i);	/* Do left-shift */
	    }
	    op10m_setz(d.LO);
	} else {
	    op10m_signclr(d.LO);	/* Ensure low sign is clear */
	    d.HI = op10ior(op10lsh(d.HI, i), op10lsh(d.LO, i-(HBITS*2-1)));
	    LSHIFTM(d.LO, i);		/* Shift left in place */
	}
	if (r) {			/* Restore original sign */
	    op10m_signset(d.HI);
	    op10m_signset(d.LO);
	} else {
	    op10m_signclr(d.HI);
	    op10m_signclr(d.LO);
	}
    } else {		/* Right shift -- bring in sign bit */
	r = wskipl(d.HI);		/* Remember sign bit */
	i = -i;				/* Get magnitude */
	if (i > (HBITS*2-1)) {
	    if (i > (HBITS*4-2)) d.LO = r ? w10mask : w10zero;
	    else d.LO = x_ash(d.HI, (HBITS*2-1)-i);
	    d.HI = r ? w10mask : w10zero;
	} else {
	    op10m_signclr(d.LO);	/* Ensure low sign is clear */
	    d.LO = op10ior(op10lsh(d.LO, -i), x_ash(d.HI, (HBITS*2-1)-i));
	    d.HI = x_ash(d.HI, -i);
	    if (r)
		op10m_signset(d.LO);	/* Set low sign if high is set */
	}
    }
    return d;
}

/* X_QASH1(q) - Special semi-arithmetic shift for quadword
**	Arith-shifts quadword left by 1, but does LSH in high word.
*/
static qw10_t x_qash1(register qw10_t q)
{
    LSHIFTM(q.D1.LO, 1);	/* Shift into sign bits */
    LSHIFTM(q.D1.HI, 1);
    LSHIFTM(q.D0.LO, 1);
    LSHIFTM(q.D0.HI, 1);
    if (wskipl(q.D1.LO)) {		/* Now carry anything necessary */
	op10m_signclr(q.D1.LO);
	op10m_iori(q.D1.HI, 1);
    }
    if (wskipl(q.D1.HI)) {
	op10m_signclr(q.D1.HI);
	op10m_iori(q.D0.LO, 1);
    }
    if (wskipl(q.D0.LO)) {
	op10m_signclr(q.D0.LO);
	op10m_iori(q.D0.HI, 1);
    }
    return q;
}


#define rothalf(hi, lo, x) /* Shift X bits of LO up into HI */ \
	((((hi) << (x)) | ((lo) >> (HBITS - (x)))) & HMASK)

/* Logical Rotate */

w10_t op10rot(register w10_t w, register h10_t h)
{
    register h10_t tmp;

    /* Canonicalize shift arg to positive value */
    if ((h = imm8op(h)) < 0)
	h = (HBITS*2) - ((-h) % (HBITS*2));	/* Right shift */
    else h %= (HBITS*2);			/* Left shift */

    if (h >= HBITS) {		/* Rotating more than one halfword? */
	tmp = LHGET(w);		/* Swap halves and reduce shifting */
	LHSET(w, RHGET(w));
	RHSET(w, tmp);
	h -= HBITS;
    }
    if (h > 0) {
	tmp = LHGET(w);			/* Save old LH */
	/* Shift bits up into LH from RH */
	LHSET(w, rothalf(LHGET(w), RHGET(w), h));
	/* Shift bits up into RH from temp */
	RHSET(w, rothalf(RHGET(w), tmp,  h));
    }
    return w;
}

/* Double Logical Rotate */

dw10_t op10rotc(register dw10_t d,
		register h10_t h)
{
    register w10_t wtmp;
    register h10_t tmp;

    /* Canonicalize shift arg to positive value */
    if ((h = imm8op(h)) < 0)
	h = (HBITS*4) - ((-h) % (HBITS*4));	/* Right shift */
    else h %= (HBITS*4);			/* Left shift */

    if (h >= HBITS*2) {		/* Rotating more than one word? */
	wtmp = d.HI;		/* Swap words and reduce shifting */
	d.HI = d.LO;
	d.LO = wtmp;
	h -= HBITS*2;
    }
    if (h >= HBITS) {		/* Rotating more than one halfword? */
	tmp = LHGET(d.HI);		/* Move halves and reduce shifting */
	LHSET(d.HI, RHGET(d.HI));
	RHSET(d.HI, LHGET(d.LO));
	LHSET(d.LO, RHGET(d.LO));
	RHSET(d.LO, tmp);
	h -= HBITS;
    }
    if (h > 0) {
	tmp = LHGET(d.HI);		/* Save high halfwd */
	/* Shift into LH from RH */
	LHSET(d.HI, rothalf(LHGET(d.HI), RHGET(d.HI), h));
	/* Shift into RH from lo LH */
	RHSET(d.HI, rothalf(RHGET(d.HI), LHGET(d.LO), h));
	/* Shift into lo LH fm lo RH */
	LHSET(d.LO, rothalf(LHGET(d.LO), RHGET(d.LO), h));
	/* Shift into lo RH fm temp */
	RHSET(d.LO, rothalf(RHGET(d.LO), tmp, h));
    }
    return d;
}

/* CIRC - Special ITS instruction!
**	Like ROTC except bits are shifted in and out of 2nd AC in the
**	opposite direction, thus CIRC 1,44 leaves in AC 2 the reverse of
**	what AC 1 contained.
** (Alan: CIRC behaves exactly like ROTC except you reverse AC+1 before and after the
	operation.  The behavior when count > 36. may well differ between
	KA/KI/KL/KS, but will be consistent with LSH, ROT, etc.)
*/
static h10_t hrev(h10_t);

/* Double Logical Circulate */

dw10_t op10circ(register dw10_t d,
		register h10_t h)
{
    register w10_t wtmp;
    register h10_t tmp;

    /* Canonicalize shift arg to positive value */
    if ((h = imm8op(h)) < 0)
	h = (HBITS*4) - ((-h) % (HBITS*4));	/* Right shift */
    else h %= (HBITS*4);			/* Left shift */

    /* Reduce to halfword case by circ'ing halfword chunks */
    if (h >= HBITS*3) {		/* Circ 54+n where (0 <= n < 18) */
	tmp = RHGET(d.HI);
	RHSET(d.HI, LHGET(d.HI));
	LHSET(d.HI, hrev((h10_t) LHGET(d.LO)));
	LHSET(d.LO, RHGET(d.LO));
	RHSET(d.LO, hrev(tmp));
	h -= HBITS*3;
    } else if (h >= HBITS*2) {	/* Circ 36+n where (0 <= n < 18) */
	wtmp = d.HI;
	LHSET(d.HI, hrev((h10_t) RHGET(d.LO)));
	RHSET(d.HI, hrev((h10_t) LHGET(d.LO)));
	LHSET(d.LO, hrev((h10_t) RHGET(wtmp)));
	RHSET(d.LO, hrev((h10_t) LHGET(wtmp)));
	h -= HBITS*2;
    } else if (h >= HBITS) {	/* Circ 18+n where (0 <= n < 18) */
	tmp = LHGET(d.HI);
	LHSET(d.HI, RHGET(d.HI));
	RHSET(d.HI, hrev((h10_t) RHGET(d.LO)));
	RHSET(d.LO, LHGET(d.LO));
	LHSET(d.LO, hrev(tmp));
	h -= HBITS;
    }

    if (h > 0) {		/* Circ n where (0 <= n < 18) */
	tmp = LHGET(d.HI);			/* Save high halfwd */
						/* Shift into      from */
	LHSET(d.HI, rothalf(LHGET(d.HI), RHGET(d.HI), h));
						/* hi LH      hi RH */
	RHSET(d.HI, rothalf(RHGET(d.HI), hrev((h10_t) RHGET(d.LO)), h));
						/* hi RH  rev lo RH */
	RHSET(d.LO, rothalf(LHGET(d.LO), RHGET(d.LO), 18-h));
						/* lo RH      lo LH */
	LHSET(d.LO, rothalf(hrev(tmp), LHGET(d.LO), 18-h));
						/* lo LH  rev hi LH */
    }
    return d;
}

/* HREV - reverse bits in halfword */
static h10_t hrev(register h10_t h)
{
    register h10_t r;
    register int i;

    if (!h) return 0;
    for (r = 0, i = 18; --i >= 0; ) {
	r = (r << 1) | (h & 01);	/* Transfer low bit to high */
	if (!(h >>= 1))			/* Check for early escape */
	    return r << i;
    }
    return r;
}
#undef rothalf

/* Fixed-point PC Flag facilities */

/*
	The algorithm used here to set the carry and overflow flags
is basically table-driven by the three sign bits of the two operands and
their result.
	Subtraction uses a different table; although one can get a correct
value simply by negating the second operand and adding it in, that procedure
does not set the flags properly.

	A + B = S	C0 C1	OVFL + TRAP1
	---------	-----	------------
[0]	0  0  0		-  -
[1]	0  0  1		-  X	AR0+TR1
[2]	0  1  0		X  X
[3]	0  1  1		-  -
[4]	1  0  0		X  X
[5]	1  0  1		-  -
[6]	1  1  0		X  -	AR0+TR1
[7]	1  1  1		X  X

	A - B = S	C0 C1	OVFL + TRAP1
	---------	-----	------------
[0]	0  0  0		X  X
[1]	0  0  1		-  -
[2]	0  1  0		-  -
[3]	0  1  1		-  X	AR0+TR1
[4]	1  0  0		X  -	AR0+TR1
[5]	1  0  1		X  X
[6]	1  1  0		X  X
[7]	1  1  1		-  -

*/


/* Define common macro for setting result flags of add operations.
** Requires that the original operands still be available.
*/
#define ADDFLAGS(sum, a, b) \
    if (wskipl(a)) {					\
	if (wskipl(b)) {			/* If both summands neg */\
	    if (wskipl(sum))			/* Always have a carry */\
		OP10_PCFSET(PCF_CR0+PCF_CR1);		\
	    else				/* Sum pos, overflowed */\
		OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_CR0);	\
	} else if (!wskipl(sum))		/* Mixed signs, so expect 1 */\
	    OP10_PCFSET(PCF_CR0+PCF_CR1);	/* Sum pos, carried out */\
    } else {						\
	if (!wskipl(b)) {			/* If both summands pos */\
	    if (wskipl(sum))			/* Sum shd be positive too */\
		OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_CR1);	\
	} else if (!wskipl(sum))		/* Mixed signs, so expect 1 */\
	    OP10_PCFSET(PCF_CR0+PCF_CR1);	/* Sum pos, carried out */\
    }

/* Define common macro for setting result flags of subtract operations.
** Requires that the original operands still be available.
*/
#define SUBFLAGS(sum, a, b) \
    if (!wskipl(a)) {				/* Check 0-3? */\
	if (!wskipl(b)) {			/* Check 0-1? */\
	    if (!wskipl(sum))			/* */\
		OP10_PCFSET(PCF_CR0+PCF_CR1);	/* [0] */\
	} else {				/* [1] is nop */\
	    if (wskipl(sum))			/* */\
		OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_CR1);	/* [3] */\
	}					/* [2] is nop */\
    } else {					/* Check 4-7 */\
	if (!wskipl(b)) {			/* Check 4-5 */\
	    if (!wskipl(sum))			/* [4] or [5] */\
		OP10_PCFSET(PCF_CR0+PCF_ARO+PCF_TR1);	/* [4] */\
	    else				\
		OP10_PCFSET(PCF_CR0+PCF_CR1);	/* [5] */\
	} else {				\
	    if (!wskipl(sum))			/* [6] or [7] */\
		OP10_PCFSET(PCF_CR0+PCF_CR1);	/* [6] */\
	}					/* [7] is nop */\
    }


/* Fixed-point Single-precision arithmetic */


w10_t op10add(register w10_t a, register w10_t b)
{
#if IFFLAGS
    register w10_t sum;

    sum = a;
    op10m_add(sum, b);		/* Add the words */
    ADDFLAGS(sum, a, b)		/* Invoke macro to set PC flags */
    return sum;
#else
    op10m_add(a,b);
    return a;
#endif
}

w10_t op10inc(register w10_t a)
{
    op10mf_inc(a);
    return a;
}

w10_t op10sub(register w10_t a, register w10_t b)
{
#if IFFLAGS
    register w10_t sum;

    sum = a;
    op10m_sub(sum, b);		/* Subtract the words */
    SUBFLAGS(sum, a, b)		/* Invoke macro to set PC flags */
    return sum;
#else
    op10m_sub(a,b);
    return a;
#endif
}

w10_t op10imul(register w10_t a, register w10_t b)
{
    register dw10_t d;
    register int sign;

    if (sign = wskipl(a))	/* Make args positive, remember signs */
	op10m_movn(a);
    if (wskipl(b)) {
	sign = !sign;
	op10m_movn(b);
#if IFFLAGS
	/* Check for screw case */
	if (wskipl(b) && wskipl(a)) {	/* If max neg * max neg, */
	    OP10_PCFSET(PCF_ARO+PCF_TR1);
	    return b;			/* Return low word (max neg) */
	}
#endif
    }
    d = op10xmul(a, b);		/* Do unsigned multiply, get double product */

    /* For IMUL, result must fit into one word, or overflow is set.
    ** Can check most efficiently here, while value is still positive,
    ** for a high-word of 0 - this will win 99.99% of the time.
    ** Just need to be careful of the 1 ? 0 case, which will fit if negated.
    */
#if IFFLAGS
    if (wskipn(d.HI)
      && (!sign || op10m_cain(d.HI, 1) || wskipn(d.LO)))
	OP10_PCFSET(PCF_ARO+PCF_TR1);
#endif

    if (sign)			/* Negate result if necessary */
	op10m_dneg(d);		/* Yup, negate in place, fixed-point style */
    return d.LO;		/* Return low-order result */
}

dw10_t op10mul(register w10_t a, register w10_t b)
{
    register dw10_t d;
    register int sign;

    if (sign = wskipl(a))
	op10m_movn(a);
    if (wskipl(b)) {
	sign = !sign;
	op10m_movn(b);
#if IFFLAGS
	if (wskipl(b) && wskipl(a)) {
	    OP10_PCFSET(PCF_ARO+PCF_TR1);
	    return dw10maxneg;		/* Return max double negative */
	}
#endif
    }
    d = op10xmul(a, b);		/* Do the multiply */
    if (sign)			/* Negate result if necessary */
	op10m_dneg(d);		/* Yup, negate in place, fixed-point style */
    return d;
}

#define DIGBTS 16
#define DIGMSK (((uint32)1<<DIGBTS)-1)
#define NDIGS(a) (((a)+DIGBTS-1)/DIGBTS)

/* OP10XMUL - Special unsigned multiplication of 36-bit integers.
**	Returns a 71-bit product, with the low order sign bit clear.
*/
dw10_t op10xmul(register w10_t a, register w10_t b)
{
    register int ai, bi;
    register dw10_t d;
    uint32 av[NDIGS(36)], bv[NDIGS(36)], pv[2*NDIGS(36)];

    /* First get args into vectors of N-bit digit integers */
#define setvec(w,v) \
	v[0] = LHGET(w) >> (HBITS-(36-(2*DIGBTS))); \
	v[1] = ((LHGET(w) << (18-DIGBTS)) & DIGMSK) | (RHGET(w) >> DIGBTS); \
	v[2] = RHGET(w) & DIGMSK; \
	DEBUGPRF(("Array: %lo %lo %lo\n", (long)v[0], (long)v[1], (long)v[2]));
    setvec(a, av);
    setvec(b, bv);
/*    pv[0] = pv[1] = pv[2] = */ pv[3] = pv[4] = pv[5] = 0;

    /* Now multiply the vectors together.  Because we have 32-bit unsigned
    ** integers, we can hold the result of multiplying 2 16-bit integers.
    */
    for (ai = NDIGS(36); --ai >= 0;) {
	register uint32 hicarry = 0;
	for (bi = NDIGS(36); --bi >= 0;) {
#if 0
	    register uint32 prod, high;
	    prod = av[ai] * bv[bi];
	    high = prod >> DIGBTS;
	    prod &= DIGMSK;
	    DEBUGPRF(("%d*b%d: %lo,%lo + %lo + p%d => ",
				ai,bi,high,prod, hicarry, ai+bi+1));
	    hicarry += prod + pv[ai+bi+1];
	    pv[ai+bi+1] = hicarry & DIGMSK;
	    hicarry = (hicarry>>DIGBTS) + high;
	    DEBUGPRF(("p%d: %lo, carry %lo\n", ai+bi+1, pv[ai+bi+1], hicarry));
#else
	    hicarry += av[ai] * bv[bi] + pv[ai+bi+1];
	    pv[ai+bi+1] = hicarry & DIGMSK;
	    hicarry >>= DIGBTS;
#endif
	}
	pv[ai] = hicarry;
    }

    /* Now put together doubleword from 71-bit product vector.
    ** Have to leave low-order sign bit clear.
    */
    DEBUGPRF(("ArOut: %lo %lo %lo %lo %lo %lo\n",
	(long)pv[0], (long)pv[1], (long)pv[2],
	(long)pv[3], (long)pv[4], (long)pv[5]));

#if DIGBTS==16
    RHSET(d.LO, (pv[4] & ((1<<2)-1))<< (18-2) |  pv[5] );
    LHSET(d.LO, (pv[3] & ((1<<3)-1))<< (17-3) | (pv[4] >> 2) );
    RHSET(d.HI, (pv[2] & ((1<<5)-1))<< (18-5) | (pv[3] >> 3) );
    LHSET(d.HI, (pv[1] & ((1<<7)-1))<< (18-7) | (pv[2] >> 5) );
#elif DIGBTS==15
    RHSET(d.LO, (pv[4] & ((1<< 3)-1))<< (18-3)  |  pv[5] );
    LHSET(d.LO, (pv[3] & ((1<< 5)-1))<< (17-5)  | (pv[4] >> 3) );
    RHSET(d.HI, (pv[2] & ((1<< 8)-1))<< (18-8)  | (pv[3] >> 5) );
    LHSET(d.HI, (pv[1] & ((1<<11)-1))<< (18-11) | (pv[2] >> 8) );
#endif
    return d;
}

/* IDIV and DIV.  Note that if these routines are used for instruction
**	emulation, some provision must be made for aborting the instruction
**	if No Divide is set -- so that I/DIVM and I/DIVB won't clobber the
**	memory operand.
** X(I)DIV exists for this purpose.  Unlike the usual OP10x routines, its first
**	arg is a pointer (through which it fetches operand and stores value)
**	and it returns a zero int if the division failed.
** IDIV may need to do a special check for 1<<35 / -1 depending on how
**	accurately we want to emulate specific processors.
**	The KS10 result is 1<<35 with no error; others produce overflow.
**	On KA, KI, and some KLs (1<<35)/1 also overflows.  See proc ref man.
*/
dw10_t op10idiv(w10_t a, w10_t b)
{
    dw10_t d;
    d.LO = a;
    d.HI = wskipl(a) ? w10mask : w10zero;
    if (!x_div(&d, b)) {
	d.HI = a;
	op10m_setz(d.LO);
#ifdef OP10_IDIVI_OK
	if (!OP10_IDIVI_OK(b))
#endif
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_DIV);
    }
    return d;
}

/* Note returns doubleword! High is quotient */
dw10_t op10div(
	       dw10_t d,	/* Dividend */
	       register w10_t w) /* Divisor */
{
    if (!x_div(&d, w))
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_DIV);
    return d;
}

int op10xidiv(dw10_t *ad,
	      w10_t a,
	      w10_t b)
{
    register int r;

    ad->LO = a;
    ad->HI = wskipl(a) ? w10mask : w10zero;
    if (r = x_div(ad, b))
	return r;	/* Success */

    /* Failure is only possible if A is 1<<35 and B is 1, 0, or -1.
    ** B=0 always causes No Divide.
    ** B=1  is OK on the KS and multi-section KLs.
    ** B=-1 is OK on the KS only.
    */
#ifdef OP10_IDIVI_OK
    if (OP10_IDIVI_OK(b)) {
	ad->HI = a;
	ad->LO = w10zero;
	return 1;
    }
#endif /* OP10_IDIVI_OK */
    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_DIV);
    return 0;
}

int op10xdiv(dw10_t *ad,	/* Dividend */
	     register w10_t w)	/* Divisor */
{
    register int r;
    return (r = x_div(ad, w)) ? r :
#if IFFLAGS
			(OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_DIV), 0);
#else
			(0);
#endif
}

/* X_DIV - Internal single-precision division.
**
** NOTE: There is a discrepancy between the PRM and what an actual KLX does.
**	It appears that if the resulting quotient would be the largest
**	negative number (400000,,0) then certain operands are permitted
**	without causing a no-divide condition.
**
**	Specifically, args like [1 ? 0] / [-1] result in SETZ even
**	though this is explicitly prohibited by PRM p.2-14 which says
**	"If the high order word of the magnitude of the double length
**	number in AC,AC+1 is greater than or *EQUAL TO* the magnitude of
**	the operand specified by M, set Trap 1, ..."
**
**	This sort of makes sense as the result is arithmetically correct;
**	the max neg # is one less in magnitude than the max positive # so
**	this one case is allowed to escape.
**
** The crucial code is in ddivstep where it checks for the above condition.
** Unfortunately at that point it doesn't know anything about the ultimate
** sign of the result and so can't tell whether to special-case equality
** or not.  At least, it can't without the stupid extra argument passed
** as the "qsign" parameter.
**
** FDVR, which uses ddivstep also, doesn't seem to suffer from this
** inconsistency.  DDIV, which uses qdivstep but otherwise is conceptually
** identical, doesn't either -- that is, it follows the PRM.  Go figure.
*/
/* Note returns doubleword! High is quotient */
static int x_div(
		 dw10_t *ad,	/* Dividend */
	 register w10_t w)	/* Divisor */
{
    dw10_t d;
    int numsign, densign;

    d = *ad;			/* Don't touch original arg in case we fail */

    /* Make args positive.  Note that if either arg is the max neg number
    ** the result will be that same number.  Don't bother checking here
    ** as the initial dividend/divisor comparison test should catch them.
    */
    if (numsign = wskipl(d.HI))
	op10m_dmovn(d);		/* Make positive (unless 1<<35) */
    if (densign = wskipl(w))
	op10m_movn(w);		/* Get absolute value, ignore max neg */

    /* Check for zero divisor, otherwise initial test could pass if
    ** dividend high word was 1<<35.  Sigh.
    */
    if (!wskipn(w))
	return 0;		/* Fail, no divide */

    /* Quick method - see if can do division in native arithmetic */
    if (!wskipn(d.HI) && !(LHGET(d.LO)&(~(HMASK>>4)))
      && !(LHGET(w)&(~(HMASK>>4)))) {
	register uint32 num, den;
	num = op10wtos(d.LO);
	den = op10wtos(w);
	d.HI = op10utow(num/den);
	d.LO = op10utow(num%den);
	DEBUGPRF(("idivquick: %o/%o Q: %lo,%lo R: %lo,%lo\n",
		numsign, densign, DBGV_W(d.HI), DBGV_W(d.LO)));
    } else {
	/* Ugh, hack big numbers */
	/* Divide them, get double result */
	if (!ddivstep(&d, w, 35, numsign != densign))
	    return 0;			/* Fail, no divide */

	/* Fix up remainder */
	if (wskipl(d.LO))		/* If final remainder was negative */
	    op10m_add(d.LO, w);		/* add divisor back in */
	DEBUGPRF(("idivslow: %o/%o Q: %lo,%lo R: %lo,%lo\n",
		numsign, densign, DBGV_W(d.HI), DBGV_W(d.LO)));
    }

    /* Unsigned division done, now fix up results */
    if (numsign)		/* Remainder has same sign as dividend */
	op10m_movn(d.LO);	/* Store remainder */

    /* Now fix up quotient if dividend & divisor had different sign */
    if (numsign != densign)
	op10m_movn(d.HI);
    *ad = d;
    return 1;
}

static int
ddivstep(dw10_t *ad,
	 register w10_t w,
	 register int nmagbits,
	 int qsign)	/* TRUE if quotient will ultimately be negative */
{
    register dw10_t d;
    register w10_t q;
    register int qbit;

    d = *ad;

    DEBUGPRF(("ddivstep(%d) D: %lo,%lo,%lo,%lo W: %lo,%lo qsign %s\n",
		nmagbits, DBGV_D(d), DBGV_W(w), (qsign ? "TRUE" : "FALSE")));

    /* Apply initial test - high order wd must be less than abs val of 
    ** divisor, else overflow & no divide.
    */
    op10m_sub(d.HI, w);		/* Subtract divisor from high word */
    if (!wskipl(d.HI)) {	/* See if divisor was larger or equal */
	/* Ugh, enter Land of the Special Case.
	** According to the PRM this should always cause No Divide, but
	** on the KLX at least, this is okay provided the final quotient
	** has the maximum negative value (400000,,0).
	** This is the only reason for "qsign".  Sigh.
	** Note that for this case, the remainder will always be the low
	** word as-is (with sign cleared, so fixup done by x_div will
	** do the right thing).
	*/
	if (qsign && !wskipn(d.HI)) {
	    op10m_signclr(d.LO);	/* Flush low sign for test & rem */
	    /* Now one last check - can't have any more quotient bits,
	    ** so make sure rest of dividend is less than divisor, i.e. it
	    ** all goes into the remainder.
	    */
	    if (op10m_ucmpl(d.LO, w)) {
		LHSET(d.HI, H10SIGN);	/* Make quotient be max neg #! */
		*ad = d;
		return 1;		/* And return success! */
	    }
	}
	return 0;		/* Overflow, no divide */
    }

    /* First quotient bit (future sign bit) is always 0 at this point. */
    qbit = 0;			/* Set first quotient bit */
    op10m_setz(q);		/* Clear current & future signs */
    LSHIFTM(d.LO, 1);		/* Squeeze out low sign bit */

    while (--nmagbits >= 0) {		/* Loop N times for N magnitude bits */
	/* Do a LSHC d,1 to shift new bit into dividend */
	LSHIFTM(d.HI, 1);
	if (wskipl(d.LO))		/* If high bit set in low wd, */
	    op10m_iori(d.HI, 1);	/* shift into low bit of high */
	LSHIFTM(d.LO, 1);

	if (qbit) op10m_sub(d.HI, w);	/* Sub or add, per previous qbit */
	else op10m_add(d.HI, w);
	qbit = (wskipl(d.HI) ? 0 : 1);	/* Get new quotient bit */
	LSHIFTM(q, 1);			/* Shift quotient over, add bit */
	op10m_iori(q, qbit);
	DEBUGPRF(("ddivstep %d: D: %lo,%lo,%lo,%lo Q: %lo,%lo \n",
		nmagbits, DBGV_D(d), DBGV_W(q)));
#if 0
	/* Minor speed hack; works, but worth the overhead of zero-testing? */
	if (!wskipn(d.HI) && !wskipn(d.LO)) {
	    /* Special copout if hit zero early */
	    LSHIFTM(d.HI, nmagbits);
	    LSHIFTM(q, nmagbits);
	    break;
	}
#endif
    }

    ad->LO = d.HI;	/* Return quotient & remainder in right places */
    ad->HI = q;
    return 1;
}

/* Fixed-point Double-precision arithmetic */

dw10_t op10dadd(register dw10_t da, register dw10_t db)
{
#if IFFLAGS
    register w10_t sum;

    op10m_signclr(da.LO);	/* Must ignore signs of low words */
    op10m_signclr(db.LO);
    op10m_add(da.LO, db.LO);	/* Add low B to low A */

    /* The following is a near-duplicate of the code for op10add.
    ** It has to be replicated, instead of just calling op10add, in order to
    ** properly handle the carry from the low word addition without
    ** losing track of what the flags should be.
    */
    sum = da.HI;		/* Set up high sum */
    if (wskipl(da.LO))		/* If low sign set, carry to high sum! */
	op10m_inc(sum);
    op10m_add(sum, db.HI);	/* Add high words */

    ADDFLAGS(sum, da.HI, db.HI)	/* Check for and set PC flags */

    da.HI = sum;		/* OK, put sum back in place */
    if (wskipl(da.HI))		/* Copy sign bit into low word of result */
	op10m_signset(da.LO);
    else op10m_signclr(da.LO);
    return da;
}

static dw10_t
x_dadd(register dw10_t da, register dw10_t db)
{
#endif /* IFFLAGS */

    op10m_signclr(da.LO);	/* Must ignore signs of low words */
    op10m_signclr(db.LO);
    op10m_add(da.LO, db.LO);	/* Add low B to low A */
    if (wskipl(da.LO))		/* If sign set, carry to high */
	op10m_inc(da.HI);
    op10m_add(da.HI, db.HI);	/* Add high words */
    if (wskipl(da.HI))		/* Copy sign bit into low word of result */
	op10m_signset(da.LO);
    else op10m_signclr(da.LO);
    return da;
}

/* OP10DINC - Increment double-prec fixed-point. */

dw10_t op10dinc(register dw10_t d)
{
    op10m_inc(d.LO);
    if (!wmagskipn(d.LO)) {	/* Check for carry to high word */
	op10m_inc(d.HI);			/* Yep, bump high word */
	LHSET(d.LO, (LHGET(d.HI) & HSIGN));	/* Copy sign bit to low */
    }
    return d;
}

/* OP10DFINC - Increment double-prec floating-point (internal)
*/
static dw10_t op10dfinc(register dw10_t d)
{
    op10m_inc(d.LO);
    if (!wmagskipn(d.LO)) {	/* Check for carry to high word */
	op10m_signclr(d.LO);	/* Yep, clear low sign */
	op10m_inc(d.HI);	/* And, bump high word */
    }
    return d;
}


dw10_t op10dsub(register dw10_t da, register dw10_t db)
{
#if IFFLAGS
    register w10_t sum;

#if 1
    op10m_signclr(da.LO);	/* Must ignore signs of low words */
    op10m_signclr(db.LO);
    op10m_sub(da.LO, db.LO);	/* Sub low B from low A */
    sum = da.HI;		/* Set up high result word */
    if (wskipl(da.LO))		/* If low sign set, carry from high sum! */
	op10m_dec(sum);
    op10m_sub(sum, db.HI);	/* Subtract high words */

    SUBFLAGS(sum, da.HI, db.HI)	/* Check for and set PC flags */

    da.HI = sum;		/* OK, put sum back in place */
    if (wskipl(da.HI))		/* Copy sign bit into low word of result */
	op10m_signset(da.LO);
    else op10m_signclr(da.LO);
    return da;

#else
    op10m_movn(db.LO);		/* Negate low word */
    if (wmagskipn(db.LO))	/* See if stuff, ignore sign bit */
	op10m_setcm(db.HI);	/* Low word has stuff, no carry into high */
    else op10m_movn(db.HI);	/* Low word clear, so carry into high word */
    return op10dadd(da, db);	/* Now add da + (-db) */
#endif
}

static dw10_t
x_dsub(register dw10_t da, register dw10_t db)
{
#endif /* IFFLAGS */
    op10m_dmovn(db);		/* Negate 2nd operand in place */
    return x_dadd(da, db);	/* and add in */
}


qw10_t op10dmul(register dw10_t da, register dw10_t db)
{
    register int sign;

    if (sign = wskipl(da.HI))
	op10m_dneg(da);
    if (wskipl(db.HI)) {
	sign = !sign;
	op10m_dneg(db);
#if IFFLAGS
	if (wskipl(db.HI) && wskipl(da.HI)) {
	    OP10_PCFSET(PCF_ARO+PCF_TR1);
	    return qw10maxneg;		/* Return max neg quadword */
	}
#endif
    }
    return sign ? x_qneg(x_dmul(da, db)) : x_dmul(da, db);
}

/* Unsigned double-int multiply */

static qw10_t x_dmul(register dw10_t da, register dw10_t db)
{
    register qw10_t q;		/* Quad-length result */
    register int ai, bi;
    uint32 av[5], bv[5], pv[10];

    /* First get 71-bit args into vectors of 16-bit integers.
    ** The main screw is ignoring the sign bit of low-order word.
    */
#define dsetvec(dw,v) \
	v[0] =  LHGET(dw.HI) >> 11;		/* High 7 bits of LH */\
	v[1] = (LHGET(dw.HI) & ((1<<11)-1))<<5 | (RHGET(dw.HI) >> 13); \
	v[2] = (RHGET(dw.HI) & ((1<<13)-1))<<3 | ((LHGET(dw.LO) >> 14) & 07); \
	v[3] = (LHGET(dw.LO) & ((1<<14)-1))<<2 | (RHGET(dw.LO) >> 16); \
	v[4] =  RHGET(dw.LO) & (((h10_t)1<<16)-1); \
	DEBUGPRF(("Array: %lo %lo %lo %lo %lo\n", \
		(long)v[0],(long)v[1],(long)v[2],(long)v[3],(long)v[4]));

    dsetvec(da, av);
    dsetvec(db, bv);
    for (ai = 10; --ai >= 0; )		/* Clear product vector */
	pv[ai] = 0;

    /* Now multiply the vectors together.  Because we have 32-bit unsigned
    ** integers, we can hold the result of multiplying 2 16-bit integers.
    */
    for (ai = 5; --ai >= 0;) {
	register uint32 hicarry = 0;
	for (bi = 5; --bi >= 0;) {
#if 1
	    register uint32 prod, high;
	    prod = av[ai] * bv[bi];
	    high = prod >> DIGBTS;
	    prod &= DIGMSK;
	    DEBUGPRF(("%d*b%d: %lo,%lo + %lo + p%d => ",
		ai,bi, (long)high, (long)prod, (long)hicarry, ai+bi+1));
	    hicarry += prod + pv[ai+bi+1];
	    pv[ai+bi+1] = hicarry & DIGMSK;
	    hicarry = (hicarry>>DIGBTS) + high;
	    DEBUGPRF(("p%d: %lo, carry %lo\n",
			ai+bi+1, (long)pv[ai+bi+1], (long)hicarry));
#else
	    hicarry += av[ai] * bv[bi] + pv[ai+bi+1];
	    pv[ai+bi+1] = hicarry & DIGMSK;
	    hicarry >>= DIGBTS;
#endif
	}
	pv[ai] = hicarry;
    }

    /* Now put together quadword from low 140+1 bits of product vector.
    ** Screw, again, is skipping over sign bits of low-order words.
    */
    DEBUGPRF(("ArOut: %lo %lo %lo %lo %lo %lo %lo %lo %lo %lo\n",
		(long)pv[0], (long)pv[1], (long)pv[2], (long)pv[3],
		(long)pv[4], (long)pv[5], (long)pv[6], (long)pv[7],
		(long)pv[8], (long)pv[9]));

    LHSET(q.D0.HI, (pv[1] & ((1<<13)-1))<< (18-13) | (pv[2] >> (16-(18-13))) );
    RHSET(q.D0.HI, (pv[2] & ((1<<11)-1))<< (18-11) | (pv[3] >> (16-(18-11))) );
    LHSET(q.D0.LO, (pv[3] & ((1<< 9)-1))<< (18-10) | (pv[4] >> (16-(18-10))) );
    RHSET(q.D0.LO, (pv[4] & ((1<< 8)-1))<< (18- 8) | (pv[5] >> (16-(18- 8))) );
    LHSET(q.D1.HI, (pv[5] & ((1<< 6)-1))<< (18- 7) | (pv[6] >> (16-(18- 7))) );
    RHSET(q.D1.HI, (pv[6] & ((1<< 5)-1))<< (18- 5) | (pv[7] >> (16-(18- 5))) );
    LHSET(q.D1.LO, (pv[7] & ((1<< 3)-1))<< (18- 4) | (pv[8] >> (16-(18- 4))) );
    RHSET(q.D1.LO, (pv[8] & ((1<< 2)-1))<< (18- 2) |  pv[9] );

    /* Must duplicate sign bit in all low-order words. */
    if (wskipl(q.D0.HI)) {
	op10m_signset(q.D0.LO);
	op10m_signset(q.D1.HI);
	op10m_signset(q.D1.LO);
    }
    return q;
}

/* Note returns quadword! High double is quotient */

qw10_t op10ddiv(
		qw10_t qw,	/* Dividend */
	register dw10_t d)	/* Divisor */
{
    qw10_t origq;
    int numsign, densign;

    origq = qw;			/* Save original arg in case of error */

    /* Make args positive.  Note that if either arg is the max neg number
    ** the result will be that same number.  Don't bother checking here
    ** as the initial dividend/divisor comparison test should catch them.
    */
    if (numsign = wskipl(qw.D0.HI))
	qw = x_qneg(qw);
    if (densign = wskipl(d.HI))
	op10m_dmovn(d);

    /* Check for zero divisor */
    if (!dlomagskipn(d)) {
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_DIV);
	return origq;		/* Overflow, no divide */
    }

    /* See if can do simpler division */
    if (!dlomagskipn(qw.D0) && !wmagskipn(qw.D1.HI) && !wskipn(d.HI)) {
	op10m_signclr(qw.D1.HI);	/* Ensure sign is off in low words */
	op10m_signclr(d.LO);
	DEBUGPRF(("ddivquick: %o/%o N: %lo,%lo,%lo,%lo D: %lo,%lo\n",
		numsign, densign, DBGV_D(qw.D1), DBGV_W(d.LO)));

	if (!x_div(&qw.D1, d.LO)) {	/* Do single-prec division */
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_DIV);
	    return origq;		/* Overflow, no divide */
	}
	qw.D0.LO = qw.D1.HI;		/* Quotient here, qw.D0.HI already 0 */
	op10m_setz(qw.D1.HI);		/* Remainder already in qw.D1.LO */
    } else {
	/* Ugh, must hack big numbers */
	if (!qdivstep(&qw, d, 70)) {		/* Divide, get 70 bits */
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_DIV);
	    return origq;		/* Overflow, no divide */
	}
	/* Fix up remainder.  Always has same sign as dividend. */
	if (wskipl(qw.D1.HI))		/* If final result was negative */
	    qw.D1 = x_dadd(qw.D1, d);	/* add divisor back in */
    }
    DEBUGPRF(("ddiv: %o/%o Q: %lo,%lo,%lo,%lo R: %lo,%lo,%lo,%lo\n",
	numsign, densign, DBGV_D(qw.D0), DBGV_D(qw.D1)));

    /* Unsigned division done, now fix up results */
    if (numsign)		/* Remainder has same sign as dividend */
	op10m_dneg(qw.D1);	/* Negate in place, fixed-point style */

    /* Now fix up quotient if dividend & divisor had different sign */
    if (numsign != densign)
	op10m_dneg(qw.D0);	/* Negate in place, fixed-point style */

    return qw;
}

/* General Floating-Point arithmetic defs
*/

/* PDP-10 Single-precision hardware floating point has the format:
**	1 sign bit
**	8 exponent bits, excess-0200, ones-complement if negative
**	27 fraction bits, twos-complement if negative
**
** PDP-10 Double-precision hardware floating point has the format:
**	1 sign bit
**	8 exponent bits, excess-0200, ones-complement if negative
**	27 fraction bits, twos-complement if negative
** plus	1 ignored low-word sign bit (set 0 by all operations)
**	35 additional fraction bits (62 total)
*/


/* Single-prec Format Exponent macros (also applicable to D format)
**	These treat the exponent as a quantity from the LH of high word.
*/
#define SFEBITS 9
#define SFEMASK (((h10_t)1<<SFEBITS)-1)	/* Right-justified mask */
#define SFESIGN (H10SIGN >> (H10BITS-SFEBITS)) /* Right-justified sign */
#define SFEMAGS  (SFESIGN-1)			/* Right-justified mag bits */
#define SFELHF   (SFEMASK << (H10BITS-SFEBITS))	/* Exponent field in LH */
#define SFEGET(w) ((LHGET(w) >> (H10BITS-SFEBITS)) & SFEMASK)
#define SFEEXCESS (SFESIGN>>1)			/* The N in Excess-N */

/* Single-prec Format Fraction macros (also applicable to D format)
**	These treat the fraction as a quantity from the LH of high word.
*/
#define SFFBITS (H10BITS-SFEBITS)		/* LH fraction bits */
#define SFFMASK (((h10_t)1<<SFFBITS)-1)
#define SFFHIBIT (H10SIGN>>SFEBITS)		/* Top bit (not a sign) */
#define SFFMAGS (SFFMASK>>1)			/* Remaining mag bits */

/* Macro to canonicalize exponent and fraction.
**	exp - must be an "int" variable
**	w - first word of float/double/gdouble
*/
#define SFSETUP(exp, w) \
    exp = SFEGET(w);		/* Get sign and exponent */\
    if (exp & SFESIGN) {	/* Negative? */\
	exp = exp ^ SFEMASK;	/* If so, get ones-complement, and */\
	op10m_tlo(w, SFELHF);	/* propagate sign thru exp */\
    } else			/* Else mask sign+exp out of fract */\
	op10m_tlz(w, SFELHF)

#define SF_POSSETUP(sign, exp, w) \
    exp = SFEGET(w);		/* Get sign and exponent */\
    if (sign = (exp & SFESIGN)) {	/* Negative? */\
	exp = exp ^ SFEMASK;	/* If so, get ones-complement, and */\
	op10m_tlo(w, SFELHF);	/* propagate sign thru exp */\
	op10m_movn(w);		/* and make positive */\
    } else			/* Else mask sign+exp out of fract */\
	op10m_tlz(w, SFELHF)

/* After SF_POSSETUP, can test for normalization with:
**	op10m_tlnn(w, SFFHIBIT|(SFFHIBIT<<1))
** If SFFHIBIT<<1 is set, fraction was negative 400,,0.
** Else if SFFHIBIT is set, fraction was either positive
**	or a negative less than 400,,0.
** If neither bit is set, fraction is unnormalized.
**	Just remember zero is valid if exponent also zero.
*/

#define SF_DPOSSETUP(sign, exp, d) \
    exp = SFEGET((d).HI);		/* Get sign and exponent */\
    if (sign = (exp & SFESIGN)) {	/* Negative? */\
	exp = exp ^ SFEMASK;		/* If so, get ones-complement, and */\
	op10m_tlo((d).HI, SFELHF);	/* propagate sign thru exp */\
	op10m_dmovn(d);			/* and make positive */\
	if (!op10m_tlnn((d).HI, SFFMASK))	/* If carried out of fract */\
	    ++exp, LHSET((d).HI, SFFHIBIT);	/* adjust it */\
    } else				/* Else mask sign+exp out of fract */\
	op10m_tlz((d).HI, SFELHF)


#if OP10_GFMT

/* PDP-10 G Format double-precision floating point has the format:
**	1 sign bit
**	11 exponent bits, excess-02000, ones-complement if negative
**	24 fraction bits, twos-complement if negative
** plus	1 ignored low-word sign bit (set 0 by all operations)
**	35 additional fraction bits (59 total)
**
** Basically identical to Double (D-format) but with 3 bits taken from
** the precision and given to the exponent.
*/

/* G Format Exponent macros
**	These treat the exponent as a quantity from the LH of high word.
*/
#define GFEBITS 12
#define GFEMASK  (((h10_t)1<<GFEBITS)-1)	/* Right-justified mask */
#define GFESIGN  (H10SIGN >> (H10BITS-GFEBITS))	/* Right-justified sign */
#define GFEMAGS  (GFESIGN-1)			/* Right-justified mag bits */
#define GFELHF   (GFEMASK << (H10BITS-GFEBITS))	/* Exponent field in LH */
#define GFEGET(w) ((LHGET(w) >> (H10BITS-GFEBITS)) & GFEMASK)
#define GFEEXCESS (GFESIGN>>1)			/* The N in Excess-N */

/* G Format Fraction macros.
**	These treat the fraction as a quantity from the LH of high word.
*/
#define GFFBITS (H10BITS-GFEBITS)		/* LH fraction bits */
#define GFFMASK (((h10_t)1<<GFFBITS)-1)
#define GFFHIBIT (H10SIGN>>GFEBITS)		/* Top bit (not a sign) */
#define GFFMAGS (GFFMASK>>1)			/* Remaining mag bits */

/* Macro to canonicalize exponent and fraction.
**	exp - must be an "int" variable
**	w - first word of float/double/gdouble
*/
#define GFSETUP(exp, w) \
    exp = GFEGET(w);		/* Get sign and exponent */\
    if (exp & GFESIGN) {	/* Negative? */\
	exp = exp ^ GFEMASK;	/* If so, get ones-complement, and */\
	op10m_tlo(w, GFELHF);	/* propagate sign thru exp */\
    } else			/* Else mask sign+exp out of fract */\
	op10m_tlz(w, GFELHF)


/* Macro to insert exponent into high word of G-format double.
*/
#define GFEXPSET(w, exp) \
    if (wskipl(w))	/* If sign set, use ones-complement of exp */\
	LHSET(w, (LHGET(w) ^ ((h10_t)(exp & GFEMAGS) << (H10BITS-GFEBITS))));\
    else								\
	LHSET(w, (LHGET(w) | ((h10_t)(exp & GFEMAGS) << (H10BITS-GFEBITS))))

#define GF_DPOSSETUP(sign, exp, d) \
    exp = GFEGET((d).HI);		/* Get sign and exponent */\
    if (sign = (exp & GFESIGN)) {	/* Negative? */\
	exp = exp ^ GFEMASK;		/* If so, get ones-complement, and */\
	op10m_tlo((d).HI, GFELHF);	/* propagate sign thru exp */\
	op10m_dmovn(d);			/* and make positive */\
	if (!op10m_tlnn((d).HI, GFFMASK))	/* If carried out of fract */\
	    ++exp, LHSET((d).HI, GFFHIBIT);	/* adjust it */\
    } else				/* Else mask sign+exp out of fract */\
	op10m_tlz((d).HI, GFELHF)


#endif /* OP10_GFMT */

/* Floating-point Single-precision arithmetic */

/* Internal flags for ffadr(), selected so most common case (fadr)
** has flags of 0.
*/
#define FADF_NONORM	04	/* Set to skip normalization */
#define FADF_NEGB	02	/* Set to negate B before adding */
#define FADF_NORND	01	/* Set to skip rounding */

w10_t op10fadr(register w10_t a, register w10_t b)
{	return ffadr(a, b, 0);		/* No negate, Round */
}

w10_t op10fad(register w10_t a, register w10_t b)
{	return ffadr(a, b, FADF_NORND);	/* No negate, no round */
}

w10_t op10fsbr(register w10_t a, register w10_t b)
{	return ffadr(a, b, FADF_NEGB);	/* Negate, Round */
}

w10_t op10fsb(register w10_t a, register w10_t b)
{	return ffadr(a, b, FADF_NEGB|FADF_NORND);	/* Negate, no round */
}

/* Common floating add subroutine shared by above external routines. */

static w10_t ffadr(register w10_t a, register w10_t b, int flg)
{
    register dw10_t d;
    register int i, expa, expb;
    register uint18 rbit;		/* May be 1<<17 */

    SFSETUP(expa, a);		/* Set up exponent and fraction for A */
    SFSETUP(expb, b);		/* Ditto for B */

    if (flg & FADF_NEGB)	/* If actually subtracting B, */
	op10m_movn(b);		/* OK to negate fraction now */

    /* Now see which exponent is smaller; get larger in A, smaller in B.
    ** This needs to be done even if one of the args is completely zero, so
    ** that the other is set up in A for the normalization code.
    */
    if ((i = expa - expb) < 0) {
	w10_t tmp;			/* Swap args */
	int etmp;
	etmp = expb, tmp = b;
	expb = expa, b = a;
	expa = etmp, a = tmp;
	i = -i;
    }

    /* Now right-shift B by the number of bits in i, then add to A. */
    d.HI = x_ash(b, -i);		/* Shift B and put in hi word */
    op10m_add(d.HI, a);			/* then add A to it */
    d.LO = x_ash(b, 35-i);		/* Set 2nd word (low of B) */

    DEBUGPRF(("FADR preN: (%o) %lo,%lo,%lo,%lo\n", expa, DBGV_D(d)));

    /* D now has double-length sum.
    ** The high word sign of A is correct; that of low word may not be.
    ** Check high 10 bits to determine normalization step.
    */
    rbit = 0;
    switch (LHGET(d.HI) >> 8) {
	case 00004>>2:		/* Normalized positive fraction */
	    DEBUGPRF(("WINNING POSITIVE\n"));
	    rbit = (LHGET(d.LO) & (HSIGN>>1));
	    break;

	case 00010>>2:		/* Overflow of pos fraction */
	case 00014>>2:
	    DEBUGPRF(("OVERFLOW POSITIVE\n"));
	    rbit = RHGET(d.HI) & 01;
	    d.HI = x_ash(d.HI, -1);	/* Normalize by shifting right */
	    ++expa;			/* And adding one to exponent */
	    break;

	case 07760>>2:		/* Overflow of neg fraction, maybe zero mag */
	case 07764>>2:		/* Overflow of neg fraction, nonzero mag */
	    DEBUGPRF(("OVERFLOW NEGATIVE\n"));
	    d = x_ashc(d, -1);	/* First normalize into 0777 case, */
	    ++expa;		/* then fall thru below to check! */

	case 07770>>2:		/* Normalized negative fraction (probably) */
	    if (!(flg & FADF_NORND)) {
		rbit = LHGET(d.LO) & (HSIGN>>1);
		if (rbit && !op10m_txnn(d.LO, (HMASK>>2), HMASK))
		    rbit = 0;		/* No round if no other bits */
	    }
	    if (op10m_txnn(d.HI, 0777, HMASK) || rbit) {
		DEBUGPRF(("WINNING NEGATIVE\n"));
		break;
	    }
	    DEBUGPRF(("OVERFLOW NEG, ZEROMAG -2^28 case\n"));
	    ++expa;
	    op10m_tlo(d.HI, 0400);	/* Effect ASHC -1 */
	    break;

	/* Anything else is an unnormalized result and must be handled
	** using slow, fully general case.
	*/
	case 0:
	    if (!wskipn(d.HI) && !wmagskipn(d.LO)) {
		DEBUGPRF(("FADR: S: ZERO RESULT!\n"));
		return w10zero;
	    }
	default:
	    if (flg & FADF_NONORM)	/* May not want normalization */
		break;

	    /* Must mess with low word, sigh.  Ensure it has correct sign bit,
	    ** then find the first fraction bit in the 2-word sum.
	    ** Normalization shifts until B0 != B9 OR (B9==1 and B10-35 == 0).
	    ** The latter case must be checked for if negative.
	    */
	    if (wskipl(d.HI)) {
		op10m_signset(d.LO);
		if ((i = op10ffo(op10setcm(d.HI))) >= 36)
		    i = op10ffo(op10setcm(d.LO)) + 35;
	    } else {
		op10m_signclr(d.LO);
		if ((i = op10ffo(d.HI)) >= 36)
		    i = op10ffo(d.LO) + 35;
	    }
	    DEBUGPRF(("RENORMALIZATION: %d", i-9));
	    /* i now has position of first fraction bit; move it into place. */
	    d = x_ashc(d, (i -= 9));		/* Derive shift and do it */
	    expa -= i;				/* Adjust exponent */
	    rbit = LHGET(d.LO) & (HSIGN>>1);	/* Find roundup bit */
	    if (wskipl(d.HI)) {			/* Special checks for neg */
		DEBUGPRF((" checking neg"));
		if (!op10m_txnn(d.HI, 0777, HMASK)) {
		    /* If fraction is all zero, must be the ugly
		    ** special case of negative 2^28 (777000,,0)
		    */
		    DEBUGPRF((" - special -2^28 case"));
		    ++expa;
		    op10m_tlo(d.HI, 0400);	/* Effect ASHC -1 */
		    rbit = 0;
		} else if (rbit && !op10m_txnn(d.LO, (HMASK>>2), HMASK)) {
		    DEBUGPRF((" - no lower-round bits, rbit=0"));
		    rbit = 0;			/* No round if no other bits */
		}
	    }
	    DEBUGPRF(("\n"));
	    break;
    }

    /* At this point, fraction in d.HI is normalized, and rbit is set if a
    ** roundup increment must be done (which in rare cases can trigger
    ** a single-shift renormalization).
    */
    if (!(flg & FADF_NORND) && rbit) {
	DEBUGPRF(("ROUNDUP DONE"));
	op10m_inc(d.HI);		/* Add 1, then check for overflow */
	if (wskipl(d.HI)) {		/* Negatives are messy */
	    if (LHGET(d.HI) & 0400) {	/* If first fract bit is "wrong" */
		if ((LHGET(d.HI) & 0377) || RHGET(d.HI)) {	/* See if really OK */
		    --expa, d.HI = x_ash(d.HI, 1);	/* Nope, shift up */
		    DEBUGPRF((" - NEG RENORMALIZED!"));
		}
	    }
	} else if (!(LHGET(d.HI) & 0400)) {
	    ++expa, LHSET(d.HI, 0400);	/* Positive overflow, fix up */
	    DEBUGPRF((" - POS RENORMALIZED!"));
	}
	DEBUGPRF(("\n"));
    }
    DEBUGPRF(("FADR: S: (%o) %lo,%lo\n", expa, DBGV_W(d.HI)));

    /* Ta da!  Put exponent back in.  We can check now for overflow or
    ** underflow.
    */
    if (wskipl(d.HI))	/* If sign set, put ones-complement of exp */
	LHSET(d.HI, (LHGET(d.HI) ^ ((h10_t)(expa & 0377) << 9)));
    else
	LHSET(d.HI, (LHGET(d.HI) | ((h10_t)(expa & 0377) << 9)));

#if IFFLAGS
    if (expa < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (expa > SFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif
    return d.HI;
}



w10_t op10fmp(register w10_t a, register w10_t b)
{	return ffmpr(a, b, 0);		/* Do without rounding */
}

w10_t op10fmpr(register w10_t a, register w10_t b)
{	return ffmpr(a, b, 1);		/* Do with rounding */
}

static w10_t ffmpr(register w10_t a, register w10_t b, int dornd)
{
    register dw10_t d;
    register int sign, exp, i, expb;

    /* Make operands positive and remember sign of result */
    SF_POSSETUP(sign, exp, a);		/* Get positive exp and fract */

    /* Repeat same steps for other arg */
    expb = SFEGET(b);			/* Get sign and exponent */
    if (expb & SFESIGN) {		/* Negative? */
	sign = !sign;			/* Set 0 if both args negative */
	expb = expb ^ SFEMASK;		/* If so, get ones-complement */
	op10m_tlo(b, SFELHF);		/* propagate sign thru exp */
	op10m_movn(b);			/* Then make positive */
    } else
	op10m_tlz(b, SFELHF);		/* Else mask sign+exp out of fract */

    /* Add exponents together */
    exp += expb - 0200;

    /* Now multiply the 27-bit magnitude numbers (maybe 28-bit if neg) */
    d = op10xmul(a, b);		/* Multiply the numbers, get double result */

    /* The high 17 bits of product should be empty (1 sign, 16 magnitude bits)
    ** If everything was normalized, the next bit (start of fraction) should
    ** be 1 and we just shift everything left 8 bits to get in place.
    ** If not, must look for first bit of fraction.
    */
    DEBUGPRF(("FMP preN: (%d.) %lo,%lo,%lo,%lo\n", exp, DBGV_D(d)));
    if (LHGET(d.HI) == 01) i = 8;
    else {
	/* High fraction bit not set... see if result was 0 */
	if (!wskipn(d.HI) && !wskipn(d.LO))
	    return d.HI;		/* Zero product, return 0.0 */
	i = adffo(d);			/* Find first fract bit in double */
	i -= SFEBITS;			/* Adjust for place shifting to */
    }

    /* OK, now normalize. */
    d = x_ashc(d, i);			/* Do ASHC bringing in zeros */
    exp -= i - 8;
    DEBUGPRF(("FMP pstN: (%d.) %lo,%lo,%lo,%lo\n", exp, DBGV_D(d)));

    /* Now do rounding.  See if bit lower than LSB is set; if so, add one
    ** to LSB.  This might possibly carry all the way to MSB, so be ready
    ** to re-normalize.
    */
    if (dornd) {
	if (op10m_tlnn(d.LO, (HSIGN>>1))) {	/* Roundup bit there? */
	    op10m_inc(d.HI);
	    if (op10m_tlnn(d.HI, SFELHF)) {	/* If carried out of fract */
		LHSET(d.HI, SFFHIBIT);	/* must be 01000, shift it back */
		++exp;			/* and adjust exponent */
	    }
	}
    }

    /* Fraction all done, just add exponent back now */
    op10m_tlo(d.HI, ((h10_t)(exp & SFEMAGS) << SFFBITS));
#if IFFLAGS
    if (exp < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (exp > SFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif

    /* Convert back to negative result if necessary.  Note special hackery
    ** if not rounding, which merely gives ones-complement if any bits were
    ** set in low-order part.  This is needed to reproduce hardware, even
    ** though it results in a less precise answer!
    */
    if (sign) {
	if (dornd || !wskipn(d.LO))
	    op10m_movn(d.HI);
	else
	    op10m_setcm(d.HI);
    }
    return d.HI;
}


/* FDV - slightly different calling sequence needed because like
**	other divides, operation can be aborted without changing
**	any operands.
**
** OP10FDV{R} are provided for backward compatibility with optest prog.
*/

w10_t op10fdv(w10_t a, w10_t b)
{
    (void) op10xfdv(&a, b, 0);	/* Do without rounding */
    return a;			/* Return whatever result was */
}				/* (no change if no divide) */

w10_t op10fdvr(w10_t a, w10_t b)
{
    (void) op10xfdv(&a, b, 1);	/* Do with rounding */
    return a;			/* Return whatever result was */
}				/* (no change if no divide) */


int op10xfdv(w10_t *aw,
	     register w10_t b,
	     int dornd)
{
    register w10_t a = *aw;
    dw10_t d;
    register int sign, exp, i, expb;

    /* Make operands positive and remember sign of result */

    if (sign = wskipl(a)) {	/* Negative? */
	op10m_movn(a);		/* Then make positive */
	if (wskipl(a)) {		/* If numerator is max neg value, */
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	    return 0;			/* PDP-10 always chokes, sigh */
	}
    }
    exp = SFEGET(a);			/* Get positive exponent */
    op10m_tlz(a, SFELHF);		/* Mask sign+exp out of fract */

    /* Repeat for other arg, but allow max neg number */
    expb = SFEGET(b);			/* Get sign and exponent */
    if (expb & SFESIGN) {		/* Negative? */
	sign = !sign;			/* Set 0 if both args negative */
	expb = expb ^ SFEMASK;		/* If so, get ones-complement */
	op10m_tlo(b, SFELHF);		/* propagate sign thru exp */
	op10m_movn(b);			/* Then make positive */
    } else
	op10m_tlz(b, SFELHF);		/* Else mask sign+exp out of fract */


    /* Subtract exponents.  The +1 compensates for pre-multiply of divisor.
    */
    exp += 0200 + 1 - expb;
    DEBUGPRF(("FDVR (%o) N: %lo,%lo D: %lo,%lo\n",
		exp, DBGV_W(a), DBGV_W(b)));

    /* Get two 27-bit (or 28 if SETZ) magnitude numbers and divide them.
    ** In order to ensure that the division can succeed, we pre-multiply the
    ** divisor by 2.  This is compensated for by adjusting the final exponent.
    ** This can fail and cause No Divide (etc) if operands aren't normalized.
    */
    b = x_ash(b, 1);		/* Make divisor be bigger */
    if (op10m_camge(a, b)) {
	DEBUGPRF(("FDVR: Fail, N: %lo,%lo >= D: %lo,%lo\n",
			DBGV_W(a), DBGV_W(b)));

	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	return 0;		/* Fail, return original dividend */
    }

    d.HI = a;			/* Set up dividend */
    op10m_setz(d.LO);
    DEBUGPRF(("FDVR: N: %lo,%lo|%lo,%lo  D: %lo,%lo\n", DBGV_D(d), DBGV_W(b)));

    /* Used to be always 30; maybe revert if needed to duplicate
    ** unnormalized behavior??
    */
    i = (dornd ? 29 : 28);
    if (!ddivstep(&d, b, i, 0)) {	/* Divide them, get double result */
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	return 0;			/* Fail */
    }
    exp += 35-i;			/* Adjust for partial division */

    DEBUGPRF(("FDV preN: Q: %lo,%lo R: %lo,%lo  E: %o\n", DBGV_D(d), exp));

    if (!wskipn(d.HI)) {		/* Check for zero result */
	*aw = d.HI;
	return TRUE;
    }

    /* To compensate for possible unnormalized operands and roundup carry,
    ** we do a completely general normalization and round.  Note we've already
    ** tested for 0 so bits definitely exist.
    */
    i = op10ffo(d.HI) - 9;		/* Find first one-bit, derive shift */
    if (dornd && i < 0) {		/* If right-shifting, */
	/* Note max bits needed are 9, so "int" is big enough */
	register int bit = 1 << (-i-1);	/* Find last bit that will vanish */
	if (RHGET(d.HI) & bit) { 		/* If need to round up, */
	    static w10_t wrnd;
	    LRHSET(wrnd, 0, bit);
	    op10m_add(d.HI, wrnd);		/* do it */
	    DEBUGPRF(("ROUNDED UP!\n"));
	    if (i != (op10ffo(d.HI)-9)) {	/* See if carried */
		DEBUGPRF(("CARRIED!\n"));
		if (wskipl(d.HI)) {		/* Yes, into sign bit? */
		    DEBUGPRF(("CARRIED INTO SIGN!\n"));
		    d.HI = x_ash(d.HI, -1);	/* Yes, yecch!  Special */
		    op10m_signclr(d.HI);	/* handling to undo damage */
		    --exp;
		} else --i;		/* Carried but compensation trivial */
	    }
	}
    }
    if (i != 8)
	DEBUGPRF(("ABNORMALIZATION: %d\n", i));

    d.HI = x_ash(d.HI, i);		/* Shift the fraction into place */
    exp -= i + 8;			/* Adjust exponent for shift */

    DEBUGPRF(("FDV pstN: Q: %lo,%lo R: %lo,%lo  E: %o\n", DBGV_D(d), exp));

    /* Fraction all done, just add exponent back now */
    LHSET(d.HI, (LHGET(d.HI) | ((h10_t)(exp & SFEMAGS) << SFFBITS)));
#if IFFLAGS
    if (exp < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (exp > SFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif

    /* Ugly hackery - see PRM 2-23 "Note" for FDV.
    ** For FDV on KL+KS a negative Q is represented by twos-complement only if
    ** remainder is 0, else ones-complement (similar to the way FMP works!).
    ** On KA+KI it is always twos-complement...
    */
    if (sign) {
#if OP10_FDV_KL		/* KL or KS */
	if (!dornd) {
	    /* Fix up remainder for testing */
	    if (wskipl(d.LO))		/* If negative, */
		op10m_add(d.LO, b);	/* add divisor back in */
	    if (wskipn(d.LO))
		op10m_setcm(d.HI);
	    else
		op10m_movn(d.HI);
	} else
#endif
	    op10m_movn(d.HI);
    }
    *aw = d.HI;
    return TRUE;
}

#if 0	/* Old version */

static w10_t ffdvr(a, b, dornd)
register w10_t a, b;
{
    dw10_t d;
    register int sign, exp, i, expb;

    /* Make operands positive and remember sign of result */
    d.HI = a;				/* Preserve A in case fail */
    if (sign = wskipl(d.HI)) {	/* Negative? */
	op10m_movn(d.HI);		/* Then make positive */
	if (wskipl(d.HI)) {		/* If numerator is max neg value, */
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	    return a;			/* PDP-10 always chokes, sigh */
	}
    }
    exp = SFEGET(d.HI);			/* Get positive exponent */
    op10m_tlz(d.HI, SFELHF);		/* Mask sign+exp out of fract */

    /* Repeat for other arg, but allow max neg number */
    expb = SFEGET(b);			/* Get sign and exponent */
    if (expb & SFESIGN) {		/* Negative? */
	sign = !sign;			/* Set 0 if both args negative */
	expb = expb ^ SFEMASK;		/* If so, get ones-complement */
	op10m_tlo(b, SFELHF);		/* propagate sign thru exp */
	op10m_movn(b);			/* Then make positive */
    } else
	op10m_tlz(b, SFELHF);		/* Else mask sign+exp out of fract */


    /* Subtract exponents.  The +1 compensates for pre-multiply of divisor.
    */
    exp += 0200 + 1 - expb;
    DEBUGPRF(("FDVR (%o) N: %lo,%lo D: %lo,%lo\n",
		exp, DBGV_W(d.HI), DBGV_W(b)));

    /* Get two 27-bit (or 28 if SETZ) magnitude numbers and divide them.
    ** In order to ensure that the division can succeed, we pre-multiply the
    ** divisor by 2.  This is compensated for by adjusting the final exponent.
    ** This can fail and cause No Divide (etc) if operands aren't normalized.
    */
    b = x_ash(b, 1);		/* Make divisor be bigger */
    if (op10m_camge(d.HI, b)) {
	DEBUGPRF(("FDVR: Fail, N: %lo,%lo >= D: %lo,%lo\n",
			DBGV_W(d.HI), DBGV_W(b)));

	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	return a;		/* Fail, return original dividend */
    }

    op10m_setz(d.LO);		/* Set up dividend */
    DEBUGPRF(("FDVR: N: %lo,%lo|%lo,%lo  D: %lo,%lo\n", DBGV_D(d), DBGV_W(b)));

    if (!ddivstep(&d, b, 30, 0)) {	/* Divide them, get double result */
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	return a;
    }
    exp += 35-30;			/* Adjust for partial division */

    DEBUGPRF(("FDVR: Q: %lo,%lo R: %lo,%lo  E: %o\n", DBGV_D(d), exp));

    if (!wskipn(d.HI))		/* Check for zero result */
	return d.HI;

    /* To compensate for possible unnormalized operands and roundup carry,
    ** we do a completely general normalization and round.  Note we've already
    ** tested for 0 so bits definitely exist.
    */
    i = op10ffo(d.HI) - 9;		/* Find first one-bit, derive shift */
    if (dornd && i < 0) {		/* If right-shifting, */
	/* Note max bits needed are 9, so "int" is big enough */
	register int bit = 1 << (-i-1);	/* Find last bit that will vanish */
	if (RHGET(d.HI) & bit) { 		/* If need to round up, */
	    static w10_t wrnd;
	    LRHSET(wrnd, 0, bit);
	    op10m_add(d.HI, wrnd);		/* do it */
	    DEBUGPRF(("ROUNDED UP!\n"));
	    if (i != (op10ffo(d.HI)-9)) {	/* See if carried */
		DEBUGPRF(("CARRIED!\n"));
		if (wskipl(d.HI)) {		/* Yes, into sign bit? */
		    DEBUGPRF(("CARRIED INTO SIGN!\n"));
		    d.HI = x_ash(d.HI, -1);	/* Yes, yecch!  Special */
		    op10m_signclr(d.HI);	/* handling to undo damage */
		    --exp;
		} else --i;		/* Carried but compensation trivial */
	    }
	}
    }
    if (i != 8)
	DEBUGPRF(("ABNORMALIZATION: %d\n", i));

    d.HI = x_ash(d.HI, i);		/* Shift the fraction into place */
    exp -= i + 8;			/* Adjust exponent for shift */

    /* Fraction all done, just add exponent back now */
    LHSET(d.HI, (LHGET(d.HI) | ((h10_t)(exp & SFEMAGS) << SFFBITS)));
#if IFFLAGS
    if (exp < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (exp > SFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif

    if (sign)
	op10m_movn(d.HI);
    return d.HI;
}

#endif /* 0 - old version */

/* Single-precision Floating-point Conversions */

/* SFNORM - Returns normalized single-precision floating point, given
**	full-word integer fraction and exponent to use if fraction were
**	already in proper position (first significant bit at B9).
**	Sets flags if exponent overflows.
*/
static w10_t sfnorm(register int exp,
		    register w10_t w,	/* Contains signed integer fraction */
		    int rnd)		/* TRUE to round result */
{
    register int i, rbit;	/* rbit only needs 9 bits, so "int" safe */

    /* Normalizing a negative number has a lot of gross hair associated
    ** with it.  By far the simplest thing to do is convert all negative
    ** numbers to positive values and only negate back at the end of
    ** normalization.
    ** Skeptical?  See the old code here, included for posterity.  Good luck.
    */
#if 0 /* Old code */
#  define bmask(n) (((h10_t)1<<(n))-1)	/* Mask of N right-justified bits */
	/* Alternative def:
		define bmask(n) (01777777>>(n))
	  Then refs change from
		bmask(18-i)	=>	bmask(i+1)
		bmask(36-i)	=>	bmask(i-17)
	  Only used in sfnorm()...
	*/
    if (wskipl(w)) {		/* Find 1st significant bit */
	i = op10ffo(op10setcm(w));
	/* Messy - must see if any lower bits are set, and if not,
	** adjust position back to compensate for 2s-complement form.
	*/
	if (i < 17) {			/* Check LH and RH both? */
	    if (!RHGET(w) && !(LHGET(w) & bmask(18-i))) --i;
	} else if (i < 36) {		/* Check only RH */
	    if (!(RHGET(w) & bmask(36-i))) --i;
	} else i = 35;			/* w is -1 */

    } else i = op10ffo(w);
    DEBUGPRF(("SFNORM: %lo,%lo i=%d\n", DBGV_W(w), i));
    if (i >= 36)
	return w10zero;

    /* i now has position of first fraction bit; move it into place.
    ** For a negative number, the roundup bit is 1 only if at least
    ** one lower bit is set.
    */
    i -= 9;			/* Make relative to proper position, B9 */
    if (rnd && i < 0) {		/* If rounding OK and losing bits */
	rbit = 1 << ((-i)-1);	/* Get bitmask for last bit lost */
	if (wskipl(w))		/* If negative, must have at least one 1-bit */
	    rbit = RHGET(w) & (rbit-1);	/* in bits lost off right  */
	else rbit = RHGET(w) & rbit;
    } else rbit = 0;
    if (i)
	w = x_ash(w, i);	/* Derive shift and do it */
    exp -= i;			/* Adjust exponent */

    /* At this point, number in is normalized, and rbit is set if a
    ** roundup increment must be done (which in rare cases can trigger
    ** a single-shift renormalization).
    */
    if (rbit) {
	DEBUGPRF(("ROUNDUP DONE"));
	op10m_inc(w);			/* Add then check for overflow */
	if (wskipl(w)) {		/* Negatives are messy */
	    if (LHGET(w) & 0400) {	/* If first fract bit is "wrong" */
		if ((LHGET(w) & 0377) || RHGET(w)) {	/* See if really OK */
		    --exp, w = x_ash(w, 1);	/* Nope, shift up */
		    DEBUGPRF((" - NEG RENORMALIZED!"));
		}
	    }
	} else if (!(LHGET(w) & 0400)) {
	    ++exp, LHSET(w, 0400);		/* Positive overflow, fix up */
	    DEBUGPRF((" - POS RENORMALIZED!"));
	}
	DEBUGPRF(("\n"));
    }

    /* Ta da!  Put exponent back in.  We can check now for overflow or
    ** underflow.
    */
    if (wskipl(w))	/* If sign set, put ones-complement of exp */
	LHSET(w, (LHGET(w) ^ ((h10_t)(exp & 0377) << 9)));
    else
	LHSET(w, (LHGET(w) | ((h10_t)(exp & 0377) << 9)));

#else 	/* End old code, start new code */

    /* New regime.  Check for negative fraction and convert to use
    ** positive value instead.  The fraction is always assumed to be signed.
    **
    ** Note that the maximum negative number cannot be made positive,
    ** but this is compensated for by remembering the sign and using
    ** unsigned operations.
    */
    register int sign;

    if (sign = wskipl(w))
	op10m_movn(w);		/* Negate the fraction, max neg is OK. */
    i = op10ffo(w);

    DEBUGPRF(("SFNORM: %lo,%lo i=%d exp=%o\n", DBGV_W(w), i, exp));

    if (i >= 36)
	return w10zero;		/* No magnitude bits, so just return 0 */

    /* i now has position of first fraction bit; move it into place.
    */
    i -= 9;			/* Make relative to proper position, B9 */
    exp -= i;			/* Adjust exponent */
    if (rnd && i < 0)		/* If rounding OK and losing bits */
	rbit = RHGET(w)		/* Find value of last bit that will be lost */
		& (1 << ((-i)-1));
    else rbit = 0;

    if (i < 0)			/* Now do logical (unsigned) shift! */
	RSHIFTM(w, -i);
    else if (i > 0)
	LSHIFTM(w, i);

    DEBUGPRF(("SFNORMED: %lo,%lo i=%d exp=%o rbit=%o\n",
				DBGV_W(w), i, exp, rbit));

    /* At this point, number in w is normalized, and rbit is set if a
    ** roundup increment must be done.  In the rare case that this
    ** increment overflows, its value will be 1000,,0 and will trigger
    ** a single-shift renormalization.
    */
    if (rbit) {
	op10m_inc(w);			/* Add, then check for overflow */
	if (LHGET(w) > SFFMASK) {
	    ++exp;			/* Positive overflow, fix up */
	    LHSET(w, 0400);
	    DEBUGPRF(("ROUNDUP - RENORMED!\n"));
	} else
	    DEBUGPRF(("ROUNDUP\n"));
    }

    /* Fraction is normalized with exponent bits clear, can just OR
    ** exponent back in, and negate if needed to finish off result.
    */
    op10m_tlo(w, ((h10_t)(exp & SFEMAGS) << SFFBITS));
    if (sign)
	op10m_movn(w);

    /* Ta da! Check before returning for overflow or underflow.
    ** Note that the negation, if one, can never cause either since a
    ** positive number can always be negated.
    */
#endif /* New code */

#if IFFLAGS
    if (exp < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (exp > SFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif

    return w;
}

w10_t
op10fltr(register w10_t w)
{
    return sfnorm(SFEEXCESS+27, w, 1);	/* Normalize and round */
}

/* FSC
**	Testing weirdness: when AC has a single-bit fraction (either
**	positive or complemented) and E is 400000, then *only* Floating
**	Underflow is set on a real DEC-2065.  The proc ref man claims
**	Overflow and Floating Overflow should always also be set, and
**	that's what this code does.
*/
w10_t
op10fsc(register w10_t w, register h10_t h)
{
    register int exp;

    SFSETUP(exp, w);		/* Set up exponent and fraction */
    exp += imm8op(h);		/* Scale per immediate arg */
    return sfnorm(exp, w, 0);	/* Put float together, don't round */
}

w10_t op10fix(w10_t w)
{
    op10xfix(&w, 0);		/* Fix, no round */
    return w;
}

w10_t op10fixr(w10_t w)
{
    op10xfix(&w, 1);		/* Fix, with round */
    return w;
}

/* Single Float->Fix conversion.
**	Note that rounding does NOT work the same way as for floats.
**	See PRM 2-28.
**	+(1.5) is rounded to +2, but -(1.5) is rounded to -1!
**
** If the first bit of the fraction is 1 for a positive number,
**	then 1 is added to the integral part (thus increasing its magnitude).
** If the first bit of the fraction is 1 for a negative number,
**	then 1 is added to the integral part (thus lowering its magnitude).
**
**	This produces asymmetric results due to twos-complement nature.
**	Apparently this was done for simplicity or to follow the "Algol
**	standard".
*/
int op10xfix(w10_t *aw, int dornd)
{
    register w10_t w = *aw;

#if 1
    register int i, sign;
    register uint18 rbit;

    SF_POSSETUP(sign, i, w);	/* Set up exponent, get positive fract */

    /* Apply exponent size test now */
    if (i > (SFEEXCESS+35)) {		/* If integer wd be too big, */
	OP10_PCFSET(PCF_ARO+PCF_TR1);	/* set flags */
	return 0;			/* and return arg unchanged */
    }
    
    /* OK to proceed!
    ** Fraction should now be a 27-bit positive integer.  It will be a 28-bit
    ** value of 1000,,0 only if given an unnormalized negative float with
    ** no one-bits in the fraction part.
    */
    i -= (SFEEXCESS+27);	/* Get # bits to shift fraction */
    if (i <= -28) {		/* If shifting it to oblivion */
	*aw = w10zero;		/* just return zero, skip rigmarole */
	return 1;
    }

    /* Determine whether to round.
    **	Note special test for negative values, needed in order to
    **	emulate peculiar way FIXR does rounding.
    */
    if (dornd && i < 0) {	/* If rounding OK and losing bits */
	if (i < -H10BITS) {	/* Find value of last bit that will be lost */
	    rbit = LHGET(w)	/* Bit comes from LH */
		    & (1 << ((-i)-(H10BITS+1)));
	} else {
	    rbit = RHGET(w)	/* Bit comes from RH */
		    & (((uint18)1) << ((-i)-1));
	}
	if (rbit && sign	/* Special test for negative .5 */
	  && !op10m_tdnn(w, w10loobits[(-i)-1]))
	    rbit = 0;		/* No round if so */
    } else rbit = 0;

    if (i < 0) {		/* Now do logical (unsigned) shift! */
	RSHIFTM(w, -i);
	if (rbit)		/* Need any rounding? */
	    op10m_inc(w);	/* Yep, trivial */
    } else if (i > 0)
	LSHIFTM(w, i);

    if (sign)			/* Now negate back if necessary */
	op10m_movn(w);

#else	/* New regime! */

    register int i;
    register uint18 rbit;

    SFSETUP(i, w);		/* Set up exponent, get signed fract */

    /* Apply exponent size test now */
    if (i > (SFEEXCESS+35)) {		/* If integer wd be too big, */
	OP10_PCFSET(PCF_ARO+PCF_TR1);	/* set flags */
	return 0;			/* and return arg unchanged */
    }
    
    /* OK to proceed!
    ** Fraction should now be a 27-bit positive or negative integer.
    ** It will be a 28-bit value of -1000,,0 only if given an unnormalized
    ** negative float with no one-bits in the fraction part.
    */
    i -= (SFEEXCESS+27);	/* Get # bits to shift fraction */
    if (i <= -28) {		/* If shifting it to oblivion */
	*aw = w10zero;		/* just return zero, skip rigmarole */
	return 1;
    }

    /* Determine whether to round.
    */
    if (dornd && i < 0) {	/* If rounding OK and losing bits */
	if (i < -H10BITS) {	/* Find value of last bit that will be lost */
	    rbit = LHGET(w)	/* Bit comes from LH */
		    & (1 << ((-i)-(H10BITS+1)));
	} else {
	    rbit = RHGET(w)	/* Bit comes from RH */
		    & (((uint18)1) << ((-i)-1));
	}
    } else rbit = 0;

    if (i < 0) {		/* Now shift! */
	if (wskipl(w))
	    NASH_RSHIFTM(w, -i); /* Negative, must do ASH */
	else
	    RSHIFTM(w, -i);	/* Positive, fast logical is OK */
	if (rbit)		/* Need any rounding? */
	    op10m_inc(w);	/* Yep, trivial */
    } else if (i > 0)
	LSHIFTM(w, i);		/* Left-shift can always be logical */
#endif

    *aw = w;			/* Return result */
    return 1;			/* Say succeeded */
}

/* Floating-point Double-precision arithmetic.  REAL moby hair! */


dw10_t op10dfad(dw10_t a, dw10_t b)
{
    return dfad(a, b, 0);
}

dw10_t op10dfsb(dw10_t a, dw10_t b)
{
    return dfad(a, b, 1);
}

/* Common DFAD/DFSB routine.

	The PRM says addition is done in a "triple-length" sum.  It's
unclear exactly how many magnitude bits this is.
	For KLH10_CPU_KLX at least, 3 words of 27+35+35 bits are not enough;
the following computation failed when using just one extra word for
the "lost" roundup bits, because the low bit of the 2nd operand was
being shifted out of the 3rd word.  Even though these operands are
not normalized, ideally the emulation should still produce the
same results.
	
2065: 	dfad   700000,0,0,0    776000,0,0,1  =>  677400,0,0,0
Emul:	dfad   700000,0,0,0    776000,0,0,1  =>  677377,777777,377777,777777

	So, for the time being this code uses 4 words, with the lowest
word used only to hold shifted-out bits for possible rounding computations.

*/

static dw10_t dfad(register dw10_t a,
		   register dw10_t b,
		   int negb)
{
    register int i, expa, expb;
    register uint18 rbit;		/* May be 1<<17 */
    register dw10_t rd;

    SFSETUP(expa, a.HI);		/* Set up exponent and fraction */
    SFSETUP(expb, b.HI);		/* Ditto for other arg */

    if (negb)				/* If actually subtracting B, */
	op10m_dmovn(b);			/* negate double fraction now */

    /* Now see which exponent is smaller; get larger in A, smaller in B.
    ** This needs to be done even if one of the args is completely zero, so
    ** that the other is set up in A for the normalization code.
    */
    if ((i = expa - expb) < 0) {
	dw10_t tmp;			/* Swap args */
	int etmp;
	etmp = expb, tmp = b;
	expb = expa, b = a;
	expa = etmp, a = tmp;
	i = -i;
    }
    DEBUGPRF(("DFAD A: (%o) %lo,%lo,%lo,%lo  B: (%o) %lo,%lo,%lo,%lo i=%d\n",
			expa, DBGV_D(a), expb, DBGV_D(b), i));

    /* Now right-shift B by the number of bits in i, then add to A.
    */
    rd = x_ashc(b, 70-i);	/* Get bottom 2 words of a 4-word ASHC */
    b = x_ashc(b, -i);		/* Shift B for adding */
    DEBUGPRF(("DFAD B: (%o) %lo,%lo,%lo,%lo  R: %lo,%lo,%lo,%lo\n",
			expa, DBGV_D(b), DBGV_D(rd)));

    a = x_dadd(a, b);		/* Add into A */

    DEBUGPRF(("  preN: (%o) %lo,%lo,%lo,%lo  R: %lo,%lo,%lo,%lo\n",
			expa, DBGV_D(a), DBGV_D(rd)));

    /* A now has sum, with the low-order 3rd word in RD's high word.
    ** The sign of A is correct; that of RD is irrelevant.
    ** Check high 10 bits to determine normalization step.
    */
    rbit = 0;
    switch (LHGET(a.HI) >> 8) {
	case 00004>>2:		/* Normalized positive fraction */
	    DEBUGPRF(("WINNING POSITIVE\n"));
	    rbit = (LHGET(rd.HI) & (HSIGN>>1));
	    break;

	case 00010>>2:		/* Overflow of pos fraction */
	case 00014>>2:
	    DEBUGPRF(("OVERFLOW POSITIVE\n"));
	    rbit = RHGET(a.LO) & 01;
	    a = x_ashc(a, -1);		/* Normalize by shifting right */
	    ++expa;			/* And adding one to exponent */
	    break;

	case 07760>>2:		/* Overflow of neg fraction, maybe zero mag */
	case 07764>>2:		/* Overflow of neg fraction, nonzero mag */
	    DEBUGPRF(("OVERFLOW NEGATIVE\n"));
	    rbit = RHGET(a.LO) & 01;
	    a = x_ashc(a, -1);	/* First normalize into 0777 case, */
	    rd = x_ashc(rd, -1);	/* Shift RD as well, */
	    if (rbit)		/* Copy bit lost from A */
		LHSET(rd.HI, (LHGET(rd.HI) | (HSIGN>>1)));
	    else
		LHSET(rd.HI, (LHGET(rd.HI) & ~(HSIGN>>1)));
	    ++expa;		/* then fall thru below to check! */

	case 07770>>2:		/* Normalized negative fraction (probably) */
	    rbit = LHGET(rd.HI) & (HSIGN>>1);
	    DEBUGPRF(("DFAD:NEG: (%o) %lo,%lo,%lo,%lo  %lo,%lo,%lo,%lo  r=%lo\n",
		expa, DBGV_D(a), DBGV_D(rd), (long)rbit ));

	    if (!(LHGET(rd.HI)&(HMASK>>2)) && !RHGET(rd.HI)
	      && !wmagskipn(rd.LO)) {
		rbit = 0;		/* No round if no other bits */
		DEBUGPRF(("NEG - RND DROPPED - "));
	    }
	    if ((LHGET(a.HI)&0777) || RHGET(a.HI) || rbit || wmagskipn(a.LO)) {
		DEBUGPRF(("WINNING NEGATIVE\n"));
		break;
	    }
	    DEBUGPRF(("OVERFLOW NEG (ZERO MAG)\n"));
	    rbit = 0;
	    LHSET(a.HI, 0777400);	/* Pretend right-shifted 1 bit */
	    ++expa;
	    break;

	/* Anything else is an unnormalized result and must be handled
	** using fully general case.
	*/
	case 0:
	    if (!wskipn(a.HI) && !wmagskipn(a.LO) && !wmagskipn(rd.HI)) {
		DEBUGPRF(("DFAD: S: ZERO RESULT!\n"));
		return dw10zero;
	    }
	default:
	    /* Must mess with RD, sigh.  Ensure it has correct sign bit,
	    ** then find the first fraction bit in the 3-word sum.
	    */
	    b.LO = rd.HI;
	    if (wskipl(a.HI)) {
		op10m_signset(b.LO);
		if ((i = op10ffo(op10setcm(a.HI))) >= 36) {
		    if ((i = op10ffo(op10setcm(a.LO))) >= 36) {
			i = op10ffo(op10setcm(b.LO)) + 70;
		    } else i += 35;
		}
	    } else {
		op10m_signclr(b.LO);		/* Clear sign bit */
		if ((i = op10ffo(a.HI)) >= 36) {
		    if ((i = op10ffo(a.LO)) >= 36) {
			i = op10ffo(b.LO) + 70;
		    } else i += 35;
		}
	    }
	    DEBUGPRF(("RENORMALIZATION: first %c-bit at %d = ASHC %d\n",
			(wskipl(a.HI) ? '0' : '1'), i, i-9));

	    /* i now has position of first fraction bit; move it into place. */
	    b.HI = a.LO;	/* Duplicate middle word for easy shift */
	    if ((i -= 9) < 35) {	/* Right shift or small left shift */
		a = x_ashc(a, i);
		b = x_ashc(b, i);
		if (i > 0) a.LO = b.HI;
	    } else if (i < 70) {	/* One-word or larger shift */
		a = x_ashc(b, i-35);
		op10m_setz(b.LO);
	    } else {			/* Two-word or larger shift */
		a.HI = x_ash(b.LO, i-70);
		op10m_setz(a.LO);
		op10m_setz(b.LO);
	    }
	    DEBUGPRF(("DFAD: S1:(%o) %lo,%lo,%lo,%lo  R: %lo,%lo\n",
			expa-i, DBGV_D(a),
			(long)LHGET(b.LO), (long)RHGET(b.LO)));

	    /* Now find roundup bit, and do fixups for negative results */
	    rbit = LHGET(b.LO) & (HSIGN>>1);	/* Find roundup bit */
	    if (wskipl(a.HI)) {			/* Special checks for neg */
		if (!(LHGET(b.LO)&(HMASK>>2)) && !RHGET(b.LO))
		    rbit = 0;			/* No round if no other bits */
		if (LHGET(a.HI) == 0777000 && !RHGET(a.HI) && !wmagskipn(a.LO)
		  && !rbit) {
		    ++expa, a = x_ashc(a, -1);	/* Make LH 0777400 */
		    DEBUGPRF(("    NEG FIXUP\n"));
		}
	    }
	    expa -= i;			/* Adjust exponent */
	    break;
    }

    /* At this point, number in A is normalized, and rbit is set if a
    ** roundup increment must be done (which in rare cases can trigger
    ** a single-shift renormalization).
    */
    if (rbit) {
	DEBUGPRF(("ROUNDUP DONE"));
	a = op10dfinc(a);		/* Add then check for overflow */
	if (wskipl(a.HI)) {		/* Negatives are messy */
	    if (LHGET(a.HI) & 0400) {	/* If first fract bit is "wrong" */
		if ((LHGET(a.HI) & 0377) || RHGET(a.HI)	/* See if really OK */
		  || wmagskipn(a.LO)) {
		    --expa, a = x_ashc(a, 1);		/* Nope, shift up */
		    DEBUGPRF((" - NEG RENORMALIZED!"));
		}
	    }
	} else if (!(LHGET(a.HI) & 0400)) {
	    ++expa, LHSET(a.HI, 0400);		/* Positive overflow, fix up */
	    DEBUGPRF((" - POS RENORMALIZED!"));
	}
	DEBUGPRF(("\n"));
    }
    DEBUGPRF(("DFAD: S: (%o) %lo,%lo,%lo,%lo\n", expa, DBGV_D(a)));

    /* Ta da!  Put exponent back in (it may have underflowed or overflowed
    ** but nothing we can do about it).
    */
    if (wskipl(a.HI))	/* If sign set, put ones-complement of exp */
	LHSET(a.HI, (LHGET(a.HI) ^ ((h10_t)(expa & SFEMAGS) << 9)));
    else
	LHSET(a.HI, (LHGET(a.HI) | ((h10_t)(expa & SFEMAGS) << 9)));
#if IFFLAGS
    if (expa < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (expa > SFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif
    op10m_signclr(a.LO);	/* Double flt always zeros low sign */
    return a;
}


dw10_t op10dfmp(register dw10_t a, register dw10_t b)
{
    register qw10_t q;
    register int exp, expb;
    register int sign;
    int signb;

    /* Make operands positive and remember sign of result
    ** These macros guarantee a positive number of at most 62 magnitude bits,
    ** even for the max-neg-# case.
    */
    SF_DPOSSETUP(sign, exp, a);
    SF_DPOSSETUP(signb, expb, b);

    DEBUGPRF(("DFMP: A: (%o) %lo,%lo,%lo,%lo  B: (%o) %lo,%lo,%lo,%lo\n",
			exp, DBGV_D(a),	expb, DBGV_D(b)));


    sign ^= signb;		/* Set 0 if both were negative */
    exp += expb - SFEEXCESS;	/* Find resulting exponent */

    /* Have two 62-bit magnitude numbers, multiply them (124-bit product) */

    q = x_dmul(a, b);		/* Multiply the numbers, get quad result */

    DEBUGPRF(("DFMP Q: %lo,%lo,%lo,%lo|%lo,%lo,%lo,%lo exp=%o (%d.)\n",
			DBGV_Q(q), exp, exp));

    /* Since there are 124 bits of product with 35 magnitude bits in each
    ** word, the high word should contain 124-105 = 19 high mag bits.
    ** Thus the top 36-19=17 bits of product should be empty.
    ** If everything was normalized, the next bit (start of fraction, B17)
    ** should be 1 and we just shift everything left by 17-9=8 bits to get
    ** in place.
    ** If not, must normalize by looking for first bit of fraction.
    */
    if (!op10m_tlnn(q.D0.HI, 1<<(H10BITS-(1+17)))) {	/* Test B17 */
	/* No high fraction bits set... see if result was 0 */
	if (!RHGET(q.D0.HI) && !wskipn(q.D0.LO) && !dskipn(q.D1))
	    return q.D0;		/* Zero product, return 0.0 */
	do {
	    q = x_qash1(q);		/* ASH quad up one bit */
	    --exp;			/* Adjust exponent accordingly */
	} while (!op10m_tlnn(q.D0.HI, 1<<(H10BITS-(1+17))));
					/* Until get high fract bit */
    }

    /* Fraction in expected place, ASH 17-9=8 bits left to align for exp */
    q.D0 = x_ashc(q.D0, 8);		/* Do ASHC bringing in 8 zeros */
    op10m_iori(q.D0.LO,			/* Add 8 bits from 3rd wd */
		(LHGET(q.D1.HI) >> ((HBITS-1)-8)));

    DEBUGPRF(("DFMP2Q: %lo,%lo,%lo,%lo|%lo,%lo,%lo,%lo exp=%o (%d.)\n",
			DBGV_Q(q), exp, exp));

    /* Now do rounding.  See if bit lower than LSB is set (this is the next
    ** bit after the ones added from 3rd word).
    ** If so, add one to LSB.  This will overflow only if the fraction was
    ** all ones, resulting in a roundup to all zeros except the low exponent
    ** bit.  If so, re-normalize by simply resetting the LH to right value
    ** and adjusting the exponent.
    */
    if (LHGET(q.D1.HI) & (1 << ((HBITS-1)-(8+1))) ) {
	/* Special hack in attempt to emulate hardware.  Don't round if
	** result is going to be negative and there are no more low bits
	** set in the product (to the right of the round bit).
	*/
	if (sign && !op10m_tlnn(q.D1.HI, 0377)
	  && !RHGET(q.D1.HI) && !wskipn(q.D1.LO)) {
	    DEBUGPRF(("DFMP ROUNDUP NEGSKIPPED\n"));
	} else {
	    DEBUGPRF(("DFMP ROUNDUP\n"));
	    q.D0 = op10dfinc(q.D0);
	    if (op10m_tlnn(q.D0.HI, SFELHF)) {	/* If carried out of fract */
		LHSET(q.D0.HI, SFFHIBIT);	/* reset to proper value */
		++exp;			/* and adjust exponent */
		DEBUGPRF(("DFMP POS RENORM\n"));
	    }
	}
    }

    /* Fraction all done, just add exponent back now */
    LHSET(q.D0.HI, (LHGET(q.D0.HI) | ((h10_t)(exp & SFEMAGS) << SFFBITS)));
#if IFFLAGS
    if (exp < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (exp > SFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif

    if (sign)
	op10m_dmovn(q.D0);
    return q.D0;
}

/* DFDV - note no-divide case simply returns 1st arg instead of
**	failing explicitly.
*/

dw10_t op10dfdv(register dw10_t a, register dw10_t b)
{
    qw10_t q;
    register int exp, i;
    register int expb, sign, signb;

#if 1	/* New code */
    /* Make operands positive and remember sign of result */
    /* The KLX appears to have a peculiarly rigid filter for the
    ** args of DFDV which the other DF ops don't share.
    ** If the 2nd arg (divisor) is 0, result is always no-divide.
    ** else, if the 1st arg is 400000,,0 0,,0
    **		then the result is a No-Divide (regardless of 2nd arg).
    ** else, if the 1st arg has no fraction bits set,
    **		then the result is always 0, even though this could be
    **		an unnormalized value.  This test is NOT applied
    **		to the divisor!!!
    */
    SF_DPOSSETUP(signb, expb, b);	/* Set up 2nd arg normally */
    if (!wskipn(b.HI) && !wmagskipn(b.LO)) {	/* Check for zero divisor */
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	return a;		/* Ugh, No Divide error.  Return original op */
    }

    /* Now set up A.  Note a different mechanism is used in order to emulate
    ** the way that the KL apparently chops off all exponent bits when
    ** extracting the fraction bits, thus turning an unnormalized negative
    ** number of (e.g.) 777000,,0 into 0,,0 instead of 1000,,0.
    **/
    q.D0 = a;			/* Set up A in Q (preserves orig value) */
    if (wskipl(q.D0.HI)) {
	op10m_dmovn(q.D0);
	if (wskipl(q.D0.HI)) {
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	    return a;		/* Ugh, No Divide error.  Return original op */
	}
	sign = !signb;
    } else sign = signb;
    exp = SFEGET(q.D0.HI);	/* Get sign and exponent */
    op10m_tlz(q.D0.HI, SFELHF);	/* Chop out exponent bits */
    if (!wskipn(q.D0.HI) && !wmagskipn(q.D0.LO))
	return dw10zero;	/* No fraction bits, return 0 */

    DEBUGPRF(("DFDV: A: (%o) %lo,%lo,%lo,%lo  B: (%o) %lo,%lo,%lo,%lo\n",
	exp, DBGV_D(q.D0), expb, DBGV_D(b)));

    exp += SFEEXCESS - expb;	/* Do division on exponents */

    /* Now guaranteed to have two positive 62-bit magnitude numbers. */


#else /* end new, start old */

    /* Make operands positive and remember sign of result */
    q.D0 = a;
    if (sign = wskipl(a.HI)) {
	op10m_dmovn(q.D0);

	/* KLX always fails immediately with No Divide if dividend is
	** 400000,,0 ? 0,,0
	*/
	if (wskipl(a.HI)) {
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	    return a;		/* Ugh, No Divide error.  Return original op */
	}
    }
    if (wskipl(b.HI)) {
	op10m_dmovn(b);
	sign = !sign;		/* Invert sense of sign flag */
    }

    /* Subtract exponents.
    */
    exp = 0200 + (LHGET(q.D0.HI) >> 9) - (LHGET(b.HI) >> 9);

    /* Get two 62-bit magnitude numbers and divide them.
    ** In order to ensure that the division can succeed, we pre-multiply the
    ** divisor by 2.  This is compensated for by adjusting the final exponent.
    ** This can fail with No Divide if the operands aren't normalized.
    */
    /* Mask out everything but fraction */
    LHSET(q.D0.HI, (LHGET(q.D0.HI) & 0777));
    LHSET(b.HI, (LHGET(b.HI) & 0777));
#endif	/* Old code */

    /* Divide them.
    ** Note that floating fractions, if normalized, will always have their most
    ** significant bit in the same place, which confuses the ordinary division
    ** test.
    ** So, in order to ensure that the division can succeed, we pre-multiply
    ** the divisor by 2 (to make it a 63-bit value).
    ** This is compensated for by adjusting the final exponent.
    ** This can fail with No Divide if the operands aren't normalized.
    **
    ** Note that the udcmpl test relies on both args being positive and the
    ** low sign bit being clear.
    */
    op10m_signclr(q.D0.LO);		/* Ensure valid for udcmpl test */
    op10m_signclr(b.LO);
    if (!op10m_udcmpl(q.D0, b)) {	/* If q.D0 >= b */
	DEBUGPRF(("DFDV Divisor pre-multiply!\n"));
	b = x_ashc(b, 1);		/* Make divisor be bigger */
	++exp;
	if (!op10m_udcmpl(q.D0, b)) {	/* If q.D0 >= b */
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	    return a;		/* Ugh, No Divide error.  Return original op */
	}
    }
    q.D1 = dw10zero;		/* Set up dividend */
#if 0
    DEBUGPRF(("DFDV: N: %lo,%lo,%lo,%lo|%lo,%lo,%lo,%lo  D: %lo,%lo,%lo,%lo\n",
			DBGV_Q(q), DBGV_D(b)));
#endif

    if (!qdivstep(&q, b, 63)) {
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	return a;
    }
    exp += 70-63;		/* Compensate for dividing only 63 steps */

    DEBUGPRF(("DFDV: Q: %lo,%lo,%lo,%lo R: %lo,%lo,%lo,%lo  E: %o\n",
			DBGV_Q(q), exp));

    if (!dskipn(q.D0))		/* Check for zero result */
	return q.D0;

    /* To compensate for possible unnormalized operands and roundup carry,
    ** we do a completely general normalization and round.  Note we've already
    ** tested for 0 so bits definitely exist.
    */
    i = adffo(q.D0) - 9;		/* Find first one-bit, derive shift */
    if (i < 0) {			/* If right-shifting, */
	/* bit only needs 9 bits, so "int" is safe */
	register int bit = 1 << (-i-1);	/* Find last bit that will vanish */
	if (RHGET(q.D0.LO) & bit) { 		/* If need to round up, */
	    static dw10_t dwrnd;
	    LRHSET(dwrnd.HI, 0, 0);
	    LRHSET(dwrnd.LO, 0, bit);
	    q.D0 = x_dadd(q.D0, dwrnd);	/* do it */
	    DEBUGPRF(("ROUNDED UP!  Bit %o\n", bit));
	    if (i != (adffo(q.D0)-9)) {		/* See if carried */
		DEBUGPRF(("CARRIED!\n"));
		if (wskipl(q.D0.HI)) {	/* Yes, into sign bit? */
		    DEBUGPRF(("CARRIED INTO SIGN!\n"));
		    q.D0 = x_ashc(q.D0, -1);	/* Yes, yecch!  Special */
		    op10m_signclr(q.D0.HI);	/* handling to undo damage */
		    --exp;
		} else --i;		/* Carried but compensation trivial */
	    }
	}
    }
    DEBUGPRF(("%sNORMALIZATION: %d\n", ((i != -8) ? "AB" : ""), i));

    q.D0 = x_ashc(q.D0, i);		/* Shift the fraction into place */
    exp -= i + 8;			/* Adjust exponent for shift */

    /* Fraction all done, just add exponent back now */
    LHSET(q.D0.HI, (LHGET(q.D0.HI) | ((h10_t)(exp & 0377) << 9)));
#if IFFLAGS
    if (exp < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (exp > 0377) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif

    if (sign)
	op10m_dmovn(q.D0);
    return q.D0;
}

static int qdivstep(qw10_t *aq,
		    register dw10_t d,
		    register int nmagbits)
{
    /*register*/ qw10_t qw;
    dw10_t quot;
    register int qbit;

    /* Apply initial test - high order double must be less than abs val of 
    ** divisor, else overflow & no divide.
    */
    qw = *aq;			/* Get quadword */
    qw.D0 = x_dsub(qw.D0, d);
    qbit = (wskipl(qw.D0.HI) ? 0 : 1);	/* Get first quotient bit */
    if (qbit)
	return 0;		/* Overflow, no divide */

    DEBUGPRF(("qdivstep: %lo,%lo,%lo,%lo|%lo,%lo,%lo,%lo  %lo,%lo,%lo,%lo  nmags=%d, qbit=%o\n",
		DBGV_Q(qw), DBGV_D(d), nmagbits, qbit));


    /* First quotient bit (future sign bit) is always 0.
    ** Because quotient is always positive, its left-shift in the loop
    ** is simple enough to do inline.
    */
    op10m_setz(quot.HI);		/* Clear current & future signs */
    op10m_setz(quot.LO);		/* Clear low word in case small N */
    while (--nmagbits >= 0) {		/* Loop N times for N magnitude bits */
	qw = x_qash1(qw);		/* Shift new bit into dividend */
	qw.D0 = (qbit ? x_dsub(qw.D0, d) : x_dadd(qw.D0, d));
	qbit = (wskipl(qw.D0.HI) ? 0 : 1);	/* Get new quotient bit */

	/* Shift quotient left by 1 */
	LSHIFTM(quot.HI, 1);		/* Shift both words left 1 */
	LSHIFTM(quot.LO, 1);
	if (wskipl(quot.LO))		/* If low sign now set, */
	    op10m_iori(quot.HI, 1);	/* that becomes low bit of high wd */

	op10m_iori(quot.LO, qbit);	/* Add qbit to right end of quotient */

	DEBUGPRF(("qdivstep %d: D: %lo,%lo,%lo,%lo|%lo,%lo,%lo,%lo  Q: %lo,%lo,%lo,%lo \n",
		nmagbits, DBGV_Q(qw), DBGV_D(quot)));
    }
    op10m_signclr(quot.LO);	/* Ensure low sign clear, compensate for
				** fact weren't doing true ASH in loop.
				*/

    qw.D1 = qw.D0;	/* Return quotient & remainder in right places */
    qw.D0 = quot;
    *aq = qw;
    return 1;
}

#if OP10_GFMT

/* Floating-point G-format Double-precision arithmetic.  More moby hair!
*/

static dw10_t x_gfad(dw10_t, dw10_t, int);	/* Aux for GFAD/GFSB */
static dw10_t gfnorm(int, dw10_t, int);
static int aqfffrac(dw10_t, dw10_t);


dw10_t op10gfad(dw10_t a, dw10_t b)
{
    return x_gfad(a, b, 0);
}

dw10_t op10gfsb(dw10_t a, dw10_t b)
{
    op10m_dmovn(b);		/* This appears to work better?!?! */
    return x_gfad(a, b, 0);
}

/* Common GFAD/GFSB routine.

	You'd expect G format to behave more or less the same as D
format.  Well, welcome to reality.  The G format implementation is
different in subtle ways which mainly reveal themselves in the handling
of unnormalized operands.  This may be because it was a late addition
to the KL and must depend completely on ucode software rather than any
inherent hardware features (such as those that might exist for handling
F or D formats).

	In particular, it appears that if the exponent difference
between the operands is larger than 72, the KL simply normalizes the
largest operand and returns that, completely ignoring the other value.

	Also, seems to have same problem as DFAD/DFSB with regard to hanging
on to the lower bits for full normalization and rounding.
	From empirical testing it appears that 2 additional bits are
needed from a "4th word", as for DFAD (36*3 + 3 real bits, or
24+35+35+2 = 96 total magnitude bits).  Simplest to just use a full
extra 35-bit word; HOWEVER, its peculiar behavior in certain cases
suggests that lower bits beyond those 2 are not used in normalization,
yet are still used in rounding -- that is, if any "lost" bits are set, the
higher bits are immediately rounded prior to the normalization shift.
That's not quite right, but they definitely have SOME effect.

	Face it, the only way to emulate this miserable lossage
accurately would be to understand and do exactly what the KL ucode
does.  But is that really a good idea?  I mean, it just affects the low
bit, usually with unnormalized operands, and it would be much simpler
to just use better algorithms - which would return maximally accurate
results even with unnormalized operands, instead of spazzing in
mysterious ways.

*/

static dw10_t x_gfad(register dw10_t a, register dw10_t b, int negb)
{
    register int i, expa, expb;
    register uint18 rbit;		/* May be 1<<17 */
    register dw10_t rd;

    GFSETUP(expa, a.HI);	/* Set up exp & fract for A */
    GFSETUP(expb, b.HI);	/* Set up exp & fract for B */

    if (negb)			/* If actually subtracting B, */
	op10m_dmovn(b);		/* negate fraction now */

    /* Now see which exponent is smaller; get larger in A, smaller in B.
    ** This needs to be done even if one of the args is completely zero, so
    ** that the other is set up in A for the normalization code.
    */
    if ((i = expa - expb) < 0) {
	dw10_t tmp;			/* Swap args */
	int etmp;
	etmp = expb, tmp = b;
	expb = expa, b = a;
	expa = etmp, a = tmp;
	i = -i;
    }

    /* Special KL G format crock. */
    if (i > 72) {
	DEBUGPRF(("GFAD i>72 NORM: i=%d (%o-%o) A: %lo,%lo,%lo,%lo\n",
			    i, expa, expb, DBGV_D(a)));
	return gfnorm(expa, a, 1);
    }

    /* Now right-shift B by the number of bits in i, then add to A. */
    rd = x_ashc(b, 70-i);	/* Get bottom 2 words of a 4-word ASHC */
    b = x_ashc(b, -i);		/* Shift B for adding */
    DEBUGPRF(("GFAD A: (%o) %lo,%lo,%lo,%lo  B: (%o) %lo,%lo,%lo,%lo i=%d\n",
		expa, DBGV_D(a), expb, DBGV_D(b), i)); 

    a = x_dadd(a, b);		/* Add into A */

    DEBUGPRF(("  preN: (%o) %lo,%lo,%lo,%lo  R: %lo,%lo,%lo,%lo\n",
		expa, DBGV_D(a), DBGV_D(rd))); 

    /* A now has sum, with the low-order 3rd word in RD's high word.
    ** The sign of A is correct; that of RD is irrelevant.
    ** Check high bits to determine normalization step, by looking
    ** at the bits of the exponent field PLUS the high fraction bit.
    */
    rbit = 0;
    switch (LHGET(a.HI) >> (H10BITS-(GFEBITS+1))) {
	case 000004>>2:		/* Normalized positive fraction */
	    DEBUGPRF(("WINNING POSITIVE\n"));
	    rbit = (LHGET(rd.HI) & (HSIGN>>1));
	    break;

	case 000010>>2:		/* Overflow of pos fraction */
	case 000014>>2:
	    DEBUGPRF(("OVERFLOW POSITIVE\n"));
	    rbit = RHGET(a.LO) & 01;
	    a = x_ashc(a, -1);		/* Normalize by shifting right */
	    ++expa;			/* And adding one to exponent */
	    break;

	case 077760>>2:		/* Overflow of neg fraction, maybe zero mag */
	case 077764>>2:		/* Overflow of neg fraction, nonzero mag */
	    DEBUGPRF(("OVERFLOW NEGATIVE\n"));
	    rbit = RHGET(a.LO) & 01;
	    a = x_ashc(a, -1);	/* First normalize into 077770 case, */
	    rd = x_ashc(rd, -1);	/* Shift RD as well, */
	    if (rbit)		/* Copy bit lost from A */
		LHSET(rd.HI, (LHGET(rd.HI) | (HSIGN>>1)));
	    else
		LHSET(rd.HI, (LHGET(rd.HI) & ~(HSIGN>>1)));
	    ++expa;		/* then fall thru below to check! */

	case 077770>>2:		/* Normalized negative fraction (probably) */
	    rbit = LHGET(rd.HI) & (HSIGN>>1);
	    DEBUGPRF(("GFAD:NEG: (%o) %lo,%lo,%lo,%lo  %lo,%lo,%lo,%lo  r=%lo\n",
		expa, DBGV_D(a), DBGV_D(rd), (long)rbit ));

	    if (rbit && !(LHGET(rd.HI)&(HMASK>>2)) && !RHGET(rd.HI)
	       && !wmagskipn(rd.LO)) {
		rbit = 0;		/* No round if no other bits */
		DEBUGPRF(("NEG - RND DROPPED - "));
	    }
	    if (op10m_tlnn(a.HI,GFFMASK) || RHGET(a.HI)
	      || rbit || wmagskipn(a.LO)) {
		DEBUGPRF(("WINNING NEGATIVE\n"));
		break;
	    }
	    DEBUGPRF(("OVERFLOW NEG (ZERO MAG)\n"));
	    rbit = 0;
	    LHSET(a.HI, 0777740);	/* Pretend right-shifted 1 bit */
	    ++expa;
	    break;

	/* Anything else is an unnormalized result and must be handled
	** using fully general case.
	*/
	case 0:
	    if (!wskipn(a.HI) && !wmagskipn(a.LO) && !wmagskipn(rd.HI)) {
		DEBUGPRF(("GFAD: S: ZERO RESULT!\n"));
		return dw10zero;
	    }
	default:
	    /* Must mess with RD, sigh.  Ensure it has correct sign bit,
	    ** then find the first fraction bit in the 3 (now 4) word sum.
	    */
	    if (wskipl(a.HI)) {
		op10m_signset(rd.HI);
		if ((i = op10ffo(op10setcm(a.HI))) >= 36) {
		    if ((i = op10ffo(op10setcm(a.LO))) >= 36) {
		/*	op10m_signset(rd.HI); */
			if ((i = op10ffo(op10setcm(rd.HI))) >= 36) {
			    op10m_signset(rd.LO);
			    i = op10ffo(op10setcm(rd.LO)) + 70 + 35;
			} else i += 70;
		    } else i += 35;
		}
	    } else {
		op10m_signclr(rd.HI);		/* Clear sign bit */
		if ((i = op10ffo(a.HI)) >= 36) {
		    if ((i = op10ffo(a.LO)) >= 36) {
		/*	op10m_signclr(rd.HI); */
			if ((i = op10ffo(rd.HI)) >= 36) {
			    op10m_signclr(rd.LO);
			    i = op10ffo(rd.LO) + 70 + 35;
			} else i += 70;
		    } else i += 35;
		}
	    }

	    /* i now has position of first fraction bit; move it into place. */
	    i -= GFEBITS;		/* Find actual shift to perform */
	    DEBUGPRF(("RENORMALIZATION: %d ", i));

	    /* G format hack, check for shift of 72.  This catches
	    ** the case of a negative number with all ones, as well as
	    ** limiting the positive all-zero case more severely than the
	    ** "case 0:" test above.
	    */
	    if (i > 72) {
		DEBUGPRF(("- i>72, RET 0\n"));
		return dw10zero;
	    }

#if 1	/* Test to see if just 2 extra bits are sufficient */

	    if (i > 35 && (LHGET(rd.LO) & 0300000)) {
		a = x_ashc(a, 2);
		op10m_tro(a.LO, (LHGET(rd.HI) >> (18-3)) & 03);
		rd = x_ashc(rd, 2);
		i -= 2;
		expa -= 2;

#if 1
	    /* This code implements "immediate rounding" such that if
	    ** normalization would use any bits from rd.LO (specifically
	    ** the low 33 bits), they are zapped (as above) but
	    ** the high bit is used for an immediate roundup.
	    **
	    ** Try using this only for negative numbers (gah)
	    */
		if (wskipl(a.HI) && (i > 35) && (LHGET(rd.LO)&(H10SIGN>>1))) {
		    DEBUGPRF((" - IMM ROUNDUP"));
		    op10m_inc(rd.HI);		/* Round up */
		    if (!wmagskipn(rd.HI))	/* Propagate */
			a = op10dfinc(a);
		    /* Ugh, must recompute i in case roundup changed it */
		    i = aqfffrac(a, rd) - GFEBITS;
		}
		op10m_setz(rd.LO);
#endif

	    }
#endif
#if 0	/* Test to see if extra 35 bits are sufficient.
	** Answer: Yes, but extra precision prevents emulating KL lossage.
	*/
	    if (i >= 35) {
		a = x_ashc(a, 35);
		a.LO = rd.HI;
		rd.HI = rd.LO;
		i -= 35;
		expa -= 35;
	    }
#endif

	    b.LO = rd.HI;	/* Set up for old code */
	    b.HI = a.LO;	/* Duplicate middle word for easy shift */
	    if (i < 35) {	/* Right shift or small left shift */
		a = x_ashc(a, i);
		b = x_ashc(b, i);
		if (i > 0) a.LO = b.HI;
	    } else if (i < 70) {	/* One-word or larger shift */
		a = x_ashc(b, i-35);
		op10m_setz(b.LO);
	    } else {			/* Two-word or larger shift */
		a.HI = x_ash(b.LO, i-70);
		op10m_setz(a.LO);
		op10m_setz(b.LO);
	    }

	    DEBUGPRF((" PostN S: %lo,%lo,%lo,%lo  R: %lo,%lo",
				DBGV_D(a), DBGV_W(b.LO)));

	    /* Now find roundup bit, and do fixups for negative results */
	    rbit = LHGET(b.LO) & (HSIGN>>1);	/* Find roundup bit */
	    if (wskipl(a.HI)) {			/* Special checks for neg */
		if (!(LHGET(b.LO)&(HMASK>>2)) && !RHGET(b.LO)) {
		    rbit = 0;			/* No round if no other bits */
		    DEBUGPRF((" - NEG RND DROPPED"));
		}
		if (LHGET(a.HI) == 0777700 && !RHGET(a.HI) && !wmagskipn(a.LO)
		  && !rbit) {
		    LHSET(a.HI, 0777740);	/* Effect an ASHC -1 */
		    ++expa;			/* Compensate for right-shft */
		    DEBUGPRF((" - NEG ZERO-MAG FIXUP"));
		}
	    }
	    DEBUGPRF(("\n"));
	    expa -= i;			/* Adjust exponent */
	    break;
    }

    /* At this point, number in A is normalized, and rbit is set if a
    ** roundup increment must be done (which in rare cases can trigger
    ** a single-shift renormalization).
    */
    if (rbit) {
	DEBUGPRF(("ROUNDUP DONE"));
	a = op10dfinc(a);		/* Add then check for overflow */
	if (wskipl(a.HI)) {		/* Negatives are messy */
	    if (op10m_tlnn(a.HI,GFFHIBIT)) {	/* If first fract bit is "wrong" */
		if (op10m_tlnn(a.HI, GFFMAGS)	/* See if really OK */
		  || RHGET(a.HI) || wmagskipn(a.LO)) {
		    --expa, a = x_ashc(a, 1);		/* Nope, shift up */
		    DEBUGPRF((" - NEG RENORMALIZED!"));
		}
	    }
	} else if (!op10m_tlnn(a.HI, GFFHIBIT)) {
	    ++expa, LHSET(a.HI, GFFHIBIT);	/* Positive overflow, fix up */
	    DEBUGPRF((" - POS RENORMALIZED!"));
	}
	DEBUGPRF(("\n"));
    }

    DEBUGPRF(("GFAD: S: (%o) %lo,%lo,%lo,%lo\n", expa, DBGV_D(a)));

    /* Ta da!  Put exponent back in (it may have underflowed or overflowed
    ** but nothing we can do about it).
    */
    if (wskipl(a.HI))	/* If sign set, put ones-complement of exp */
	 LHSET(a.HI, (LHGET(a.HI) ^ ((h10_t)(expa & GFEMAGS) << GFFBITS)));
    else LHSET(a.HI, (LHGET(a.HI) | ((h10_t)(expa & GFEMAGS) << GFFBITS)));
#if IFFLAGS
    if (expa < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (expa > GFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif
    op10m_signclr(a.LO);	/* Double flt always zeros low sign */
    return a;
}

#if 1	/* Only needed for special KL-lossage-emulation experiment */

static int aqfffrac(register dw10_t a, register dw10_t rd)
{
    register int i;

    if (wskipl(a.HI)) {
	op10m_signset(rd.HI);
	if ((i = op10ffo(op10setcm(a.HI))) >= 36) {
	    if ((i = op10ffo(op10setcm(a.LO))) >= 36) {
	/*	op10m_signset(rd.HI); */
		if ((i = op10ffo(op10setcm(rd.HI))) >= 36) {
		    op10m_signset(rd.LO);
		    i = op10ffo(op10setcm(rd.LO)) + 70 + 35;
		} else i += 70;
	    } else i += 35;
	}
    } else {
	op10m_signclr(rd.HI);		/* Clear sign bit */
	if ((i = op10ffo(a.HI)) >= 36) {
	    if ((i = op10ffo(a.LO)) >= 36) {
	/*	op10m_signclr(rd.HI); */
		if ((i = op10ffo(rd.HI)) >= 36) {
		    op10m_signclr(rd.LO);
		    i = op10ffo(rd.LO) + 70 + 35;
		} else i += 70;
	    } else i += 35;
	}
    }
    return i;
}
#endif /* special hack */


/* GFMP on the real KL has some problems with failure to set the floating
    underflow flag in some situations.
    For example, this test:
VERIFY ERROR: "gfmp"
  Input: gfmp 40,0,0,2 777740,0,0,0 => 600037,777777,377777,777776  [440000]
  E Ver: gfmp 40,0,0,2 777740,0,0,0 => 600037,777777,377777,777776  [440100]

	The emulation result (floating underflow) is correct.
*/

dw10_t op10gfmp(register dw10_t a, register dw10_t b)
{
    register qw10_t q;
    register int exp, expb;
    register int sign;
    int signb;

    /* Make operands positive and remember sign of result */
    GF_DPOSSETUP(sign, exp, a);
    GF_DPOSSETUP(signb, expb, b);

    DEBUGPRF(("GFMP: A: (%o) %lo,%lo,%lo,%lo  B: (%o) %lo,%lo,%lo,%lo\n",
			exp, DBGV_D(a), expb, DBGV_D(b)));

    sign ^= signb;		/* Set 0 if both were negative */
    exp += expb - GFEEXCESS;	/* Add exponents to effect multiply */

    /* Now have two 59-bit magnitude positive numbers.  Multiply them
    ** to get 118 bits of product in 4 words.
    */
    q = x_dmul(a, b);		/* Multiply the numbers, get quad result */

    DEBUGPRF(("GFMP Q: %lo,%lo,%lo,%lo|%lo,%lo,%lo,%lo exp=%o (%d.)\n",
			DBGV_Q(q), exp, exp));

    /* Since there are 118 bits of product with 35 magnitude bits in each
    ** word, the high word should contain 118-105 = 13 high mag bits.
    ** Thus the top 36-13=23 bits of product should be empty, and if
    ** everything was normalized, the next bit (B23, start of fraction)
    ** should be 1 and we just shift everything left by 23-12=11 bits to get
    ** things into place.
    ** If not, must normalize by looking for first bit of fraction,
    ** which is a painfully slow process.
    */
    if (!op10m_trnn(q.D0.HI, 1<<(H10BITS-(1+23-H10BITS)))) {	/* Test B23 */
	/* No high fraction bits set... see if result was 0 */
	if (!RHGET(q.D0.HI) && !wskipn(q.D0.LO) && !dskipn(q.D1))
	    return q.D0;		/* Zero product, return 0.0 */
	do {
	    q = x_qash1(q);		/* ASH quad up one bit */
	    --exp;			/* Adjust exponent accordingly */
	} while (!op10m_trnn(q.D0.HI, 1<<(H10BITS-(1+23-H10BITS))));
					/* Until get high fract bit */
    }

    /* Fraction in expected place, ASH 23-12=11 bits left to align for exp */
    q.D0 = x_ashc(q.D0, 11);		/* Do ASHC bringing in 11 zeros */
    op10m_iori(q.D0.LO,			/* Add 11 bits from 3rd wd */
		(LHGET(q.D1.HI) >> ((HBITS-1)-11)));

    /* Now do rounding.  See if bit lower than LSB is set (this is the next
    ** bit after the ones added from 3rd word).
    ** If so, add one to LSB.  This will overflow only if the fraction was
    ** all ones, resulting in a roundup to all zeros except the low exponent
    ** bit.  If so, re-normalize by simply resetting the LH to right value
    ** and adjusting the exponent.
    */
    if (LHGET(q.D1.HI) & (1 << ((HBITS-1)-(11+1))) ) {
	q.D0 = op10dfinc(q.D0);
	if (op10m_tlnn(q.D0.HI, GFELHF)) {	/* If carried out of fract */
	    LHSET(q.D0.HI, GFFHIBIT);	/* reset to proper value */
	    ++exp;			/* and adjust exponent */
	}
    }

    /* Fraction all done, just add exponent back now */
    LHSET(q.D0.HI, (LHGET(q.D0.HI) | ((h10_t)(exp & GFEMAGS) << GFFBITS)));
#if IFFLAGS
    if (exp < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (exp > GFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif

    if (sign)
	op10m_dmovn(q.D0);
    return q.D0;
}

/* GFDV - note no-divide case simply returns 1st arg instead of
**	failing explicitly.
*/

dw10_t op10gfdv(register dw10_t a, register dw10_t b)
{
    qw10_t q;
    register int exp;
    register int expb, sign, signb;
    dw10_t origa;

    /* Make operands positive and remember sign of result */
    /* The odd code here is similar to that for DFDV for the same reasons.
    */
    GF_DPOSSETUP(signb, expb, b);	/* Set up 2nd arg normally */
    if (!wskipn(b.HI) && !wmagskipn(b.LO)) {	/* Check for zero divisor */
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	return a;		/* Ugh, No Divide error.  Return original op */
    }

    /* Now set up A.  Note a different mechanism is used in order to emulate
    ** the way that the KL apparently chops off all exponent bits when
    ** extracting the fraction bits, thus turning an unnormalized negative
    ** number of (e.g.) 777000,,0 into 0,,0 instead of 1000,,0.
    **/
    origa = a;			/* Save original value of A */
    if (wskipl(a.HI)) {
	op10m_dmovn(a);
	if (wskipl(a.HI)) {
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	    return origa;	/* Ugh, No Divide error.  Return original op */
	}
	sign = !signb;
    } else sign = signb;
    exp = GFEGET(a.HI);		/* Get sign and exponent */
    op10m_tlz(a.HI, GFELHF);	/* Chop out exponent bits */
    if (!wskipn(a.HI) && !wmagskipn(a.LO))
	return dw10zero;	/* No fraction bits, return 0 */


    DEBUGPRF(("GFDV: A: (%o) %lo,%lo,%lo,%lo  B: (%o) %lo,%lo,%lo,%lo\n",
			exp, DBGV_D(a), expb, DBGV_D(b)));

    exp += GFEEXCESS - expb;	/* Do division on exponents */

    /* For G-Format, 59-bit fracts */

    /* Now have two 59-bit magnitude numbers; calculate A/B.
    ** In order to ensure that the division can succeed, we pre-multiply the
    ** divisor by 2.  This is compensated for by adjusting the final exponent.
    ** This can fail with No Divide if the operands aren't normalized.
    **
    ** Note that the udcmpl test relies on both args being positive and the
    ** low sign bit being clear.
    */
    op10m_signclr(a.LO);		/* Ensure valid for udcmpl test */
    op10m_signclr(b.LO);
    if (!op10m_udcmpl(a, b)) {		/* If a >= b */
	DEBUGPRF(("GFDV Divisor pre-multiply!\n"));
	b = x_ashc(b, 1);		/* Make divisor be bigger */
	++exp;
	if (!op10m_udcmpl(a, b)) {	/* If a >= b */
	    OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	    return origa;	/* Ugh, No Divide error.  Return original op */
	}
    }
    q.D0 = a;
    q.D1 = dw10zero;		/* Set up dividend */

    if (!qdivstep(&q, b, 59+1)) {
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_DIV);
	return a;
    }

    /* Because the fractions are aligned, the result should be aligned as
    ** well.  However, because of the extra divide step needed to compute
    ** the rounding bit, the result is left-shifted 1 bit greater than
    ** it should be.  Adjust for this by decrementing the exponent.
    */
    --exp;

    DEBUGPRF(("GFDV: Q: %lo,%lo,%lo,%lo R: %lo,%lo,%lo,%lo  exp=%o (%d.)\n",
			DBGV_Q(q), exp, exp));

    if (!dskipn(q.D0))		/* Check for zero result */
	return q.D0;

    /* To compensate for possible unnormalized operands and roundup carry,
    ** we do a completely general normalization and round.  Note we've already
    ** tested for 0 so bits definitely exist.
    ** If an exponent under/overflow happens, GFNORM will catch it.
    */
    if (sign) {
	a = gfnorm(exp, q.D0, 1);	/* Normalize and round */
	op10m_dmovn(a);
	return a;
    }
    return gfnorm(exp, q.D0, 1);	/* Normalize and round */
}
#endif /* OP10_GFMT */

/* Floating-point Double-precision Conversions */

#if OP10_GFMT

/* GFNORM - Returns normalized G-format floating point, given
**	double-word integer fraction and exponent to use if fraction were
**	already in proper position (first significant bit at B12)
**	Sets flags if exponent overflows.
*/
static dw10_t gfnorm(register int exp,
		     register dw10_t d,	/* Contains signed integer fraction */
		     int rnd)		/* TRUE to round result */
{
    register int i, rbit;	/* rbit only needs 12 bits, so "int" safe */
    register int sign;

    /* This algorithm interprets the fraction as a signed double integer
    ** and converts it immediately to positive form, if negative.
    ** This is much easier than trying to do the right thing for negatives
    ** all through the code.  See the old code for sfnorm() to see what I mean.
    */

    /* Check for negative fraction and convert to use
    ** positive value instead.  Note that the maximum negative number
    ** cannot be made positive, but this is compensated for by
    ** adjusting its exponent and fraction.
    ** sfnorm() gets away with doing nothing because it can use a logical
    ** shift.  This won't work here because ASHC looks at the sign bit,
    ** hence the trickery here to ensure the sign is always 0.
    */
    if (sign = wskipl(d.HI)) {
	op10m_dmovn(d);		/* Negate the fraction. */
	if (wskipl(d.HI)) {		/* If still negative, */
	    ++exp;			/* adjust exponent and value */
	    LHSET(d.HI, H10SIGN>>1);	/* to get positive # */
	}
    }
    i = adffo(d);		/* Find first 1-bit in double, skip low sign */

    DEBUGPRF(("GFNORM: %lo,%lo,%lo,%lo  i=%d, exp=%o\n",
			DBGV_D(d), i, exp));
    if (i >= 71)
	return dw10zero;	/* No sig bits, force result to zero */

    /* i now has position of first fraction bit; move it into place.
    */
    i -= GFEBITS;		/* Make relative to proper position, B12 */
    exp -= i;			/* Adjust exponent */
    if (rnd && i < 0)		/* If rounding OK and losing bits */
	rbit = RHGET(d.LO)	/* Find value of last bit that will be lost */
		& (1 << ((-i)-1));
    else rbit = 0;

    if (i)
	d = x_ashc(d, i);	/* Now do double arithmetic shift! */

    DEBUGPRF(("GNORMED: %lo,%lo,%lo,%lo  exp=%d, rbit=%o\n",
			DBGV_D(d), exp, rbit));

    /* At this point, number in w is normalized, and rbit is set if a
    ** roundup increment must be done.  In the rare case that this
    ** increment overflows, its value will be 100,,0 and will trigger
    ** a single-shift renormalization.
    */
    if (rbit) {
	d = op10dfinc(d);		/* Add then check for overflow */
	/* (above call sets low sign, but it gets flushed later) */

	if (LHGET(d.HI) > GFFMASK) {
	    ++exp;			/* Positive overflow, fix up */
	    LHSET(d.HI, GFFHIBIT);
	    DEBUGPRF(("ROUNDUP - RENORMED!\n"));
	} else
	    DEBUGPRF(("ROUNDUP\n"));
    }

    /* Fraction is normalized with exponent bits clear, can just OR
    ** exponent back in, and negate if needed to finish off result.
    */
    op10m_tlo(d.HI, ((h10_t)(exp & GFEMAGS) << GFFBITS));
    if (sign)
	op10m_dmovn(d);		/* This forces low sign 0 as well */
    else
	op10m_signclr(d.LO);	/* Unsure, so force low sign 0 */

    /* Ta da! Check before returning for overflow or underflow.
    ** Note that the negation, if one, can never cause either since a
    ** positive number can always be negated.
    */

#if IFFLAGS
    if (exp < 0) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
    else if (exp > GFEMAGS) OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
#endif

    return d;
}


/* GFSC - G Format floating scale.
**	Note peculiarity wherein GFSC tests to see if the first word is zero.
**	It should actually test the MAGNITUDE of the FRACTION, meaning
**	that if the number is negative, the test is made for zero bits
**	rather than one bits.
**	However, the real machine appears to always check for zero bits
**	in the ENTIRE word, including the sign bit.  So a negative number
**	or a non-zero exponent never triggers this case, regardless of
**	whether its high word fraction is all ones or all zeroes.  Sigh.
**
** GFSC on a KL screws up on the flags when given an E of 400000 applied to
**	a fairly low exponent.  It sets Overflow and Floating Overflow,
**	but not Floating Underflow!
**	This code correctly sets Floating Underflow.
**	(Oddly enough, FSC also messes this case up, but in a different way!)
**
** A real machine also screws up peculiarly in another way:
** VERIFY ERROR: "gfsc"
**  Input: gfsc   777777,777777,777777,777777    100  =>  777200,0,20,0  [0]
**  E Ver: gfsc   777777,777777,777777,777777    100  =>  777140,0,0,0  [0]
**	Note the bogus "20" in the low LH!
**      Admittedly the argument is unnormalized, but this is still weird.
*/
dw10_t
op10gfsc(register dw10_t d,
	 register h10_t h)
{
    register int exp;

#if 1	/* See note above about this stupid test. */
    if (!wskipn(d.HI))
	return dw10zero;	/* High word, return all zeros */
#endif

    /* Canonicalize exponent and fraction */
    GFSETUP(exp, d.HI);		/* Set up exp and fraction */

    /* Do special check of 1st wd fractional part.
    ** Note that at this point the exponent has been replaced by copies of
    ** the sign bit, so the test is simplified.
    */
#if 0
    if (wskipl(d.HI)) {
	register dw10_t nd;
	nd = d;
	op10m_dmovn(nd);	/* Do a double negate, then test */
	if (!wskipn(nd.HI))
	    return dw10zero;
    } else
	if (!wskipn(d.HI))
	    return dw10zero;	/* High fraction zero, return all zeros */
#endif

    /* Scale up exponent.  Note immediate arg is taken as 11 bits, not
    ** 8 bits as for all other similar ops!  (fsc, lsh, etc)
    */
    exp += (((h)&HSIGN) ? ((int)(h)|~03777) : ((int)(h)&03777));

    return gfnorm(exp, d, 0);	/* Put G-float together, don't round */
}

/* GSNGL - G Format to Single-precision Format, with round.
**
** Note: According to PRM, AC is unchanged if the G exponent is out of
**	range to begin with.  However, this appears not to be always true, per
**	following counter-example:
** VERIFY ERROR: "gsngl"
**   Input: gsngl  157373,175573,76600,173751  =>  373731,755732  [440100]
**   E Ver: gsngl  157373,175573,76600,173751  =>  157373,175573  [440100]
**
**	Note result was clobbered.  Doesn't even have the excuse of
**	unnormalized arg for this one!
*/
/* G -> S, with round, TRUE if succeeded */
int
op10xgsngl(register dw10_t *ad)
{
    register int exp;
    register dw10_t d;

    d = *ad;			/* Fetch the double */

    /* Canonicalize exponent and fraction */
    GFSETUP(exp, d.HI);		/* Set up exp and fraction */

    /* Special check for zero, needed to avoid setting flags */
    if (!exp && !wskipn(d.HI)) {
	ad->w[0] = w10zero;
	return 1;
    }

    /* Convert to new exponent range */
    exp = (exp - GFEEXCESS) + SFEEXCESS;

    /* Now see if OK for single-prec */
    DEBUGPRF(("GSNGL: new exp=%o (%d.)\n", exp, exp));
    if (exp < 0) {
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV+PCF_FXU);
	return 0;
    } else if (exp > SFEMAGS) {
	OP10_PCFSET(PCF_ARO+PCF_TR1+PCF_FOV);
	return 0;
    }

    /* Now what?  PRM seems to imply no normalization is done, just
    ** rounding.  However, it seems best to carry out a fully general
    ** normalize and round.
    ** Need to find out what real machine does.
    */

    /* The following code uses cleverness to take advantage of
    ** existing single-prec renormalization & round.
    ** First we shift the fraction up by 10 bits, not 3, so as to take
    ** maximum advantage of any rounding from the low bits.  This is
    ** particularly necessary for negative numbers, otherwise it's possible
    ** for sfnorm's neg->pos algorithm to get a spurious roundup for an
    ** all-ones fraction.
    ** This is then fed to sfnorm(), with the new exponent adjusted down
    ** by 10-3 to compensate for the extra bit shifts; sfnorm() knows how to
    ** do a right-shift to put things back in place.
    **  sfnorm() will set the flags if rounding overflows, just as PRM says.
    */
    d = x_ashc(d, 3+7);		/* Shift double fraction 10 places left */

    ad->w[0] = sfnorm(exp-7, d.HI, 1);	/* Return result in 1st wd of dbl */
    return 1;
}

/* S -> G */

dw10_t
op10gdble(register w10_t w)
{
    register dw10_t d;
    register int exp;
    register int sign;

    SF_POSSETUP(sign, exp, w);	/* Set up sign, exponent and fraction */
    d.HI = w;			/* Set high word */
    op10m_setz(d.LO);		/* Clear low */

    /* See if fraction is already normalized, to save the time of
    ** a full renormalization call.
    */
    if (op10m_tlnn(d.HI, SFFHIBIT)) {	/* Normalized?  Most common case */
	d = x_ashc(d, -3);			/* Right-shift 3 bits */
	exp = (exp - SFEEXCESS) + GFEEXCESS;	/* Convert up to G exponent */
	LHSET(d.HI, (LHGET(d.HI) | ((h10_t)(exp) << (H10BITS-GFEBITS))));
    } else if (!wskipn(d.HI)) {		/* Zero?  Next most common */
	return d;
    } else {
	/* Either unnormalized operand, or fraction is negative 400,,0.
	** Invoke slow but fully general normalization.
	*/
	d = gfnorm(((exp - SFEEXCESS) + GFEEXCESS - 3), d, 0);
    }
    if (sign)			/* If necessary, convert back */
	op10m_dmovn(d);		/* Note low sign always left 0 */
    return d;
}

/* Dfix -> G, with round */

dw10_t
op10dgfltr(register dw10_t d)
{
    /* Double already in arg, just need to provide normalizer with
    ** right exponent.  Binary point is GFFBITS+H10BITS+35 away from
    ** normal position.
    ** Note also that gfnorm() correctly handles max-neg-# case!
    */
    return gfnorm(GFEEXCESS+GFFBITS+H10BITS+35, d, 1);
}

/* fix -> G, with round (not really) */

dw10_t
op10gfltr(register w10_t w)
{
    register dw10_t d;

    d.HI = w;			/* Set up double-word */
    op10m_setz(d.LO);

    /* Provide exponent & normalize it.
    ** Binary point of the single-prec fix is to right of 1st word.
    ** This is GFFBITS+H10BITS (ie all fraction bits in 1st word)
    ** away from desired position, so add those to exponent.
    */
    return gfnorm(GFEEXCESS+GFFBITS+H10BITS, d, 0);
}


/* The following 2 functions are placeholders for GFIX{R} and GDFIX{R},
** which on a real KL are actually simulated in the monitor.
** These functions HAVE NOT BEEN VERIFIED yet!
*/

/* G -> Dfix, TRUE if succeeded */

int
op10xgdfix(register dw10_t *ad, int rnd)
{
    register dw10_t d = *ad;
    register int i, rbit, sign;

    GF_DPOSSETUP(sign, i, d);	/* Set up exponent, get positive fract */

    /* Apply exponent size test now */
    if (i > (GFEEXCESS+70)) {		/* If integer wd be too big, */
	OP10_PCFSET(PCF_ARO+PCF_TR1);	/* set flags */
	return 0;			/* and return arg unchanged */
    }
    
    /* OK to proceed!
    ** Fraction is normally a 24+35= 59-bit positive integer.  It will be a
    ** 60-bit value of 100,,0 ? 0 only if the float was negative with no
    ** one-bits in the fraction part (which is an unnormalized value).
    */
    i -= (GFEEXCESS+59);	/* Get # bits to shift fraction */
    if (i <= -60) {		/* If shifting it to oblivion */
	*ad = dw10zero;		/* just return zero, skip rigmarole */
	return 1;
    }

    if (rnd && i < 0) {		/* If rounding OK and losing bits */
	/* Find value of last bit that will be lost */
	if (i <= -35) {		/* Bit in high word? */
	    if (i <= -(35+H10BITS))	/* Bit comes from High LH? */
		rbit = LHGET(d.HI)
		    & (1 << ((-i)-(35+H10BITS+1)));
	    else		/* Bit comes from High RH */
		rbit = RHGET(d.HI)
		    & (1 << ((-i)-(35+1)));
	} else {		/* Bit in low word (shift -35 < i < 0) */
	    if (i <= -H10BITS)	/* Bit comes from Low LH? */
		rbit = LHGET(d.LO)
		    & (1 << ((-i)-(H10BITS+1)));
	    else		/* Bit comes from Low RH */
		rbit = RHGET(d.LO)
		    & (1 << ((-i)-1));
	}
    } else rbit = 0;

    if (i)
	d = x_ashc(d, i);	/* Do the shift */

    if (rbit) {			/* Need any rounding? */
	op10m_inc(d.LO);	/* Yep, trivial */
	if (!wmagskipn(d.LO))	/* Check for carry to high */
	    op10m_inc(d.HI);	/* Will never overflow this one */
    }

    /* Since returning a double fix, make sure result conforms to fixed-point
    ** sign convention (low sign is copy of high's)
    */
    if (sign)
	op10m_dneg(d);		/* Negate (sets low sign) */
    else
	op10m_signclr(d.LO);	/* Else clear low sign */

    *ad = d;			/* Return result */
    return 1;			/* Say succeeded */
}


/* G -> fix, TRUE if succeeded */

int
op10xgfix(register dw10_t *ad, int rnd)
{
    /* See code for op10xgdfix() above (and op10xfix()) for comments.
    ** This is a concise version.
    */
    register dw10_t d = *ad;
    register int i, rbit, sign;

    GF_DPOSSETUP(sign, i, d);	/* Set up exponent, get positive fract */
    if (i > (GFEEXCESS+35)) {		/* If integer wd be too big, */
	OP10_PCFSET(PCF_ARO+PCF_TR1);	/* set flags */
	return 0;			/* and return arg unchanged */
    }
    i -= (GFEEXCESS+35);	/* Get # bits to shift fraction */
    if (i <= -36) {		/* If shifting it to oblivion */
	ad->HI = w10zero;	/* just return zero, skip rigmarole */
	return 1;
    }

    if (rnd && i < 0) {		/* If rounding OK and losing bits */
	/* Find value of last bit that will be lost */
	if (i <= -35) {		/* Bit in high word? */
	    if (i <= -(35+H10BITS))	/* Bit comes from High LH? */
		rbit = LHGET(d.HI)
		    & (1 << ((-i)-(35+H10BITS+1)));
	    else		/* Bit comes from High RH */
		rbit = RHGET(d.HI)
		    & (1 << ((-i)-(35+1)));
	} else {		/* Bit in low word (shift -35 < i < 0) */
	    if (i <= -H10BITS)	/* Bit comes from Low LH? */
		rbit = LHGET(d.LO)
		    & (1 << ((-i)-(H10BITS+1)));
	    else		/* Bit comes from Low RH */
		rbit = RHGET(d.LO)
		    & (1 << ((-i)-1));
	}
    } else rbit = 0;

    if (i)
	d = x_ashc(d, i);	/* Do the shift */

    if (rbit) {			/* Need any rounding? */
	op10m_inc(d.HI);	/* Yep, trivial */
	if (!wmagskipn(d.HI)) {	/* Check for carry out */
	    if (!sign) {	/* Ugh.  OK if result will be negative, */
		OP10_PCFSET(PCF_ARO+PCF_TR1);	/* else set flags */
		return 0;			/* and return arg unchanged */
	    }
	}
    }

    /* Since returning a double fix, make sure result conforms to fixed-point
    ** sign convention (low sign is copy of high's)
    */
    if (sign)
	op10m_movn(d.HI);	/* Negate (sets low sign) */

    ad->HI = d.HI;		/* Return result */
    return 1;			/* Say succeeded */
}

#endif /* OP10_GFMT */

/* Old obsolete operations for PDP-6, KA-10, KI-10.
*/

/* DFN - Double Floating Negate
**	No flags are set.  "Negating a correctly formatted floating point
**	number cannot cause overflow".
*/
dw10_t op10dfn(register dw10_t d)
{
    register h10_t h;

    op10m_setcm(d.HI);			/* Complement high order word */
    h = (LHGET(d.LO) ^ HMASK) & (HMASK>>9);	/* Ditto low but clear sign&exp */
    RHSET(d.LO, (-RHGET(d.LO)) & HMASK);	/* Negate low RH */
    if (!RHGET(d.LO))
	/* If RH clear, fix up LH, but avoid sign & exp bits!! */
	if (!(h = (h+1) & (HMASK>>9))) {	/* If need more, */
	    op10m_inc(d.HI);			/* increment hi wd */
	}
    /* Combine low-order LH back in, must leave high 9 bits alone.
    ** The high 9 bits of H are already clear.
    */
    LHSET(d.LO, ((LHGET(d.LO) & ~(HMASK>>9)) | h));
    return d;
}

/* UFA - Unnormalized Floating Add
**	Similar enough to FAD that can use the same code!  ffadr()
**	may normalize 1 bit if an overflow happens, which is the same
**	thing that UFA does.
*/
w10_t op10ufa(register w10_t a, register w10_t b)
{
    return ffadr(a, b, FADF_NONORM|FADF_NORND);	/* No norm, no round */
}


/* Long form single-precision (immediate mode of FAD, FSB, FMP, FDV) */

static dw10_t dtokad(register dw10_t d)	/* (double) -> (KA double) */
{
    register int exp;

    op10m_signclr(d.LO);		/* Sign shd be 0, but make sure */
    RSHIFTM(d.LO, 8);			/* Make room for exponent */
    if (wskipn(d.LO)) {
	exp = SFEGET(d.HI);		/* Get sign and exponent */
	if (exp & SFESIGN) {		/* Negative? */
	    exp = exp ^ SFEMASK;	/* If so, get ones-complement */
	}
	if (exp < 27)			/* If exponent too small, */
	    op10m_setz(d.LO);		/* just clear 2nd word */
	else				/* Ah, put 2nd exp into low wd */
	    LHSET(d.LO, (LHGET(d.LO) | ((h10_t)(exp-27) << SFFBITS)));
    }
    return d;
}

dw10_t op10fadl(register w10_t a, register w10_t b)
{
    return dtokad(dfad(op10fdble(a), op10fdble(b), 0));
}

dw10_t op10fsbl(register w10_t a, register w10_t b)
{
    return dtokad(dfad(op10fdble(a), op10fdble(b), 1));
}

/* FMPL - This implementation differs when the exponent overflows.
**	Hardware clears the low word; this code still attempts to set
**	the exponent.
*/
dw10_t op10fmpl(register w10_t a, register w10_t b)
{
    return dtokad(op10dfmp(op10fdble(a), op10fdble(b)));
}

/* FDVL - I'm not positive this emulates the hardware precisely for
**	abnormal cases.  In particular, no-divide probably fails.
*/
dw10_t op10fdvl(register dw10_t d,
		register w10_t w)
{
    /* Convert from KA software double to hardware double */
    LSHIFTM(d.LO, 8);
    op10m_signclr(d.LO);		/* Ensure sign flushed */
    return dtokad(op10dfdv(d, op10fdble(w)));
}

/* OP10FDBLE - Converts Float to Double Float.
**	Auxiliary for Long-form instructions above, also a KCC C simop
**	in its own right.
*/
/* (double)(float)w */
dw10_t op10fdble(register w10_t w)
{
    dw10_t d;
    d.HI = w;
    op10m_setz(d.LO);
    return d;
}

/* Floating-point conversions to and from native platform format.
*/
#if 0	/* Not needed currently - eventually for KCC only? */

#include <math.h>		/* For native assistance */

/* OP10_DPUTF - Native double float to single-word PDP-10 float
*/
w10_t op10_dputf(dbl)
double dbl;
{
    register w10_t w;
    register uint32 imant;
    int exp, neg = 0;

    if (dbl <= 0) {	/* Get magnitude for easier hacking */
	if (dbl == 0)
	    return w10zero;
	neg++;
	dbl = -dbl;
	if (dbl < 0)
	    return w10maxneg;
    }
    /* Split value into normalized exp and fraction, then
    ** floating-shift fraction so we have 27 binary digits when the double
    ** is converted into an integer.
    */
    imant = (int32)ldexp(frexp(dbl, &exp), 27);	/* Get 27-bit integ fraction */

    /* Form high half from exponent and high fraction.  It is perfectly
    ** possible for the exponent to overflow -- tough.
    */
    LHSET(w, (((h10_t)((exp+0200)&0377)<<9) | ((imant>>18)&0777)));
    RHSET(w, (imant & HMASK));
    if (neg)
	op10m_movn(w);
    return w;
}

/* OP10_FGETD - Single-word PDP-10 float to native double float
*/
double op10_fgetd(w)
register w10_t w;
{
    register uint32 imant;
    register int exp;
    register double dbl;
    int neg = 0;

    if (wskipl(w)) {
	neg++;
	op10m_movn(w);
    }
    imant = op10wtos(w) & (((uint32)1<<27)-1);	/* Get 27-bit fraction */
    exp = (LHGET(w) >> 9) & 0377;		/* Get 8-bit exponent */
    if (!exp && !imant) return 0.0;
    dbl = ldexp((double)imant, (exp-0200)-27);	/* Put it together */
    return neg ? -dbl : dbl;
}


/* OP10_DPUTD - Native double float to double-word PDP-10 double
*/
static dw10_t op10_dputd(dbl)
double dbl;
{
#if 0
#endif
}

/* OP10_DGETD - Double-word PDP-10 doublefloat to native double float
*/
static double op10_dgetd(d)
register dw10_t d;
{
#if 0
    register uint32 imant;
    register int exp;
    register double dbl;
    int neg = 0;

    if (wskipl(w)) {
	neg++;
	op10m_movn(w);
    }
    imant = op10wtos(w) & (((uint32)1<<27)-1);	/* Get 27-bit fraction */
    exp = (LHGET(w) >> 9) & 0377;		/* Get 8-bit exponent */
    if (!exp && !imant) return 0.0;
    dbl = ldexp((double)imant, (exp-0200)-27);	/* Put it together */
    return neg ? -dbl : dbl;
#endif
}
#endif /* 0 */

/* Simulated Ops for KCC C operations (not actual PDP-10 ops) */

#if OP10_KCCOPS

/* Comparisons and testing */

/* C simop: Gives result < = > 0 for A-B */

int op10cmp(register w10_t a, register w10_t b)
{
    return op10m_caml(a,b) ? -1 : op10m_camn(a,b);
}

/* C simop: Gives result < = > 0 for unsigned A-B */

int op10ucmp(register w10_t a, register w10_t b)
{
    return op10m_ucmpl(a,b) ? -1 : op10m_camn(a,b);
}

/* C simop: Gives result < = > 0 for W */

int op10tst(register w10_t w)
{
    return op10m_skipl(w) ? -1 : op10m_skipn(w);
}

/* C simop: Gives result < = > 0 for DW A-B */
**	Double compare is a bit weird since sign of low word must be ignored.
*/
int op10dcmp(register dw10_t da, register dw10_t db)
{
    register int i;
    if (i = op10cmp(da.HI, db.HI))
	return i;
    op10m_signclr(da.LO);	/* Clear sign bits of low words */
    op10m_signclr(db.LO);
    return op10cmp(da.LO, db.LO);
}

/* C simop: Gives result < = > 0 for doubleword */

int op10dtst(register dw10_t d)
{
    return op10m_signtst(d.HI) ? -1 : dskipn(d);
}

/* C simop: Double fix negate */

dw10_t op10dneg(register dw10_t d)
{
    op10m_dneg(d);
    return d;
}


/* C simop: Unsigned Multiply */

w10_t op10uimul(w10_t a, w10_t b)
{
    dw10_t d;
    d = op10xmul(a, b);
    if (RHGET(d.HI) & 01)	/* If low bit of high wd is set, */
	op10m_signset(d.LO);	/* set high bit of low wd */
    return d.LO;
}

/* C simop: Unsigned Divide */

dw10_t op10uidiv(register w10_t a, register w10_t b)
{
    dw10_t db;
    qw10_t q;

    if (!wskipl(a) && !wskipl(b))
	return op10idiv(a, b);

    /* Punt for now by doing double fix division */
    q.D0 = dw10zero;
    q.D1.HI = wskipl(a) ? w10one : w10zero;
    q.D1.LO = a;
    db.HI = wskipl(b) ? w10one : w10zero;
    db.LO = b;
    q = op10ddiv(q, db);
    db.HI = q.D0.LO;
    db.LO = q.D1.LO;
    return db;
}

/* C simop: (float)(unsigned)w */

w10_t op10uflt(register w10_t w)
{
    if (wskipl(w))
	return sfnorm(0200+27+1, op10lsh(w, -1), 1);
     return sfnorm(0200+27, w, 1);	/* Normalize and round */
}

/* C simop: (int)(double)d */

w10_t op10dfix(register dw10_t d)
{
    register int exp, sign;

    SF_DPOSSETUP(sign, exp, d);	/* Extract parts and set up positive fract */
    d = x_ashc(d, exp - 0233);
    if (sign)
	op10m_dmovn(d);
    return d.HI;
}

/* C simop: (double)(int)w */

dw10_t op10dflt(register w10_t w)
{
    register dw10_t d;
    d.HI = w;
    op10m_setz(d.LO);
    d = x_ashc(d, -8);		/* Make room for exponent */
    LHSET(d.HI, (LHGET(d.HI) ^ 0243000));	/* Set proper exponent */
    return dfad(d, dw10zero, 0);	/* Normalize and return */
}

/* C simop: (double)(unsigned)w */

dw10_t op10udflt(register w10_t w)
{
    register dw10_t d;
    d.HI = w;
    op10m_setz(d.LO);
    d = op10lshc(d, -9);	/* Make room for exponent */
    RSHIFTM(d.LO, 1);		/* Get out of way of sign bit */
    LHSET(d.HI, (LHGET(d.HI) ^ 0244000));	/* Set proper exponent */
    return dfad(d, dw10zero, 0);	/* Normalize and return */
}

/* C simop: (float)(double)d */

w10_t op10dsngl(register dw10_t d)
{
    int sign;
    if (sign = wskipl(d.HI))
	op10m_dmovn(d);
    if (LHGET(d.LO) & (HSIGN>>1)) {
	if (RHGET(d.HI) & 01) {
	    w10_t one;
	    LRHSET(one, (LHGET(d.HI) & 0777000), 01);
	    d.HI = ffadr(d.HI, one, FADF_NORND);	/* No rounding */
	} else
	    op10m_iori(d.HI, 1);
    }
    if (sign)
	op10m_movn(d.HI);
    return d.HI;
}

#endif /* OP10_KCCOPS */
