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