|
152
|
1 /* Implementation of the degree trignometric functions COSD, SIND, TAND.
|
|
|
2 Copyright (C) 2020 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 SIND (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 (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);
|
|
|
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
|
|
|
237 RETTYPE
|
|
|
238 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 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 (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);
|
|
|
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
|
|
|
385 RETTYPE
|
|
|
386 TAND (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 (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);
|
|
|
486
|
|
|
487 #else
|
|
|
488 ERROR_RETURN(tand, KIND, x);
|
|
|
489 #endif // ENABLE_TAND
|
|
|
490 }
|
|
|
491 #endif // TAND
|
|
|
492
|
|
|
493 /* vim: set ft=c: */
|
