1 /* Implementation of the degree trignometric functions COSD, SIND, TAND.
2 Copyright (C) 2020-2024 Free Software Foundation, Inc.
3 Contributed by Steven G. Kargl <kargl@gcc.gnu.org>
5 This file is part of the GNU Fortran runtime library (libgfortran).
7 Libgfortran is free software; you can redistribute it and/or
8 modify it under the terms of the GNU General Public
9 License as published by the Free Software Foundation; either
10 version 3 of the License, or (at your option) any later version.
12 Libgfortran is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
17 Under Section 7 of GPL version 3, you are granted additional
18 permissions described in the GCC Runtime Library Exception, version
19 3.1, as published by the Free Software Foundation.
21 You should have received a copy of the GNU General Public License and
22 a copy of the GCC Runtime Library Exception along with this program;
23 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
24 <http://www.gnu.org/licenses/>. */
26 #include "libgfortran.h"
30 /* Body of library functions which are cannot be implemented on the current
31 * platform because it lacks a capability, such as an underlying trigonometric
32 * function (sin, cos, tan) or C99 floating-point function (fabs, fmod). */
33 #define STRINGIFY_EXPAND(x) #x
34 #define ERROR_RETURN(f, k, x) runtime_error (#f " is unavailable for" \
35 " REAL(KIND=" STRINGIFY_EXPAND(k) ") because the system math library" \
36 " lacks support for it"); \
40 For real x, let {x}_P or x_P be the closest representible number in the
41 floating point representation which uses P binary bits of fractional
42 precision (with IEEE rounding semantics).
44 Similarly, let f_P(x) be shorthand for {f(x)}_P.
46 Let ulp_P(x) be the unit of least precision for x: in other words the
47 maximal value of |a_P - b_P| where a_P <= x <= b_P and a_P != b_P.
49 Let x ~= y <-> | x - y | < ulp_P(x - y).
51 Let deg(x) be the value of x radians in degrees.
53 Values for each precision P were selected as follows.
56 COSD_SMALL = 2**{-N} such that for all x <= COSD_SMALL:
58 * cos(deg(x)) ~= 1, or equivalently:
60 | 1 - cos(deg(x)) | < ulp_P(1).
62 Unfortunately for SIND (and therefore TAND) a similar relation is only
63 possible for REAL(4) and REAL(8). With REAL(10) and REAL(16), enough
64 precision is available such that sin_P(x) != x_P for some x less than any
65 value. (There are values where this equality holds, but the distance has
68 For REAL(4) and REAL(8), we can select SIND_SMALL such that:
70 * sin(deg(x)) ~= deg(x), or equivalently:
72 | deg(x) - sin(deg(x)) | < ulp_P(deg(x)).
76 #ifdef HAVE_GFC_REAL_4
78 /* Build _gfortran_sind_r4, _gfortran_cosd_r4, and _gfortran_tand_r4 */
81 #define TINY 0x1.p-100 /* ~= 7.889e-31 */
82 #define COSD_SMALL 0x1.p-7 /* = 7.8125e-3 */
83 #define SIND_SMALL 0x1.p-5 /* = 3.125e-2 */
84 #define COSD30 8.66025388e-01
85 #define PIO180H 1.74560547e-02 /* high 12 bits. */
86 #define PIO180L -2.76216747e-06 /* Next 24 bits. */
88 #if defined(HAVE_FABSF) && defined(HAVE_FMODF) && defined(HAVE_COPYSIGNF)
102 #endif /* HAVE_FABSF && HAVE_FMODF && HAVE_COPYSIGNF */
104 #ifdef GFC_REAL_4_INFINITY
105 #define HAVE_INFINITY_KIND
108 #include "trigd_lib.inc"
120 #undef HAVE_INFINITY_KIND
122 #endif /* HAVE_GFC_REAL_4... */
125 #ifdef HAVE_GFC_REAL_8
127 /* Build _gfortran_sind_r8, _gfortran_cosd_r8, and _gfortran_tand_r8 */
130 #define TINY 0x1.p-1000 /* ~= 9.33e-302 (min exp -1074) */
131 #define COSD_SMALL 0x1.p-21 /* ~= 4.768e-7 */
132 #define SIND_SMALL 0x1.p-19 /* ~= 9.537e-7 */
133 #define COSD30 8.6602540378443860e-01
134 #define PIO180H 1.7453283071517944e-02 /* high 21 bits. */
135 #define PIO180L 9.4484253514332993e-09 /* Next 53 bits. */
137 #if defined(HAVE_FABS) && defined(HAVE_FMOD) && defined(HAVE_COPYSIGN)
151 #endif /* HAVE_FABS && HAVE_FMOD && HAVE_COPYSIGN */
153 #ifdef GFC_REAL_8_INFINITY
154 #define HAVE_INFINITY_KIND
157 #include "trigd_lib.inc"
169 #undef HAVE_INFINITY_KIND
171 #endif /* HAVE_GFC_REAL_8... */
174 #ifdef HAVE_GFC_REAL_10
176 /* Build _gfortran_sind_r10, _gfortran_cosd_r10, and _gfortran_tand_r10 */
179 #define TINY 0x1.p-16400 /* ~= 1.28e-4937 (min exp -16494) */
180 #define COSD_SMALL 0x1.p-26 /* ~= 1.490e-8 */
181 #undef SIND_SMALL /* not precise */
182 #define COSD30 8.66025403784438646787e-01
183 #define PIO180H 1.74532925229868851602e-02 /* high 32 bits */
184 #define PIO180L -3.04358939097084072823e-12 /* Next 64 bits */
186 #if defined(HAVE_FABSL) && defined(HAVE_FMODL) && defined(HAVE_COPYSIGNL)
200 #endif /* HAVE_FABSL && HAVE_FMODL && HAVE_COPYSIGNL */
202 #ifdef GFC_REAL_10_INFINITY
203 #define HAVE_INFINITY_KIND
206 #include "trigd_lib.inc"
218 #undef HAVE_INFINITY_KIND
220 #endif /* HAVE_GFC_REAL_10 */
223 #ifdef HAVE_GFC_REAL_16
225 /* Build _gfortran_sind_r16, _gfortran_cosd_r16, and _gfortran_tand_r16 */
228 #define TINY 0x1.p-16400 /* ~= 1.28e-4937 */
229 #undef SIND_SMALL /* not precise */
231 #if GFC_REAL_16_DIGITS == 64
232 /* 80 bit precision, use constants from REAL(10). */
233 #define COSD_SMALL 0x1.p-26 /* ~= 1.490e-8 */
234 #define COSD30 8.66025403784438646787e-01
235 #define PIO180H 1.74532925229868851602e-02 /* high 32 bits */
236 #define PIO180L -3.04358939097084072823e-12 /* Next 64 bits */
239 /* Proper float128 precision. */
240 #define COSD_SMALL 0x1.p-51 /* ~= 4.441e-16 */
241 #define COSD30 8.66025403784438646763723170752936183e-01
242 #define PIO180H 1.74532925199433197605003442731685936e-02
243 #define PIO180L -2.39912634365882824665106671063098954e-17
246 #ifdef GFC_REAL_16_IS_LONG_DOUBLE
248 #if defined(HAVE_FABSL) && defined(HAVE_FMODL) && defined(HAVE_COPYSIGNL)
262 #endif /* HAVE_FABSL && HAVE_FMODL && HAVE_COPYSIGNL */
264 #elif defined(GFC_REAL_16_USE_IEC_60559)
266 #if defined(HAVE_FABSF128) && defined(HAVE_FMODF128) && defined(HAVE_COPYSIGNF128)
280 #endif /* HAVE_FABSF128 && HAVE_FMODF128 && HAVE_COPYSIGNF128 */
284 /* libquadmath: HAVE_*Q are never defined. They must be available. */
289 #endif /* GFC_REAL_16_IS_LONG_DOUBLE */
291 #ifdef GFC_REAL_16_INFINITY
292 #define HAVE_INFINITY_KIND
295 #include "trigd_lib.inc"
307 #undef HAVE_INFINITY_KIND
309 #endif /* HAVE_GFC_REAL_16 */
311 #ifdef HAVE_GFC_REAL_17
313 /* Build _gfortran_sind_r17, _gfortran_cosd_r17, and _gfortran_tand_r17 */
316 #define TINY 0x1.p-16400 /* ~= 1.28e-4937 */
317 #undef SIND_SMALL /* not precise */
319 /* Proper float128 precision. */
320 #define COSD_SMALL 0x1.p-51 /* ~= 4.441e-16 */
321 #define COSD30 8.66025403784438646763723170752936183e-01
322 #define PIO180H 1.74532925199433197605003442731685936e-02
323 #define PIO180L -2.39912634365882824665106671063098954e-17
325 /* libquadmath or glibc 2.32+: HAVE_*Q are never defined. They must be available. */
330 #ifdef GFC_REAL_17_INFINITY
331 #define HAVE_INFINITY_KIND
335 #define COPYSIGN __copysignieee128
336 #define FMOD __fmodieee128
337 #define FABS __fabsieee128
338 #define FMA __fmaieee128
339 #define SIN __sinieee128
340 #define COS __cosieee128
341 #define TAN __tanieee128
344 #include "trigd_lib.inc"
356 #undef HAVE_INFINITY_KIND
358 #endif /* HAVE_GFC_REAL_17 */