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>
4 and Fritz Reese <foreese@gcc.gnu.org>
5
6This file is part of the GNU Fortran runtime library (libgfortran).
7
8Libgfortran is free software; you can redistribute it and/or
9modify it under the terms of the GNU General Public
10License as published by the Free Software Foundation; either
11version 3 of the License, or (at your option) any later version.
12
13Libgfortran is distributed in the hope that it will be useful,
14but WITHOUT ANY WARRANTY; without even the implied warranty of
15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16GNU General Public License for more details.
17
18Under Section 7 of GPL version 3, you are granted additional
19permissions described in the GCC Runtime Library Exception, version
203.1, as published by the Free Software Foundation.
21
22You should have received a copy of the GNU General Public License and
23a copy of the GCC Runtime Library Exception along with this program;
24see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
25<http://www.gnu.org/licenses/>. */
26
27
28/*
29
30This file is included from both the FE and the runtime library code.
31Operations are generalized using GMP/MPFR functions. When included from
32libgfortran, these should be overridden using macros which will use native
33operations conforming to the same API. From the FE, the GMP/MPFR functions can
34be used as-is.
35
36The following macros are used and must be defined, unless listed as [optional]:
37
38FTYPE
39 Type name for the real-valued parameter.
40 Variables of this type are constructed/destroyed using mpfr_init()
41 and mpfr_clear.
42
43RETTYPE
44 Return type of the functions.
45
46RETURN(x)
47 Insert code to return a value.
48 The parameter x is the result variable, which was also the input parameter.
49
50ITYPE
51 Type name for integer types.
52
53SIND, COSD, TRIGD
54 Names for the degree-valued trig functions defined by this module.
55
56ENABLE_SIND, ENABLE_COSD, ENABLE_TAND
57 Whether the degree-valued trig functions can be enabled.
58
59ERROR_RETURN(f, k, x)
60 If ENABLE_<xxx>D is not defined, this is substituted to assert an
61 error condition for function f, kind k, and parameter x.
62 The function argument is one of {sind, cosd, tand}.
63
64ISFINITE(x)
65 Whether x is a regular number or zero (not inf or NaN).
66
67D2R(x)
68 Convert x from radians to degrees.
69
70SET_COSD30(x)
71 Set x to COSD(30), or equivalently, SIND(60).
72
73TINY_LITERAL [optional]
74 Value subtracted from 1 to cause raise INEXACT for COSD(x) for x << 1.
75 If not set, COSD(x) for x <= COSD_SMALL_LITERAL simply returns 1.
76
77COSD_SMALL_LITERAL [optional]
78 Value such that x <= COSD_SMALL_LITERAL implies COSD(x) = 1 to within the
79 precision of FTYPE. If not set, this condition is not checked.
80
81SIND_SMALL_LITERAL [optional]
82 Value such that x <= SIND_SMALL_LITERAL implies SIND(x) = D2R(x) to within
83 the precision of FTYPE. If not set, this condition is not checked.
84
85*/
86
87
88#ifdef SIND
89/* Compute sind(x) = sin(x * pi / 180). */
90
91RETTYPE
92SIND (FTYPE x)
93{
94#ifdef ENABLE_SIND
95 if (ISFINITE (x))
96 {
97 FTYPE s, one;
98
99 /* sin(-x) = - sin(x). */
100 mpfr_init (s);
101 mpfr_init_set_ui (one, 1, GFC_RND_MODE);
102 mpfr_copysign (s, one, x, GFC_RND_MODE);
103 mpfr_clear (one);
104
105#ifdef SIND_SMALL_LITERAL
106 /* sin(x) = x as x -> 0; but only for some precisions. */
107 FTYPE ax;
108 mpfr_init (ax);
109 mpfr_abs (ax, x, GFC_RND_MODE);
110 if (mpfr_cmp_ld (ax, SIND_SMALL_LITERAL) < 0)
111 {
112 D2R (x);
113 mpfr_clear (ax);
114 return x;
115 }
116
117 mpfr_swap (x, ax);
118 mpfr_clear (ax);
119
120#else
121 mpfr_abs (x, x, GFC_RND_MODE);
122#endif /* SIND_SMALL_LITERAL */
123
124 /* Reduce angle to x in [0,360]. */
125 FTYPE period;
126 mpfr_init_set_ui (period, 360, GFC_RND_MODE);
127 mpfr_fmod (x, x, period, GFC_RND_MODE);
128 mpfr_clear (period);
129
130 /* Special cases with exact results. */
131 ITYPE n;
132 mpz_init (n);
133 if (mpfr_get_z (z: n, f: x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 30))
134 {
135 /* Flip sign for odd n*pi (x is % 360 so this is only for 180).
136 This respects sgn(sin(x)) = sgn(d/dx sin(x)) = sgn(cos(x)). */
137 if (mpz_divisible_ui_p (n, 180))
138 {
139 mpfr_set_ui (x, 0, GFC_RND_MODE);
140 if (mpz_cmp_ui (n, 180) == 0)
141 mpfr_neg (s, s, GFC_RND_MODE);
142 }
143 else if (mpz_divisible_ui_p (n, 90))
144 mpfr_set_si (x, (mpz_cmp_ui (n, 90) == 0 ? 1 : -1), GFC_RND_MODE);
145 else if (mpz_divisible_ui_p (n, 60))
146 {
147 SET_COSD30 (x);
148 if (mpz_cmp_ui (n, 180) >= 0)
149 mpfr_neg (x, x, GFC_RND_MODE);
150 }
151 else
152 mpfr_set_ld (x, (mpz_cmp_ui (n, 180) < 0 ? 0.5L : -0.5L),
153 GFC_RND_MODE);
154 }
155
156 /* Fold [0,360] into the range [0,45], and compute either SIN() or
157 COS() depending on symmetry of shifting into the [0,45] range. */
158 else
159 {
160 bool fold_cos = false;
161 if (mpfr_cmp_ui (x, 180) <= 0)
162 {
163 if (mpfr_cmp_ui (x, 90) <= 0)
164 {
165 if (mpfr_cmp_ui (x, 45) > 0)
166 {
167 /* x = COS(D2R(90 - x)) */
168 mpfr_ui_sub (x, 90, x, GFC_RND_MODE);
169 fold_cos = true;
170 }
171 }
172 else
173 {
174 if (mpfr_cmp_ui (x, 135) <= 0)
175 {
176 mpfr_sub_ui (x, x, 90, GFC_RND_MODE);
177 fold_cos = true;
178 }
179 else
180 mpfr_ui_sub (x, 180, x, GFC_RND_MODE);
181 }
182 }
183
184 else if (mpfr_cmp_ui (x, 270) <= 0)
185 {
186 if (mpfr_cmp_ui (x, 225) <= 0)
187 mpfr_sub_ui (x, x, 180, GFC_RND_MODE);
188 else
189 {
190 mpfr_ui_sub (x, 270, x, GFC_RND_MODE);
191 fold_cos = true;
192 }
193 mpfr_neg (s, s, GFC_RND_MODE);
194 }
195
196 else
197 {
198 if (mpfr_cmp_ui (x, 315) <= 0)
199 {
200 mpfr_sub_ui (x, x, 270, GFC_RND_MODE);
201 fold_cos = true;
202 }
203 else
204 mpfr_ui_sub (x, 360, x, GFC_RND_MODE);
205 mpfr_neg (s, s, GFC_RND_MODE);
206 }
207
208 D2R (x);
209
210 if (fold_cos)
211 mpfr_cos (x, x, GFC_RND_MODE);
212 else
213 mpfr_sin (x, x, GFC_RND_MODE);
214 }
215
216 mpfr_mul (x, x, s, GFC_RND_MODE);
217 mpz_clear (n);
218 mpfr_clear (s);
219 }
220
221 /* Return NaN for +-Inf and NaN and raise exception. */
222 else
223 mpfr_sub (x, x, x, GFC_RND_MODE);
224
225 RETURN (x);
226
227#else
228 ERROR_RETURN(sind, KIND, x);
229#endif // ENABLE_SIND
230}
231#endif // SIND
232
233
234#ifdef COSD
235/* Compute cosd(x) = cos(x * pi / 180). */
236
237RETTYPE
238COSD (FTYPE x)
239{
240#ifdef ENABLE_COSD
241#if defined(TINY_LITERAL) && defined(COSD_SMALL_LITERAL)
242 static const volatile FTYPE tiny = TINY_LITERAL;
243#endif
244
245 if (ISFINITE (x))
246 {
247#ifdef COSD_SMALL_LITERAL
248 FTYPE ax;
249 mpfr_init (ax);
250
251 mpfr_abs (ax, x, GFC_RND_MODE);
252 /* No spurious underflows!. In radians, cos(x) = 1-x*x/2 as x -> 0. */
253 if (mpfr_cmp_ld (ax, COSD_SMALL_LITERAL) <= 0)
254 {
255 mpfr_set_ui (x, 1, GFC_RND_MODE);
256#ifdef TINY_LITERAL
257 /* Cause INEXACT. */
258 if (!mpfr_zero_p (ax))
259 mpfr_sub_d (x, x, tiny, GFC_RND_MODE);
260#endif
261
262 mpfr_clear (ax);
263 return x;
264 }
265
266 mpfr_swap (x, ax);
267 mpfr_clear (ax);
268#else
269 mpfr_abs (x, x, GFC_RND_MODE);
270#endif /* COSD_SMALL_LITERAL */
271
272 /* Reduce angle to ax in [0,360]. */
273 FTYPE period;
274 mpfr_init_set_ui (period, 360, GFC_RND_MODE);
275 mpfr_fmod (x, x, period, GFC_RND_MODE);
276 mpfr_clear (period);
277
278 /* Special cases with exact results.
279 Return negative zero for cosd(270) for consistency with libm cos(). */
280 ITYPE n;
281 mpz_init (n);
282 if (mpfr_get_z (z: n, f: x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 30))
283 {
284 if (mpz_divisible_ui_p (n, 180))
285 mpfr_set_si (x, (mpz_cmp_ui (n, 180) == 0 ? -1 : 1),
286 GFC_RND_MODE);
287 else if (mpz_divisible_ui_p (n, 90))
288 mpfr_set_zero (x, 0);
289 else if (mpz_divisible_ui_p (n, 60))
290 {
291 mpfr_set_ld (x, 0.5, GFC_RND_MODE);
292 if (mpz_cmp_ui (n, 60) != 0 && mpz_cmp_ui (n, 300) != 0)
293 mpfr_neg (x, x, GFC_RND_MODE);
294 }
295 else
296 {
297 SET_COSD30 (x);
298 if (mpz_cmp_ui (n, 30) != 0 && mpz_cmp_ui (n, 330) != 0)
299 mpfr_neg (x, x, GFC_RND_MODE);
300 }
301 }
302
303 /* Fold [0,360] into the range [0,45], and compute either SIN() or
304 COS() depending on symmetry of shifting into the [0,45] range. */
305 else
306 {
307 bool neg = false;
308 bool fold_sin = false;
309 if (mpfr_cmp_ui (x, 180) <= 0)
310 {
311 if (mpfr_cmp_ui (x, 90) <= 0)
312 {
313 if (mpfr_cmp_ui (x, 45) > 0)
314 {
315 mpfr_ui_sub (x, 90, x, GFC_RND_MODE);
316 fold_sin = true;
317 }
318 }
319 else
320 {
321 if (mpfr_cmp_ui (x, 135) <= 0)
322 {
323 mpfr_sub_ui (x, x, 90, GFC_RND_MODE);
324 fold_sin = true;
325 }
326 else
327 mpfr_ui_sub (x, 180, x, GFC_RND_MODE);
328 neg = true;
329 }
330 }
331
332 else if (mpfr_cmp_ui (x, 270) <= 0)
333 {
334 if (mpfr_cmp_ui (x, 225) <= 0)
335 mpfr_sub_ui (x, x, 180, GFC_RND_MODE);
336 else
337 {
338 mpfr_ui_sub (x, 270, x, GFC_RND_MODE);
339 fold_sin = true;
340 }
341 neg = true;
342 }
343
344 else
345 {
346 if (mpfr_cmp_ui (x, 315) <= 0)
347 {
348 mpfr_sub_ui (x, x, 270, GFC_RND_MODE);
349 fold_sin = true;
350 }
351 else
352 mpfr_ui_sub (x, 360, x, GFC_RND_MODE);
353 }
354
355 D2R (x);
356
357 if (fold_sin)
358 mpfr_sin (x, x, GFC_RND_MODE);
359 else
360 mpfr_cos (x, x, GFC_RND_MODE);
361
362 if (neg)
363 mpfr_neg (x, x, GFC_RND_MODE);
364 }
365
366 mpz_clear (n);
367 }
368
369 /* Return NaN for +-Inf and NaN and raise exception. */
370 else
371 mpfr_sub (x, x, x, GFC_RND_MODE);
372
373 RETURN (x);
374
375#else
376 ERROR_RETURN(cosd, KIND, x);
377#endif // ENABLE_COSD
378}
379#endif // COSD
380
381
382#ifdef TAND
383/* Compute tand(x) = tan(x * pi / 180). */
384
385RETTYPE
386TAND (FTYPE x)
387{
388#ifdef ENABLE_TAND
389 if (ISFINITE (x))
390 {
391 FTYPE s, one;
392
393 /* tan(-x) = - tan(x). */
394 mpfr_init (s);
395 mpfr_init_set_ui (one, 1, GFC_RND_MODE);
396 mpfr_copysign (s, one, x, GFC_RND_MODE);
397 mpfr_clear (one);
398
399#ifdef SIND_SMALL_LITERAL
400 /* tan(x) = x as x -> 0; but only for some precisions. */
401 FTYPE ax;
402 mpfr_init (ax);
403 mpfr_abs (ax, x, GFC_RND_MODE);
404 if (mpfr_cmp_ld (ax, SIND_SMALL_LITERAL) < 0)
405 {
406 D2R (x);
407 mpfr_clear (ax);
408 return x;
409 }
410
411 mpfr_swap (x, ax);
412 mpfr_clear (ax);
413
414#else
415 mpfr_abs (x, x, GFC_RND_MODE);
416#endif /* SIND_SMALL_LITERAL */
417
418 /* Reduce angle to x in [0,360]. */
419 FTYPE period;
420 mpfr_init_set_ui (period, 360, GFC_RND_MODE);
421 mpfr_fmod (x, x, period, GFC_RND_MODE);
422 mpfr_clear (period);
423
424 /* Special cases with exact results. */
425 ITYPE n;
426 mpz_init (n);
427 if (mpfr_get_z (z: n, f: x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 45))
428 {
429 if (mpz_divisible_ui_p (n, 180))
430 mpfr_set_zero (x, 0);
431
432 /* Though mathematically NaN is more appropriate for tan(n*90),
433 returning +/-Inf offers the advantage that 1/tan(n*90) returns 0,
434 which is mathematically sound. In fact we rely on this behavior
435 to implement COTAND(x) = 1 / TAND(x).
436 */
437 else if (mpz_divisible_ui_p (n, 90))
438 mpfr_set_inf (x, mpz_cmp_ui (n, 90) == 0 ? 0 : 1);
439
440 else
441 {
442 mpfr_set_ui (x, 1, GFC_RND_MODE);
443 if (mpz_cmp_ui (n, 45) != 0 && mpz_cmp_ui (n, 225) != 0)
444 mpfr_neg (x, x, GFC_RND_MODE);
445 }
446 }
447
448 else
449 {
450 /* Fold [0,360] into the range [0,90], and compute TAN(). */
451 if (mpfr_cmp_ui (x, 180) <= 0)
452 {
453 if (mpfr_cmp_ui (x, 90) > 0)
454 {
455 mpfr_ui_sub (x, 180, x, GFC_RND_MODE);
456 mpfr_neg (s, s, GFC_RND_MODE);
457 }
458 }
459 else
460 {
461 if (mpfr_cmp_ui (x, 270) <= 0)
462 {
463 mpfr_sub_ui (x, x, 180, GFC_RND_MODE);
464 }
465 else
466 {
467 mpfr_ui_sub (x, 360, x, GFC_RND_MODE);
468 mpfr_neg (s, s, GFC_RND_MODE);
469 }
470 }
471
472 D2R (x);
473 mpfr_tan (x, x, GFC_RND_MODE);
474 }
475
476 mpfr_mul (x, x, s, GFC_RND_MODE);
477 mpz_clear (n);
478 mpfr_clear (s);
479 }
480
481 /* Return NaN for +-Inf and NaN and raise exception. */
482 else
483 mpfr_sub (x, x, x, GFC_RND_MODE);
484
485 RETURN (x);
486
487#else
488 ERROR_RETURN(tand, KIND, x);
489#endif // ENABLE_TAND
490}
491#endif // TAND
492
493/* vim: set ft=c: */
494

source code of libgfortran/intrinsics/trigd.inc