1/* $NetBSD: n_atan2.S,v 1.5 2002/02/24 01:06:21 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 79ENTRY(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 97xnot1: 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 112begin: 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 119L1: 120 .word L20-L1 121 .word L20-L1 122 .word L30-L1 123 .word L40-L1 124 .word L40-L1 125L10: 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 138L20: 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 144L30: 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 154L50: 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 162L40: 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) 168L60: 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 176L80: 177 movw -8(%fp),%r2 178 bneq pim 179 bisw2 -4(%fp),%r0 # return sign(y)*%r0 180 ret 181L90: # x >= 2**25 182 movq %r6,%r0 183 jbr L80 184pim: 185 subd3 %r0,$0x68c2a2210fda4149,%r0 # pi-t 186 bisw2 -4(%fp),%r0 187 ret 188yeq0: 189 movw -8(%fp),%r2 190 beql zero # if sign(x)=1 return pi 191 movq $0x68c2a2210fda4149,%r0 # pi=3.1415926535897932270d1 192 ret 193zero: 194 clrq %r0 # return 0 195 ret 196pio2: 197 movq $0x68c2a2210fda40c9,%r0 # pi/2=1.5707963267948966135d1 198 bisw2 -4(%fp),%r0 # return sign(y)*pi/2 199 ret 200resop: 201 movq $0x8000,%r0 # propagate the reserved operand 202 ret 203 204 _ALIGN_TEXT 205ptable: 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