|
|
1.1 root 1: /*
2: * Copyright (c) 1985 Regents of the University of California.
3: *
4: * Use and reproduction of this software are granted in accordance with
5: * the terms and conditions specified in the Berkeley Software License
6: * Agreement (in particular, this entails acknowledgement of the programs'
7: * source, and inclusion of this notice) with the additional understanding
8: * that all recipients should regard themselves as participants in an
9: * ongoing research project and hence should feel obligated to report
10: * their experiences (good or bad) with these elementary function codes,
11: * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
12: *
13: * @(#)support.s 1.3 (Berkeley) 8/21/85
14: *
15: * copysign(x,y),
16: * logb(x),
17: * scalb(x,N),
18: * finite(x),
19: * drem(x,y),
20: * Coded in vax assembly language by K.C. Ng, 3/14/85.
21: * Revised by K.C. Ng on 4/9/85.
22: */
23:
24: /*
25: * double copysign(x,y)
26: * double x,y;
27: */
28: .globl _copysign
29: .text
30: .align 1
31: _copysign:
32: .word 0x4
33: movq 4(ap),r0 # load x into r0
34: bicw3 $0x807f,r0,r2 # mask off the exponent of x
35: beql Lz # if zero or reserved op then return x
36: bicw3 $0x7fff,12(ap),r2 # copy the sign bit of y into r2
37: bicw2 $0x8000,r0 # replace x by |x|
38: bisw2 r2,r0 # copy the sign bit of y to x
39: Lz: ret
40:
41: /*
42: * double logb(x)
43: * double x;
44: */
45: .globl _logb
46: .text
47: .align 1
48: _logb:
49: .word 0x0
50: bicl3 $0xffff807f,4(ap),r0 # mask off the exponent of x
51: beql Ln
52: ashl $-7,r0,r0 # get the bias exponent
53: subl2 $129,r0 # get the unbias exponent
54: cvtld r0,r0 # return the answer in double
55: ret
56: Ln: movq 4(ap),r0 # r0:1 = x (zero or reserved op)
57: bneq 1f # simply return if reserved op
58: movq $0x0000fe00ffffcfff,r0 # -2147483647.0
59: 1: ret
60:
61: /*
62: * long finite(x)
63: * double x;
64: */
65: .globl _finite
66: .text
67: .align 1
68: _finite:
69: .word 0x0000
70: bicw3 $0x7f,4(ap),r0 # mask off the mantissa
71: cmpw r0,$0x8000 # to see if x is the reserved op
72: beql 1f # if so, return FALSE (0)
73: movl $1,r0 # else return TRUE (1)
74: ret
75: 1: clrl r0
76: ret
77:
78: /*
79: * double scalb(x,N)
80: * double x; int N;
81: */
82: .globl _scalb
83: .set ERANGE,34
84: .text
85: .align 1
86: _scalb:
87: .word 0xc
88: movq 4(ap),r0
89: bicl3 $0xffff807f,r0,r3
90: beql ret1 # 0 or reserved operand
91: movl 12(ap),r2
92: cmpl r2,$0x12c
93: bgeq ovfl
94: cmpl r2,$-0x12c
95: bleq unfl
96: ashl $7,r2,r2
97: addl2 r2,r3
98: bleq unfl
99: cmpl r3,$0x8000
100: bgeq ovfl
101: addl2 r2,r0
102: ret
103: ovfl: pushl $ERANGE
104: calls $1,_infnan # if it returns
105: bicw3 $0x7fff,4(ap),r2 # get the sign of input arg
106: bisw2 r2,r0 # re-attach the sign to r0/1
107: ret
108: unfl: movq $0,r0
109: ret1: ret
110:
111: /*
112: * DREM(X,Y)
113: * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
114: * DOUBLE PRECISION (VAX D format 56 bits)
115: * CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/8/85.
116: */
117: .globl _drem
118: .set EDOM,33
119: .text
120: .align 1
121: _drem:
122: .word 0xffc
123: subl2 $12,sp
124: movq 4(ap),r0 #r0=x
125: movq 12(ap),r2 #r2=y
126: jeql Rop #if y=0 then generate reserved op fault
127: bicw3 $0x007f,r0,r4 #check if x is Rop
128: cmpw r4,$0x8000
129: jeql Ret #if x is Rop then return Rop
130: bicl3 $0x007f,r2,r4 #check if y is Rop
131: cmpw r4,$0x8000
132: jeql Ret #if y is Rop then return Rop
133: bicw2 $0x8000,r2 #y := |y|
134: movw $0,-4(fp) #-4(fp) = nx := 0
135: cmpw r2,$0x1c80 #yexp ? 57
136: bgtr C1 #if yexp > 57 goto C1
137: addw2 $0x1c80,r2 #scale up y by 2**57
138: movw $0x1c80,-4(fp) #nx := 57 (exponent field)
139: C1:
140: movw -4(fp),-8(fp) #-8(fp) = nf := nx
141: bicw3 $0x7fff,r0,-12(fp) #-12(fp) = sign of x
142: bicw2 $0x8000,r0 #x := |x|
143: movq r2,r10 #y1 := y
144: bicl2 $0xffff07ff,r11 #clear the last 27 bits of y1
145: loop:
146: cmpd r0,r2 #x ? y
147: bleq E1 #if x <= y goto E1
148: /* begin argument reduction */
149: movq r2,r4 #t =y
150: movq r10,r6 #t1=y1
151: bicw3 $0x807f,r0,r8 #xexp= exponent of x
152: bicw3 $0x807f,r2,r9 #yexp= exponent fo y
153: subw2 r9,r8 #xexp-yexp
154: subw2 $0x0c80,r8 #k=xexp-yexp-25(exponent bit field)
155: blss C2 #if k<0 goto C2
156: addw2 r8,r4 #t +=k
157: addw2 r8,r6 #t1+=k, scale up t and t1
158: C2:
159: divd3 r4,r0,r8 #x/t
160: cvtdl r8,r8 #n=[x/t] truncated
161: cvtld r8,r8 #float(n)
162: subd2 r6,r4 #t:=t-t1
163: muld2 r8,r4 #n*(t-t1)
164: muld2 r8,r6 #n*t1
165: subd2 r6,r0 #x-n*t1
166: subd2 r4,r0 #(x-n*t1)-n*(t-t1)
167: brb loop
168: E1:
169: movw -4(fp),r6 #r6=nx
170: beql C3 #if nx=0 goto C3
171: addw2 r6,r0 #x:=x*2**57 scale up x by nx
172: movw $0,-4(fp) #clear nx
173: brb loop
174: C3:
175: movq r2,r4 #r4 = y
176: subw2 $0x80,r4 #r4 = y/2
177: cmpd r0,r4 #x:y/2
178: blss E2 #if x < y/2 goto E2
179: bgtr C4 #if x > y/2 goto C4
180: cvtdl r8,r8 #ifix(float(n))
181: blbc r8,E2 #if the last bit is zero, goto E2
182: C4:
183: subd2 r2,r0 #x-y
184: E2:
185: xorw2 -12(fp),r0 #x^sign (exclusive or)
186: movw -8(fp),r6 #r6=nf
187: bicw3 $0x807f,r0,r8 #r8=exponent of x
188: bicw2 $0x7f80,r0 #clear the exponent of x
189: subw2 r6,r8 #r8=xexp-nf
190: bgtr C5 #if xexp-nf is positive goto C5
191: movw $0,r8 #clear r8
192: movq $0,r0 #x underflow to zero
193: C5:
194: bisw2 r8,r0 #put r8 into x's exponent field
195: ret
196: Rop: #Reserved operand
197: pushl $EDOM
198: calls $1,_infnan #generate reserved op fault
199: ret
200: Ret:
201: movq $0x8000,r0 #propagate reserved op
202: ret
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.