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