xref: /original-bsd/lib/libm/national/support.s (revision 2301fdfb) 
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
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.
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;
21;	@(#)support.s	5.3 (Berkeley) 06/30/88
22;
23
24; IEEE recommended functions
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
40	andd	0x80000000,r1
41	cmpqd	0,r1
42	beq	end
43	negl	f0,f0
44end:	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 <= x*2**-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
72L1:	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;
78L2:
79	movl	4(sp),f0	; logb(+inf)=+inf, logb(NaN)=NaN
80	ret	0
81;
82; x is 0 or subnormal
83;
84L3:
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;
102L5:
103	movl	0f1.0e300,f0
104	mull	0f-1.0e300,f0	; multiply two big numbers to get -INF
105	ret	0
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)
129L1:	addl	f2,f0		; f0 = x + s
130	subl	f2,f0		; f0 = f0 - s
131	ret	0
132itself:	movl	4(sp),f0
133	ret	0
134L52:	.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 x*2**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 x*2**N
199end:	ret	0
200inof:	bcs	underflow	; negative int overflow if Carry bit is set
201overflow:
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
209underflow:
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
219zero:	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
226L19:		; 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
242end:	ret	0
243L30:	.double	0x0,0x3bf00000	; floating point 2**-64
244L32:	.double	0x0,0x43f00000	; floating point 2**64
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;
301L10:
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
319LOOP:
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
3341:
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
350fnad:
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
360final:
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
3711:
372	subl	f2,f0		; x = x - y
3732:
374	cmpqd	0,-40(fp)
375	beq	3f
376	divl	-12(fp),f0	; scale down the answer
3773:
378	movl	f0,tos
379	xord	r6,tos
380	movl	tos,f0
381	movd	-16(fp),r0
382	lfsr	r0		; restore the fsr
383
384end:	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;
391L2:
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
399L4:	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;
404L3:
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;
413L5:
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)
423L6:
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;
433L1:
434	movl	8(fp),f0
435	subl	f0,f0		; if x is INF, then INF-INF is NaN
436	ret	0
437L64:	.double 0x0,0x43f00000	; L64 = 2**64
438LTINY:	.double 0x0,0x04000000	; LTINY = 2**-959
439