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