![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to atan2.s CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSD / usr.lib / libm / VAX|
Return to atan2.s CVS log | Up to [CSRG BSD Unix] / 43BSD / usr.lib / libm / VAX |
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:
14: # @(#)atan2.s 1.2 (Berkeley) 8/21/85
15:
16: # ATAN2(Y,X)
17: # RETURN ARG (X+iY)
18: # VAX D FORMAT (56 BITS PRECISION)
19: # CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/16/85;
20: #
21: #
22: # Method :
23: # 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
24: # 2. Reduce x to positive by (if x and y are unexceptional):
25: # ARG (x+iy) = arctan(y/x) ... if x > 0,
26: # ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
27: # 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument
28: # is further reduced to one of the following intervals and the
29: # arctangent of y/x is evaluated by the corresponding formula:
30: #
31: # [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
32: # [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) )
33: # [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) )
34: # [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) )
35: # [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y )
36: #
37: # Special cases:
38: # Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y).
39: #
40: # ARG( NAN , (anything) ) is NaN;
41: # ARG( (anything), NaN ) is NaN;
42: # ARG(+(anything but NaN), +-0) is +-0 ;
43: # ARG(-(anything but NaN), +-0) is +-PI ;
44: # ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2;
45: # ARG( +INF,+-(anything but INF and NaN) ) is +-0 ;
46: # ARG( -INF,+-(anything but INF and NaN) ) is +-PI;
47: # ARG( +INF,+-INF ) is +-PI/4 ;
48: # ARG( -INF,+-INF ) is +-3PI/4;
49: # ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
50: #
51: # Accuracy:
52: # atan2(y,x) returns the exact ARG(x+iy) nearly rounded.
53: #
54: .text
55: .align 1
56: .globl _atan2
57: _atan2 :
58: .word 0x0ff4
59: movq 4(ap),r2 # r2 = y
60: movq 12(ap),r4 # r4 = x
61: bicw3 $0x7f,r2,r0
62: bicw3 $0x7f,r4,r1
63: cmpw r0,$0x8000 # y is the reserved operand
64: jeql resop
65: cmpw r1,$0x8000 # x is the reserved operand
66: jeql resop
67: subl2 $8,sp
68: bicw3 $0x7fff,r2,-4(fp) # copy y sign bit to -4(fp)
69: bicw3 $0x7fff,r4,-8(fp) # copy x sign bit to -8(fp)
70: cmpd r4,$0x4080 # x = 1.0 ?
71: bneq xnot1
72: movq r2,r0
73: bicw2 $0x8000,r0 # t = |y|
74: movq r0,r2 # y = |y|
75: brb begin
76: xnot1:
77: bicw3 $0x807f,r2,r11 # yexp
78: jeql yeq0 # if y=0 goto yeq0
79: bicw3 $0x807f,r4,r10 # xexp
80: jeql pio2 # if x=0 goto pio2
81: subw2 r10,r11 # k = yexp - xexp
82: cmpw r11,$0x2000 # k >= 64 (exp) ?
83: jgeq pio2 # atan2 = +-pi/2
84: divd3 r4,r2,r0 # t = y/x never overflow
85: bicw2 $0x8000,r0 # t > 0
86: bicw2 $0xff80,r2 # clear the exponent of y
87: bicw2 $0xff80,r4 # clear the exponent of x
88: bisw2 $0x4080,r2 # normalize y to [1,2)
89: bisw2 $0x4080,r4 # normalize x to [1,2)
90: subw2 r11,r4 # scale x so that yexp-xexp=k
91: begin:
92: cmpw r0,$0x411c # t : 39/16
93: jgeq L50
94: addl3 $0x180,r0,r10 # 8*t
95: cvtrfl r10,r10 # [8*t] rounded to int
96: ashl $-1,r10,r10 # [8*t]/2
97: casel r10,$0,$4
98: L1:
99: .word L20-L1
100: .word L20-L1
101: .word L30-L1
102: .word L40-L1
103: .word L40-L1
104: L10:
105: movq $0xb4d9940f985e407b,r6 # Hi=.98279372324732906796d0
106: movq $0x21b1879a3bc2a2fc,r8 # Lo=-.17092002525602665777d-17
107: subd3 r4,r2,r0 # y-x
108: addw2 $0x80,r0 # 2(y-x)
109: subd2 r4,r0 # 2(y-x)-x
110: addw2 $0x80,r4 # 2x
111: movq r2,r10
112: addw2 $0x80,r10 # 2y
113: addd2 r10,r2 # 3y
114: addd2 r4,r2 # 3y+2x
115: divd2 r2,r0 # (2y-3x)/(2x+3y)
116: brw L60
117: L20:
118: cmpw r0,$0x3280 # t : 2**(-28)
119: jlss L80
120: clrq r6 # Hi=r6=0, Lo=r8=0
121: clrq r8
122: brw L60
123: L30:
124: movq $0xda7b2b0d63383fed,r6 # Hi=.46364760900080611433d0
125: movq $0xf0ea17b2bf912295,r8 # Lo=.10147340032515978826d-17
126: movq r2,r0
127: addw2 $0x80,r0 # 2y
128: subd2 r4,r0 # 2y-x
129: addw2 $0x80,r4 # 2x
130: addd2 r2,r4 # 2x+y
131: divd2 r4,r0 # (2y-x)/(2x+y)
132: brb L60
133: L50:
134: movq $0x68c2a2210fda40c9,r6 # Hi=1.5707963267948966135d1
135: movq $0x06e0145c26332326,r8 # Lo=.22517417741562176079d-17
136: cmpw r0,$0x5100 # y : 2**57
137: bgeq L90
138: divd3 r2,r4,r0
139: bisw2 $0x8000,r0 # -x/y
140: brb L60
141: L40:
142: movq $0x68c2a2210fda4049,r6 # Hi=.78539816339744830676d0
143: movq $0x06e0145c263322a6,r8 # Lo=.11258708870781088040d-17
144: subd3 r4,r2,r0 # y-x
145: addd2 r4,r2 # y+x
146: divd2 r2,r0 # (y-x)/(y+x)
147: L60:
148: movq r0,r10
149: muld2 r0,r0
150: polyd r0,$12,ptable
151: muld2 r10,r0
152: subd2 r0,r8
153: addd3 r8,r10,r0
154: addd2 r6,r0
155: L80:
156: movw -8(fp),r2
157: bneq pim
158: bisw2 -4(fp),r0 # return sign(y)*r0
159: ret
160: L90: # x >= 2**25
161: movq r6,r0
162: brb L80
163: pim:
164: subd3 r0,$0x68c2a2210fda4149,r0 # pi-t
165: bisw2 -4(fp),r0
166: ret
167: yeq0:
168: movw -8(fp),r2
169: beql zero # if sign(x)=1 return pi
170: movq $0x68c2a2210fda4149,r0 # pi=3.1415926535897932270d1
171: ret
172: zero:
173: clrq r0 # return 0
174: ret
175: pio2:
176: movq $0x68c2a2210fda40c9,r0 # pi/2=1.5707963267948966135d1
177: bisw2 -4(fp),r0 # return sign(y)*pi/2
178: ret
179: resop:
180: movq $0x8000,r0 # propagate the reserved operand
181: ret
182: .align 2
183: ptable:
184: .quad 0xb50f5ce96e7abd60
185: .quad 0x51e44a42c1073e02
186: .quad 0x3487e3289643be35
187: .quad 0xdb62066dffba3e54
188: .quad 0xcf8e2d5199abbe70
189: .quad 0x26f39cb884883e88
190: .quad 0x135117d18998be9d
191: .quad 0x602ce9742e883eba
192: .quad 0xa35ad0be8e38bee3
193: .quad 0xffac922249243f12
194: .quad 0x7f14ccccccccbf4c
195: .quad 0xaa8faaaaaaaa3faa
196: .quad 0x0000000000000000
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.