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