n_atan2.S revision 1.6.22.1
1 /* n_atan2.S,v 1.6 2003/08/07 16:44:44 agc Exp */ 2 /* 3 * Copyright (c) 1985, 1993 4 * The Regents of the University of California. All rights reserved. 5 * 6 * Redistribution and use in source and binary forms, with or without 7 * modification, are permitted provided that the following conditions 8 * are met: 9 * 1. Redistributions of source code must retain the above copyright 10 * notice, this list of conditions and the following disclaimer. 11 * 2. Redistributions in binary form must reproduce the above copyright 12 * notice, this list of conditions and the following disclaimer in the 13 * documentation and/or other materials provided with the distribution. 14 * 3. Neither the name of the University nor the names of its contributors 15 * may be used to endorse or promote products derived from this software 16 * without specific prior written permission. 17 * 18 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 19 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 20 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 21 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 22 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 23 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 24 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 25 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 26 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 27 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 28 * SUCH DAMAGE. 29 * 30 * @(#)atan2.s 8.1 (Berkeley) 6/4/93 31 */ 32 33 #include <machine/asm.h> 34 35 /* 36 * ATAN2(Y,X) 37 * RETURN ARG (X+iY) 38 * VAX D FORMAT (56 BITS PRECISION) 39 * CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/16/85; 40 * 41 * 42 * Method : 43 * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). 44 * 2. Reduce x to positive by (if x and y are unexceptional): 45 * ARG (x+iy) = arctan(y/x) ... if x > 0, 46 * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, 47 * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument 48 * is further reduced to one of the following intervals and the 49 * arctangent of y/x is evaluated by the corresponding formula: 50 * 51 * [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) 52 * [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) ) 53 * [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) ) 54 * [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) ) 55 * [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y ) 56 * 57 * Special cases: 58 * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y). 59 * 60 * ARG( NAN , (anything) ) is NaN; 61 * ARG( (anything), NaN ) is NaN; 62 * ARG(+(anything but NaN), +-0) is +-0 ; 63 * ARG(-(anything but NaN), +-0) is +-PI ; 64 * ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2; 65 * ARG( +INF,+-(anything but INF and NaN) ) is +-0 ; 66 * ARG( -INF,+-(anything but INF and NaN) ) is +-PI; 67 * ARG( +INF,+-INF ) is +-PI/4 ; 68 * ARG( -INF,+-INF ) is +-3PI/4; 69 * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2; 70 * 71 * Accuracy: 72 * atan2(y,x) returns the exact ARG(x+iy) nearly rounded. 73 */ 74 75 #ifdef WEAK_ALIAS 76 WEAK_ALIAS(atan2f, _atan2f) 77 #endif 78 79 ENTRY(_atan2f, 0) 80 cvtfd 4(%ap),-(%sp) 81 calls $2,_C_LABEL(_atan2) 82 cvtdf %r0,%r0 83 ret 84 85 #ifdef WEAK_ALIAS 86 WEAK_ALIAS(atan2, _atan2) 87 #endif 88 89 ENTRY(_atan2, 0x0fc0) 90 movq 4(%ap),%r2 # %r2 = y 91 movq 12(%ap),%r4 # %r4 = x 92 bicw3 $0x7f,%r2,%r0 93 bicw3 $0x7f,%r4,%r1 94 cmpw %r0,$0x8000 # y is the reserved operand 95 jeql resop 96 cmpw %r1,$0x8000 # x is the reserved operand 97 jeql resop 98 subl2 $8,%sp 99 bicw3 $0x7fff,%r2,-4(%fp) # copy y sign bit to -4(%fp) 100 bicw3 $0x7fff,%r4,-8(%fp) # copy x sign bit to -8(%fp) 101 cmpd %r4,$0x4080 # x = 1.0 ? 102 bneq xnot1 103 movq %r2,%r0 104 bicw2 $0x8000,%r0 # t = |y| 105 movq %r0,%r2 # y = |y| 106 jbr begin 107 xnot1: 108 bicw3 $0x807f,%r2,%r11 # yexp 109 jeql yeq0 # if y=0 goto yeq0 110 bicw3 $0x807f,%r4,%r10 # xexp 111 jeql pio2 # if x=0 goto pio2 112 subw2 %r10,%r11 # k = yexp - xexp 113 cmpw %r11,$0x2000 # k >= 64 (exp) ? 114 jgeq pio2 # atan2 = +-pi/2 115 divd3 %r4,%r2,%r0 # t = y/x never overflow 116 bicw2 $0x8000,%r0 # t > 0 117 bicw2 $0xff80,%r2 # clear the exponent of y 118 bicw2 $0xff80,%r4 # clear the exponent of x 119 bisw2 $0x4080,%r2 # normalize y to [1,2) 120 bisw2 $0x4080,%r4 # normalize x to [1,2) 121 subw2 %r11,%r4 # scale x so that yexp-xexp=k 122 begin: 123 cmpw %r0,$0x411c # t : 39/16 124 jgeq L50 125 addl3 $0x180,%r0,%r10 # 8*t 126 cvtrfl %r10,%r10 # [8*t] rounded to int 127 ashl $-1,%r10,%r10 # [8*t]/2 128 casel %r10,$0,$4 129 L1: 130 .word L20-L1 131 .word L20-L1 132 .word L30-L1 133 .word L40-L1 134 .word L40-L1 135 L10: 136 movq $0xb4d9940f985e407b,%r6 # Hi=.98279372324732906796d0 137 movq $0x21b1879a3bc2a2fc,%r8 # Lo=-.17092002525602665777d-17 138 subd3 %r4,%r2,%r0 # y-x 139 addw2 $0x80,%r0 # 2(y-x) 140 subd2 %r4,%r0 # 2(y-x)-x 141 addw2 $0x80,%r4 # 2x 142 movq %r2,%r10 143 addw2 $0x80,%r10 # 2y 144 addd2 %r10,%r2 # 3y 145 addd2 %r4,%r2 # 3y+2x 146 divd2 %r2,%r0 # (2y-3x)/(2x+3y) 147 jbr L60 148 L20: 149 cmpw %r0,$0x3280 # t : 2**(-28) 150 jlss L80 151 clrq %r6 # Hi=%r6=0, Lo=%r8=0 152 clrq %r8 153 jbr L60 154 L30: 155 movq $0xda7b2b0d63383fed,%r6 # Hi=.46364760900080611433d0 156 movq $0xf0ea17b2bf912295,%r8 # Lo=.10147340032515978826d-17 157 movq %r2,%r0 158 addw2 $0x80,%r0 # 2y 159 subd2 %r4,%r0 # 2y-x 160 addw2 $0x80,%r4 # 2x 161 addd2 %r2,%r4 # 2x+y 162 divd2 %r4,%r0 # (2y-x)/(2x+y) 163 jbr L60 164 L50: 165 movq $0x68c2a2210fda40c9,%r6 # Hi=1.5707963267948966135d1 166 movq $0x06e0145c26332326,%r8 # Lo=.22517417741562176079d-17 167 cmpw %r0,$0x5100 # y : 2**57 168 bgeq L90 169 divd3 %r2,%r4,%r0 170 bisw2 $0x8000,%r0 # -x/y 171 jbr L60 172 L40: 173 movq $0x68c2a2210fda4049,%r6 # Hi=.78539816339744830676d0 174 movq $0x06e0145c263322a6,%r8 # Lo=.11258708870781088040d-17 175 subd3 %r4,%r2,%r0 # y-x 176 addd2 %r4,%r2 # y+x 177 divd2 %r2,%r0 # (y-x)/(y+x) 178 L60: 179 movq %r0,%r10 180 muld2 %r0,%r0 181 polyd %r0,$12,ptable 182 muld2 %r10,%r0 183 subd2 %r0,%r8 184 addd3 %r8,%r10,%r0 185 addd2 %r6,%r0 186 L80: 187 movw -8(%fp),%r2 188 bneq pim 189 bisw2 -4(%fp),%r0 # return sign(y)*%r0 190 ret 191 L90: # x >= 2**25 192 movq %r6,%r0 193 jbr L80 194 pim: 195 subd3 %r0,$0x68c2a2210fda4149,%r0 # pi-t 196 bisw2 -4(%fp),%r0 197 ret 198 yeq0: 199 movw -8(%fp),%r2 200 beql zero # if sign(x)=1 return pi 201 movq $0x68c2a2210fda4149,%r0 # pi=3.1415926535897932270d1 202 ret 203 zero: 204 clrq %r0 # return 0 205 ret 206 pio2: 207 movq $0x68c2a2210fda40c9,%r0 # pi/2=1.5707963267948966135d1 208 bisw2 -4(%fp),%r0 # return sign(y)*pi/2 209 ret 210 resop: 211 movq $0x8000,%r0 # propagate the reserved operand 212 ret 213 214 _ALIGN_TEXT 215 ptable: 216 .quad 0xb50f5ce96e7abd60 217 .quad 0x51e44a42c1073e02 218 .quad 0x3487e3289643be35 219 .quad 0xdb62066dffba3e54 220 .quad 0xcf8e2d5199abbe70 221 .quad 0x26f39cb884883e88 222 .quad 0x135117d18998be9d 223 .quad 0x602ce9742e883eba 224 .quad 0xa35ad0be8e38bee3 225 .quad 0xffac922249243f12 226 .quad 0x7f14ccccccccbf4c 227 .quad 0xaa8faaaaaaaa3faa 228 .quad 0x0000000000000000 229
Indexes created Thu Oct 02 07:10:07 GMT 2025