[unix-history] / usr / src / lib / libm / national / support.s

; Copyright (c) 1985 Regents of the University of California.
; All rights reserved.
;
; Redistribution and use in source and binary forms are permitted
; provided that the above copyright notice and this paragraph are
; duplicated in all such forms and that any documentation,
; advertising materials, and other materials related to such
; distribution and use acknowledge that the software was developed
; by the University of California, Berkeley.  The name of the
; University may not be used to endorse or promote products derived
; from this software without specific prior written permission.
; THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
; IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
; WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
;
; All recipients should regard themselves as participants in an ongoing
; research project and hence should feel obligated to report their
; experiences (good or bad) with these elementary function codes, using
; the sendbug(8) program, to the authors.
;
;	@(#)support.s	5.3 (Berkeley) %G%
;

; IEEE recommended functions
; 
; double copysign(x,y)
; double x,y;
; IEEE 754 recommended function, return x*sign(y)
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
;
	.vers	2
	.text
	.align	2
	.globl	_copysign
_copysign:
	movl	4(sp),f0
	movd	8(sp),r0
	movd	16(sp),r1
	xord	r0,r1
	andd	0x80000000,r1
	cmpqd	0,r1
	beq	end
	negl	f0,f0
end:	ret	0

; 
; double logb(x)
; double x;
; IEEE p854 recommended function, return the exponent of x (return float(N) 
; such that 1 <= x*2**-N < 2, even for subnormal number.
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; Note: subnormal number (if implemented) will be taken care of. 
;
	.vers	2
	.text
	.align	2
	.globl	_logb
_logb:
;
; extract the exponent of x
; glossaries:	r0 = high part of x
;		r1 = unbias exponent of x
;		r2 = 20 (first exponent bit position)
;
	movd	8(sp),r0
	movd	20,r2
	extd	r2,r0,r1,11	; extract the exponent of x
	cmpqd	0,r1		; if exponent bits = 0, goto L3
	beq	L3
	cmpd	0x7ff,r1
	beq	L2		; if exponent bits = 0x7ff, goto L2
L1:	subd	1023,r1		; unbias the exponent
	movdl	r1,f0		; convert the exponent to floating value
	ret	0
;
; x is INF or NaN, simply return x
;
L2:
	movl	4(sp),f0	; logb(+inf)=+inf, logb(NaN)=NaN
	ret	0
;
; x is 0 or subnormal
;
L3:
	movl	4(sp),f0
	cmpl	0f0,f0
	beq	L5		; x is 0 , goto L5 (return -inf)
;
; Now x is subnormal
;
	mull	L64,f0		; scale up f0 with 2**64
	movl	f0,tos
	movd	tos,r0
	movd	tos,r0		; now r0 = new high part of x
	extd	r2,r0,r1,11	; extract the exponent of x to r1
	subd	1087,r1		; unbias the exponent with correction 
	movdl	r1,f0		; convert the exponent to floating value
	ret	0
;
; x is 0, return logb(0)= -INF
;
L5:
	movl	0f1.0e300,f0
	mull	0f-1.0e300,f0	; multiply two big numbers to get -INF
	ret	0
; 
; double rint(x)
; double x;
; ... delivers integer nearest x in direction of prevailing rounding
; ... mode
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; Note: subnormal number (if implemented) will be taken care of. 
;
	.vers	2
	.text
	.align	2
	.globl	_rint
_rint:
;
	movd	8(sp),r0
	movd	20,r2
	extd	r2,r0,r1,11	; extract the exponent of x
	cmpd	0x433,r1
	ble	itself
	movl	L52,f2		; f2 = L = 2**52
	cmpqd	0,r0
	ble	L1
	negl	f2,f2		; f2 = s = copysign(L,x)
L1:	addl	f2,f0		; f0 = x + s
	subl	f2,f0		; f0 = f0 - s
	ret	0
itself:	movl	4(sp),f0
	ret	0
L52:	.double	0x0,0x43300000	; L52=2**52
; 
; int finite(x)
; double x;
; IEEE 754 recommended function, return 0 if x is NaN or INF, else 0
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
;
	.vers	2
	.text
	.align	2
	.globl	_finite
_finite:
	movd	4(sp),r1
	andd	0x800fffff,r1
	cmpd	0x7ff00000,r1
	sned	r0		; r0=0 if exponent(x) = 0x7ff
	ret	0
; 
; double scalb(x,N)
; double x; int N;
; IEEE 754 recommended function, return x*2**N by adjusting 
; exponent of x.
; Coded by K.C. Ng in National 32k assembler, 11/9/85. 
; Note: subnormal number (if implemented) will be taken care of 
;
	.vers	2
	.text
	.align	2
	.globl	_scalb
_scalb:
;
; if x=0 return 0
;
	movl	4(sp),f0
	cmpl	0f0,f0
	beq	end		; scalb(0,N) is x itself
;
; extract the exponent of x
; glossaries:	r0 = high part of x, 
;		r1 = unbias exponent of x,
;		r2 = 20 (first exponent bit position).
;
	movd	8(sp),r0	; r0 = high part of x
	movd	20,r2		; r2 = 20
	extd	r2,r0,r1,11	; extract the exponent of x in r1
	cmpd	0x7ff,r1	
;
; if exponent of x is 0x7ff, then x is NaN or INF; simply return x  
;
	beq	end		
	cmpqd	0,r1
;
; if exponent of x is zero, then x is subnormal; goto L19
;
	beq	L19		
	addd	12(sp),r1	; r1 = (exponent of x) + N
	bfs	inof		; if integer overflows, goto inof
	cmpqd	0,r1		; if new exponent <= 0, goto underflow
	bge	underflow
	cmpd	2047,r1		; if new exponent >= 2047 goto overflow
	ble	overflow
	insd	r2,r1,r0,11	; insert the new exponent 
	movd	r0,tos
	movd	8(sp),tos
	movl	tos,f0		; return x*2**N
end:	ret	0
inof:	bcs	underflow	; negative int overflow if Carry bit is set
overflow:
	andd	0x80000000,r0	; keep the sign of x
	ord	0x7fe00000,r0	; set x to a huge number
	movd	r0,tos
	movqd	0,tos
	movl	tos,f0
	mull	0f1.0e300,f0	; multiply two huge number to get overflow
	ret	0
underflow:
	addd	64,r1		; add 64 to exonent to see if it is subnormal
	cmpqd	0,r1
	bge	zero		; underflow to zero
	insd	r2,r1,r0,11	; insert the new exponent 
	movd	r0,tos
	movd	8(sp),tos
	movl	tos,f0
	mull	L30,f0		; readjust x by multiply it with 2**-64
	ret	0
zero:	andd	0x80000000,r0	; keep the sign of x
	ord	0x00100000,r0	; set x to a tiny number
	movd	r0,tos
	movqd	0,tos
	movl	tos,f0
	mull	0f1.0e-300,f0	; underflow to 0  by multipling two tiny nos.
	ret	0
L19:		; subnormal number
	mull	L32,f0		; scale up x by 2**64
	movl	f0,tos
	movd	tos,r0
	movd	tos,r0		; get the high part of new x
	extd	r2,r0,r1,11	; extract the exponent of x in r1
	addd	12(sp),r1	; exponent of x + N
	subd	64,r1		; adjust it by subtracting 64
	cmpqd	0,r1
	bge	underflow
	cmpd	2047,r1
	ble	overflow
	insd	r2,r1,r0,11	; insert back the incremented exponent 
	movd	r0,tos
	movd	8(sp),tos
	movl	tos,f0
end:	ret	0
L30:	.double	0x0,0x3bf00000	; floating point 2**-64
L32:	.double	0x0,0x43f00000	; floating point 2**64
; 
; double drem(x,y)
; double x,y;
; IEEE double remainder function, return x-n*y, where n=x/y rounded to 
; nearest integer (half way case goes to even). Result exact.
; Coded by K.C. Ng in National 32k assembly, 11/19/85.
;
	.vers	2
	.text
	.align	2
	.globl	_drem
_drem:
;
; glossaries:	
;		r2 = high part of x
;		r3 = exponent of x
;		r4 = high part of y
;		r5 = exponent of y
;		r6 = sign of x
;		r7 = constant 0x7ff00000
;
;  16(fp) : y
;   8(fp) : x
; -12(fp) : adjustment on y when y is subnormal
; -16(fp) : fsr
; -20(fp) : nx
; -28(fp) : t
; -36(fp) : t1
; -40(fp) : nf
; 
;
	enter	[r3,r4,r5,r6,r7],40
	movl	f6,tos
	movl	f4,tos
	movl	0f0,-12(fp)
	movd	0,-20(fp)
	movd	0,-40(fp)
	movd	0x7ff00000,r7	; initialize r7=0x7ff00000
	movd	12(fp),r2	; r2 = high(x)
	movd	r2,r3
	andd	r7,r3		; r3 = xexp
	cmpd	r7,r3
; if x is NaN or INF goto L1
	beq	L1		
	movd	20(fp),r4
	bicd	[31],r4		; r4 = high part of |y|
	movd	r4,20(fp)	; y = |y|
	movd	r4,r5
	andd	r7,r5		; r5 = yexp
	cmpd	r7,r5
	beq	L2		; if y is NaN or INF goto L2
	cmpd	0x04000000,r5	; 
	bgt	L3		; if y is tiny goto L3
;
; now y != 0 , x is finite
;
L10:
	movd	r2,r6
	andd	0x80000000,r6	; r6 = sign(x)
	bicd	[31],r2		; x <- |x|
	sfsr	r1
	movd	r1,-16(fp)	; save fsr in -16(fp)
	bicd	[5],r1
	lfsr	r1		; disable inexact interupt
	movd	16(fp),r0	; r0 = low part of y
	movd	r0,r1		; r1 = r0 = low part of y
	andd	0xf8000000,r1	; mask off the lsb 27 bits of y

	movd	r2,12(fp)	; update x to |x|
	movd	r0,-28(fp)	; 
	movd	r4,-24(fp)	; t  = y
	movd	r4,-32(fp)	; 
	movd	r1,-36(fp)	; t1 = y with trialing 27 zeros
	movd	0x01900000,r1	; r1 = 25 in exponent field
LOOP:
	movl	8(fp),f0	; f0 = x
	movl	16(fp),f2	; f2 = y
	cmpl	f0,f2
	ble	fnad		; goto fnad (final adjustment) if x <= y
	movd	r4,-32(fp)
	movd	r3,r0
	subd	r5,r0		; xexp - yexp
	subd	r1,r0		; r0 = xexp - yexp - m25
	cmpqd	0,r0		; r0 > 0 ?
	bge	1f
	addd	r4,r0		; scale up (high) y
	movd	r0,-24(fp)	; scale up t
	movl	-28(fp),f2	; t
	movd	r0,-32(fp)	; scale up t1
1:
	movl	-36(fp),f4	; t1
	movl	f0,f6
	divl	f2,f6		; f6 = x/t
	floorld	f6,r0		; r0 = [x/t]
	movdl	r0,f6		; f6 = n
	subl	f4,f2		; t = t - t1 (tail of t1)
	mull	f6,f4		; f4 = n*t1	...exact
	subl	f4,f0		; x = x - n*t1
	mull	f6,f2		; n*(t-t1)	...exact
	subl	f2,f0		; x = x - n*(t-t1)
; update xexp
	movl	f0,8(fp)
	movd	12(fp),r3	
	andd	r7,r3
	jump	LOOP
fnad:
	cmpqd	0,-20(fp)	; 0 = nx?
	beq	final
	mull	-12(fp),8(fp)	; scale up x the same amount as y
	movd	0,-20(fp)
	movd	12(fp),r2
	movd	r2,r3
	andd	r7,r3		; update exponent of x
	jump	LOOP

final:
	movl	16(fp),f2	; f2 = y (f0=x, r0=n)
	subd	0x100000,r4	; high y /2
	movd	r4,-24(fp)
	movl	-28(fp),f4	; f4 = y/2
	cmpl	f0,f4		; x > y/2 ?
	bgt	1f
	bne	2f
	andd	1,r0		; n is odd or even
	cmpqd	0,r0
	beq	2f
1:
	subl	f2,f0		; x = x - y
2:
	cmpqd	0,-40(fp)
	beq	3f
	divl	-12(fp),f0	; scale down the answer
3:
	movl	f0,tos
	xord	r6,tos
	movl	tos,f0
	movd	-16(fp),r0
	lfsr	r0		; restore the fsr

end:	movl	tos,f4
	movl	tos,f6
	exit	[r3,r4,r5,r6,r7]
	ret	0
;
; y is NaN or INF
;
L2:	
	movd	16(fp),r0	; r0 = low part of y
	andd	0xfffff,r4	; r4 = high part of y & 0x000fffff
	ord	r4,r0
	cmpqd	0,r0
	beq	L4
	movl	16(fp),f0	; y is NaN, return y
	jump	end
L4:	movl	8(fp),f0	; y is inf, return x
	jump	end
;
; exponent of y is less than 64, y may be zero or subnormal
;
L3:
	movl	16(fp),f0
	cmpl	0f0,f0
	bne	L5
	divl	f0,f0		; y is 0, return NaN by doing 0/0
	jump	end
;
; subnormal y or tiny y
;
L5:	
	movd	0x04000000,-20(fp)	; nx = 64 in exponent field
	movl	L64,f2
	movl	f2,-12(fp)
	mull	f2,f0
	cmpl	f0,LTINY
	bgt	L6
	mull	f2,f0
	addd	0x04000000,-20(fp)	; nx = nx + 64 in exponent field
	mull	f2,-12(fp)
L6:
	movd	-20(fp),-40(fp)
	movl	f0,16(fp)
	movd	20(fp),r4
	movd	r4,r5
	andd	r7,r5		; exponent of new y
	jump	L10
;
; x is NaN or INF, return x-x
;
L1:
	movl	8(fp),f0
	subl	f0,f0		; if x is INF, then INF-INF is NaN
	ret	0
L64:	.double 0x0,0x43f00000	; L64 = 2**64
LTINY:	.double 0x0,0x04000000	; LTINY = 2**-959
Commit	Line	Data
bc0f2b0e KB	1	; Copyright (c) 1985 Regents of the University of California.
	2	; All rights reserved.
	3	;
	4	; Redistribution and use in source and binary forms are permitted
a399f6c8 KB	5	; provided that the above copyright notice and this paragraph are
	6	; duplicated in all such forms and that any documentation,
	7	; advertising materials, and other materials related to such
	8	; distribution and use acknowledge that the software was developed
	9	; by the University of California, Berkeley. The name of the
	10	; University may not be used to endorse or promote products derived
	11	; from this software without specific prior written permission.
	12	; THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
	13	; IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
	14	; WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
bc0f2b0e KB	15	;
	16	; All recipients should regard themselves as participants in an ongoing
	17	; research project and hence should feel obligated to report their
	18	; experiences (good or bad) with these elementary function codes, using
	19	; the sendbug(8) program, to the authors.
	20	;
a399f6c8	21	; @(#)support.s 5.3 (Berkeley) %G%
bc0f2b0e KB	22	;
bc0f2b0e KB	23
ab5cf571	24	; IEEE recommended functions
ab5cf571 GK	25	;
	26	; double copysign(x,y)
	27	; double x,y;
	28	; IEEE 754 recommended function, return x*sign(y)
	29	; Coded by K.C. Ng in National 32k assembler, 11/9/85.
	30	;
	31	.vers 2
	32	.text
	33	.align 2
	34	.globl _copysign
	35	_copysign:
	36	movl 4(sp),f0
	37	movd 8(sp),r0
	38	movd 16(sp),r1
	39	xord r0,r1
599b6e68	40	andd 0x80000000,r1
ab5cf571 GK	41	cmpqd 0,r1
	42	beq end
	43	negl f0,f0
	44	end: ret 0
	45
	46	;
	47	; double logb(x)
	48	; double x;
	49	; IEEE p854 recommended function, return the exponent of x (return float(N)
	50	; such that 1 <= x2*-N < 2, even for subnormal number.
	51	; Coded by K.C. Ng in National 32k assembler, 11/9/85.
	52	; Note: subnormal number (if implemented) will be taken care of.
	53	;
	54	.vers 2
	55	.text
	56	.align 2
	57	.globl _logb
	58	_logb:
	59	;
	60	; extract the exponent of x
	61	; glossaries: r0 = high part of x
	62	; r1 = unbias exponent of x
	63	; r2 = 20 (first exponent bit position)
	64	;
	65	movd 8(sp),r0
	66	movd 20,r2
	67	extd r2,r0,r1,11 ; extract the exponent of x
	68	cmpqd 0,r1 ; if exponent bits = 0, goto L3
	69	beq L3
	70	cmpd 0x7ff,r1
	71	beq L2 ; if exponent bits = 0x7ff, goto L2
	72	L1: subd 1023,r1 ; unbias the exponent
	73	movdl r1,f0 ; convert the exponent to floating value
	74	ret 0
	75	;
	76	; x is INF or NaN, simply return x
	77	;
	78	L2:
	79	movl 4(sp),f0 ; logb(+inf)=+inf, logb(NaN)=NaN
	80	ret 0
	81	;
	82	; x is 0 or subnormal
	83	;
	84	L3:
	85	movl 4(sp),f0
	86	cmpl 0f0,f0
	87	beq L5 ; x is 0 , goto L5 (return -inf)
	88	;
	89	; Now x is subnormal
	90	;
	91	mull L64,f0 ; scale up f0 with 2**64
	92	movl f0,tos
	93	movd tos,r0
	94	movd tos,r0 ; now r0 = new high part of x
	95	extd r2,r0,r1,11 ; extract the exponent of x to r1
	96	subd 1087,r1 ; unbias the exponent with correction
	97	movdl r1,f0 ; convert the exponent to floating value
	98	ret 0
	99	;
	100	; x is 0, return logb(0)= -INF
	101	;
	102	L5:
	103	movl 0f1.0e300,f0
	104	mull 0f-1.0e300,f0 ; multiply two big numbers to get -INF
105	ret 0
ab5cf571 GK	106	;
	107	; double rint(x)
	108	; double x;
	109	; ... delivers integer nearest x in direction of prevailing rounding
	110	; ... mode
	111	; Coded by K.C. Ng in National 32k assembler, 11/9/85.
	112	; Note: subnormal number (if implemented) will be taken care of.
	113	;
	114	.vers 2
	115	.text
	116	.align 2
	117	.globl _rint
	118	_rint:
	119	;
	120	movd 8(sp),r0
	121	movd 20,r2
	122	extd r2,r0,r1,11 ; extract the exponent of x
	123	cmpd 0x433,r1
	124	ble itself
	125	movl L52,f2 ; f2 = L = 2**52
	126	cmpqd 0,r0
	127	ble L1
	128	negl f2,f2 ; f2 = s = copysign(L,x)
	129	L1: addl f2,f0 ; f0 = x + s
	130	subl f2,f0 ; f0 = f0 - s
	131	ret 0
	132	itself: movl 4(sp),f0
	133	ret 0
	134	L52: .double 0x0,0x43300000 ; L52=2**52
	135	;
	136	; int finite(x)
	137	; double x;
	138	; IEEE 754 recommended function, return 0 if x is NaN or INF, else 0
	139	; Coded by K.C. Ng in National 32k assembler, 11/9/85.
	140	;
	141	.vers 2
	142	.text
	143	.align 2
	144	.globl _finite
	145	_finite:
	146	movd 4(sp),r1
	147	andd 0x800fffff,r1
	148	cmpd 0x7ff00000,r1
	149	sned r0 ; r0=0 if exponent(x) = 0x7ff
	150	ret 0
	151	;
	152	; double scalb(x,N)
	153	; double x; int N;
	154	; IEEE 754 recommended function, return x2*N by adjusting
	155	; exponent of x.
	156	; Coded by K.C. Ng in National 32k assembler, 11/9/85.
	157	; Note: subnormal number (if implemented) will be taken care of
	158	;
	159	.vers 2
	160	.text
	161	.align 2
	162	.globl _scalb
	163	_scalb:
	164	;
	165	; if x=0 return 0
	166	;
	167	movl 4(sp),f0
	168	cmpl 0f0,f0
	169	beq end ; scalb(0,N) is x itself
170	;
171	; extract the exponent of x
172	; glossaries: r0 = high part of x,
173	; r1 = unbias exponent of x,
174	; r2 = 20 (first exponent bit position).
175	;
176	movd 8(sp),r0 ; r0 = high part of x
177	movd 20,r2 ; r2 = 20
178	extd r2,r0,r1,11 ; extract the exponent of x in r1
179	cmpd 0x7ff,r1
180	;
181	; if exponent of x is 0x7ff, then x is NaN or INF; simply return x
182	;
183	beq end
184	cmpqd 0,r1
185	;
186	; if exponent of x is zero, then x is subnormal; goto L19
187	;
188	beq L19
189	addd 12(sp),r1 ; r1 = (exponent of x) + N
190	bfs inof ; if integer overflows, goto inof
191	cmpqd 0,r1 ; if new exponent <= 0, goto underflow
192	bge underflow
193	cmpd 2047,r1 ; if new exponent >= 2047 goto overflow
194	ble overflow
195	insd r2,r1,r0,11 ; insert the new exponent
196	movd r0,tos
197	movd 8(sp),tos
198	movl tos,f0 ; return x2*N
199	end: ret 0
200	inof: bcs underflow ; negative int overflow if Carry bit is set
201	overflow:
202	andd 0x80000000,r0 ; keep the sign of x
203	ord 0x7fe00000,r0 ; set x to a huge number
204	movd r0,tos
205	movqd 0,tos
206	movl tos,f0
207	mull 0f1.0e300,f0 ; multiply two huge number to get overflow
208	ret 0
209	underflow:
210	addd 64,r1 ; add 64 to exonent to see if it is subnormal
211	cmpqd 0,r1
212	bge zero ; underflow to zero
213	insd r2,r1,r0,11 ; insert the new exponent
214	movd r0,tos
215	movd 8(sp),tos
216	movl tos,f0
217	mull L30,f0 ; readjust x by multiply it with 2**-64
218	ret 0
219	zero: andd 0x80000000,r0 ; keep the sign of x
220	ord 0x00100000,r0 ; set x to a tiny number
221	movd r0,tos
222	movqd 0,tos
223	movl tos,f0
224	mull 0f1.0e-300,f0 ; underflow to 0 by multipling two tiny nos.
225	ret 0
226	L19: ; subnormal number
227	mull L32,f0 ; scale up x by 2**64
228	movl f0,tos
229	movd tos,r0
230	movd tos,r0 ; get the high part of new x
231	extd r2,r0,r1,11 ; extract the exponent of x in r1
232	addd 12(sp),r1 ; exponent of x + N
233	subd 64,r1 ; adjust it by subtracting 64
234	cmpqd 0,r1
235	bge underflow
236	cmpd 2047,r1
237	ble overflow
238	insd r2,r1,r0,11 ; insert back the incremented exponent
239	movd r0,tos
240	movd 8(sp),tos
241	movl tos,f0
242	end: ret 0
243	L30: .double 0x0,0x3bf00000 ; floating point 2**-64
244	L32: .double 0x0,0x43f00000 ; floating point 2**64
dc4bce59 ZAL	245	;
	246	; double drem(x,y)
	247	; double x,y;
	248	; IEEE double remainder function, return x-n*y, where n=x/y rounded to
	249	; nearest integer (half way case goes to even). Result exact.
	250	; Coded by K.C. Ng in National 32k assembly, 11/19/85.
	251	;
	252	.vers 2
	253	.text
	254	.align 2
	255	.globl _drem
	256	_drem:
	257	;
	258	; glossaries:
	259	; r2 = high part of x
	260	; r3 = exponent of x
	261	; r4 = high part of y
	262	; r5 = exponent of y
	263	; r6 = sign of x
	264	; r7 = constant 0x7ff00000
	265	;
	266	; 16(fp) : y
	267	; 8(fp) : x
	268	; -12(fp) : adjustment on y when y is subnormal
	269	; -16(fp) : fsr
	270	; -20(fp) : nx
	271	; -28(fp) : t
	272	; -36(fp) : t1
	273	; -40(fp) : nf
	274	;
	275	;
	276	enter [r3,r4,r5,r6,r7],40
	277	movl f6,tos
	278	movl f4,tos
	279	movl 0f0,-12(fp)
	280	movd 0,-20(fp)
	281	movd 0,-40(fp)
	282	movd 0x7ff00000,r7 ; initialize r7=0x7ff00000
	283	movd 12(fp),r2 ; r2 = high(x)
	284	movd r2,r3
	285	andd r7,r3 ; r3 = xexp
	286	cmpd r7,r3
	287	; if x is NaN or INF goto L1
	288	beq L1
	289	movd 20(fp),r4
	290	bicd [31],r4 ; r4 = high part of \|y\|
	291	movd r4,20(fp) ; y = \|y\|
	292	movd r4,r5
	293	andd r7,r5 ; r5 = yexp
	294	cmpd r7,r5
	295	beq L2 ; if y is NaN or INF goto L2
	296	cmpd 0x04000000,r5 ;
	297	bgt L3 ; if y is tiny goto L3
	298	;
	299	; now y != 0 , x is finite
	300	;
	301	L10:
	302	movd r2,r6
	303	andd 0x80000000,r6 ; r6 = sign(x)
	304	bicd [31],r2 ; x <- \|x\|
	305	sfsr r1
	306	movd r1,-16(fp) ; save fsr in -16(fp)
	307	bicd [5],r1
	308	lfsr r1 ; disable inexact interupt
309	movd 16(fp),r0 ; r0 = low part of y
310	movd r0,r1 ; r1 = r0 = low part of y
311	andd 0xf8000000,r1 ; mask off the lsb 27 bits of y
312
313	movd r2,12(fp) ; update x to \|x\|
314	movd r0,-28(fp) ;
315	movd r4,-24(fp) ; t = y
316	movd r4,-32(fp) ;
317	movd r1,-36(fp) ; t1 = y with trialing 27 zeros
318	movd 0x01900000,r1 ; r1 = 25 in exponent field
319	LOOP:
320	movl 8(fp),f0 ; f0 = x
321	movl 16(fp),f2 ; f2 = y
322	cmpl f0,f2
323	ble fnad ; goto fnad (final adjustment) if x <= y
324	movd r4,-32(fp)
325	movd r3,r0
326	subd r5,r0 ; xexp - yexp
327	subd r1,r0 ; r0 = xexp - yexp - m25
328	cmpqd 0,r0 ; r0 > 0 ?
329	bge 1f
330	addd r4,r0 ; scale up (high) y
331	movd r0,-24(fp) ; scale up t
332	movl -28(fp),f2 ; t
333	movd r0,-32(fp) ; scale up t1
334	1:
335	movl -36(fp),f4 ; t1
336	movl f0,f6
337	divl f2,f6 ; f6 = x/t
338	floorld f6,r0 ; r0 = [x/t]
339	movdl r0,f6 ; f6 = n
340	subl f4,f2 ; t = t - t1 (tail of t1)
341	mull f6,f4 ; f4 = n*t1 ...exact
342	subl f4,f0 ; x = x - n*t1
343	mull f6,f2 ; n*(t-t1) ...exact
344	subl f2,f0 ; x = x - n*(t-t1)
345	; update xexp
346	movl f0,8(fp)
347	movd 12(fp),r3
348	andd r7,r3
349	jump LOOP
350	fnad:
351	cmpqd 0,-20(fp) ; 0 = nx?
352	beq final
353	mull -12(fp),8(fp) ; scale up x the same amount as y
354	movd 0,-20(fp)
355	movd 12(fp),r2
356	movd r2,r3
357	andd r7,r3 ; update exponent of x
358	jump LOOP
359
360	final:
361	movl 16(fp),f2 ; f2 = y (f0=x, r0=n)
362	subd 0x100000,r4 ; high y /2
363	movd r4,-24(fp)
364	movl -28(fp),f4 ; f4 = y/2
365	cmpl f0,f4 ; x > y/2 ?
366	bgt 1f
367	bne 2f
368	andd 1,r0 ; n is odd or even
369	cmpqd 0,r0
370	beq 2f
371	1:
372	subl f2,f0 ; x = x - y
373	2:
374	cmpqd 0,-40(fp)
375	beq 3f
376	divl -12(fp),f0 ; scale down the answer
377	3:
378	movl f0,tos
379	xord r6,tos
380	movl tos,f0
381	movd -16(fp),r0
382	lfsr r0 ; restore the fsr
383
384	end: movl tos,f4
385	movl tos,f6
386	exit [r3,r4,r5,r6,r7]
387	ret 0
388	;
389	; y is NaN or INF
390	;
391	L2:
392	movd 16(fp),r0 ; r0 = low part of y
393	andd 0xfffff,r4 ; r4 = high part of y & 0x000fffff
394	ord r4,r0
395	cmpqd 0,r0
396	beq L4
397	movl 16(fp),f0 ; y is NaN, return y
398	jump end
399	L4: movl 8(fp),f0 ; y is inf, return x
400	jump end
401	;
402	; exponent of y is less than 64, y may be zero or subnormal
403	;
404	L3:
405	movl 16(fp),f0
406	cmpl 0f0,f0
407	bne L5
408	divl f0,f0 ; y is 0, return NaN by doing 0/0
409	jump end
410	;
411	; subnormal y or tiny y
412	;
413	L5:
414	movd 0x04000000,-20(fp) ; nx = 64 in exponent field
415	movl L64,f2
416	movl f2,-12(fp)
417	mull f2,f0
418	cmpl f0,LTINY
419	bgt L6
420	mull f2,f0
421	addd 0x04000000,-20(fp) ; nx = nx + 64 in exponent field
422	mull f2,-12(fp)
423	L6:
424	movd -20(fp),-40(fp)
425	movl f0,16(fp)
426	movd 20(fp),r4
427	movd r4,r5
428	andd r7,r5 ; exponent of new y
429	jump L10
430	;
431	; x is NaN or INF, return x-x
432	;
433	L1:
434	movl 8(fp),f0
435	subl f0,f0 ; if x is INF, then INF-INF is NaN
436	ret 0
437	L64: .double 0x0,0x43f00000 ; L64 = 2**64
438	LTINY: .double 0x0,0x04000000 ; LTINY = 2**-959