code:
/*	$OpenBSD: impys.S,v 1.12 2011/04/16 20:52:12 deraadt Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)impys.s: Revision: 1.11.88.1 Date: 93/12/07 15:06:28 */

#include <machine/asm.h>
#include <machine/frame.h>

;****************************************************************************
;
;Implement an integer multiply routine for 32-bit operands and 64-bit product
;with operand values of zero (multiplicand only) and -2**31 treated specially.
;The algorithm uses the absolute value of the multiplier, four bits at a time,
;from right to left, to generate partial product.  Execution speed is more
;important than program size in this implementation.
;
;****************************************************************************
;
; Definitions - General registers
;
gr0	.reg		%r0		; General register zero
pu	.reg		%r3		; upper part of product
pl	.reg		%r4		; lower part of product
op2	.reg		%r4		; multiplier
op1	.reg		%r5		; multiplicand
cnt	.reg		%r6		; count in multiply
brindex	.reg		%r7		; index into the br. table
sign	.reg		%r8		; sign of product
pc	.reg		%r9		; carry bit of product, = 00...01
pm	.reg		%r10		; value of -1 used in shifting

;*****************************************************************************
	.text

LEAF_ENTRY(s_xmpy)
	stws,ma		pu,4(sp)		; save registers on stack
	stws,ma		pl,4(sp)		; save registers on stack
	stws,ma		op1,4(sp)		; save registers on stack
	stws,ma		cnt,4(sp)		; save registers on stack
	stws,ma		brindex,4(sp)		; save registers on stack
	stws,ma		sign,4(sp)		; save registers on stack
	stws,ma		pc,4(sp)		; save registers on stack
	stws,ma		pm,4(sp)		; save registers on stack
;
;   Start multiply process
;
	ldws		0(arg1),op2		; get multiplier
	ldws		0(arg0),op1		; get multiplicand
	addi		-1,gr0,pm		; initialize pm to 111...1
	comb,<		op2,gr0,mpyb		; br. if multiplier < 0
	xor		op2,op1,sign		; sign(0) = sign of product
mpy1	comb,<		op1,gr0,mpya		; br. if multiplicand < 0
	addi		0,gr0,pu		; clear product
	addib,=		0,op1,fini0		; op1 = 0, product = 0
mpy2	addi		1,gr0,pc		; initialize pc to 00...01
	movib,tr	8,cnt,mloop		; set count for mpy loop
	extru		op2,31,4,brindex	; 4 bits as index into table
;
	.align		8
;
	b		sh4c			; br. if sign overflow
sh4n	shd		pu,pl,4,pl		; shift product right 4 bits
	addib,<=	-1,cnt,mulend		; reduce count by 1, exit if
	extru		pu,27,28,pu		;   <= zero
;
mloop	blr		brindex,gr0		; br. into table
						;   entries of 2 words
	extru		op2,27,4,brindex	; next 4 bits into index
;
;
;	branch table for the multiplication process with four multiplier bits
;
mtable						; two words per entry
;
; ----	bits = 0000 ---- shift product 4 bits -------------------------------
;
	b		sh4n+4			; just shift partial
	shd		pu,pl,4,pl		;   product right 4 bits
;
;  ----	bits = 0001 ---- add op1, then shift 4 bits
;
	addb,tr		op1,pu,sh4n+4		; add op1 to product, to shift
	shd		pu,pl,4,pl		;   product right 4 bits
;
;  ----	bits = 0010 ---- add op1, add op1, then shift 4 bits
;
	addb,tr		op1,pu,sh4n		; add 2*op1, to shift
	addb,uv		op1,pu,sh4c		;   product right 4 bits
;
;  ---- bits = 0011 ---- add op1, add 2*op1, shift 4 bits
;
	addb,tr		op1,pu,sh4n-4		; add op1 & 2*op1, shift
	sh1add,nsv	op1,pu,pu			;   product right 4 bits
;
;  ----	bits = 0100 ---- shift 2, add op1, shift 2
;
	b		sh2sa
	shd		pu,pl,2,pl		; shift product 2 bits
;
;  ----	bits = 0101 ---- add op1, shift 2, add op1, and shift 2 again
;
	addb,tr		op1,pu,sh2us		; add op1 to product
	shd		pu,pl,2,pl		; shift 2 bits
;
;  ----	bits = 0110 ---- add op1, add op1, shift 2, add op1, and shift 2 again
;
	addb,tr		op1,pu,sh2c		; add 2*op1, to shift 2 bits
	addb,nuv	op1,pu,sh2us		; br. if not overflow
;
;  ----	bits = 0111 ---- subtract op1, shift 3, add op1, and shift 1
;
	b		sh3s
	sub		pu,op1,pu		; subtract op1, br. to sh3s

;
;  ----	bits = 1000 ---- shift 3, add op1, shift 1
;
	b		sh3sa
	shd		pu,pl,3,pl		; shift product right 3 bits
;
;  ----	bits = 1001 ---- add op1, shift 3, add op1, shift 1
;
	addb,tr		op1,pu,sh3us		; add op1, to shift 3, add op1,
	shd		pu,pl,3,pl		;   and shift 1
;
;  ----	bits = 1010 ---- add op1, add op1, shift 3, add op1, shift 1
;
	addb,tr		op1,pu,sh3c		; add 2*op1, to shift 3 bits
	addb,nuv	op1,pu,sh3us		;   br. if no overflow
;
;  ----	bits = 1011 ---- add -op1, shift 2, add -op1, shift 2, inc. next index
;
	addib,tr	1,brindex,sh2s		; add 1 to index, subtract op1,
	sub		pu,op1,pu		;   shift 2 with minus sign
;
;  ----	bits = 1100 ---- shift 2, subtract op1, shift 2, increment next index
;
	addib,tr	1,brindex,sh2sb		; add 1 to index, to shift
	shd		pu,pl,2,pl		; shift right 2 bits signed
;
;  ----	bits = 1101 ---- add op1, shift 2, add -op1, shift 2
;
	addb,tr		op1,pu,sh2ns		; add op1, to shift 2
	shd		pu,pl,2,pl		;   right 2 unsigned, etc.
;
;  ----	bits = 1110 ---- shift 1 signed, add -op1, shift 3 signed
;
	addib,tr	1,brindex,sh1sa		; add 1 to index, to shift
	shd		pu,pl,1,pl		; shift 1 bit
;
;  ----	bits = 1111 ---- add -op1, shift 4 signed
;
	addib,tr	1,brindex,sh4s		; add 1 to index, subtract op1,
	sub		pu,op1,pu		;   to shift 4 signed

;
;  ----	bits = 10000 ---- shift 4 signed
;
	addib,tr	1,brindex,sh4s+4	; add 1 to index
	shd		pu,pl,4,pl		; shift 4 signed
;
;  ---- end of table ---------------------------------------------------------
;
sh4s	shd		pu,pl,4,pl
	addib,tr	-1,cnt,mloop		; loop (count > 0 always here)
	shd		pm,pu,4,pu		; shift 4, minus signed
;
sh4c	addib,>		-1,cnt,mloop		; decrement count, loop if > 0
	shd		pc,pu,4,pu		; shift 4 with overflow
	b		signs			; end of multiply
	bb,>=,n		sign,0,fini		; test sign of procduct
;
mpyb	add,=		op2,op2,gr0		; if <> 0, back to main sect.
	b		mpy1
	sub		0,op2,op2		; op2 = |multiplier|
	add,>=		op1,gr0,gr0		; if op1 < 0, invert sign,
	xor		pm,sign,sign		;   for correct result
;
;	special case for multiplier = -2**31, op1 = signed multiplicand
;		or multiplicand = -2**31, op1 = signed multiplier
;
	shd		op1,0,1,pl		; shift op1 left 31 bits
mmax	extrs		op1,30,31,pu
	b		signs			; negate product (if needed)
	bb,>=,n		sign,0,fini		; test sign of product
;
mpya	add,=		op1,op1,gr0		; op1 = -2**31, special case
	b		mpy2
	sub		0,op1,op1		; op1 = |multiplicand|
	add,>=		op2,gr0,gr0		; if op2 < 0, invert sign,
	xor		pm,sign,sign		;   for correct result
	movb,tr		op2,op1,mmax		; use op2 as multiplicand
	shd		op1,0,1,pl		; shift it left 31 bits
;
sh3c	shd		pu,pl,3,pl		; shift product 3 bits
	shd		pc,pu,3,pu		; shift 3 signed
	addb,tr		op1,pu,sh1		; add op1, to shift 1 bit
	shd		pu,pl,1,pl
;
sh3us	extru		pu,28,29,pu		; shift 3 unsigned
	addb,tr		op1,pu,sh1		; add op1, to shift 1 bit
	shd		pu,pl,1,pl
;
sh3sa	extrs		pu,28,29,pu		; shift 3 signed
	addb,tr		op1,pu,sh1		; add op1, to shift 1 bit
	shd		pu,pl,1,pl
;
sh3s	shd		pu,pl,3,pl		; shift 3 minus signed
	shd		pm,pu,3,pu
	addb,tr		op1,pu,sh1		; add op1, to shift 1 bit
	shd		pu,pl,1,pl
;
sh1	addib,>		-1,cnt,mloop		; loop if count > 0
	extru		pu,30,31,pu
	b		signs			; end of multiply
	bb,>=,n		sign,0,fini		; test sign of product
;
sh2ns	addib,tr	1,brindex,sh2sb+4	; increment index
	extru		pu,29,30,pu		; shift unsigned
;
sh2s	shd		pu,pl,2,pl		; shift with minus sign
	shd		pm,pu,2,pu		;
	sub		pu,op1,pu		; subtract op1
	shd		pu,pl,2,pl		; shift with minus sign
	addib,tr	-1,cnt,mloop		; decrement count, loop
	shd		pm,pu,2,pu		; shift with minus sign
						; count never reaches 0 here
;
sh2sb	extrs		pu,29,30,pu		; shift 2 signed
	sub		pu,op1,pu		; subtract op1 from product
	shd		pu,pl,2,pl		; shift with minus sign
	addib,tr	-1,cnt,mloop		; decrement count, loop
	shd		pm,pu,2,pu		; shift with minus sign
						; count never reaches 0 here
;
sh1sa	extrs		pu,30,31,pu		;   signed
	sub		pu,op1,pu		; subtract op1 from product
	shd		pu,pl,3,pl		; shift 3 with minus sign
	addib,tr	-1,cnt,mloop		; dec. count, to loop
	shd		pm,pu,3,pu		; count never reaches 0 here
;
fini0	movib,tr,n	0,pl,fini		; product = 0 as op1 = 0
;
sh2us	extru		pu,29,30,pu		; shift 2 unsigned
	addb,tr		op1,pu,sh2a		; add op1
	shd		pu,pl,2,pl		; shift 2 bits
;
sh2c	shd		pu,pl,2,pl
	shd		pc,pu,2,pu		; shift with carry
	addb,tr		op1,pu,sh2a		; add op1 to product
	shd		pu,pl,2,pl		; br. to sh2 to shift pu
;
sh2sa	extrs		pu,29,30,pu		; shift with sign
	addb,tr		op1,pu,sh2a		; add op1 to product
	shd		pu,pl,2,pl		; br. to sh2 to shift pu
;
sh2a	addib,>		-1,cnt,mloop		; loop if count > 0
	extru		pu,29,30,pu
;
mulend	bb,>=,n		sign,0,fini		; test sign of product
signs	sub		0,pl,pl			; negate product if sign
	subb		0,pu,pu			;   is negative
;
;	finish
;
fini	stws		pu,0(arg2)		; save high part of result
	stws		pl,4(arg2)		; save low part of result

	ldws,mb		-4(sp),pm		; restore registers
	ldws,mb		-4(sp),pc		; restore registers
	ldws,mb		-4(sp),sign		; restore registers
	ldws,mb		-4(sp),brindex		; restore registers
	ldws,mb		-4(sp),cnt		; restore registers
	ldws,mb		-4(sp),op1		; restore registers
	ldws,mb		-4(sp),pl		; restore registers
	bv		0(rp)			; return
	ldws,mb		-4(sp),pu		; restore registers
EXIT(s_xmpy)

	.end


code:
/*	$OpenBSD: impyu.S,v 1.12 2011/04/16 20:52:12 deraadt Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)impyu.s: Revision: 1.11.88.1 Date: 93/12/07 15:06:31 */

#include <machine/asm.h>
#include <machine/frame.h>

;****************************************************************************
;
;Implement an integer multiply routine for 32-bit operands and 64-bit product
; with operand values of zero (multiplicand only) and 2**32reated specially.
; The algorithm uses the multiplier, four bits at a time, from right to left,
; to generate partial product.  Execution speed is more important than program
; size in this implementation.
;
;******************************************************************************
;
; Definitions - General registers
;
gr0	.reg		%r0		; General register zero
pu	.reg		%r3		; upper part of product
pl	.reg		%r4		; lower part of product
op2	.reg		%r4		; multiplier
op1	.reg		%r5		; multiplicand
cnt	.reg		%r6		; count in multiply
brindex	.reg		%r7		; index into the br. table
saveop2	.reg		%r8		; save op2 if high bit of multiplicand
					; is set
pc	.reg		%r9		; carry bit of product, = 00...01
pm	.reg		%r10		; value of -1 used in shifting
temp	.reg		%r6

;****************************************************************************
	.text
LEAF_ENTRY(u_xmpy)
	stws,ma		pu,4(sp)		; save registers on stack
	stws,ma		pl,4(sp)		; save registers on stack
	stws,ma		op1,4(sp)		; save registers on stack
	stws,ma		cnt,4(sp)		; save registers on stack
	stws,ma		brindex,4(sp)		; save registers on stack
	stws,ma		saveop2,4(sp)		; save registers on stack
	stws,ma		pc,4(sp)		; save registers on stack
	stws,ma		pm,4(sp)		; save registers on stack
;
;   Start multiply process
;
	ldws		0(arg0),op1		; get multiplicand
	ldws		0(arg1),op2		; get multiplier
	addib,=		0,op1,fini0		; op1 = 0, product = 0
	addi		0,gr0,pu		; clear product
	bb,>=		op1,0,mpy1		; test msb of multiplicand
	addi		0,gr0,saveop2		; clear saveop2
;
; msb of multiplicand is set so will save multiplier for a final
; addition into the result
;
	extru,=		op1,31,31,op1		; clear msb of multiplicand
	b		mpy1			; if op1 < 2**32, start multiply
	add		op2,gr0,saveop2		;   save op2 in saveop2
	shd		gr0,op2,1,pu		; shift op2 left 31 for result
	b		fini			; go to finish
	shd		op2,gr0,1,pl
;
mpy1	addi		-1,gr0,pm		; initialize pm to 111...1
	addi		1,gr0,pc		; initialize pc to 00...01
	movib,tr	8,cnt,mloop		; set count for mpy loop
	extru		op2,31,4,brindex	; 4 bits as index into table
;
	.align		8
;
	b		sh4c			; br. if sign overflow
sh4n	shd		pu,pl,4,pl		; shift product right 4 bits
	addib,<=	-1,cnt,mulend		; reduce count by 1, exit if
	extru		pu,27,28,pu		;   <= zero
;
mloop	blr		brindex,gr0		; br. into table
						;   entries of 2 words
	extru		op2,27,4,brindex	; next 4 bits into index
;
;
;	branch table for the multiplication process with four multiplier bits
;
mtable						; two words per entry
;
; ----	bits = 0000 ---- shift product 4 bits -------------------------------
;
	b		sh4n+4			; just shift partial
	shd		pu,pl,4,pl		;   product right 4 bits
;
;  ----	bits = 0001 ---- add op1, then shift 4 bits
;
	addb,tr		op1,pu,sh4n+4		; add op1 to product, to shift
	shd		pu,pl,4,pl		;   product right 4 bits
;
;  ----	bits = 0010 ---- add op1, add op1, then shift 4 bits
;
	addb,tr		op1,pu,sh4n		; add 2*op1, to shift
	addb,uv		op1,pu,sh4c		;   product right 4 bits
;
;  ---- bits = 0011 ---- add op1, add 2*op1, shift 4 bits
;
	addb,tr		op1,pu,sh4n-4		; add op1 & 2*op1, shift
	sh1add,nuv	op1,pu,pu			;   product right 4 bits
;
;  ----	bits = 0100 ---- shift 2, add op1, shift 2
;
	b		sh2sa
	shd		pu,pl,2,pl		; shift product 2 bits
;
;  ----	bits = 0101 ---- add op1, shift 2, add op1, and shift 2 again
;
	addb,tr		op1,pu,sh2us		; add op1 to product
	shd		pu,pl,2,pl		; shift 2 bits
;
;  ----	bits = 0110 ---- add op1, add op1, shift 2, add op1, and shift 2 again
;
	addb,tr		op1,pu,sh2c		; add 2*op1, to shift 2 bits
	addb,nuv	op1,pu,sh2us		; br. if not overflow
;
;  ----	bits = 0111 ---- subtract op1, shift 3, add op1, and shift 1
;
	b		sh3s
	sub		pu,op1,pu		; subtract op1, br. to sh3s

;
;  ----	bits = 1000 ---- shift 3, add op1, shift 1
;
	b		sh3sa
	shd		pu,pl,3,pl		; shift product right 3 bits
;
;  ----	bits = 1001 ---- add op1, shift 3, add op1, shift 1
;
	addb,tr		op1,pu,sh3us		; add op1, to shift 3, add op1,
	shd		pu,pl,3,pl		;   and shift 1
;
;  ----	bits = 1010 ---- add op1, add op1, shift 3, add op1, shift 1
;
	addb,tr		op1,pu,sh3c		; add 2*op1, to shift 3 bits
	addb,nuv	op1,pu,sh3us		;   br. if no overflow
;
;  ----	bits = 1011 ---- add -op1, shift 2, add -op1, shift 2, inc. next index
;
	addib,tr	1,brindex,sh2s		; add 1 to index, subtract op1,
	sub		pu,op1,pu		;   shift 2 with minus sign
;
;  ----	bits = 1100 ---- shift 2, subtract op1, shift 2, increment next index
;
	addib,tr	1,brindex,sh2sb		; add 1 to index, to shift
	shd		pu,pl,2,pl		; shift right 2 bits signed
;
;  ----	bits = 1101 ---- add op1, shift 2, add -op1, shift 2
;
	addb,tr		op1,pu,sh2ns		; add op1, to shift 2
	shd		pu,pl,2,pl		;   right 2 unsigned, etc.
;
;  ----	bits = 1110 ---- shift 1 signed, add -op1, shift 3 signed
;
	addib,tr	1,brindex,sh1sa		; add 1 to index, to shift
	shd		pu,pl,1,pl		; shift 1 bit
;
;  ----	bits = 1111 ---- add -op1, shift 4 signed
;
	addib,tr	1,brindex,sh4s		; add 1 to index, subtract op1,
	sub		pu,op1,pu		;   to shift 4 signed

;
;  ----	bits = 10000 ---- shift 4 signed
;
	addib,tr	1,brindex,sh4s+4		; add 1 to index
	shd		pu,pl,4,pl		; shift 4 signed
;
;  ---- end of table ---------------------------------------------------------
;
sh4s	shd		pu,pl,4,pl
	addib,>		-1,cnt,mloop		; decrement count, loop if > 0
	shd		pm,pu,4,pu		; shift 4, minus signed
	addb,tr		op1,pu,lastadd		; do one more add, then finish
	addb,=,n	saveop2,gr0,fini	; check saveop2
;
sh4c	addib,>		-1,cnt,mloop		; decrement count, loop if > 0
	shd		pc,pu,4,pu		; shift 4 with overflow
	b		lastadd			; end of multiply
	addb,=,n	saveop2,gr0,fini	; check saveop2
;
sh3c	shd		pu,pl,3,pl		; shift product 3 bits
	shd		pc,pu,3,pu		; shift 3 signed
	addb,tr		op1,pu,sh1		; add op1, to shift 1 bit
	shd		pu,pl,1,pl
;
sh3us	extru		pu,28,29,pu		; shift 3 unsigned
	addb,tr		op1,pu,sh1		; add op1, to shift 1 bit
	shd		pu,pl,1,pl
;
sh3sa	extrs		pu,28,29,pu		; shift 3 signed
	addb,tr		op1,pu,sh1		; add op1, to shift 1 bit
	shd		pu,pl,1,pl
;
sh3s	shd		pu,pl,3,pl		; shift 3 minus signed
	shd		pm,pu,3,pu
	addb,tr		op1,pu,sh1		; add op1, to shift 1 bit
	shd		pu,pl,1,pl
;
sh1	addib,>		-1,cnt,mloop		; loop if count > 0
	extru		pu,30,31,pu
	b		lastadd			; end of multiply
	addb,=,n	saveop2,gr0,fini	; check saveop2
;
sh2ns	addib,tr	1,brindex,sh2sb+4	; increment index
	extru		pu,29,30,pu		; shift unsigned
;
sh2s	shd		pu,pl,2,pl		; shift with minus sign
	shd		pm,pu,2,pu		;
	sub		pu,op1,pu		; subtract op1
	shd		pu,pl,2,pl		; shift with minus sign
	addib,>		-1,cnt,mloop		; decrement count, loop if > 0
	shd		pm,pu,2,pu		; shift with minus sign
	addb,tr		op1,pu,lastadd		; do one more add, then finish
	addb,=,n	saveop2,gr0,fini	; check saveop2
;
sh2sb	extrs		pu,29,30,pu		; shift 2 signed
	sub		pu,op1,pu		; subtract op1 from product
	shd		pu,pl,2,pl		; shift with minus sign
	addib,>		-1,cnt,mloop		; decrement count, loop if > 0
	shd		pm,pu,2,pu		; shift with minus sign
	addb,tr		op1,pu,lastadd		; do one more add, then finish
	addb,=,n	saveop2,gr0,fini	; check saveop2
;
sh1sa	extrs		pu,30,31,pu		;   signed
	sub		pu,op1,pu		; subtract op1 from product
	shd		pu,pl,3,pl		; shift 3 with minus sign
	addib,>		-1,cnt,mloop		; decrement count, loop if >0
	shd		pm,pu,3,pu
	addb,tr		op1,pu,lastadd		; do one more add, then finish
	addb,=,n	saveop2,gr0,fini	; check saveop2
;
fini0	movib,tr	0,pl,fini		; product = 0 as op1 = 0
	stws		pu,0(arg2)		; save high part of result
;
sh2us	extru		pu,29,30,pu		; shift 2 unsigned
	addb,tr		op1,pu,sh2a		; add op1
	shd		pu,pl,2,pl		; shift 2 bits
;
sh2c	shd		pu,pl,2,pl
	shd		pc,pu,2,pu		; shift with carry
	addb,tr		op1,pu,sh2a		; add op1 to product
	shd		pu,pl,2,pl		; br. to sh2 to shift pu
;
sh2sa	extrs		pu,29,30,pu		; shift with sign
	addb,tr		op1,pu,sh2a		; add op1 to product
	shd		pu,pl,2,pl		; br. to sh2 to shift pu
;
sh2a	addib,>		-1,cnt,mloop		; loop if count > 0
	extru		pu,29,30,pu
;
mulend	addb,=,n	saveop2,gr0,fini	; check saveop2
lastadd	shd		saveop2,gr0,1,temp	;  if saveop2 <> 0, shift it
	shd		gr0,saveop2,1,saveop2	;  left 31 and add to result
	add		pl,temp,pl
	addc		pu,saveop2,pu
;
;	finish
;
fini	stws		pu,0(arg2)		; save high part of result
	stws		pl,4(arg2)		; save low part of result

	ldws,mb		-4(sp),pm		; restore registers
	ldws,mb		-4(sp),pc		; restore registers
	ldws,mb		-4(sp),saveop2		; restore registers
	ldws,mb		-4(sp),brindex		; restore registers
	ldws,mb		-4(sp),cnt		; restore registers
	ldws,mb		-4(sp),op1		; restore registers
	ldws,mb		-4(sp),pl		; restore registers
	bv		0(rp)			; return
	ldws,mb		-4(sp),pu		; restore registers
EXIT(u_xmpy)

	.end

license: unknown


code:
/*	$OpenBSD: md.h,v 1.7 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)md.h: Revision: 1.8.88.1 Date: 93/12/07 15:06:33 */

#include <sys/cdefs.h>

/*****************************************************************
 * Muliply/Divide SFU Internal State                             *
 *****************************************************************/
struct mdsfu_register {
    int rslt_hi,
	rslt_lo,
	ovflow;
    };

#define result_hi result->rslt_hi
#define result_lo result->rslt_lo
#define overflow result->ovflow

/*
 *  Constants
 */

#undef FALSE
#define FALSE 0

#undef TRUE
#define TRUE (!FALSE)

#undef NIL
#define NIL 0

#define WORD_LEN 32
#define BIT0 1<<31
#define BIT28 0x8
#define BIT29 0x4
#define BIT30 0x2
#define BIT31 0x1

	/* Simply copy the arguments to the emulated copies of the registers */
#define	mdrr(reg1,reg2,result)	{result_hi = reg1;result_lo = reg2;}

/*
 *  Structures
 */

struct md_state {
	int resulthi,	/* high word of result */
	    resultlo;	/* low word of result */
};

void divsfm(int, int, struct mdsfu_register *);
void divsfr(int, int, struct mdsfu_register *);
void divsim(int, int, struct mdsfu_register *);
void divsir(int, int, struct mdsfu_register *);

void divu(int, int, int, struct mdsfu_register *);
void divufr(unsigned int, unsigned int, struct mdsfu_register *);
void divuir(unsigned int, unsigned int, struct mdsfu_register *);

void mpyaccs(int, int, struct mdsfu_register *);
void mpyaccu(unsigned int, unsigned int, struct mdsfu_register *);
void mpys(int, int, struct mdsfu_register *);
void mpyscv(int, int, struct mdsfu_register *);
void mpyu(unsigned int, unsigned int, struct mdsfu_register *);
void mpyucv(unsigned int, unsigned int, struct mdsfu_register *);

int s_xmpy(int *, int *, struct mdsfu_register *);
int u_xmpy(int *, int *, struct mdsfu_register *);


code:
/*	$OpenBSD: mpyaccs.c,v 1.6 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)mpyaccs.c: Revision: 1.6.88.1 Date: 93/12/07 15:06:39 */

#include "md.h"

void
mpyaccs(opnd1,opnd2,result)
	int opnd1, opnd2;
	struct mdsfu_register *result;
{
	struct mdsfu_register temp;
	int carry, sign;

	s_xmpy(&opnd1,&opnd2,&temp);

	/* get result of low word add, and check for carry out */
	if ((result_lo += (unsigned)temp.rslt_lo) < (unsigned)temp.rslt_lo)
		carry = 1;
	else
		carry = 0;

	/* get result of high word add, and determine overflow status */
	sign = result_hi ^ temp.rslt_hi;
	result_hi += temp.rslt_hi + carry;
	if (sign >= 0 && (temp.rslt_hi ^ result_hi) < 0)
		overflow = TRUE;
}


code:
/*	$OpenBSD: mpyaccu.c,v 1.7 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)mpyaccu.c: Revision: 1.6.88.1 Date: 93/12/07 15:06:41 */

#include "md.h"

void
mpyaccu(opnd1,opnd2,result)
	unsigned int opnd1, opnd2;
	struct mdsfu_register *result;
{
	struct mdsfu_register temp;
	int carry;

	u_xmpy(&opnd1,&opnd2,&temp);

	/* get result of low word add, and check for carry out */
	if ((result_lo += (unsigned)temp.rslt_lo) < (unsigned)temp.rslt_lo)
		carry = 1;
	else
		carry = 0;

	/* get result of high word add, and determine overflow status */
	if ((result_hi += (unsigned)temp.rslt_hi + carry) <
	    (unsigned)temp.rslt_hi)
		overflow = TRUE;
}


code:
/*	$OpenBSD: mpys.c,v 1.7 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)mpys.c: Revision: 1.6.88.1 Date: 93/12/07 15:06:43 */

#include "md.h"

void
mpys(opnd1,opnd2,result)
	int opnd1, opnd2;
	struct mdsfu_register *result;
{
	s_xmpy(&opnd1,&opnd2,result);

	/* determine overflow status */
	if ((result_hi == 0 && result_lo >= 0) ||
	    (result_hi == -1 && result_lo < 0))
		overflow = FALSE;
	else
		overflow = TRUE;
}


code:
/*	$OpenBSD: mpyscv.c,v 1.7 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)mpyscv.c: Revision: 1.6.88.1 Date: 93/12/07 15:06:45 */

#include "md.h"

void
mpyscv(opnd1,opnd2,result)
	int opnd1, opnd2;
	struct mdsfu_register *result;
{
	s_xmpy(&opnd1,&opnd2,result);
	overflow = FALSE;
}


code:
/*	$OpenBSD: mpyu.c,v 1.6 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)mpyu.c: Revision: 1.6.88.1 Date: 93/12/07 15:06:48 */

#include "md.h"

void
mpyu(opnd1,opnd2,result)
	unsigned int opnd1, opnd2;
	struct mdsfu_register *result;
{
	u_xmpy(&opnd1,&opnd2,result);

	/* determine overflow status */
	if (result_hi)
		overflow = TRUE;
	else
		overflow = FALSE;
	return;
}


code:
/*	$OpenBSD: mpyucv.c,v 1.6 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)mpyucv.c: Revision: 1.6.88.1 Date: 93/12/07 15:06:50 */

#include "md.h"

void
mpyucv(opnd1,opnd2,result)
	unsigned int opnd1, opnd2;
	struct mdsfu_register *result;
{
	u_xmpy(&opnd1,&opnd2,result);
	overflow = FALSE;
}


code:
/*	$OpenBSD: quad_float.h,v 1.8 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)quad_float.h: Revision: 1.7.88.1 Date: 93/12/07 15:06:55 */

/******************************
 *  Quad precision functions  *
 ******************************/

/* 32-bit word grabing functions */
#define Quad_firstword(value) Qallp1(value)
#define Quad_secondword(value) Qallp2(value)
#define Quad_thirdword(value)  Qallp3(value)
#define Quad_fourthword(value)  Qallp4(value)

/* This magnitude comparison uses the signless first words and
 * the regular part2 words.  The comparison is graphically:
 *
 *       1st greater?  ----------->|
 *                                 |
 *       1st less?-----------------+------->|
 *                                 |        |
 *       2nd greater?------------->|        |
 *                                 |        |
 *       2nd less?-----------------+------->|
 *                                 |        |
 *       3rd greater?------------->|        |
 *                                 |        |
 *       3rd less?-----------------+------->|
 *                                 |        |
 *       4th greater or equal?---->|        |
 *                                 |        |
 *                               False     True
 */
#define Quad_ismagnitudeless(leftp3,leftp4,rightp1,rightp2,rightp3,rightp4,signlessleft,signlessright) \
/*  Quad_floating_point left, right;          *				\
 *  unsigned int signlessleft, signlessright; */			\
      ( signlessleft<=signlessright &&					\
       (signlessleft<signlessright || (Qallp2(leftp2)<=Qallp2(rightp2) && \
	(Qallp2(leftp2)<Qallp2(rightp2) || (Qallp3(leftp3)<=Qallp3(rightp3) && \
	 (Qallp3(leftp3)<Qallp3(rightp3) || Qallp4(leftp4)<Qallp4(rightp4)))))))

#define Quad_xor_to_intp1(leftp1,rightp1,result)		\
    /* quad_floating_point left, right;				\
     * unsigned int result; */					\
    result = Qallp1(leftp1) XOR Qallp1(rightp1);

#define Quad_xor_from_intp1(leftp1,rightp1,result)		\
    /* quad_floating_point right, result;			\
     * unsigned int left; */					\
    Qallp1(resultp1) = left XOR Qallp1(rightp1)

#define Quad_swap_lower(leftp1,leftp2,leftp3,leftp4,rightp1,rightp2,rightp3,rightp4)  \
    /* quad_floating_point left, right; */			\
    Qallp2(leftp2)  = Qallp2(leftp2) XOR Qallp2(rightp2)	\
    Qallp2(rightp2) = Qallp2(leftp2) XOR Qallp2(rightp2)	\
    Qallp2(leftp2)  = Qallp2(leftp2) XOR Qallp2(rightp2)	\
    Qallp3(leftp3)  = Qallp3(leftp3) XOR Qallp3(rightp3)	\
    Qallp3(rightp3) = Qallp3(leftp3) XOR Qallp3(rightp3)	\
    Qallp3(leftp3)  = Qallp3(leftp3) XOR Qallp3(rightp3)	\
    Qallp4(leftp4)  = Qallp4(leftp4) XOR Qallp4(rightp4)	\
    Qallp4(rightp4) = Qallp4(leftp4) XOR Qallp4(rightp4)	\
    Qallp4(leftp4)  = Qallp4(leftp4) XOR Qallp4(rightp4)

int quad_frnd(quad_floating_point *, quad_floating_point *, quad_floating_point *, unsigned int *);


code:
/*	$OpenBSD: sfadd.c,v 1.6 2002/11/29 09:27:34 deraadt Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)sfadd.c: Revision: 2.7.88.1 Date: 93/12/07 15:07:00 */

#include "float.h"
#include "sgl_float.h"

/*
 * Single_add: add two single precision values.
 */
int
sgl_fadd(leftptr, rightptr, dstptr, status)
    sgl_floating_point *leftptr, *rightptr, *dstptr;
    unsigned int *status;
{
    register unsigned int left, right, result, extent;
    register unsigned int signless_upper_left, signless_upper_right, save;

    register int result_exponent, right_exponent, diff_exponent;
    register int sign_save, jumpsize;
    register int inexact = FALSE;
    register int underflowtrap;

    /* Create local copies of the numbers */
    left = *leftptr;
    right = *rightptr;

    /* A zero "save" helps discover equal operands (for later),	*
     * and is used in swapping operands (if needed).		*/
    Sgl_xortointp1(left,right,/*to*/save);

    /*
     * check first operand for NaN's or infinity
     */
    if ((result_exponent = Sgl_exponent(left)) == SGL_INFINITY_EXPONENT)
	{
	if (Sgl_iszero_mantissa(left))
	    {
	    if (Sgl_isnotnan(right))
		{
		if (Sgl_isinfinity(right) && save!=0)
		    {
		    /*
		     * invalid since operands are opposite signed infinity's
		     */
		    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		    Set_invalidflag();
		    Sgl_makequietnan(result);
		    *dstptr = result;
		    return(NOEXCEPTION);
		    }
		/*
		 * return infinity
		 */
		*dstptr = left;
		return(NOEXCEPTION);
		}
	    }
	else
	    {
	    /*
	     * is NaN; signaling or quiet?
	     */
	    if (Sgl_isone_signaling(left))
		{
		/* trap if INVALIDTRAP enabled */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		/* make NaN quiet */
		Set_invalidflag();
		Sgl_set_quiet(left);
		}
	    /*
	     * is second operand a signaling NaN?
	     */
	    else if (Sgl_is_signalingnan(right))
		{
		/* trap if INVALIDTRAP enabled */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		/* make NaN quiet */
		Set_invalidflag();
		Sgl_set_quiet(right);
		*dstptr = right;
		return(NOEXCEPTION);
		}
	    /*
	     * return quiet NaN
	     */
	    *dstptr = left;
	    return(NOEXCEPTION);
	    }
	} /* End left NaN or Infinity processing */
    /*
     * check second operand for NaN's or infinity
     */
    if (Sgl_isinfinity_exponent(right))
	{
	if (Sgl_iszero_mantissa(right))
	    {
	    /* return infinity */
	    *dstptr = right;
	    return(NOEXCEPTION);
	    }
	/*
	 * is NaN; signaling or quiet?
	 */
	if (Sgl_isone_signaling(right))
	    {
	    /* trap if INVALIDTRAP enabled */
	    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
	    /* make NaN quiet */
	    Set_invalidflag();
	    Sgl_set_quiet(right);
	    }
	/*
	 * return quiet NaN
	 */
	*dstptr = right;
	return(NOEXCEPTION);
	} /* End right NaN or Infinity processing */

    /* Invariant: Must be dealing with finite numbers */

    /* Compare operands by removing the sign */
    Sgl_copytoint_exponentmantissa(left,signless_upper_left);
    Sgl_copytoint_exponentmantissa(right,signless_upper_right);

    /* sign difference selects add or sub operation. */
    if(Sgl_ismagnitudeless(signless_upper_left,signless_upper_right))
	{
	/* Set the left operand to the larger one by XOR swap *
	 *  First finish the first word using "save"	  */
	Sgl_xorfromintp1(save,right,/*to*/right);
	Sgl_xorfromintp1(save,left,/*to*/left);
	result_exponent = Sgl_exponent(left);
	}
    /* Invariant:  left is not smaller than right. */

    if((right_exponent = Sgl_exponent(right)) == 0)
	{
	/* Denormalized operands.  First look for zeroes */
	if(Sgl_iszero_mantissa(right))
	    {
	    /* right is zero */
	    if(Sgl_iszero_exponentmantissa(left))
		{
		/* Both operands are zeros */
		if(Is_rounding_mode(ROUNDMINUS))
		    {
		    Sgl_or_signs(left,/*with*/right);
		    }
		else
		    {
		    Sgl_and_signs(left,/*with*/right);
		    }
		}
	    else
		{
		/* Left is not a zero and must be the result.  Trapped
		 * underflows are signaled if left is denormalized.  Result
		 * is always exact. */
		if( (result_exponent == 0) && Is_underflowtrap_enabled() )
		    {
		    /* need to normalize results mantissa */
		    sign_save = Sgl_signextendedsign(left);
		    Sgl_leftshiftby1(left);
		    Sgl_normalize(left,result_exponent);
		    Sgl_set_sign(left,/*using*/sign_save);
		    Sgl_setwrapped_exponent(left,result_exponent,unfl);
		    *dstptr = left;
		    return(UNDERFLOWEXCEPTION);
		    }
		}
	    *dstptr = left;
	    return(NOEXCEPTION);
	    }

	/* Neither are zeroes */
	Sgl_clear_sign(right);	/* Exponent is already cleared */
	if(result_exponent == 0 )
	    {
	    /* Both operands are denormalized.  The result must be exact
	     * and is simply calculated.  A sum could become normalized and a
	     * difference could cancel to a true zero. */
	    if( (/*signed*/int) save < 0 )
		{
		Sgl_subtract(left,/*minus*/right,/*into*/result);
		if(Sgl_iszero_mantissa(result))
		    {
		    if(Is_rounding_mode(ROUNDMINUS))
			{
			Sgl_setone_sign(result);
			}
		    else
			{
			Sgl_setzero_sign(result);
			}
		    *dstptr = result;
		    return(NOEXCEPTION);
		    }
		}
	    else
		{
		Sgl_addition(left,right,/*into*/result);
		if(Sgl_isone_hidden(result))
		    {
		    *dstptr = result;
		    return(NOEXCEPTION);
		    }
		}
	    if(Is_underflowtrap_enabled())
		{
		/* need to normalize result */
		sign_save = Sgl_signextendedsign(result);
		Sgl_leftshiftby1(result);
		Sgl_normalize(result,result_exponent);
		Sgl_set_sign(result,/*using*/sign_save);
		Sgl_setwrapped_exponent(result,result_exponent,unfl);
		*dstptr = result;
		return(UNDERFLOWEXCEPTION);
		}
	    *dstptr = result;
	    return(NOEXCEPTION);
	    }
	right_exponent = 1;	/* Set exponent to reflect different bias
				 * with denomalized numbers. */
	}
    else
	{
	Sgl_clear_signexponent_set_hidden(right);
	}
    Sgl_clear_exponent_set_hidden(left);
    diff_exponent = result_exponent - right_exponent;

    /*
     * Special case alignment of operands that would force alignment
     * beyond the extent of the extension.  A further optimization
     * could special case this but only reduces the path length for this
     * infrequent case.
     */
    if(diff_exponent > SGL_THRESHOLD)
	{
	diff_exponent = SGL_THRESHOLD;
	}

    /* Align right operand by shifting to right */
    Sgl_right_align(/*operand*/right,/*shifted by*/diff_exponent,
     /*and lower to*/extent);

    /* Treat sum and difference of the operands separately. */
    if( (/*signed*/int) save < 0 )
	{
	/*
	 * Difference of the two operands.  Their can be no overflow.  A
	 * borrow can occur out of the hidden bit and force a post
	 * normalization phase.
	 */
	Sgl_subtract_withextension(left,/*minus*/right,/*with*/extent,/*into*/result);
	if(Sgl_iszero_hidden(result))
	    {
	    /* Handle normalization */
	    /* A straight forward algorithm would now shift the result
	     * and extension left until the hidden bit becomes one.  Not
	     * all of the extension bits need participate in the shift.
	     * Only the two most significant bits (round and guard) are
	     * needed.  If only a single shift is needed then the guard
	     * bit becomes a significant low order bit and the extension
	     * must participate in the rounding.  If more than a single
	     * shift is needed, then all bits to the right of the guard
	     * bit are zeros, and the guard bit may or may not be zero. */
	    sign_save = Sgl_signextendedsign(result);
	    Sgl_leftshiftby1_withextent(result,extent,result);

	    /* Need to check for a zero result.  The sign and exponent
	     * fields have already been zeroed.  The more efficient test
	     * of the full object can be used.
	     */
	    if(Sgl_iszero(result))
		/* Must have been "x-x" or "x+(-x)". */
		{
		if(Is_rounding_mode(ROUNDMINUS)) Sgl_setone_sign(result);
		*dstptr = result;
		return(NOEXCEPTION);
		}
	    result_exponent--;
	    /* Look to see if normalization is finished. */
	    if(Sgl_isone_hidden(result))
		{
		if(result_exponent==0)
		    {
		    /* Denormalized, exponent should be zero.  Left operand *
		     * was normalized, so extent (guard, round) was zero    */
		    goto underflow;
		    }
		else
		    {
		    /* No further normalization is needed. */
		    Sgl_set_sign(result,/*using*/sign_save);
		    Ext_leftshiftby1(extent);
		    goto round;
		    }
		}

	    /* Check for denormalized, exponent should be zero.  Left    *
	     * operand was normalized, so extent (guard, round) was zero */
	    if(!(underflowtrap = Is_underflowtrap_enabled()) &&
	       result_exponent==0) goto underflow;

	    /* Shift extension to complete one bit of normalization and
	     * update exponent. */
	    Ext_leftshiftby1(extent);

	    /* Discover first one bit to determine shift amount.  Use a
	     * modified binary search.  We have already shifted the result
	     * one position right and still not found a one so the remainder
	     * of the extension must be zero and simplifies rounding. */
	    /* Scan bytes */
	    while(Sgl_iszero_hiddenhigh7mantissa(result))
		{
		Sgl_leftshiftby8(result);
		if((result_exponent -= 8) <= 0  && !underflowtrap)
		    goto underflow;
		}
	    /* Now narrow it down to the nibble */
	    if(Sgl_iszero_hiddenhigh3mantissa(result))
		{
		/* The lower nibble contains the normalizing one */
		Sgl_leftshiftby4(result);
		if((result_exponent -= 4) <= 0 && !underflowtrap)
		    goto underflow;
		}
	    /* Select case were first bit is set (already normalized)
	     * otherwise select the proper shift. */
	    if((jumpsize = Sgl_hiddenhigh3mantissa(result)) > 7)
		{
		/* Already normalized */
		if(result_exponent <= 0) goto underflow;
		Sgl_set_sign(result,/*using*/sign_save);
		Sgl_set_exponent(result,/*using*/result_exponent);
		*dstptr = result;
		return(NOEXCEPTION);
		}
	    Sgl_sethigh4bits(result,/*using*/sign_save);
	    switch(jumpsize)
		{
		case 1:
		    {
		    Sgl_leftshiftby3(result);
		    result_exponent -= 3;
		    break;
		    }
		case 2:
		case 3:
		    {
		    Sgl_leftshiftby2(result);
		    result_exponent -= 2;
		    break;
		    }
		case 4:
		case 5:
		case 6:
		case 7:
		    {
		    Sgl_leftshiftby1(result);
		    result_exponent -= 1;
		    break;
		    }
		}
	    if(result_exponent > 0)
		{
		Sgl_set_exponent(result,/*using*/result_exponent);
		*dstptr = result;
		return(NOEXCEPTION); /* Sign bit is already set */
		}
	    /* Fixup potential underflows */
	  underflow:
	    if(Is_underflowtrap_enabled())
		{
		Sgl_set_sign(result,sign_save);
		Sgl_setwrapped_exponent(result,result_exponent,unfl);
		*dstptr = result;
		/* inexact = FALSE; */
		return(UNDERFLOWEXCEPTION);
		}
	    /*
	     * Since we cannot get an inexact denormalized result,
	     * we can now return.
	     */
	    Sgl_right_align(result,/*by*/(1-result_exponent),extent);
	    Sgl_clear_signexponent(result);
	    Sgl_set_sign(result,sign_save);
	    *dstptr = result;
	    return(NOEXCEPTION);
	    } /* end if(hidden...)... */
	/* Fall through and round */
	} /* end if(save < 0)... */
    else
	{
	/* Add magnitudes */
	Sgl_addition(left,right,/*to*/result);
	if(Sgl_isone_hiddenoverflow(result))
	    {
	    /* Prenormalization required. */
	    Sgl_rightshiftby1_withextent(result,extent,extent);
	    Sgl_arithrightshiftby1(result);
	    result_exponent++;
	    } /* end if hiddenoverflow... */
	} /* end else ...add magnitudes... */

    /* Round the result.  If the extension is all zeros,then the result is
     * exact.  Otherwise round in the correct direction.  No underflow is
     * possible. If a postnormalization is necessary, then the mantissa is
     * all zeros so no shift is needed. */
  round:
    if(Ext_isnotzero(extent))
	{
	inexact = TRUE;
	switch(Rounding_mode())
	    {
	    case ROUNDNEAREST: /* The default. */
	    if(Ext_isone_sign(extent))
		{
		/* at least 1/2 ulp */
		if(Ext_isnotzero_lower(extent)  ||
		  Sgl_isone_lowmantissa(result))
		    {
		    /* either exactly half way and odd or more than 1/2ulp */
		    Sgl_increment(result);
		    }
		}
	    break;

	    case ROUNDPLUS:
	    if(Sgl_iszero_sign(result))
		{
		/* Round up positive results */
		Sgl_increment(result);
		}
	    break;

	    case ROUNDMINUS:
	    if(Sgl_isone_sign(result))
		{
		/* Round down negative results */
		Sgl_increment(result);
		}

	    case ROUNDZERO:;
	    /* truncate is simple */
	    } /* end switch... */
	if(Sgl_isone_hiddenoverflow(result)) result_exponent++;
	}
    if(result_exponent == SGL_INFINITY_EXPONENT)
	{
	/* Overflow */
	if(Is_overflowtrap_enabled())
	    {
	    Sgl_setwrapped_exponent(result,result_exponent,ovfl);
	    *dstptr = result;
	    if (inexact) {
		if (Is_inexacttrap_enabled())
		    return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
		else Set_inexactflag();
	    }
	    return(OVERFLOWEXCEPTION);
	    }
	else
	    {
	    Set_overflowflag();
	    inexact = TRUE;
	    Sgl_setoverflow(result);
	    }
	}
    else Sgl_set_exponent(result,result_exponent);
    *dstptr = result;
    if(inexact) {
	if(Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
	else Set_inexactflag();
    }
    return(NOEXCEPTION);
}


code:
/*	$OpenBSD: sfcmp.c,v 1.6 2002/09/20 19:26:59 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)sfcmp.c: Revision: 1.7.88.2 Date: 93/12/08 13:28:43 */

#include "float.h"
#include "sgl_float.h"

/*
 * sgl_cmp: compare two values
 */
int
sgl_fcmp(leftptr, rightptr, cond, status)
    sgl_floating_point *leftptr, *rightptr;
    unsigned int cond; /* The predicate to be tested */
    unsigned int *status;
{
    register unsigned int left, right;
    register int xorresult;

    /* Create local copies of the numbers */
    left = *leftptr;
    right = *rightptr;
    /*
     * Test for NaN
     */
    if(    (Sgl_exponent(left) == SGL_INFINITY_EXPONENT)
	|| (Sgl_exponent(right) == SGL_INFINITY_EXPONENT) )
	{
	/* Check if a NaN is involved.  Signal an invalid exception when
	 * comparing a signaling NaN or when comparing quiet NaNs and the
	 * low bit of the condition is set */
	if( (  (Sgl_exponent(left) == SGL_INFINITY_EXPONENT)
	    && Sgl_isnotzero_mantissa(left)
	    && (Exception(cond) || Sgl_isone_signaling(left)))
	   ||
	    (  (Sgl_exponent(right) == SGL_INFINITY_EXPONENT)
	    && Sgl_isnotzero_mantissa(right)
	    && (Exception(cond) || Sgl_isone_signaling(right)) ) )
	    {
	    if( Is_invalidtrap_enabled() ) {
		Set_status_cbit(Unordered(cond));
		return(INVALIDEXCEPTION);
	    }
	    else Set_invalidflag();
	    Set_status_cbit(Unordered(cond));
	    return(NOEXCEPTION);
	    }
	/* All the exceptional conditions are handled, now special case
	   NaN compares */
	else if( ((Sgl_exponent(left) == SGL_INFINITY_EXPONENT)
	    && Sgl_isnotzero_mantissa(left))
	   ||
	    ((Sgl_exponent(right) == SGL_INFINITY_EXPONENT)
	    && Sgl_isnotzero_mantissa(right)) )
	    {
	    /* NaNs always compare unordered. */
	    Set_status_cbit(Unordered(cond));
	    return(NOEXCEPTION);
	    }
	/* infinities will drop down to the normal compare mechanisms */
	}
    /* First compare for unequal signs => less or greater or
     * special equal case */
    Sgl_xortointp1(left,right,xorresult);
    if( xorresult < 0 )
	{
	/* left negative => less, left positive => greater.
	 * equal is possible if both operands are zeros. */
	if( Sgl_iszero_exponentmantissa(left)
	  && Sgl_iszero_exponentmantissa(right) )
	    {
	    Set_status_cbit(Equal(cond));
	    }
	else if( Sgl_isone_sign(left) )
	    {
	    Set_status_cbit(Lessthan(cond));
	    }
	else
	    {
	    Set_status_cbit(Greaterthan(cond));
	    }
	}
    /* Signs are the same.  Treat negative numbers separately
     * from the positives because of the reversed sense.  */
    else if( Sgl_all(left) == Sgl_all(right) )
	{
	Set_status_cbit(Equal(cond));
	}
    else if( Sgl_iszero_sign(left) )
	{
	/* Positive compare */
	if( Sgl_all(left) < Sgl_all(right) )
	    {
	    Set_status_cbit(Lessthan(cond));
	    }
	else
	    {
	    Set_status_cbit(Greaterthan(cond));
	    }
	}
    else
	{
	/* Negative compare.  Signed or unsigned compares
	 * both work the same.  That distinction is only
	 * important when the sign bits differ. */
	if( Sgl_all(left) > Sgl_all(right) )
	    {
	    Set_status_cbit(Lessthan(cond));
	    }
	else
	    {
	    Set_status_cbit(Greaterthan(cond));
	    }
	}
	return(NOEXCEPTION);
}


code:
/*	$OpenBSD: sfdiv.c,v 1.5 2002/05/07 22:19:30 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)sfdiv.c: Revision: 2.9.88.1 Date: 93/12/07 15:07:05 */

#include "float.h"
#include "sgl_float.h"

/*
 *  Single Precision Floating-point Divide
 */
int
sgl_fdiv(srcptr1,srcptr2,dstptr,status)
	sgl_floating_point *srcptr1, *srcptr2, *dstptr;
	unsigned int *status;
{
	register unsigned int opnd1, opnd2, opnd3, result;
	register int dest_exponent, count;
	register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
	int is_tiny;

	opnd1 = *srcptr1;
	opnd2 = *srcptr2;
	/*
	 * set sign bit of result
	 */
	if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);
	else Sgl_setzero(result);
	/*
	 * check first operand for NaN's or infinity
	 */
	if (Sgl_isinfinity_exponent(opnd1)) {
		if (Sgl_iszero_mantissa(opnd1)) {
			if (Sgl_isnotnan(opnd2)) {
				if (Sgl_isinfinity(opnd2)) {
					/*
					 * invalid since both operands
					 * are infinity
					 */
					if (Is_invalidtrap_enabled())
						return(INVALIDEXCEPTION);
					Set_invalidflag();
					Sgl_makequietnan(result);
					*dstptr = result;
					return(NOEXCEPTION);
				}
				/*
				 * return infinity
				 */
				Sgl_setinfinity_exponentmantissa(result);
				*dstptr = result;
				return(NOEXCEPTION);
			}
		}
		else {
			/*
			 * is NaN; signaling or quiet?
			 */
			if (Sgl_isone_signaling(opnd1)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Sgl_set_quiet(opnd1);
			}
			/*
			 * is second operand a signaling NaN?
			 */
			else if (Sgl_is_signalingnan(opnd2)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Sgl_set_quiet(opnd2);
				*dstptr = opnd2;
				return(NOEXCEPTION);
			}
			/*
			 * return quiet NaN
			 */
			*dstptr = opnd1;
			return(NOEXCEPTION);
		}
	}
	/*
	 * check second operand for NaN's or infinity
	 */
	if (Sgl_isinfinity_exponent(opnd2)) {
		if (Sgl_iszero_mantissa(opnd2)) {
			/*
			 * return zero
			 */
			Sgl_setzero_exponentmantissa(result);
			*dstptr = result;
			return(NOEXCEPTION);
		}
		/*
		 * is NaN; signaling or quiet?
		 */
		if (Sgl_isone_signaling(opnd2)) {
			/* trap if INVALIDTRAP enabled */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
			/* make NaN quiet */
			Set_invalidflag();
			Sgl_set_quiet(opnd2);
		}
		/*
		 * return quiet NaN
		 */
		*dstptr = opnd2;
		return(NOEXCEPTION);
	}
	/*
	 * check for division by zero
	 */
	if (Sgl_iszero_exponentmantissa(opnd2)) {
		if (Sgl_iszero_exponentmantissa(opnd1)) {
			/* invalid since both operands are zero */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
			Set_invalidflag();
			Sgl_makequietnan(result);
			*dstptr = result;
			return(NOEXCEPTION);
		}
		if (Is_divisionbyzerotrap_enabled())
			return(DIVISIONBYZEROEXCEPTION);
		Set_divisionbyzeroflag();
		Sgl_setinfinity_exponentmantissa(result);
		*dstptr = result;
		return(NOEXCEPTION);
	}
	/*
	 * Generate exponent
	 */
	dest_exponent = Sgl_exponent(opnd1) - Sgl_exponent(opnd2) + SGL_BIAS;

	/*
	 * Generate mantissa
	 */
	if (Sgl_isnotzero_exponent(opnd1)) {
		/* set hidden bit */
		Sgl_clear_signexponent_set_hidden(opnd1);
	}
	else {
		/* check for zero */
		if (Sgl_iszero_mantissa(opnd1)) {
			Sgl_setzero_exponentmantissa(result);
			*dstptr = result;
			return(NOEXCEPTION);
		}
		/* is denormalized; want to normalize */
		Sgl_clear_signexponent(opnd1);
		Sgl_leftshiftby1(opnd1);
		Sgl_normalize(opnd1,dest_exponent);
	}
	/* opnd2 needs to have hidden bit set with msb in hidden bit */
	if (Sgl_isnotzero_exponent(opnd2)) {
		Sgl_clear_signexponent_set_hidden(opnd2);
	}
	else {
		/* is denormalized; want to normalize */
		Sgl_clear_signexponent(opnd2);
		Sgl_leftshiftby1(opnd2);
		while(Sgl_iszero_hiddenhigh7mantissa(opnd2)) {
			Sgl_leftshiftby8(opnd2);
			dest_exponent += 8;
		}
		if(Sgl_iszero_hiddenhigh3mantissa(opnd2)) {
			Sgl_leftshiftby4(opnd2);
			dest_exponent += 4;
		}
		while(Sgl_iszero_hidden(opnd2)) {
			Sgl_leftshiftby1(opnd2);
			dest_exponent += 1;
		}
	}

	/* Divide the source mantissas */

	/*
	 * A non_restoring divide algorithm is used.
	 */
	Sgl_subtract(opnd1,opnd2,opnd1);
	Sgl_setzero(opnd3);
	for (count=1;count<=SGL_P && Sgl_all(opnd1);count++) {
		Sgl_leftshiftby1(opnd1);
		Sgl_leftshiftby1(opnd3);
		if (Sgl_iszero_sign(opnd1)) {
			Sgl_setone_lowmantissa(opnd3);
			Sgl_subtract(opnd1,opnd2,opnd1);
		}
		else Sgl_addition(opnd1,opnd2,opnd1);
	}
	if (count <= SGL_P) {
		Sgl_leftshiftby1(opnd3);
		Sgl_setone_lowmantissa(opnd3);
		Sgl_leftshift(opnd3,SGL_P-count);
		if (Sgl_iszero_hidden(opnd3)) {
			Sgl_leftshiftby1(opnd3);
			dest_exponent--;
		}
	}
	else {
		if (Sgl_iszero_hidden(opnd3)) {
			/* need to get one more bit of result */
			Sgl_leftshiftby1(opnd1);
			Sgl_leftshiftby1(opnd3);
			if (Sgl_iszero_sign(opnd1)) {
				Sgl_setone_lowmantissa(opnd3);
				Sgl_subtract(opnd1,opnd2,opnd1);
			}
			else Sgl_addition(opnd1,opnd2,opnd1);
			dest_exponent--;
		}
		if (Sgl_iszero_sign(opnd1)) guardbit = TRUE;
		stickybit = Sgl_all(opnd1);
	}
	inexact = guardbit | stickybit;

	/*
	 * round result
	 */
	if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
		Sgl_clear_signexponent(opnd3);
		switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Sgl_iszero_sign(result))
					Sgl_increment_mantissa(opnd3);
				break;
			case ROUNDMINUS:
				if (Sgl_isone_sign(result))
					Sgl_increment_mantissa(opnd3);
				break;
			case ROUNDNEAREST:
				if (guardbit &&
				    (stickybit || Sgl_isone_lowmantissa(opnd3)))
					Sgl_increment_mantissa(opnd3);
		}
		if (Sgl_isone_hidden(opnd3)) dest_exponent++;
	}
	Sgl_set_mantissa(result,opnd3);

	/*
	 * Test for overflow
	 */
	if (dest_exponent >= SGL_INFINITY_EXPONENT) {
		/* trap if OVERFLOWTRAP enabled */
		if (Is_overflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
			*dstptr = result;
			if (inexact) {
			    if (Is_inexacttrap_enabled())
				return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
			    else Set_inexactflag();
			}
			return(OVERFLOWEXCEPTION);
		}
		Set_overflowflag();
		/* set result to infinity or largest number */
		Sgl_setoverflow(result);
		inexact = TRUE;
	}
	/*
	 * Test for underflow
	 */
	else if (dest_exponent <= 0) {
		/* trap if UNDERFLOWTRAP enabled */
		if (Is_underflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Sgl_setwrapped_exponent(result,dest_exponent,unfl);
			*dstptr = result;
			if (inexact) {
			    if (Is_inexacttrap_enabled())
				return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
			    else Set_inexactflag();
			}
			return(UNDERFLOWEXCEPTION);
		}

		/* Determine if should set underflow flag */
		is_tiny = TRUE;
		if (dest_exponent == 0 && inexact) {
			switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Sgl_iszero_sign(result)) {
					Sgl_increment(opnd3);
					if (Sgl_isone_hiddenoverflow(opnd3))
						is_tiny = FALSE;
					Sgl_decrement(opnd3);
				}
				break;
			case ROUNDMINUS:
				if (Sgl_isone_sign(result)) {
					Sgl_increment(opnd3);
					if (Sgl_isone_hiddenoverflow(opnd3))
						is_tiny = FALSE;
					Sgl_decrement(opnd3);
				}
				break;
			case ROUNDNEAREST:
				if (guardbit && (stickybit ||
				    Sgl_isone_lowmantissa(opnd3))) {
					Sgl_increment(opnd3);
					if (Sgl_isone_hiddenoverflow(opnd3))
						is_tiny = FALSE;
					Sgl_decrement(opnd3);
				}
				break;
			}
		}

		/*
		 * denormalize result or set to signed zero
		 */
		stickybit = inexact;
		Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);

		/* return rounded number */
		if (inexact) {
			switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Sgl_iszero_sign(result)) {
					Sgl_increment(opnd3);
				}
				break;
			case ROUNDMINUS:
				if (Sgl_isone_sign(result))  {
					Sgl_increment(opnd3);
				}
				break;
			case ROUNDNEAREST:
				if (guardbit && (stickybit ||
				    Sgl_isone_lowmantissa(opnd3))) {
					Sgl_increment(opnd3);
				}
				break;
			}
			if (is_tiny) Set_underflowflag();
		}
		Sgl_set_exponentmantissa(result,opnd3);
	}
	else Sgl_set_exponent(result,dest_exponent);
	*dstptr = result;
	/* check for inexact */
	if (inexact) {
		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
		else  Set_inexactflag();
	}
	return(NOEXCEPTION);
}


code:
/*	$OpenBSD: sfmpy.c,v 1.5 2002/05/07 22:19:30 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)sfmpy.c: Revision: 2.9.88.1 Date: 93/12/07 15:07:07 */

#include "float.h"
#include "sgl_float.h"

/*
 *  Single Precision Floating-point Multiply
 */
int
sgl_fmpy(srcptr1,srcptr2,dstptr,status)
	sgl_floating_point *srcptr1, *srcptr2, *dstptr;
	unsigned int *status;
{
	register unsigned int opnd1, opnd2, opnd3, result;
	register int dest_exponent, count;
	register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
	int is_tiny;

	opnd1 = *srcptr1;
	opnd2 = *srcptr2;
	/*
	 * set sign bit of result
	 */
	if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);
	else Sgl_setzero(result);
	/*
	 * check first operand for NaN's or infinity
	 */
	if (Sgl_isinfinity_exponent(opnd1)) {
		if (Sgl_iszero_mantissa(opnd1)) {
			if (Sgl_isnotnan(opnd2)) {
				if (Sgl_iszero_exponentmantissa(opnd2)) {
					/*
					 * invalid since operands are infinity
					 * and zero
					 */
					if (Is_invalidtrap_enabled())
						return(INVALIDEXCEPTION);
					Set_invalidflag();
					Sgl_makequietnan(result);
					*dstptr = result;
					return(NOEXCEPTION);
				}
				/*
				 * return infinity
				 */
				Sgl_setinfinity_exponentmantissa(result);
				*dstptr = result;
				return(NOEXCEPTION);
			}
		}
		else {
			/*
			 * is NaN; signaling or quiet?
			 */
			if (Sgl_isone_signaling(opnd1)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Sgl_set_quiet(opnd1);
			}
			/*
			 * is second operand a signaling NaN?
			 */
			else if (Sgl_is_signalingnan(opnd2)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Sgl_set_quiet(opnd2);
				*dstptr = opnd2;
				return(NOEXCEPTION);
			}
			/*
			 * return quiet NaN
			 */
			*dstptr = opnd1;
			return(NOEXCEPTION);
		}
	}
	/*
	 * check second operand for NaN's or infinity
	 */
	if (Sgl_isinfinity_exponent(opnd2)) {
		if (Sgl_iszero_mantissa(opnd2)) {
			if (Sgl_iszero_exponentmantissa(opnd1)) {
				/* invalid since operands are zero & infinity */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				Set_invalidflag();
				Sgl_makequietnan(opnd2);
				*dstptr = opnd2;
				return(NOEXCEPTION);
			}
			/*
			 * return infinity
			 */
			Sgl_setinfinity_exponentmantissa(result);
			*dstptr = result;
			return(NOEXCEPTION);
		}
		/*
		 * is NaN; signaling or quiet?
		 */
		if (Sgl_isone_signaling(opnd2)) {
			/* trap if INVALIDTRAP enabled */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);

			/* make NaN quiet */
			Set_invalidflag();
			Sgl_set_quiet(opnd2);
		}
		/*
		 * return quiet NaN
		 */
		*dstptr = opnd2;
		return(NOEXCEPTION);
	}
	/*
	 * Generate exponent
	 */
	dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS;

	/*
	 * Generate mantissa
	 */
	if (Sgl_isnotzero_exponent(opnd1)) {
		/* set hidden bit */
		Sgl_clear_signexponent_set_hidden(opnd1);
	}
	else {
		/* check for zero */
		if (Sgl_iszero_mantissa(opnd1)) {
			Sgl_setzero_exponentmantissa(result);
			*dstptr = result;
			return(NOEXCEPTION);
		}
		/* is denormalized, adjust exponent */
		Sgl_clear_signexponent(opnd1);
		Sgl_leftshiftby1(opnd1);
		Sgl_normalize(opnd1,dest_exponent);
	}
	/* opnd2 needs to have hidden bit set with msb in hidden bit */
	if (Sgl_isnotzero_exponent(opnd2)) {
		Sgl_clear_signexponent_set_hidden(opnd2);
	}
	else {
		/* check for zero */
		if (Sgl_iszero_mantissa(opnd2)) {
			Sgl_setzero_exponentmantissa(result);
			*dstptr = result;
			return(NOEXCEPTION);
		}
		/* is denormalized; want to normalize */
		Sgl_clear_signexponent(opnd2);
		Sgl_leftshiftby1(opnd2);
		Sgl_normalize(opnd2,dest_exponent);
	}

	/* Multiply two source mantissas together */

	Sgl_leftshiftby4(opnd2);     /* make room for guard bits */
	Sgl_setzero(opnd3);
	/*
	 * Four bits at a time are inspected in each loop, and a
	 * simple shift and add multiply algorithm is used.
	 */
	for (count=1;count<SGL_P;count+=4) {
		stickybit |= Slow4(opnd3);
		Sgl_rightshiftby4(opnd3);
		if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3);
		if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2);
		if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1);
		if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2);
		Sgl_rightshiftby4(opnd1);
	}
	/* make sure result is left-justified */
	if (Sgl_iszero_sign(opnd3)) {
		Sgl_leftshiftby1(opnd3);
	}
	else {
		/* result mantissa >= 2. */
		dest_exponent++;
	}
	/* check for denormalized result */
	while (Sgl_iszero_sign(opnd3)) {
		Sgl_leftshiftby1(opnd3);
		dest_exponent--;
	}
	/*
	 * check for guard, sticky and inexact bits
	 */
	stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1);
	guardbit = Sbit24(opnd3);
	inexact = guardbit | stickybit;

	/* re-align mantissa */
	Sgl_rightshiftby8(opnd3);

	/*
	 * round result
	 */
	if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) {
		Sgl_clear_signexponent(opnd3);
		switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Sgl_iszero_sign(result))
					Sgl_increment(opnd3);
				break;
			case ROUNDMINUS:
				if (Sgl_isone_sign(result))
					Sgl_increment(opnd3);
				break;
			case ROUNDNEAREST:
				if (guardbit &&
				    (stickybit || Sgl_isone_lowmantissa(opnd3)))
					Sgl_increment(opnd3);
				break;
		}
		if (Sgl_isone_hidden(opnd3)) dest_exponent++;
	}
	Sgl_set_mantissa(result,opnd3);

	/*
	 * Test for overflow
	 */
	if (dest_exponent >= SGL_INFINITY_EXPONENT) {
		/* trap if OVERFLOWTRAP enabled */
		if (Is_overflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
			*dstptr = result;
			if (inexact) {
			    if (Is_inexacttrap_enabled())
				return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
			    else Set_inexactflag();
			}
			return(OVERFLOWEXCEPTION);
		}
		inexact = TRUE;
		Set_overflowflag();
		/* set result to infinity or largest number */
		Sgl_setoverflow(result);
	}
	/*
	 * Test for underflow
	 */
	else if (dest_exponent <= 0) {
		/* trap if UNDERFLOWTRAP enabled */
		if (Is_underflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Sgl_setwrapped_exponent(result,dest_exponent,unfl);
			*dstptr = result;
			if (inexact) {
			    if (Is_inexacttrap_enabled())
				return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
			    else Set_inexactflag();
			}
			return(UNDERFLOWEXCEPTION);
		}

		/* Determine if should set underflow flag */
		is_tiny = TRUE;
		if (dest_exponent == 0 && inexact) {
			switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Sgl_iszero_sign(result)) {
					Sgl_increment(opnd3);
					if (Sgl_isone_hiddenoverflow(opnd3))
						is_tiny = FALSE;
					Sgl_decrement(opnd3);
				}
				break;
			case ROUNDMINUS:
				if (Sgl_isone_sign(result)) {
					Sgl_increment(opnd3);
					if (Sgl_isone_hiddenoverflow(opnd3))
						is_tiny = FALSE;
					Sgl_decrement(opnd3);
				}
				break;
			case ROUNDNEAREST:
				if (guardbit && (stickybit ||
				    Sgl_isone_lowmantissa(opnd3))) {
					Sgl_increment(opnd3);
					if (Sgl_isone_hiddenoverflow(opnd3))
						is_tiny = FALSE;
					Sgl_decrement(opnd3);
				}
				break;
			}
		}

		/*
		 * denormalize result or set to signed zero
		 */
		stickybit = inexact;
		Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);

		/* return zero or smallest number */
		if (inexact) {
			switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Sgl_iszero_sign(result))
					Sgl_increment(opnd3);
				break;
			case ROUNDMINUS:
				if (Sgl_isone_sign(result))
					Sgl_increment(opnd3);
				break;
			case ROUNDNEAREST:
				if (guardbit && (stickybit ||
				    Sgl_isone_lowmantissa(opnd3)))
					Sgl_increment(opnd3);
				break;
			}
		if (is_tiny) Set_underflowflag();
		}
		Sgl_set_exponentmantissa(result,opnd3);
	}
	else Sgl_set_exponent(result,dest_exponent);
	*dstptr = result;

	/* check for inexact */
	if (inexact) {
		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
		else Set_inexactflag();
	}
	return(NOEXCEPTION);
}


code:
/*	$OpenBSD: sfrem.c,v 1.5 2002/05/07 22:19:30 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)sfrem.c: Revision: 1.7.88.1 Date: 93/12/07 15:07:10 */

#include "float.h"
#include "sgl_float.h"

/*
 *  Single Precision Floating-point Remainder
 */
int
sgl_frem(srcptr1,srcptr2,dstptr,status)
	sgl_floating_point *srcptr1, *srcptr2, *dstptr;
	unsigned int *status;
{
	register unsigned int opnd1, opnd2, result;
	register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
	register int roundup = FALSE;

	opnd1 = *srcptr1;
	opnd2 = *srcptr2;
	/*
	 * check first operand for NaN's or infinity
	 */
	if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) {
		if (Sgl_iszero_mantissa(opnd1)) {
			if (Sgl_isnotnan(opnd2)) {
				/* invalid since first operand is infinity */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				Set_invalidflag();
				Sgl_makequietnan(result);
				*dstptr = result;
				return(NOEXCEPTION);
			}
		}
		else {
			/*
			 * is NaN; signaling or quiet?
			 */
			if (Sgl_isone_signaling(opnd1)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Sgl_set_quiet(opnd1);
			}
			/*
			 * is second operand a signaling NaN?
			 */
			else if (Sgl_is_signalingnan(opnd2)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Sgl_set_quiet(opnd2);
				*dstptr = opnd2;
				return(NOEXCEPTION);
			}
			/*
			 * return quiet NaN
			 */
			*dstptr = opnd1;
			return(NOEXCEPTION);
		}
	}
	/*
	 * check second operand for NaN's or infinity
	 */
	if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) {
		if (Sgl_iszero_mantissa(opnd2)) {
			/*
			 * return first operand
			 */
			*dstptr = opnd1;
			return(NOEXCEPTION);
		}
		/*
		 * is NaN; signaling or quiet?
		 */
		if (Sgl_isone_signaling(opnd2)) {
			/* trap if INVALIDTRAP enabled */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
			/* make NaN quiet */
			Set_invalidflag();
			Sgl_set_quiet(opnd2);
		}
		/*
		 * return quiet NaN
		 */
		*dstptr = opnd2;
		return(NOEXCEPTION);
	}
	/*
	 * check second operand for zero
	 */
	if (Sgl_iszero_exponentmantissa(opnd2)) {
		/* invalid since second operand is zero */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		Set_invalidflag();
		Sgl_makequietnan(result);
		*dstptr = result;
		return(NOEXCEPTION);
	}

	/*
	 * get sign of result
	 */
	result = opnd1;

	/*
	 * check for denormalized operands
	 */
	if (opnd1_exponent == 0) {
		/* check for zero */
		if (Sgl_iszero_mantissa(opnd1)) {
			*dstptr = opnd1;
			return(NOEXCEPTION);
		}
		/* normalize, then continue */
		opnd1_exponent = 1;
		Sgl_normalize(opnd1,opnd1_exponent);
	}
	else {
		Sgl_clear_signexponent_set_hidden(opnd1);
	}
	if (opnd2_exponent == 0) {
		/* normalize, then continue */
		opnd2_exponent = 1;
		Sgl_normalize(opnd2,opnd2_exponent);
	}
	else {
		Sgl_clear_signexponent_set_hidden(opnd2);
	}

	/* find result exponent and divide step loop count */
	dest_exponent = opnd2_exponent - 1;
	stepcount = opnd1_exponent - opnd2_exponent;

	/*
	 * check for opnd1/opnd2 < 1
	 */
	if (stepcount < 0) {
		/*
		 * check for opnd1/opnd2 > 1/2
		 *
		 * In this case n will round to 1, so
		 *    r = opnd1 - opnd2
		 */
		if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) {
			Sgl_all(result) = ~Sgl_all(result);   /* set sign */
			/* align opnd2 with opnd1 */
			Sgl_leftshiftby1(opnd2);
			Sgl_subtract(opnd2,opnd1,opnd2);
			/* now normalize */
			while (Sgl_iszero_hidden(opnd2)) {
				Sgl_leftshiftby1(opnd2);
				dest_exponent--;
			}
			Sgl_set_exponentmantissa(result,opnd2);
			goto testforunderflow;
		}
		/*
		 * opnd1/opnd2 <= 1/2
		 *
		 * In this case n will round to zero, so
		 *    r = opnd1
		 */
		Sgl_set_exponentmantissa(result,opnd1);
		dest_exponent = opnd1_exponent;
		goto testforunderflow;
	}

	/*
	 * Generate result
	 *
	 * Do iterative subtract until remainder is less than operand 2.
	 */
	while (stepcount-- > 0 && Sgl_all(opnd1)) {
		if (Sgl_isnotlessthan(opnd1,opnd2))
			Sgl_subtract(opnd1,opnd2,opnd1);
		Sgl_leftshiftby1(opnd1);
	}
	/*
	 * Do last subtract, then determine which way to round if remainder
	 * is exactly 1/2 of opnd2
	 */
	if (Sgl_isnotlessthan(opnd1,opnd2)) {
		Sgl_subtract(opnd1,opnd2,opnd1);
		roundup = TRUE;
	}
	if (stepcount > 0 || Sgl_iszero(opnd1)) {
		/* division is exact, remainder is zero */
		Sgl_setzero_exponentmantissa(result);
		*dstptr = result;
		return(NOEXCEPTION);
	}

	/*
	 * Check for cases where opnd1/opnd2 < n
	 *
	 * In this case the result's sign will be opposite that of
	 * opnd1.  The mantissa also needs some correction.
	 */
	Sgl_leftshiftby1(opnd1);
	if (Sgl_isgreaterthan(opnd1,opnd2)) {
		Sgl_invert_sign(result);
		Sgl_subtract((opnd2<<1),opnd1,opnd1);
	}
	/* check for remainder being exactly 1/2 of opnd2 */
	else if (Sgl_isequal(opnd1,opnd2) && roundup) {
		Sgl_invert_sign(result);
	}

	/* normalize result's mantissa */
	while (Sgl_iszero_hidden(opnd1)) {
		dest_exponent--;
		Sgl_leftshiftby1(opnd1);
	}
	Sgl_set_exponentmantissa(result,opnd1);

	/*
	 * Test for underflow
	 */
    testforunderflow:
	if (dest_exponent <= 0) {
		/* trap if UNDERFLOWTRAP enabled */
		if (Is_underflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Sgl_setwrapped_exponent(result,dest_exponent,unfl);
			*dstptr = result;
			/* frem is always exact */
			return(UNDERFLOWEXCEPTION);
		}
		/*
		 * denormalize result or set to signed zero
		 */
		if (dest_exponent >= (1 - SGL_P)) {
			Sgl_rightshift_exponentmantissa(result,1-dest_exponent);
		} else {
			Sgl_setzero_exponentmantissa(result);
		}
	}
	else Sgl_set_exponent(result,dest_exponent);
	*dstptr = result;
	return(NOEXCEPTION);
}


code:
/*	$OpenBSD: sfsqrt.c,v 1.7 2003/04/10 17:27:59 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)sfsqrt.c: Revision: 1.9.88.1 Date: 93/12/07 15:07:13 */

#include "float.h"
#include "sgl_float.h"

/*
 *  Single Floating-point Square Root
 */

/*ARGSUSED*/
int
sgl_fsqrt(srcptr, null, dstptr, status)
	sgl_floating_point *srcptr, *null, *dstptr;
	unsigned int *status;
{
	register unsigned int src, result;
	register int src_exponent, newbit, sum;
	register int guardbit = FALSE, even_exponent;

	src = *srcptr;
	/*
	 * check source operand for NaN or infinity
	 */
	if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) {
		/*
		 * is signaling NaN?
		 */
		if (Sgl_isone_signaling(src)) {
			/* trap if INVALIDTRAP enabled */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
			/* make NaN quiet */
			Set_invalidflag();
			Sgl_set_quiet(src);
		}
		/*
		 * Return quiet NaN or positive infinity.
		 *  Fall thru to negative test if negative infinity.
		 */
		if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) {
			*dstptr = src;
			return(NOEXCEPTION);
		}
	}

	/*
	 * check for zero source operand
	 */
	if (Sgl_iszero_exponentmantissa(src)) {
		*dstptr = src;
		return(NOEXCEPTION);
	}

	/*
	 * check for negative source operand
	 */
	if (Sgl_isone_sign(src)) {
		/* trap if INVALIDTRAP enabled */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		/* make NaN quiet */
		Set_invalidflag();
		Sgl_makequietnan(src);
		*dstptr = src;
		return(NOEXCEPTION);
	}

	/*
	 * Generate result
	 */
	if (src_exponent > 0) {
		even_exponent = Sgl_hidden(src);
		Sgl_clear_signexponent_set_hidden(src);
	}
	else {
		/* normalize operand */
		Sgl_clear_signexponent(src);
		src_exponent++;
		Sgl_normalize(src,src_exponent);
		even_exponent = src_exponent & 1;
	}
	if (even_exponent) {
		/* exponent is even */
		/* Add comment here.  Explain why odd exponent needs correction */
		Sgl_leftshiftby1(src);
	}
	/*
	 * Add comment here.  Explain following algorithm.
	 *
	 * Trust me, it works.
	 *
	 */
	Sgl_setzero(result);
	newbit = 1 << SGL_P;
	while (newbit && Sgl_isnotzero(src)) {
		Sgl_addition(result,newbit,sum);
		if(sum <= Sgl_all(src)) {
			/* update result */
			Sgl_addition(result,(newbit<<1),result);
			Sgl_subtract(src,sum,src);
		}
		Sgl_rightshiftby1(newbit);
		Sgl_leftshiftby1(src);
	}
	/* correct exponent for pre-shift */
	if (even_exponent) {
		Sgl_rightshiftby1(result);
	}

	/* check for inexact */
	if (Sgl_isnotzero(src)) {
		if (!even_exponent & Sgl_islessthan(result,src))
			Sgl_increment(result);
		guardbit = Sgl_lowmantissa(result);
		Sgl_rightshiftby1(result);

		/*  now round result  */
		switch (Rounding_mode()) {
		case ROUNDPLUS:
		     Sgl_increment(result);
		     break;
		case ROUNDNEAREST:
		     /* stickybit is always true, so guardbit
		      * is enough to determine rounding */
		     if (guardbit) {
			Sgl_increment(result);
		     }
		     break;
		}
		/* increment result exponent by 1 if mantissa overflowed */
		if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2;

		if (Is_inexacttrap_enabled()) {
			Sgl_set_exponent(result,
			 ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
			*dstptr = result;
			return(INEXACTEXCEPTION);
		}
		else Set_inexactflag();
	}
	else {
		Sgl_rightshiftby1(result);
	}
	Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
	*dstptr = result;
	return(NOEXCEPTION);
}


code:
/*	$OpenBSD: sfsub.c,v 1.6 2002/11/29 09:27:34 deraadt Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)sfsub.c: Revision: 2.7.88.1 Date: 93/12/07 15:07:15 */

#include "float.h"
#include "sgl_float.h"

/*
 * Single_subtract: subtract two single precision values.
 */
int
sgl_fsub(leftptr, rightptr, dstptr, status)
    sgl_floating_point *leftptr, *rightptr, *dstptr;
    unsigned int *status;
{
    register unsigned int left, right, result, extent;
    register unsigned int signless_upper_left, signless_upper_right, save;

    register int result_exponent, right_exponent, diff_exponent;
    register int sign_save, jumpsize;
    register int inexact = FALSE, underflowtrap;

    /* Create local copies of the numbers */
    left = *leftptr;
    right = *rightptr;

    /* A zero "save" helps discover equal operands (for later),	*
     * and is used in swapping operands (if needed).		*/
    Sgl_xortointp1(left,right,/*to*/save);

    /*
     * check first operand for NaN's or infinity
     */
    if ((result_exponent = Sgl_exponent(left)) == SGL_INFINITY_EXPONENT)
	{
	if (Sgl_iszero_mantissa(left))
	    {
	    if (Sgl_isnotnan(right))
		{
		if (Sgl_isinfinity(right) && save==0)
		    {
		    /*
		     * invalid since operands are same signed infinity's
		     */
		    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		    Set_invalidflag();
		    Sgl_makequietnan(result);
		    *dstptr = result;
		    return(NOEXCEPTION);
		    }
		/*
		 * return infinity
		 */
		*dstptr = left;
		return(NOEXCEPTION);
		}
	    }
	else
	    {
	    /*
	     * is NaN; signaling or quiet?
	     */
	    if (Sgl_isone_signaling(left))
		{
		/* trap if INVALIDTRAP enabled */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		/* make NaN quiet */
		Set_invalidflag();
		Sgl_set_quiet(left);
		}
	    /*
	     * is second operand a signaling NaN?
	     */
	    else if (Sgl_is_signalingnan(right))
		{
		/* trap if INVALIDTRAP enabled */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		/* make NaN quiet */
		Set_invalidflag();
		Sgl_set_quiet(right);
		*dstptr = right;
		return(NOEXCEPTION);
		}
	    /*
	     * return quiet NaN
	     */
	    *dstptr = left;
	    return(NOEXCEPTION);
	    }
	} /* End left NaN or Infinity processing */
    /*
     * check second operand for NaN's or infinity
     */
    if (Sgl_isinfinity_exponent(right))
	{
	if (Sgl_iszero_mantissa(right))
	    {
	    /* return infinity */
	    Sgl_invert_sign(right);
	    *dstptr = right;
	    return(NOEXCEPTION);
	    }
	/*
	 * is NaN; signaling or quiet?
	 */
	if (Sgl_isone_signaling(right))
	    {
	    /* trap if INVALIDTRAP enabled */
	    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
	    /* make NaN quiet */
	    Set_invalidflag();
	    Sgl_set_quiet(right);
	    }
	/*
	 * return quiet NaN
	 */
	*dstptr = right;
	return(NOEXCEPTION);
	} /* End right NaN or Infinity processing */

    /* Invariant: Must be dealing with finite numbers */

    /* Compare operands by removing the sign */
    Sgl_copytoint_exponentmantissa(left,signless_upper_left);
    Sgl_copytoint_exponentmantissa(right,signless_upper_right);

    /* sign difference selects sub or add operation. */
    if(Sgl_ismagnitudeless(signless_upper_left,signless_upper_right))
	{
	/* Set the left operand to the larger one by XOR swap *
	 *  First finish the first word using "save"	  */
	Sgl_xorfromintp1(save,right,/*to*/right);
	Sgl_xorfromintp1(save,left,/*to*/left);
	result_exponent = Sgl_exponent(left);
	Sgl_invert_sign(left);
	}
    /* Invariant:  left is not smaller than right. */

    if((right_exponent = Sgl_exponent(right)) == 0)
	{
	/* Denormalized operands.  First look for zeroes */
	if(Sgl_iszero_mantissa(right))
	    {
	    /* right is zero */
	    if(Sgl_iszero_exponentmantissa(left))
		{
		/* Both operands are zeros */
		Sgl_invert_sign(right);
		if(Is_rounding_mode(ROUNDMINUS))
		    {
		    Sgl_or_signs(left,/*with*/right);
		    }
		else
		    {
		    Sgl_and_signs(left,/*with*/right);
		    }
		}
	    else
		{
		/* Left is not a zero and must be the result.  Trapped
		 * underflows are signaled if left is denormalized.  Result
		 * is always exact. */
		if( (result_exponent == 0) && Is_underflowtrap_enabled() )
		    {
		    /* need to normalize results mantissa */
		    sign_save = Sgl_signextendedsign(left);
		    Sgl_leftshiftby1(left);
		    Sgl_normalize(left,result_exponent);
		    Sgl_set_sign(left,/*using*/sign_save);
		    Sgl_setwrapped_exponent(left,result_exponent,unfl);
		    *dstptr = left;
		    /* inexact = FALSE */
		    return(UNDERFLOWEXCEPTION);
		    }
		}
	    *dstptr = left;
	    return(NOEXCEPTION);
	    }

	/* Neither are zeroes */
	Sgl_clear_sign(right);	/* Exponent is already cleared */
	if(result_exponent == 0 )
	    {
	    /* Both operands are denormalized.  The result must be exact
	     * and is simply calculated.  A sum could become normalized and a
	     * difference could cancel to a true zero. */
	    if( (/*signed*/int) save >= 0 )
		{
		Sgl_subtract(left,/*minus*/right,/*into*/result);
		if(Sgl_iszero_mantissa(result))
		    {
		    if(Is_rounding_mode(ROUNDMINUS))
			{
			Sgl_setone_sign(result);
			}
		    else
			{
			Sgl_setzero_sign(result);
			}
		    *dstptr = result;
		    return(NOEXCEPTION);
		    }
		}
	    else
		{
		Sgl_addition(left,right,/*into*/result);
		if(Sgl_isone_hidden(result))
		    {
		    *dstptr = result;
		    return(NOEXCEPTION);
		    }
		}
	    if(Is_underflowtrap_enabled())
		{
		/* need to normalize result */
		sign_save = Sgl_signextendedsign(result);
		Sgl_leftshiftby1(result);
		Sgl_normalize(result,result_exponent);
		Sgl_set_sign(result,/*using*/sign_save);
		Sgl_setwrapped_exponent(result,result_exponent,unfl);
		*dstptr = result;
		/* inexact = FALSE */
		return(UNDERFLOWEXCEPTION);
		}
	    *dstptr = result;
	    return(NOEXCEPTION);
	    }
	right_exponent = 1;	/* Set exponent to reflect different bias
				 * with denomalized numbers. */
	}
    else
	{
	Sgl_clear_signexponent_set_hidden(right);
	}
    Sgl_clear_exponent_set_hidden(left);
    diff_exponent = result_exponent - right_exponent;

    /*
     * Special case alignment of operands that would force alignment
     * beyond the extent of the extension.  A further optimization
     * could special case this but only reduces the path length for this
     * infrequent case.
     */
    if(diff_exponent > SGL_THRESHOLD)
	{
	diff_exponent = SGL_THRESHOLD;
	}

    /* Align right operand by shifting to right */
    Sgl_right_align(/*operand*/right,/*shifted by*/diff_exponent,
      /*and lower to*/extent);

    /* Treat sum and difference of the operands separately. */
    if( (/*signed*/int) save >= 0 )
	{
	/*
	 * Difference of the two operands.  Their can be no overflow.  A
	 * borrow can occur out of the hidden bit and force a post
	 * normalization phase.
	 */
	Sgl_subtract_withextension(left,/*minus*/right,/*with*/extent,/*into*/result);
	if(Sgl_iszero_hidden(result))
	    {
	    /* Handle normalization */
	    /* A straight forward algorithm would now shift the result
	     * and extension left until the hidden bit becomes one.  Not
	     * all of the extension bits need participate in the shift.
	     * Only the two most significant bits (round and guard) are
	     * needed.  If only a single shift is needed then the guard
	     * bit becomes a significant low order bit and the extension
	     * must participate in the rounding.  If more than a single
	     * shift is needed, then all bits to the right of the guard
	     * bit are zeros, and the guard bit may or may not be zero. */
	    sign_save = Sgl_signextendedsign(result);
	    Sgl_leftshiftby1_withextent(result,extent,result);

	    /* Need to check for a zero result.  The sign and exponent
	     * fields have already been zeroed.  The more efficient test
	     * of the full object can be used.
	     */
	    if(Sgl_iszero(result))
		/* Must have been "x-x" or "x+(-x)". */
		{
		if(Is_rounding_mode(ROUNDMINUS)) Sgl_setone_sign(result);
		*dstptr = result;
		return(NOEXCEPTION);
		}
	    result_exponent--;
	    /* Look to see if normalization is finished. */
	    if(Sgl_isone_hidden(result))
		{
		if(result_exponent==0)
		    {
		    /* Denormalized, exponent should be zero.  Left operand *
		     * was normalized, so extent (guard, round) was zero    */
		    goto underflow;
		    }
		else
		    {
		    /* No further normalization is needed. */
		    Sgl_set_sign(result,/*using*/sign_save);
		    Ext_leftshiftby1(extent);
		    goto round;
		    }
		}

	    /* Check for denormalized, exponent should be zero.  Left    *
	     * operand was normalized, so extent (guard, round) was zero */
	    if(!(underflowtrap = Is_underflowtrap_enabled()) &&
	       result_exponent==0) goto underflow;

	    /* Shift extension to complete one bit of normalization and
	     * update exponent. */
	    Ext_leftshiftby1(extent);

	    /* Discover first one bit to determine shift amount.  Use a
	     * modified binary search.  We have already shifted the result
	     * one position right and still not found a one so the remainder
	     * of the extension must be zero and simplifies rounding. */
	    /* Scan bytes */
	    while(Sgl_iszero_hiddenhigh7mantissa(result))
		{
		Sgl_leftshiftby8(result);
		if((result_exponent -= 8) <= 0  && !underflowtrap)
		    goto underflow;
		}
	    /* Now narrow it down to the nibble */
	    if(Sgl_iszero_hiddenhigh3mantissa(result))
		{
		/* The lower nibble contains the normalizing one */
		Sgl_leftshiftby4(result);
		if((result_exponent -= 4) <= 0 && !underflowtrap)
		    goto underflow;
		}
	    /* Select case were first bit is set (already normalized)
	     * otherwise select the proper shift. */
	    if((jumpsize = Sgl_hiddenhigh3mantissa(result)) > 7)
		{
		/* Already normalized */
		if(result_exponent <= 0) goto underflow;
		Sgl_set_sign(result,/*using*/sign_save);
		Sgl_set_exponent(result,/*using*/result_exponent);
		*dstptr = result;
		return(NOEXCEPTION);
		}
	    Sgl_sethigh4bits(result,/*using*/sign_save);
	    switch(jumpsize)
		{
		case 1:
		    {
		    Sgl_leftshiftby3(result);
		    result_exponent -= 3;
		    break;
		    }
		case 2:
		case 3:
		    {
		    Sgl_leftshiftby2(result);
		    result_exponent -= 2;
		    break;
		    }
		case 4:
		case 5:
		case 6:
		case 7:
		    {
		    Sgl_leftshiftby1(result);
		    result_exponent -= 1;
		    break;
		    }
		}
	    if(result_exponent > 0)
		{
		Sgl_set_exponent(result,/*using*/result_exponent);
		*dstptr = result;	/* Sign bit is already set */
		return(NOEXCEPTION);
		}
	    /* Fixup potential underflows */
	  underflow:
	    if(Is_underflowtrap_enabled())
		{
		Sgl_set_sign(result,sign_save);
		Sgl_setwrapped_exponent(result,result_exponent,unfl);
		*dstptr = result;
		/* inexact = FALSE */
		return(UNDERFLOWEXCEPTION);
		}
	    /*
	     * Since we cannot get an inexact denormalized result,
	     * we can now return.
	     */
	    Sgl_right_align(result,/*by*/(1-result_exponent),extent);
	    Sgl_clear_signexponent(result);
	    Sgl_set_sign(result,sign_save);
	    *dstptr = result;
	    return(NOEXCEPTION);
	    } /* end if(hidden...)... */
	/* Fall through and round */
	} /* end if(save >= 0)... */
    else
	{
	/* Add magnitudes */
	Sgl_addition(left,right,/*to*/result);
	if(Sgl_isone_hiddenoverflow(result))
	    {
	    /* Prenormalization required. */
	    Sgl_rightshiftby1_withextent(result,extent,extent);
	    Sgl_arithrightshiftby1(result);
	    result_exponent++;
	    } /* end if hiddenoverflow... */
	} /* end else ...sub magnitudes... */

    /* Round the result.  If the extension is all zeros,then the result is
     * exact.  Otherwise round in the correct direction.  No underflow is
     * possible. If a postnormalization is necessary, then the mantissa is
     * all zeros so no shift is needed. */
  round:
    if(Ext_isnotzero(extent))
	{
	inexact = TRUE;
	switch(Rounding_mode())
	    {
	    case ROUNDNEAREST: /* The default. */
	    if(Ext_isone_sign(extent))
		{
		/* at least 1/2 ulp */
		if(Ext_isnotzero_lower(extent)  ||
		  Sgl_isone_lowmantissa(result))
		    {
		    /* either exactly half way and odd or more than 1/2ulp */
		    Sgl_increment(result);
		    }
		}
	    break;

	    case ROUNDPLUS:
	    if(Sgl_iszero_sign(result))
		{
		/* Round up positive results */
		Sgl_increment(result);
		}
	    break;

	    case ROUNDMINUS:
	    if(Sgl_isone_sign(result))
		{
		/* Round down negative results */
		Sgl_increment(result);
		}

	    case ROUNDZERO:;
	    /* truncate is simple */
	    } /* end switch... */
	if(Sgl_isone_hiddenoverflow(result)) result_exponent++;
	}
    if(result_exponent == SGL_INFINITY_EXPONENT)
	{
	/* Overflow */
	if(Is_overflowtrap_enabled())
	    {
	    Sgl_setwrapped_exponent(result,result_exponent,ovfl);
	    *dstptr = result;
	    if (inexact) {
		if (Is_inexacttrap_enabled())
		    return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
		else Set_inexactflag();
	    }
	    return(OVERFLOWEXCEPTION);
	    }
	else
	    {
	    Set_overflowflag();
	    inexact = TRUE;
	    Sgl_setoverflow(result);
	    }
	}
    else Sgl_set_exponent(result,result_exponent);
    *dstptr = result;
    if(inexact) {
	if(Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
	else Set_inexactflag();
    }
    return(NOEXCEPTION);
}