Branch data Line data Source code
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 : :
6 : : This file is part of the GNU Fortran runtime library (libgfortran).
7 : :
8 : : Libgfortran is free software; you can redistribute it and/or
9 : : modify it under the terms of the GNU General Public
10 : : License as published by the Free Software Foundation; either
11 : : version 3 of the License, or (at your option) any later version.
12 : :
13 : : Libgfortran is distributed in the hope that it will be useful,
14 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
15 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 : : GNU General Public License for more details.
17 : :
18 : : Under Section 7 of GPL version 3, you are granted additional
19 : : permissions described in the GCC Runtime Library Exception, version
20 : : 3.1, as published by the Free Software Foundation.
21 : :
22 : : You should have received a copy of the GNU General Public License and
23 : : a copy of the GCC Runtime Library Exception along with this program;
24 : : see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
25 : : <http://www.gnu.org/licenses/>. */
26 : :
27 : :
28 : : /*
29 : :
30 : : This file is included from both the FE and the runtime library code.
31 : : Operations are generalized using GMP/MPFR functions. When included from
32 : : libgfortran, these should be overridden using macros which will use native
33 : : operations conforming to the same API. From the FE, the GMP/MPFR functions can
34 : : be used as-is.
35 : :
36 : : The following macros are used and must be defined, unless listed as [optional]:
37 : :
38 : : FTYPE
39 : : Type name for the real-valued parameter.
40 : : Variables of this type are constructed/destroyed using mpfr_init()
41 : : and mpfr_clear.
42 : :
43 : : RETTYPE
44 : : Return type of the functions.
45 : :
46 : : RETURN(x)
47 : : Insert code to return a value.
48 : : The parameter x is the result variable, which was also the input parameter.
49 : :
50 : : ITYPE
51 : : Type name for integer types.
52 : :
53 : : SIND, COSD, TRIGD
54 : : Names for the degree-valued trig functions defined by this module.
55 : :
56 : : ENABLE_SIND, ENABLE_COSD, ENABLE_TAND
57 : : Whether the degree-valued trig functions can be enabled.
58 : :
59 : : ERROR_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 : :
64 : : ISFINITE(x)
65 : : Whether x is a regular number or zero (not inf or NaN).
66 : :
67 : : D2R(x)
68 : : Convert x from radians to degrees.
69 : :
70 : : SET_COSD30(x)
71 : : Set x to COSD(30), or equivalently, SIND(60).
72 : :
73 : : TINY_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 : :
77 : : COSD_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 : :
81 : : SIND_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 : :
91 : : RETTYPE
92 : 25 : SIND (FTYPE x)
93 : : {
94 : : #ifdef ENABLE_SIND
95 : 25 : if (ISFINITE (x))
96 : : {
97 : 25 : FTYPE s, one;
98 : :
99 : : /* sin(-x) = - sin(x). */
100 : 25 : mpfr_init (s);
101 : 25 : mpfr_init_set_ui (one, 1, GFC_RND_MODE);
102 : 25 : mpfr_copysign (s, one, x, GFC_RND_MODE);
103 : 25 : 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 : 25 : mpfr_abs (x, x, GFC_RND_MODE);
122 : : #endif /* SIND_SMALL_LITERAL */
123 : :
124 : : /* Reduce angle to x in [0,360]. */
125 : 25 : FTYPE period;
126 : 25 : mpfr_init_set_ui (period, 360, GFC_RND_MODE);
127 : 25 : mpfr_fmod (x, x, period, GFC_RND_MODE);
128 : 25 : mpfr_clear (period);
129 : :
130 : : /* Special cases with exact results. */
131 : 25 : ITYPE n;
132 : 25 : mpz_init (n);
133 : 25 : if (mpfr_get_z (n, 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 : 25 : if (mpz_divisible_ui_p (n, 180))
138 : : {
139 : 1 : mpfr_set_ui (x, 0, GFC_RND_MODE);
140 : 1 : if (mpz_cmp_ui (n, 180) == 0)
141 : 0 : mpfr_neg (s, s, GFC_RND_MODE);
142 : : }
143 : 24 : else if (mpz_divisible_ui_p (n, 90))
144 : 0 : mpfr_set_si (x, (mpz_cmp_ui (n, 90) == 0 ? 1 : -1), GFC_RND_MODE);
145 : 24 : else if (mpz_divisible_ui_p (n, 60))
146 : : {
147 : 24 : SET_COSD30 (x);
148 : 24 : if (mpz_cmp_ui (n, 180) >= 0)
149 : 0 : mpfr_neg (x, x, GFC_RND_MODE);
150 : : }
151 : : else
152 : 0 : 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 : 0 : bool fold_cos = false;
161 : 0 : if (mpfr_cmp_ui (x, 180) <= 0)
162 : : {
163 : 0 : if (mpfr_cmp_ui (x, 90) <= 0)
164 : : {
165 : 0 : if (mpfr_cmp_ui (x, 45) > 0)
166 : : {
167 : : /* x = COS(D2R(90 - x)) */
168 : 0 : mpfr_ui_sub (x, 90, x, GFC_RND_MODE);
169 : 0 : fold_cos = true;
170 : : }
171 : : }
172 : : else
173 : : {
174 : 0 : if (mpfr_cmp_ui (x, 135) <= 0)
175 : : {
176 : 0 : mpfr_sub_ui (x, x, 90, GFC_RND_MODE);
177 : 0 : fold_cos = true;
178 : : }
179 : : else
180 : 0 : mpfr_ui_sub (x, 180, x, GFC_RND_MODE);
181 : : }
182 : : }
183 : :
184 : 0 : else if (mpfr_cmp_ui (x, 270) <= 0)
185 : : {
186 : 0 : if (mpfr_cmp_ui (x, 225) <= 0)
187 : 0 : mpfr_sub_ui (x, x, 180, GFC_RND_MODE);
188 : : else
189 : : {
190 : 0 : mpfr_ui_sub (x, 270, x, GFC_RND_MODE);
191 : 0 : fold_cos = true;
192 : : }
193 : 0 : mpfr_neg (s, s, GFC_RND_MODE);
194 : : }
195 : :
196 : : else
197 : : {
198 : 0 : if (mpfr_cmp_ui (x, 315) <= 0)
199 : : {
200 : 0 : mpfr_sub_ui (x, x, 270, GFC_RND_MODE);
201 : 0 : fold_cos = true;
202 : : }
203 : : else
204 : 0 : mpfr_ui_sub (x, 360, x, GFC_RND_MODE);
205 : 0 : mpfr_neg (s, s, GFC_RND_MODE);
206 : : }
207 : :
208 : 0 : D2R (x);
209 : :
210 : 0 : if (fold_cos)
211 : 0 : mpfr_cos (x, x, GFC_RND_MODE);
212 : : else
213 : 0 : mpfr_sin (x, x, GFC_RND_MODE);
214 : : }
215 : :
216 : 25 : mpfr_mul (x, x, s, GFC_RND_MODE);
217 : 25 : mpz_clear (n);
218 : 25 : mpfr_clear (s);
219 : : }
220 : :
221 : : /* Return NaN for +-Inf and NaN and raise exception. */
222 : : else
223 : 0 : mpfr_sub (x, x, x, GFC_RND_MODE);
224 : :
225 : 25 : RETURN (x);
226 : :
227 : : #else
228 : : ERROR_RETURN(sind, KIND, x);
229 : : #endif // ENABLE_SIND
230 : 25 : }
231 : : #endif // SIND
232 : :
233 : :
234 : : #ifdef COSD
235 : : /* Compute cosd(x) = cos(x * pi / 180). */
236 : :
237 : : RETTYPE
238 : 25 : COSD (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 : 25 : 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 : 25 : mpfr_abs (x, x, GFC_RND_MODE);
270 : : #endif /* COSD_SMALL_LITERAL */
271 : :
272 : : /* Reduce angle to ax in [0,360]. */
273 : 25 : FTYPE period;
274 : 25 : mpfr_init_set_ui (period, 360, GFC_RND_MODE);
275 : 25 : mpfr_fmod (x, x, period, GFC_RND_MODE);
276 : 25 : mpfr_clear (period);
277 : :
278 : : /* Special cases with exact results.
279 : : Return negative zero for cosd(270) for consistency with libm cos(). */
280 : 25 : ITYPE n;
281 : 25 : mpz_init (n);
282 : 25 : if (mpfr_get_z (n, x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 30))
283 : : {
284 : 25 : if (mpz_divisible_ui_p (n, 180))
285 : 2 : mpfr_set_si (x, (mpz_cmp_ui (n, 180) == 0 ? -1 : 1),
286 : : GFC_RND_MODE);
287 : 24 : else if (mpz_divisible_ui_p (n, 90))
288 : 0 : mpfr_set_zero (x, 0);
289 : 24 : else if (mpz_divisible_ui_p (n, 60))
290 : : {
291 : 24 : mpfr_set_ld (x, 0.5, GFC_RND_MODE);
292 : 24 : if (mpz_cmp_ui (n, 60) != 0 && mpz_cmp_ui (n, 300) != 0)
293 : 0 : mpfr_neg (x, x, GFC_RND_MODE);
294 : : }
295 : : else
296 : : {
297 : 0 : SET_COSD30 (x);
298 : 0 : if (mpz_cmp_ui (n, 30) != 0 && mpz_cmp_ui (n, 330) != 0)
299 : 0 : 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 : 0 : bool neg = false;
308 : 0 : bool fold_sin = false;
309 : 0 : if (mpfr_cmp_ui (x, 180) <= 0)
310 : : {
311 : 0 : if (mpfr_cmp_ui (x, 90) <= 0)
312 : : {
313 : 0 : if (mpfr_cmp_ui (x, 45) > 0)
314 : : {
315 : 0 : mpfr_ui_sub (x, 90, x, GFC_RND_MODE);
316 : 0 : fold_sin = true;
317 : : }
318 : : }
319 : : else
320 : : {
321 : 0 : if (mpfr_cmp_ui (x, 135) <= 0)
322 : : {
323 : 0 : mpfr_sub_ui (x, x, 90, GFC_RND_MODE);
324 : 0 : fold_sin = true;
325 : : }
326 : : else
327 : 0 : mpfr_ui_sub (x, 180, x, GFC_RND_MODE);
328 : : neg = true;
329 : : }
330 : : }
331 : :
332 : 0 : else if (mpfr_cmp_ui (x, 270) <= 0)
333 : : {
334 : 0 : if (mpfr_cmp_ui (x, 225) <= 0)
335 : 0 : mpfr_sub_ui (x, x, 180, GFC_RND_MODE);
336 : : else
337 : : {
338 : 0 : mpfr_ui_sub (x, 270, x, GFC_RND_MODE);
339 : 0 : fold_sin = true;
340 : : }
341 : : neg = true;
342 : : }
343 : :
344 : : else
345 : : {
346 : 0 : if (mpfr_cmp_ui (x, 315) <= 0)
347 : : {
348 : 0 : mpfr_sub_ui (x, x, 270, GFC_RND_MODE);
349 : 0 : fold_sin = true;
350 : : }
351 : : else
352 : 0 : mpfr_ui_sub (x, 360, x, GFC_RND_MODE);
353 : : }
354 : :
355 : 0 : D2R (x);
356 : :
357 : 0 : if (fold_sin)
358 : 0 : mpfr_sin (x, x, GFC_RND_MODE);
359 : : else
360 : 0 : mpfr_cos (x, x, GFC_RND_MODE);
361 : :
362 : 0 : if (neg)
363 : 0 : mpfr_neg (x, x, GFC_RND_MODE);
364 : : }
365 : :
366 : 25 : mpz_clear (n);
367 : : }
368 : :
369 : : /* Return NaN for +-Inf and NaN and raise exception. */
370 : : else
371 : 0 : mpfr_sub (x, x, x, GFC_RND_MODE);
372 : :
373 : 25 : RETURN (x);
374 : :
375 : : #else
376 : : ERROR_RETURN(cosd, KIND, x);
377 : : #endif // ENABLE_COSD
378 : 25 : }
379 : : #endif // COSD
380 : :
381 : :
382 : : #ifdef TAND
383 : : /* Compute tand(x) = tan(x * pi / 180). */
384 : :
385 : : RETTYPE
386 : 50 : TAND (FTYPE x)
387 : : {
388 : : #ifdef ENABLE_TAND
389 : 50 : if (ISFINITE (x))
390 : : {
391 : 50 : FTYPE s, one;
392 : :
393 : : /* tan(-x) = - tan(x). */
394 : 50 : mpfr_init (s);
395 : 50 : mpfr_init_set_ui (one, 1, GFC_RND_MODE);
396 : 50 : mpfr_copysign (s, one, x, GFC_RND_MODE);
397 : 50 : 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 : 50 : mpfr_abs (x, x, GFC_RND_MODE);
416 : : #endif /* SIND_SMALL_LITERAL */
417 : :
418 : : /* Reduce angle to x in [0,360]. */
419 : 50 : FTYPE period;
420 : 50 : mpfr_init_set_ui (period, 360, GFC_RND_MODE);
421 : 50 : mpfr_fmod (x, x, period, GFC_RND_MODE);
422 : 50 : mpfr_clear (period);
423 : :
424 : : /* Special cases with exact results. */
425 : 50 : ITYPE n;
426 : 50 : mpz_init (n);
427 : 50 : if (mpfr_get_z (n, x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 45))
428 : : {
429 : 2 : if (mpz_divisible_ui_p (n, 180))
430 : 2 : 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 : 0 : else if (mpz_divisible_ui_p (n, 90))
438 : 0 : mpfr_set_inf (x, mpz_cmp_ui (n, 90) == 0 ? 0 : 1);
439 : :
440 : : else
441 : : {
442 : 0 : mpfr_set_ui (x, 1, GFC_RND_MODE);
443 : 0 : if (mpz_cmp_ui (n, 45) != 0 && mpz_cmp_ui (n, 225) != 0)
444 : 0 : 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 : 48 : if (mpfr_cmp_ui (x, 180) <= 0)
452 : : {
453 : 48 : if (mpfr_cmp_ui (x, 90) > 0)
454 : : {
455 : 24 : mpfr_ui_sub (x, 180, x, GFC_RND_MODE);
456 : 24 : mpfr_neg (s, s, GFC_RND_MODE);
457 : : }
458 : : }
459 : : else
460 : : {
461 : 0 : if (mpfr_cmp_ui (x, 270) <= 0)
462 : : {
463 : 0 : mpfr_sub_ui (x, x, 180, GFC_RND_MODE);
464 : : }
465 : : else
466 : : {
467 : 0 : mpfr_ui_sub (x, 360, x, GFC_RND_MODE);
468 : 0 : mpfr_neg (s, s, GFC_RND_MODE);
469 : : }
470 : : }
471 : :
472 : 48 : D2R (x);
473 : 48 : mpfr_tan (x, x, GFC_RND_MODE);
474 : : }
475 : :
476 : 50 : mpfr_mul (x, x, s, GFC_RND_MODE);
477 : 50 : mpz_clear (n);
478 : 50 : mpfr_clear (s);
479 : : }
480 : :
481 : : /* Return NaN for +-Inf and NaN and raise exception. */
482 : : else
483 : 0 : mpfr_sub (x, x, x, GFC_RND_MODE);
484 : :
485 : 50 : RETURN (x);
486 : :
487 : : #else
488 : : ERROR_RETURN(tand, KIND, x);
489 : : #endif // ENABLE_TAND
490 : 50 : }
491 : : #endif // TAND
492 : :
493 : : /* vim: set ft=c: */
|
