|
|
1.1 root 1: /*
2: * Copyright (c) 1987 Regents of the University of California.
3: * All rights reserved.
4: *
5: * Redistribution and use in source and binary forms are permitted
6: * provided that the above copyright notice and this paragraph are
7: * duplicated in all such forms and that any documentation,
8: * advertising materials, and other materials related to such
9: * distribution and use acknowledge that the software was developed
10: * by the University of California, Berkeley. The name of the
11: * University may not be used to endorse or promote products derived
12: * from this software without specific prior written permission.
13: * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
14: * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
15: * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
16: *
17: * All recipients should regard themselves as participants in an ongoing
18: * research project and hence should feel obligated to report their
19: * experiences (good or bad) with these elementary function codes, using
20: * the sendbug(8) program, to the authors.
21: */
22: .data
23: .align 2
24: _sccsid:
25: .asciz "@(#)support.s 5.4 (ucb.elefunt) 6/30/88"
26: /*
27: * copysign(x,y),
28: * logb(x),
29: * scalb(x,N),
30: * finite(x),
31: * drem(x,y),
32: * Coded in vax assembly language by K. C. Ng 4/9/85.
33: * Re-coded in tahoe assembly language by Z. Alex Liu 7/13/87.
34: */
35: /*
36: * double copysign(x,y)
37: * double x,y;
38: */
39: .globl _copysign
40: .text
41: .align 2
42: _copysign:
43: .word 0x0004 # save r2
44: movl 8(fp),r1
45: movl 4(fp),r0 # r0:r1 = x
46: andl3 $0x7f800000,r0,r2 # r2 = biased exponent of x
47: beql 1f # if 0 or reserved op then return x
48: andl3 $0x80000000,12(fp),r2 # r2 = sign bit of y at bit-31
49: andl2 $0x7fffffff,r0 # replace x by |x|
50: orl2 r2,r0 # copy the sign bit of y to x
51: 1: ret
52: /*
53: * double logb(x)
54: * double x;
55: */
56: .globl _logb
57: .text
58: .align 2
59: _logb:
60: .word 0x0000 # save nothing
61: andl3 $0x7f800000,4(fp),r0 # r0[b23:b30] = biased exponent of x
62: beql 1f
63: shrl $23,r0,r0 # r0[b0:b7] = biased exponent of x
64: subl2 $129,r0 # r0 = unbiased exponent of x
65: cvld r0 # acc = unbiased exponent of x (double)
66: std r0 # r0 = unbiased exponent of x (double)
67: ret
68: 1: movl 8(fp),r1 # 8(fp) must be moved first
69: movl 4(fp),r0 # r0:r1 = x (zero or reserved op)
70: blss 2f # simply return if reserved op
71: movl $0xfe000000,r1
72: movl $0xcfffffff,r0 # -2147483647.0
73: 2: ret
74: /*
75: * long finite(x)
76: * double x;
77: */
78: .globl _finite
79: .text
80: .align 2
81: _finite:
82: .word 0x0000 # save nothing
83: andl3 $0xff800000,4(fp),r0 # r0 = sign of x & its biased exponent
84: cmpl r0,$0x80000000 # is x a reserved op?
85: beql 1f # if so, return FALSE (0)
86: movl $1,r0 # else return TRUE (1)
87: ret
88: 1: clrl r0
89: ret
90: /*
91: * double scalb(x,N)
92: * double x; int N;
93: */
94: .globl _scalb
95: .set ERANGE,34
96: .text
97: .align 2
98: _scalb:
99: .word 0x000c # save r2-r3
100: movl 8(fp),r1
101: movl 4(fp),r0 # r0:r1 = x (-128 <= Ex <= 126)
102: andl3 $0x7f800000,r0,r3 # r3[b23:b30] = biased exponent of x
103: beql 1f # is x a 0 or a reserved operand?
104: movl 12(fp),r2 # r2 = N
105: cmpl r2,$0xff # if N >= 255
106: bgeq 2f # then the result must overflow
107: cmpl r2,$-0xff # if N <= -255
108: bleq 3f # then the result must underflow
109: shrl $23,r3,r3 # r3[b0:b7] = biased exponent of x
110: addl2 r2,r3 # r3 = biased exponent of the result
111: bleq 3f # if <= 0 then the result underflows
112: cmpl r3,$0x100 # if >= 256 then the result overflows
113: bgeq 2f
114: shll $23,r3,r3 # r3[b23:b30] = biased exponent of res.
115: andl2 $0x807fffff,r0
116: orl2 r3,r0 # r0:r1 = x*2^N
117: 1: ret
118: 2: pushl $ERANGE # if the result would overflow
119: callf $8,_infnan # and _infnan returns
120: andl3 $0x80000000,4(fp),r2 # get the sign of input arg
121: orl2 r2,r0 # re-attach the sign to r0:r1
122: ret
123: 3: clrl r1 # if the result would underflow
124: clrl r0 # then return 0
125: ret
126: /*
127: * double drem(x,y)
128: * double x,y;
129: * Returns x-n*y where n=[x/y] rounded (to even in the half way case).
130: */
131: .globl _drem
132: .set EDOM,33
133: .text
134: .align 2
135: _drem:
136: .word 0x1ffc # save r2-r12
137: movl 16(fp),r3
138: movl 12(fp),r2 # r2:r3 = y
139: movl 8(fp),r1
140: movl 4(fp),r0 # r0:r1 = x
141: andl3 $0xff800000,r0,r4
142: cmpl r4,$0x80000000 # is x a reserved operand?
143: beql 1f # if yes then propagate x and return
144: andl3 $0xff800000,r2,r4
145: cmpl r4,$0x80000000 # is y a reserved operand?
146: bneq 2f
147: movl r3,r1
148: movl r2,r0 # if yes then propagate y and return
149: 1: ret
150:
151: 2: tstl r4 # is y a 0?
152: bneq 3f
153: pushl $EDOM # if so then generate reserved op fault
154: callf $8,_infnan
155: ret
156:
157: 3: andl2 $0x7fffffff,r2 # r2:r3 = y <- |y|
158: clrl r12 # r12 = nx := 0
159: cmpl r2,$0x1c800000 # Ey ? 57
160: bgtr 4f # if Ey > 57 goto 4
161: addl2 $0x1c800000,r2 # scale up y by 2**57
162: movl $0x1c800000,r12 # r12[b23:b30] = nx = 57
163: 4: pushl r12 # pushed onto stack: nf := nx
164: andl3 $0x80000000,r0,-(sp) # pushed onto stack: sign of x
165: andl2 $0x7fffffff,r0 # r0:r1 = x <- |x|
166: movl r3,r11 # r10:r11 = y1 = y w/ last 27 bits 0
167: andl3 $0xf8000000,r10,r11 # clear last 27 bits of y1
168:
169: Loop: cmpd2 r0,r2 # x ? y
170: bleq 6f # if x <= y goto 6
171: /* # begin argument reduction */
172: movl r3,r5
173: movl r2,r4 # r4:r5 = t = y
174: movl r11,r7
175: movl r10,r6 # r6:r7 = t1 = y1
176: andl3 $0x7f800000,r0,r8 # r8[b23:b30] = Ex:biased exponent of x
177: andl3 $0x7f800000,r2,r9 # r9[b23:b30] = Ey:biased exponent of y
178: subl2 r9,r8 # r8[b23:b30] = Ex-Ey
179: subl2 $0x0c800000,r8 # r8[b23:b30] = k = Ex-Ey-25
180: blss 5f # if k < 0 goto 5
181: addl2 r8,r4 # t += k
182: addl2 r8,r6 # t1 += k, scale up t and t1
183: 5: ldd r0 # acc = x
184: divd r4 # acc = x/t
185: cvdl r8 # r8 = n = [x/t] truncated
186: cvld r8 # acc = dble(n)
187: std r8 # r8:r9 = dble(n)
188: ldd r4 # acc = t
189: subd r6 # acc = t-t1
190: muld r8 # acc = n*(t-t1)
191: std r4 # r4:r5 = n*(t-t1)
192: ldd r6 # acc = t1
193: muld r8 # acc = n*t1
194: subd r0 # acc = n*t1-x
195: negd # acc = x-n*t1
196: subd r4 # acc = (x-n*t1)-n*(t-t1)
197: std r0 # r0:r1 = (x-n*t1)-n*(t-t1)
198: brb Loop
199:
200: 6: movl r12,r6 # r6 = nx
201: beql 7f # if nx == 0 goto 7
202: addl2 r6,r0 # x <- x*2**57:scale x up by nx
203: clrl r12 # clear nx
204: brb Loop
205:
206: 7: movl r3,r5
207: movl r2,r4 # r4:r5 = y
208: subl2 $0x800000,r4 # r4:r5 = y/2
209: cmpd2 r0,r4 # x ? y/2
210: blss 9f # if x < y/2 goto 9
211: bgtr 8f # if x > y/2 goto 8
212: ldd r8 # acc = dble(n)
213: cvdl r8 # r8 = ifix(dble(n))
214: bbc $0,r8,9f # if the last bit is zero, goto 9
215: 8: ldd r0 # acc = x
216: subd r2 # acc = x-y
217: std r0 # r0:r1 = x-y
218: 9: xorl2 (sp)+,r0 # x^sign (exclusive or)
219: movl (sp)+,r6 # r6 = nf
220: andl3 $0x7f800000,r0,r8 # r8 = biased exponent of x
221: andl2 $0x807fffff,r0 # r0 = x w/ exponent zapped
222: subl2 r6,r8 # r8 = Ex-nf
223: bgtr 0f # if Ex-nf > 0 goto 0
224: clrl r8 # clear r8
225: clrl r0
226: clrl r1 # x underflows to zero
227: 0: orl2 r8,r0 # put r8 into x's exponent field
228: ret
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.