src/libm/e_atan2.c
author Ryan C. Gordon <icculus@icculus.org>
Thu, 21 Apr 2016 03:16:44 -0400
changeset 11729 d1ce8396c356
parent 11685 1788f503254d
permissions -rw-r--r--
Initial shot at a renderer target for Apple's Metal API.

This isn't complete, but is enough to run testsprite2. It's currently
Mac-only; with a little work to figure out how to properly glue in a Metal
layer to a UIView, this will likely work on iOS, too.

This is only wired up to the configure script right now, and disabled by
default. CMake and Xcode still need their bits filled in as appropriate.
     1 /*
     2  * ====================================================
     3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     4  *
     5  * Developed at SunPro, a Sun Microsystems, Inc. business.
     6  * Permission to use, copy, modify, and distribute this
     7  * software is freely granted, provided that this notice
     8  * is preserved.
     9  * ====================================================
    10  */
    11 
    12 /* __ieee754_atan2(y,x)
    13  * Method :
    14  *	1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
    15  *	2. Reduce x to positive by (if x and y are unexceptional):
    16  *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
    17  *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
    18  *
    19  * Special cases:
    20  *
    21  *	ATAN2((anything), NaN ) is NaN;
    22  *	ATAN2(NAN , (anything) ) is NaN;
    23  *	ATAN2(+-0, +(anything but NaN)) is +-0  ;
    24  *	ATAN2(+-0, -(anything but NaN)) is +-pi ;
    25  *	ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
    26  *	ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
    27  *	ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
    28  *	ATAN2(+-INF,+INF ) is +-pi/4 ;
    29  *	ATAN2(+-INF,-INF ) is +-3pi/4;
    30  *	ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
    31  *
    32  * Constants:
    33  * The hexadecimal values are the intended ones for the following
    34  * constants. The decimal values may be used, provided that the
    35  * compiler will convert from decimal to binary accurately enough
    36  * to produce the hexadecimal values shown.
    37  */
    38 
    39 #include "math_libm.h"
    40 #include "math_private.h"
    41 
    42 static const double
    43 tiny  = 1.0e-300,
    44 zero  = 0.0,
    45 pi_o_4  = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
    46 pi_o_2  = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
    47 pi      = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
    48 pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
    49 
    50 double attribute_hidden __ieee754_atan2(double y, double x)
    51 {
    52 	double z;
    53 	int32_t k,m,hx,hy,ix,iy;
    54 	u_int32_t lx,ly;
    55 
    56 	EXTRACT_WORDS(hx,lx,x);
    57 	ix = hx&0x7fffffff;
    58 	EXTRACT_WORDS(hy,ly,y);
    59 	iy = hy&0x7fffffff;
    60 	if(((ix|((lx|-(int32_t)lx)>>31))>0x7ff00000)||
    61 	   ((iy|((ly|-(int32_t)ly)>>31))>0x7ff00000))	/* x or y is NaN */
    62 	   return x+y;
    63 	if(((hx-0x3ff00000)|lx)==0) return atan(y);   /* x=1.0 */
    64 	m = ((hy>>31)&1)|((hx>>30)&2);	/* 2*sign(x)+sign(y) */
    65 
    66     /* when y = 0 */
    67 	if((iy|ly)==0) {
    68 	    switch(m) {
    69 		case 0:
    70 		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
    71 		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
    72 		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
    73 	    }
    74 	}
    75     /* when x = 0 */
    76 	if((ix|lx)==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;
    77 
    78     /* when x is INF */
    79 	if(ix==0x7ff00000) {
    80 	    if(iy==0x7ff00000) {
    81 		switch(m) {
    82 		    case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
    83 		    case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
    84 		    case 2: return  3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
    85 		    case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
    86 		}
    87 	    } else {
    88 		switch(m) {
    89 		    case 0: return  zero  ;	/* atan(+...,+INF) */
    90 		    case 1: return -zero  ;	/* atan(-...,+INF) */
    91 		    case 2: return  pi+tiny  ;	/* atan(+...,-INF) */
    92 		    case 3: return -pi-tiny  ;	/* atan(-...,-INF) */
    93 		}
    94 	    }
    95 	}
    96     /* when y is INF */
    97 	if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
    98 
    99     /* compute y/x */
   100 	k = (iy-ix)>>20;
   101 	if(k > 60) z=pi_o_2+0.5*pi_lo; 	/* |y/x| >  2**60 */
   102 	else if(hx<0&&k<-60) z=0.0; 	/* |y|/x < -2**60 */
   103 	else z=atan(fabs(y/x));		/* safe to do y/x */
   104 	switch (m) {
   105 	    case 0: return       z  ;	/* atan(+,+) */
   106 	    case 1: {
   107 	    	      u_int32_t zh;
   108 		      GET_HIGH_WORD(zh,z);
   109 		      SET_HIGH_WORD(z,zh ^ 0x80000000);
   110 		    }
   111 		    return       z  ;	/* atan(-,+) */
   112 	    case 2: return  pi-(z-pi_lo);/* atan(+,-) */
   113 	    default: /* case 3 */
   114 	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
   115 	}
   116 }
   117 
   118 /*
   119  * wrapper atan2(y,x)
   120  */
   121 #ifndef _IEEE_LIBM
   122 double atan2(double y, double x)
   123 {
   124 	double z = __ieee754_atan2(y, x);
   125 	if (_LIB_VERSION == _IEEE_ || isnan(x) || isnan(y))
   126 		return z;
   127 	if (x == 0.0 && y == 0.0)
   128 		return __kernel_standard(y,x,3); /* atan2(+-0,+-0) */
   129 	return z;
   130 }
   131 #else
   132 strong_alias(__ieee754_atan2, atan2)
   133 #endif
   134 libm_hidden_def(atan2)