Mercurial > hg > CbC > CbC_gcc

comparison libgfortran/intrinsics/trigd.inc @ 152:2b5abeee2509

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update gcc11
author anatofuz
date Mon, 25 May 2020 07:50:57 +0900
parents
children
comparison
equal deleted inserted replaced
145:1830386684a0 152:2b5abeee2509
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: */