src/libm/s_atan.c
author Ryan C. Gordon <icculus@icculus.org>
Thu, 21 Apr 2016 03:16:44 -0400
changeset 11729 d1ce8396c356
parent 11683 48bcba563d9c
child 12420 4a6c91d9cc33
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.
slouken@4873
     1
/*
slouken@4873
     2
 * ====================================================
slouken@4873
     3
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
slouken@4873
     4
 *
slouken@4873
     5
 * Developed at SunPro, a Sun Microsystems, Inc. business.
slouken@4873
     6
 * Permission to use, copy, modify, and distribute this
slouken@4873
     7
 * software is freely granted, provided that this notice
slouken@4873
     8
 * is preserved.
slouken@4873
     9
 * ====================================================
slouken@4873
    10
 */
slouken@4873
    11
slouken@4873
    12
/* atan(x)
slouken@4873
    13
 * Method
slouken@4873
    14
 *   1. Reduce x to positive by atan(x) = -atan(-x).
slouken@4873
    15
 *   2. According to the integer k=4t+0.25 chopped, t=x, the argument
slouken@4873
    16
 *      is further reduced to one of the following intervals and the
slouken@4873
    17
 *      arctangent of t is evaluated by the corresponding formula:
slouken@4873
    18
 *
slouken@4873
    19
 *      [0,7/16]      atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
slouken@4873
    20
 *      [7/16,11/16]  atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
slouken@4873
    21
 *      [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
slouken@4873
    22
 *      [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
slouken@4873
    23
 *      [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )
slouken@4873
    24
 *
slouken@4873
    25
 * Constants:
slouken@4873
    26
 * The hexadecimal values are the intended ones for the following
slouken@4873
    27
 * constants. The decimal values may be used, provided that the
slouken@4873
    28
 * compiler will convert from decimal to binary accurately enough
slouken@4873
    29
 * to produce the hexadecimal values shown.
slouken@4873
    30
 */
slouken@4873
    31
slouken@6044
    32
#include "math_libm.h"
slouken@4873
    33
#include "math_private.h"
slouken@4873
    34
slouken@4873
    35
static const double atanhi[] = {
slouken@4873
    36
  4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
slouken@4873
    37
  7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
slouken@4873
    38
  9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
slouken@4873
    39
  1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
slouken@4873
    40
};
slouken@4873
    41
slouken@4873
    42
static const double atanlo[] = {
slouken@4873
    43
  2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
slouken@4873
    44
  3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
slouken@4873
    45
  1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
slouken@4873
    46
  6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
slouken@4873
    47
};
slouken@4873
    48
slouken@4873
    49
static const double aT[] = {
slouken@4873
    50
  3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
slouken@4873
    51
 -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
slouken@4873
    52
  1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
slouken@4873
    53
 -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
slouken@4873
    54
  9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
slouken@4873
    55
 -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
slouken@4873
    56
  6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
slouken@4873
    57
 -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
slouken@4873
    58
  4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
slouken@4873
    59
 -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
slouken@4873
    60
  1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
slouken@4873
    61
};
slouken@4873
    62
slouken@4873
    63
static const double
slouken@4873
    64
one   = 1.0,
slouken@4873
    65
huge   = 1.0e300;
slouken@4873
    66
slouken@4873
    67
double atan(double x)
slouken@4873
    68
{
slouken@4873
    69
	double w,s1,s2,z;
slouken@4873
    70
	int32_t ix,hx,id;
slouken@4873
    71
slouken@4873
    72
	GET_HIGH_WORD(hx,x);
slouken@4873
    73
	ix = hx&0x7fffffff;
slouken@4873
    74
	if(ix>=0x44100000) {	/* if |x| >= 2^66 */
slouken@4873
    75
	    u_int32_t low;
slouken@4873
    76
	    GET_LOW_WORD(low,x);
slouken@4873
    77
	    if(ix>0x7ff00000||
slouken@4873
    78
		(ix==0x7ff00000&&(low!=0)))
slouken@4873
    79
		return x+x;		/* NaN */
slouken@4873
    80
	    if(hx>0) return  atanhi[3]+atanlo[3];
slouken@4873
    81
	    else     return -atanhi[3]-atanlo[3];
slouken@4873
    82
	} if (ix < 0x3fdc0000) {	/* |x| < 0.4375 */
slouken@4873
    83
	    if (ix < 0x3e200000) {	/* |x| < 2^-29 */
slouken@4873
    84
		if(huge+x>one) return x;	/* raise inexact */
slouken@4873
    85
	    }
slouken@4873
    86
	    id = -1;
slouken@4873
    87
	} else {
slouken@4873
    88
	x = fabs(x);
slouken@4873
    89
	if (ix < 0x3ff30000) {		/* |x| < 1.1875 */
slouken@4873
    90
	    if (ix < 0x3fe60000) {	/* 7/16 <=|x|<11/16 */
slouken@4873
    91
		id = 0; x = (2.0*x-one)/(2.0+x);
slouken@4873
    92
	    } else {			/* 11/16<=|x|< 19/16 */
slouken@4873
    93
		id = 1; x  = (x-one)/(x+one);
slouken@4873
    94
	    }
slouken@4873
    95
	} else {
slouken@4873
    96
	    if (ix < 0x40038000) {	/* |x| < 2.4375 */
slouken@4873
    97
		id = 2; x  = (x-1.5)/(one+1.5*x);
slouken@4873
    98
	    } else {			/* 2.4375 <= |x| < 2^66 */
slouken@4873
    99
		id = 3; x  = -1.0/x;
slouken@4873
   100
	    }
slouken@4873
   101
	}}
slouken@4873
   102
    /* end of argument reduction */
slouken@4873
   103
	z = x*x;
slouken@4873
   104
	w = z*z;
slouken@4873
   105
    /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
slouken@4873
   106
	s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
slouken@4873
   107
	s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
slouken@4873
   108
	if (id<0) return x - x*(s1+s2);
slouken@4873
   109
	else {
slouken@4873
   110
	    z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
slouken@4873
   111
	    return (hx<0)? -z:z;
slouken@4873
   112
	}
slouken@4873
   113
}
slouken@4873
   114
libm_hidden_def(atan)