n_atan2.S revision 1.6
1 /* $NetBSD: 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 ENTRY(atan2, 0x0fc0) 76 movq 4(%ap),%r2 # %r2 = y 77 movq 12(%ap),%r4 # %r4 = x 78 bicw3 $0x7f,%r2,%r0 79 bicw3 $0x7f,%r4,%r1 80 cmpw %r0,$0x8000 # y is the reserved operand 81 jeql resop 82 cmpw %r1,$0x8000 # x is the reserved operand 83 jeql resop 84 subl2 $8,%sp 85 bicw3 $0x7fff,%r2,-4(%fp) # copy y sign bit to -4(%fp) 86 bicw3 $0x7fff,%r4,-8(%fp) # copy x sign bit to -8(%fp) 87 cmpd %r4,$0x4080 # x = 1.0 ? 88 bneq xnot1 89 movq %r2,%r0 90 bicw2 $0x8000,%r0 # t = |y| 91 movq %r0,%r2 # y = |y| 92 jbr begin 93 xnot1: 94 bicw3 $0x807f,%r2,%r11 # yexp 95 jeql yeq0 # if y=0 goto yeq0 96 bicw3 $0x807f,%r4,%r10 # xexp 97 jeql pio2 # if x=0 goto pio2 98 subw2 %r10,%r11 # k = yexp - xexp 99 cmpw %r11,$0x2000 # k >= 64 (exp) ? 100 jgeq pio2 # atan2 = +-pi/2 101 divd3 %r4,%r2,%r0 # t = y/x never overflow 102 bicw2 $0x8000,%r0 # t > 0 103 bicw2 $0xff80,%r2 # clear the exponent of y 104 bicw2 $0xff80,%r4 # clear the exponent of x 105 bisw2 $0x4080,%r2 # normalize y to [1,2) 106 bisw2 $0x4080,%r4 # normalize x to [1,2) 107 subw2 %r11,%r4 # scale x so that yexp-xexp=k 108 begin: 109 cmpw %r0,$0x411c # t : 39/16 110 jgeq L50 111 addl3 $0x180,%r0,%r10 # 8*t 112 cvtrfl %r10,%r10 # [8*t] rounded to int 113 ashl $-1,%r10,%r10 # [8*t]/2 114 casel %r10,$0,$4 115 L1: 116 .word L20-L1 117 .word L20-L1 118 .word L30-L1 119 .word L40-L1 120 .word L40-L1 121 L10: 122 movq $0xb4d9940f985e407b,%r6 # Hi=.98279372324732906796d0 123 movq $0x21b1879a3bc2a2fc,%r8 # Lo=-.17092002525602665777d-17 124 subd3 %r4,%r2,%r0 # y-x 125 addw2 $0x80,%r0 # 2(y-x) 126 subd2 %r4,%r0 # 2(y-x)-x 127 addw2 $0x80,%r4 # 2x 128 movq %r2,%r10 129 addw2 $0x80,%r10 # 2y 130 addd2 %r10,%r2 # 3y 131 addd2 %r4,%r2 # 3y+2x 132 divd2 %r2,%r0 # (2y-3x)/(2x+3y) 133 jbr L60 134 L20: 135 cmpw %r0,$0x3280 # t : 2**(-28) 136 jlss L80 137 clrq %r6 # Hi=%r6=0, Lo=%r8=0 138 clrq %r8 139 jbr L60 140 L30: 141 movq $0xda7b2b0d63383fed,%r6 # Hi=.46364760900080611433d0 142 movq $0xf0ea17b2bf912295,%r8 # Lo=.10147340032515978826d-17 143 movq %r2,%r0 144 addw2 $0x80,%r0 # 2y 145 subd2 %r4,%r0 # 2y-x 146 addw2 $0x80,%r4 # 2x 147 addd2 %r2,%r4 # 2x+y 148 divd2 %r4,%r0 # (2y-x)/(2x+y) 149 jbr L60 150 L50: 151 movq $0x68c2a2210fda40c9,%r6 # Hi=1.5707963267948966135d1 152 movq $0x06e0145c26332326,%r8 # Lo=.22517417741562176079d-17 153 cmpw %r0,$0x5100 # y : 2**57 154 bgeq L90 155 divd3 %r2,%r4,%r0 156 bisw2 $0x8000,%r0 # -x/y 157 jbr L60 158 L40: 159 movq $0x68c2a2210fda4049,%r6 # Hi=.78539816339744830676d0 160 movq $0x06e0145c263322a6,%r8 # Lo=.11258708870781088040d-17 161 subd3 %r4,%r2,%r0 # y-x 162 addd2 %r4,%r2 # y+x 163 divd2 %r2,%r0 # (y-x)/(y+x) 164 L60: 165 movq %r0,%r10 166 muld2 %r0,%r0 167 polyd %r0,$12,ptable 168 muld2 %r10,%r0 169 subd2 %r0,%r8 170 addd3 %r8,%r10,%r0 171 addd2 %r6,%r0 172 L80: 173 movw -8(%fp),%r2 174 bneq pim 175 bisw2 -4(%fp),%r0 # return sign(y)*%r0 176 ret 177 L90: # x >= 2**25 178 movq %r6,%r0 179 jbr L80 180 pim: 181 subd3 %r0,$0x68c2a2210fda4149,%r0 # pi-t 182 bisw2 -4(%fp),%r0 183 ret 184 yeq0: 185 movw -8(%fp),%r2 186 beql zero # if sign(x)=1 return pi 187 movq $0x68c2a2210fda4149,%r0 # pi=3.1415926535897932270d1 188 ret 189 zero: 190 clrq %r0 # return 0 191 ret 192 pio2: 193 movq $0x68c2a2210fda40c9,%r0 # pi/2=1.5707963267948966135d1 194 bisw2 -4(%fp),%r0 # return sign(y)*pi/2 195 ret 196 resop: 197 movq $0x8000,%r0 # propagate the reserved operand 198 ret 199 200 _ALIGN_TEXT 201 ptable: 202 .quad 0xb50f5ce96e7abd60 203 .quad 0x51e44a42c1073e02 204 .quad 0x3487e3289643be35 205 .quad 0xdb62066dffba3e54 206 .quad 0xcf8e2d5199abbe70 207 .quad 0x26f39cb884883e88 208 .quad 0x135117d18998be9d 209 .quad 0x602ce9742e883eba 210 .quad 0xa35ad0be8e38bee3 211 .quad 0xffac922249243f12 212 .quad 0x7f14ccccccccbf4c 213 .quad 0xaa8faaaaaaaa3faa 214 .quad 0x0000000000000000 215
Indexes created Thu Oct 02 07:10:07 GMT 2025