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