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