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