1*25c28e83SPiotr Jasiukajtis /* 2*25c28e83SPiotr Jasiukajtis * CDDL HEADER START 3*25c28e83SPiotr Jasiukajtis * 4*25c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the 5*25c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License"). 6*25c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License. 7*25c28e83SPiotr Jasiukajtis * 8*25c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9*25c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing. 10*25c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions 11*25c28e83SPiotr Jasiukajtis * and limitations under the License. 12*25c28e83SPiotr Jasiukajtis * 13*25c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each 14*25c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15*25c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the 16*25c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying 17*25c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner] 18*25c28e83SPiotr Jasiukajtis * 19*25c28e83SPiotr Jasiukajtis * CDDL HEADER END 20*25c28e83SPiotr Jasiukajtis */ 21*25c28e83SPiotr Jasiukajtis 22*25c28e83SPiotr Jasiukajtis /* 23*25c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24*25c28e83SPiotr Jasiukajtis */ 25*25c28e83SPiotr Jasiukajtis /* 26*25c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27*25c28e83SPiotr Jasiukajtis * Use is subject to license terms. 28*25c28e83SPiotr Jasiukajtis */ 29*25c28e83SPiotr Jasiukajtis 30*25c28e83SPiotr Jasiukajtis #pragma weak atan = __atan 31*25c28e83SPiotr Jasiukajtis 32*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 33*25c28e83SPiotr Jasiukajtis /* 34*25c28e83SPiotr Jasiukajtis * atan(x) 35*25c28e83SPiotr Jasiukajtis * Accurate Table look-up algorithm with polynomial approximation in 36*25c28e83SPiotr Jasiukajtis * partially product form. 37*25c28e83SPiotr Jasiukajtis * 38*25c28e83SPiotr Jasiukajtis * -- K.C. Ng, October 17, 2004 39*25c28e83SPiotr Jasiukajtis * 40*25c28e83SPiotr Jasiukajtis * Algorithm 41*25c28e83SPiotr Jasiukajtis * 42*25c28e83SPiotr Jasiukajtis * (1). Purge off Inf and NaN and 0 43*25c28e83SPiotr Jasiukajtis * (2). Reduce x to positive by atan(x) = -atan(-x). 44*25c28e83SPiotr Jasiukajtis * (3). For x <= 1/8 and let z = x*x, return 45*25c28e83SPiotr Jasiukajtis * (2.1) if x < 2^(-prec/2), atan(x) = x with inexact flag raised 46*25c28e83SPiotr Jasiukajtis * (2.2) if x < 2^(-prec/4-1), atan(x) = x+(x/3)(x*x) 47*25c28e83SPiotr Jasiukajtis * (2.3) if x < 2^(-prec/6-2), atan(x) = x+(z-5/3)(z*x/5) 48*25c28e83SPiotr Jasiukajtis * (2.4) Otherwise 49*25c28e83SPiotr Jasiukajtis * atan(x) = poly1(x) = x + A * B, 50*25c28e83SPiotr Jasiukajtis * where 51*25c28e83SPiotr Jasiukajtis * A = (p1*x*z) * (p2+z(p3+z)) 52*25c28e83SPiotr Jasiukajtis * B = (p4+z)+z*z) * (p5+z(p6+z)) 53*25c28e83SPiotr Jasiukajtis * Note: (i) domain of poly1 is [0, 1/8], (ii) remez relative 54*25c28e83SPiotr Jasiukajtis * approximation error of poly1 is bounded by 55*25c28e83SPiotr Jasiukajtis * |(atan(x)-poly1(x))/x| <= 2^-57.61 56*25c28e83SPiotr Jasiukajtis * (4). For x >= 8 then 57*25c28e83SPiotr Jasiukajtis * (3.1) if x >= 2^prec, atan(x) = atan(inf) - pio2lo 58*25c28e83SPiotr Jasiukajtis * (3.2) if x >= 2^(prec/3), atan(x) = atan(inf) - 1/x 59*25c28e83SPiotr Jasiukajtis * (3.3) if x <= 65, atan(x) = atan(inf) - poly1(1/x) 60*25c28e83SPiotr Jasiukajtis * (3.4) otherwise atan(x) = atan(inf) - poly2(1/x) 61*25c28e83SPiotr Jasiukajtis * where 62*25c28e83SPiotr Jasiukajtis * poly2(r) = (q1*r) * (q2+z(q3+z)) * (q4+z), 63*25c28e83SPiotr Jasiukajtis * its domain is [0, 0.0154]; and its remez absolute 64*25c28e83SPiotr Jasiukajtis * approximation error is bounded by 65*25c28e83SPiotr Jasiukajtis * |atan(x)-poly2(x)|<= 2^-59.45 66*25c28e83SPiotr Jasiukajtis * 67*25c28e83SPiotr Jasiukajtis * (5). Now x is in (0.125, 8). 68*25c28e83SPiotr Jasiukajtis * Recall identity 69*25c28e83SPiotr Jasiukajtis * atan(x) = atan(y) + atan((x-y)/(1+x*y)). 70*25c28e83SPiotr Jasiukajtis * Let j = (ix - 0x3fc00000) >> 16, 0 <= j < 96, where ix is the high 71*25c28e83SPiotr Jasiukajtis * part of x in IEEE double format. Then 72*25c28e83SPiotr Jasiukajtis * atan(x) = atan(y[j]) + poly2((x-y[j])/(1+x*y[j])) 73*25c28e83SPiotr Jasiukajtis * where y[j] are carefully chosen so that it matches x to around 4.5 74*25c28e83SPiotr Jasiukajtis * bits and at the same time atan(y[j]) is very close to an IEEE double 75*25c28e83SPiotr Jasiukajtis * floating point number. Calculation indicates that 76*25c28e83SPiotr Jasiukajtis * max|(x-y[j])/(1+x*y[j])| < 0.0154 77*25c28e83SPiotr Jasiukajtis * j,x 78*25c28e83SPiotr Jasiukajtis * 79*25c28e83SPiotr Jasiukajtis * Accuracy: Maximum error observed is bounded by 0.6 ulp after testing 80*25c28e83SPiotr Jasiukajtis * more than 10 million random arguments 81*25c28e83SPiotr Jasiukajtis */ 82*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 83*25c28e83SPiotr Jasiukajtis 84*25c28e83SPiotr Jasiukajtis #include "libm.h" 85*25c28e83SPiotr Jasiukajtis #include "libm_synonyms.h" 86*25c28e83SPiotr Jasiukajtis #include "libm_protos.h" 87*25c28e83SPiotr Jasiukajtis 88*25c28e83SPiotr Jasiukajtis extern const double _TBL_atan[]; 89*25c28e83SPiotr Jasiukajtis static const double g[] = { 90*25c28e83SPiotr Jasiukajtis /* one = */ 1.0, 91*25c28e83SPiotr Jasiukajtis /* p1 = */ 8.02176624254765935351230154992663301527500152588e-0002, 92*25c28e83SPiotr Jasiukajtis /* p2 = */ 1.27223421700559402580665846471674740314483642578e+0000, 93*25c28e83SPiotr Jasiukajtis /* p3 = */ -1.20606901800503640842521235754247754812240600586e+0000, 94*25c28e83SPiotr Jasiukajtis /* p4 = */ -2.36088967922325565496066701598465442657470703125e+0000, 95*25c28e83SPiotr Jasiukajtis /* p5 = */ 1.38345799501389166152875986881554126739501953125e+0000, 96*25c28e83SPiotr Jasiukajtis /* p6 = */ 1.06742368078953453469637224770849570631980895996e+0000, 97*25c28e83SPiotr Jasiukajtis /* q1 = */ -1.42796626333911796935538518482644576579332351685e-0001, 98*25c28e83SPiotr Jasiukajtis /* q2 = */ 3.51427110447873227059810477159863497078605962912e+0000, 99*25c28e83SPiotr Jasiukajtis /* q3 = */ 5.92129112708164262457444237952586263418197631836e-0001, 100*25c28e83SPiotr Jasiukajtis /* q4 = */ -1.99272234785683144409063061175402253866195678711e+0000, 101*25c28e83SPiotr Jasiukajtis /* pio2hi */ 1.570796326794896558e+00, 102*25c28e83SPiotr Jasiukajtis /* pio2lo */ 6.123233995736765886e-17, 103*25c28e83SPiotr Jasiukajtis /* t1 = */ -0.333333333333333333333333333333333, 104*25c28e83SPiotr Jasiukajtis /* t2 = */ 0.2, 105*25c28e83SPiotr Jasiukajtis /* t3 = */ -1.666666666666666666666666666666666, 106*25c28e83SPiotr Jasiukajtis }; 107*25c28e83SPiotr Jasiukajtis 108*25c28e83SPiotr Jasiukajtis #define one g[0] 109*25c28e83SPiotr Jasiukajtis #define p1 g[1] 110*25c28e83SPiotr Jasiukajtis #define p2 g[2] 111*25c28e83SPiotr Jasiukajtis #define p3 g[3] 112*25c28e83SPiotr Jasiukajtis #define p4 g[4] 113*25c28e83SPiotr Jasiukajtis #define p5 g[5] 114*25c28e83SPiotr Jasiukajtis #define p6 g[6] 115*25c28e83SPiotr Jasiukajtis #define q1 g[7] 116*25c28e83SPiotr Jasiukajtis #define q2 g[8] 117*25c28e83SPiotr Jasiukajtis #define q3 g[9] 118*25c28e83SPiotr Jasiukajtis #define q4 g[10] 119*25c28e83SPiotr Jasiukajtis #define pio2hi g[11] 120*25c28e83SPiotr Jasiukajtis #define pio2lo g[12] 121*25c28e83SPiotr Jasiukajtis #define t1 g[13] 122*25c28e83SPiotr Jasiukajtis #define t2 g[14] 123*25c28e83SPiotr Jasiukajtis #define t3 g[15] 124*25c28e83SPiotr Jasiukajtis 125*25c28e83SPiotr Jasiukajtis 126*25c28e83SPiotr Jasiukajtis double 127*25c28e83SPiotr Jasiukajtis atan(double x) { 128*25c28e83SPiotr Jasiukajtis double y, z, r, p, s; 129*25c28e83SPiotr Jasiukajtis int ix, lx, hx, j; 130*25c28e83SPiotr Jasiukajtis 131*25c28e83SPiotr Jasiukajtis hx = ((int *) &x)[HIWORD]; 132*25c28e83SPiotr Jasiukajtis lx = ((int *) &x)[LOWORD]; 133*25c28e83SPiotr Jasiukajtis ix = hx & ~0x80000000; 134*25c28e83SPiotr Jasiukajtis j = ix >> 20; 135*25c28e83SPiotr Jasiukajtis 136*25c28e83SPiotr Jasiukajtis /* for |x| < 1/8 */ 137*25c28e83SPiotr Jasiukajtis if (j < 0x3fc) { 138*25c28e83SPiotr Jasiukajtis if (j < 0x3f5) { /* when |x| < 2**(-prec/6-2) */ 139*25c28e83SPiotr Jasiukajtis if (j < 0x3e3) { /* if |x| < 2**(-prec/2-2) */ 140*25c28e83SPiotr Jasiukajtis return ((int) x == 0 ? x : one); 141*25c28e83SPiotr Jasiukajtis } 142*25c28e83SPiotr Jasiukajtis if (j < 0x3f1) { /* if |x| < 2**(-prec/4-1) */ 143*25c28e83SPiotr Jasiukajtis return (x + (x * t1) * (x * x)); 144*25c28e83SPiotr Jasiukajtis } else { /* if |x| < 2**(-prec/6-2) */ 145*25c28e83SPiotr Jasiukajtis z = x * x; 146*25c28e83SPiotr Jasiukajtis s = t2 * x; 147*25c28e83SPiotr Jasiukajtis return (x + (t3 + z) * (s * z)); 148*25c28e83SPiotr Jasiukajtis } 149*25c28e83SPiotr Jasiukajtis } 150*25c28e83SPiotr Jasiukajtis z = x * x; s = p1 * x; 151*25c28e83SPiotr Jasiukajtis return (x + ((s * z) * (p2 + z * (p3 + z))) * 152*25c28e83SPiotr Jasiukajtis (((p4 + z) + z * z) * (p5 + z * (p6 + z)))); 153*25c28e83SPiotr Jasiukajtis } 154*25c28e83SPiotr Jasiukajtis 155*25c28e83SPiotr Jasiukajtis /* for |x| >= 8.0 */ 156*25c28e83SPiotr Jasiukajtis if (j >= 0x402) { 157*25c28e83SPiotr Jasiukajtis if (j < 0x436) { 158*25c28e83SPiotr Jasiukajtis r = one / x; 159*25c28e83SPiotr Jasiukajtis if (hx >= 0) { 160*25c28e83SPiotr Jasiukajtis y = pio2hi; p = pio2lo; 161*25c28e83SPiotr Jasiukajtis } else { 162*25c28e83SPiotr Jasiukajtis y = -pio2hi; p = -pio2lo; 163*25c28e83SPiotr Jasiukajtis } 164*25c28e83SPiotr Jasiukajtis if (ix < 0x40504000) { /* x < 65 */ 165*25c28e83SPiotr Jasiukajtis z = r * r; 166*25c28e83SPiotr Jasiukajtis s = p1 * r; 167*25c28e83SPiotr Jasiukajtis return (y + ((p - r) - ((s * z) * 168*25c28e83SPiotr Jasiukajtis (p2 + z * (p3 + z))) * 169*25c28e83SPiotr Jasiukajtis (((p4 + z) + z * z) * 170*25c28e83SPiotr Jasiukajtis (p5 + z * (p6 + z))))); 171*25c28e83SPiotr Jasiukajtis } else if (j < 0x412) { 172*25c28e83SPiotr Jasiukajtis z = r * r; 173*25c28e83SPiotr Jasiukajtis return (y + (p - ((q1 * r) * (q4 + z)) * 174*25c28e83SPiotr Jasiukajtis (q2 + z * (q3 + z)))); 175*25c28e83SPiotr Jasiukajtis } else 176*25c28e83SPiotr Jasiukajtis return (y + (p - r)); 177*25c28e83SPiotr Jasiukajtis } else { 178*25c28e83SPiotr Jasiukajtis if (j >= 0x7ff) /* x is inf or NaN */ 179*25c28e83SPiotr Jasiukajtis if (((ix - 0x7ff00000) | lx) != 0) 180*25c28e83SPiotr Jasiukajtis #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN) 181*25c28e83SPiotr Jasiukajtis return (ix >= 0x7ff80000 ? x : x - x); 182*25c28e83SPiotr Jasiukajtis /* assumes sparc-like QNaN */ 183*25c28e83SPiotr Jasiukajtis #else 184*25c28e83SPiotr Jasiukajtis return (x - x); 185*25c28e83SPiotr Jasiukajtis #endif 186*25c28e83SPiotr Jasiukajtis y = -pio2lo; 187*25c28e83SPiotr Jasiukajtis return (hx >= 0 ? pio2hi - y : y - pio2hi); 188*25c28e83SPiotr Jasiukajtis } 189*25c28e83SPiotr Jasiukajtis } else { /* now x is between 1/8 and 8 */ 190*25c28e83SPiotr Jasiukajtis double *w, w0, w1, s, z; 191*25c28e83SPiotr Jasiukajtis w = (double *) _TBL_atan + (((ix - 0x3fc00000) >> 16) << 1); 192*25c28e83SPiotr Jasiukajtis w0 = (hx >= 0)? w[0] : -w[0]; 193*25c28e83SPiotr Jasiukajtis s = (x - w0) / (one + x * w0); 194*25c28e83SPiotr Jasiukajtis w1 = (hx >= 0)? w[1] : -w[1]; 195*25c28e83SPiotr Jasiukajtis z = s * s; 196*25c28e83SPiotr Jasiukajtis return (((q1 * s) * (q4 + z)) * (q2 + z * (q3 + z)) + w1); 197*25c28e83SPiotr Jasiukajtis } 198*25c28e83SPiotr Jasiukajtis } 199