1 /* $NetBSD: s_atanl.c,v 1.8 2024/07/17 12:00:48 riastradh Exp $ */
2
3 /* FreeBSD: head/lib/msun/src/s_atan.c 176451 2008-02-22 02:30:36Z das */
4 /*
5 * ====================================================
6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
7 *
8 * Developed at SunPro, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
11 * is preserved.
12 * ====================================================
13 */
14
15 #include <sys/cdefs.h>
16 __RCSID("$NetBSD: s_atanl.c,v 1.8 2024/07/17 12:00:48 riastradh Exp $");
17
18 #include "namespace.h"
19
20 #include <float.h>
21 #include <machine/ieee.h>
22
23 #include "math.h"
24 #include "math_private.h"
25
26 #ifdef __HAVE_LONG_DOUBLE
27
28 /*
29 * See comments in s_atan.c.
30 * Converted to long double by David Schultz <das@FreeBSD.ORG>.
31 */
32
33 __weak_alias(atanl, _atanl)
34
35 #if LDBL_MANT_DIG == 64
36 #include "../ld80/invtrig.h"
37 #elif LDBL_MANT_DIG == 113
38 #include "../ld128/invtrig.h"
39 #else
40 #error "Unsupported long double format"
41 #endif
42
43 #ifdef LDBL_IMPLICIT_NBIT
44 #define LDBL_NBIT 0
45 #endif
46
47 static const long double
48 one = 1.0,
49 huge = 1.0e300;
50
51 long double
atanl(long double x)52 atanl(long double x)
53 {
54 union ieee_ext_u u;
55 long double w,s1,s2,z;
56 int id;
57 int16_t expsign, expt;
58 int32_t expman;
59
60 u.extu_ld = x;
61 expsign = GET_EXPSIGN(&u);
62 expt = expsign & 0x7fff;
63 if(expt >= ATAN_CONST) { /* if |x| is large, atan(x)~=pi/2 */
64 if(expt == BIAS + LDBL_MAX_EXP &&
65 ((u.extu_frach&~LDBL_NBIT)|u.extu_fracl)!=0)
66 return x+x; /* NaN */
67 if(expsign>0) return atanhi[3]+atanlo[3];
68 else return -atanhi[3]-atanlo[3];
69 }
70 /* Extract the exponent and the first few bits of the significand. */
71 /* XXX There should be a more convenient way to do this. */
72 expman = (expt << 8) | ((u.extu_frach >> (MANH_SIZE - 9)) & 0xff);
73 if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
74 if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */
75 if(huge+x>one) return x; /* raise inexact */
76 }
77 id = -1;
78 } else {
79 x = fabsl(x);
80 if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */
81 if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <=|x|<11/16 */
82 id = 0; x = (2.0*x-one)/(2.0+x);
83 } else { /* 11/16<=|x|< 19/16 */
84 id = 1; x = (x-one)/(x+one);
85 }
86 } else {
87 if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */
88 id = 2; x = (x-1.5)/(one+1.5*x);
89 } else { /* 2.4375 <= |x| < 2^ATAN_CONST */
90 id = 3; x = -1.0/x;
91 }
92 }}
93 /* end of argument reduction */
94 z = x*x;
95 w = z*z;
96 /* break sum aT[i]z**(i+1) into odd and even poly */
97 s1 = z*T_even(w);
98 s2 = w*T_odd(w);
99 if (id<0) return x - x*(s1+s2);
100 else {
101 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
102 return (expsign<0)? -z:z;
103 }
104 }
105 #endif
106