![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to support.s CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDTahoe / usr.lib / libm / national|
Return to support.s CVS log | Up to [CSRG BSD Unix] / 43BSDTahoe / usr.lib / libm / national |
1.1 root 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) 6/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
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 <= 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
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
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 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
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
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
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.