]>
Commit | Line | Data |
---|---|---|
57391dda | 1 | /* Implementation of the degree trignometric functions COSD, SIND, TAND. |
83ffe9cd | 2 | Copyright (C) 2020-2023 Free Software Foundation, Inc. |
57391dda FR |
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 | ||
e8eecc2a | 36 | The following macros are used and must be defined, unless listed as [optional]: |
57391dda FR |
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 | ||
e8eecc2a FR |
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}. | |
57391dda | 63 | |
e8eecc2a FR |
64 | ISFINITE(x) |
65 | Whether x is a regular number or zero (not inf or NaN). | |
57391dda | 66 | |
e8eecc2a FR |
67 | D2R(x) |
68 | Convert x from radians to degrees. | |
57391dda | 69 | |
e8eecc2a FR |
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 | |
57391dda FR |
79 | precision of FTYPE. If not set, this condition is not checked. |
80 | ||
e8eecc2a FR |
81 | SIND_SMALL_LITERAL [optional] |
82 | Value such that x <= SIND_SMALL_LITERAL implies SIND(x) = D2R(x) to within | |
57391dda FR |
83 | the precision of FTYPE. If not set, this condition is not checked. |
84 | ||
57391dda FR |
85 | */ |
86 | ||
87 | ||
e8eecc2a | 88 | #ifdef SIND |
57391dda FR |
89 | /* Compute sind(x) = sin(x * pi / 180). */ |
90 | ||
91 | RETTYPE | |
92 | SIND (FTYPE x) | |
93 | { | |
e8eecc2a | 94 | #ifdef ENABLE_SIND |
57391dda FR |
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 | ||
e8eecc2a | 105 | #ifdef SIND_SMALL_LITERAL |
57391dda FR |
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); | |
e8eecc2a | 110 | if (mpfr_cmp_ld (ax, SIND_SMALL_LITERAL) < 0) |
57391dda FR |
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); | |
e8eecc2a | 122 | #endif /* SIND_SMALL_LITERAL */ |
57391dda FR |
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 (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 | 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); | |
e8eecc2a FR |
226 | |
227 | #else | |
228 | ERROR_RETURN(sind, KIND, x); | |
229 | #endif // ENABLE_SIND | |
57391dda | 230 | } |
e8eecc2a | 231 | #endif // SIND |
57391dda FR |
232 | |
233 | ||
e8eecc2a | 234 | #ifdef COSD |
57391dda FR |
235 | /* Compute cosd(x) = cos(x * pi / 180). */ |
236 | ||
237 | RETTYPE | |
238 | COSD (FTYPE x) | |
239 | { | |
e8eecc2a FR |
240 | #ifdef ENABLE_COSD |
241 | #if defined(TINY_LITERAL) && defined(COSD_SMALL_LITERAL) | |
242 | static const volatile FTYPE tiny = TINY_LITERAL; | |
57391dda FR |
243 | #endif |
244 | ||
245 | if (ISFINITE (x)) | |
246 | { | |
e8eecc2a | 247 | #ifdef COSD_SMALL_LITERAL |
57391dda FR |
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. */ | |
e8eecc2a | 253 | if (mpfr_cmp_ld (ax, COSD_SMALL_LITERAL) <= 0) |
57391dda FR |
254 | { |
255 | mpfr_set_ui (x, 1, GFC_RND_MODE); | |
e8eecc2a | 256 | #ifdef TINY_LITERAL |
57391dda FR |
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); | |
e8eecc2a | 270 | #endif /* COSD_SMALL_LITERAL */ |
57391dda FR |
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 (n, 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); | |
e8eecc2a FR |
374 | |
375 | #else | |
376 | ERROR_RETURN(cosd, KIND, x); | |
377 | #endif // ENABLE_COSD | |
57391dda | 378 | } |
e8eecc2a | 379 | #endif // COSD |
57391dda FR |
380 | |
381 | ||
e8eecc2a | 382 | #ifdef TAND |
57391dda FR |
383 | /* Compute tand(x) = tan(x * pi / 180). */ |
384 | ||
385 | RETTYPE | |
386 | TAND (FTYPE x) | |
387 | { | |
e8eecc2a | 388 | #ifdef ENABLE_TAND |
57391dda FR |
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 | ||
e8eecc2a | 399 | #ifdef SIND_SMALL_LITERAL |
57391dda FR |
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); | |
e8eecc2a | 404 | if (mpfr_cmp_ld (ax, SIND_SMALL_LITERAL) < 0) |
57391dda FR |
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); | |
e8eecc2a | 416 | #endif /* SIND_SMALL_LITERAL */ |
57391dda FR |
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 (n, 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); | |
e8eecc2a FR |
486 | |
487 | #else | |
488 | ERROR_RETURN(tand, KIND, x); | |
489 | #endif // ENABLE_TAND | |
57391dda | 490 | } |
e8eecc2a | 491 | #endif // TAND |
57391dda FR |
492 | |
493 | /* vim: set ft=c: */ |