s_cexpl.c revision 1.1 1 1.1 christos /*-
2 1.1 christos * SPDX-License-Identifier: BSD-2-Clause
3 1.1 christos *
4 1.1 christos * Copyright (c) 2011 David Schultz <das (at) FreeBSD.ORG>
5 1.1 christos * All rights reserved.
6 1.1 christos *
7 1.1 christos * Redistribution and use in source and binary forms, with or without
8 1.1 christos * modification, are permitted provided that the following conditions
9 1.1 christos * are met:
10 1.1 christos * 1. Redistributions of source code must retain the above copyright
11 1.1 christos * notice, this list of conditions and the following disclaimer.
12 1.1 christos * 2. Redistributions in binary form must reproduce the above copyright
13 1.1 christos * notice, this list of conditions and the following disclaimer in the
14 1.1 christos * documentation and/or other materials provided with the distribution.
15 1.1 christos *
16 1.1 christos * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
17 1.1 christos * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18 1.1 christos * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19 1.1 christos * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
20 1.1 christos * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
21 1.1 christos * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
22 1.1 christos * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
23 1.1 christos * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
24 1.1 christos * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
25 1.1 christos * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
26 1.1 christos * SUCH DAMAGE.
27 1.1 christos *
28 1.1 christos * src/s_cexp.c converted to long double complex by Steven G. Kargl
29 1.1 christos */
30 1.1 christos
31 1.1 christos #include <sys/cdefs.h>
32 1.1 christos #include <complex.h>
33 1.1 christos #include <float.h>
34 1.1 christos #ifdef __i386__
35 1.1 christos #include <ieeefp.h>
36 1.1 christos #endif
37 1.1 christos
38 1.1 christos #include "fpmath.h"
39 1.1 christos #include "math.h"
40 1.1 christos #include "math_private.h"
41 1.1 christos #include "k_expl.h"
42 1.1 christos
43 1.1 christos long double complex
44 1.1 christos cexpl (long double complex z)
45 1.1 christos {
46 1.1 christos long double c, exp_x, s, x, y;
47 1.1 christos uint64_t lx, ly;
48 1.1 christos uint16_t hx, hy;
49 1.1 christos
50 1.1 christos ENTERI();
51 1.1 christos
52 1.1 christos x = creall(z);
53 1.1 christos y = cimagl(z);
54 1.1 christos
55 1.1 christos EXTRACT_LDBL80_WORDS(hy, ly, y);
56 1.1 christos hy &= 0x7fff;
57 1.1 christos
58 1.1 christos /* cexp(x + I 0) = exp(x) + I 0 */
59 1.1 christos if ((hy | ly) == 0)
60 1.1 christos RETURNI(CMPLXL(expl(x), y));
61 1.1 christos EXTRACT_LDBL80_WORDS(hx, lx, x);
62 1.1 christos /* cexp(0 + I y) = cos(y) + I sin(y) */
63 1.1 christos if (((hx & 0x7fff) | lx) == 0) {
64 1.1 christos sincosl(y, &s, &c);
65 1.1 christos RETURNI(CMPLXL(c, s));
66 1.1 christos }
67 1.1 christos
68 1.1 christos if (hy >= 0x7fff) {
69 1.1 christos if ((hx & 0x7fff) < 0x7fff || ((hx & 0x7fff) == 0x7fff &&
70 1.1 christos (lx & 0x7fffffffffffffffULL) != 0)) {
71 1.1 christos /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
72 1.1 christos RETURNI(CMPLXL(y - y, y - y));
73 1.1 christos } else if (hx & 0x8000) {
74 1.1 christos /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
75 1.1 christos RETURNI(CMPLXL(0.0, 0.0));
76 1.1 christos } else {
77 1.1 christos /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
78 1.1 christos RETURNI(CMPLXL(x, y - y));
79 1.1 christos }
80 1.1 christos }
81 1.1 christos
82 1.1 christos /*
83 1.1 christos * exp_ovfl = 11356.5234062941439497
84 1.1 christos * cexp_ovfl = 22755.3287906024445633
85 1.1 christos */
86 1.1 christos if ((hx == 0x400c && lx > 0xb17217f7d1cf79acULL) ||
87 1.1 christos (hx == 0x400d && lx < 0xb1c6a8573de9768cULL)) {
88 1.1 christos /*
89 1.1 christos * x is between exp_ovfl and cexp_ovfl, so we must scale to
90 1.1 christos * avoid overflow in exp(x).
91 1.1 christos */
92 1.1 christos RETURNI(__ldexp_cexpl(z, 0));
93 1.1 christos } else {
94 1.1 christos /*
95 1.1 christos * Cases covered here:
96 1.1 christos * - x < exp_ovfl and exp(x) won't overflow (common case)
97 1.1 christos * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0
98 1.1 christos * - x = +-Inf (generated by exp())
99 1.1 christos * - x = NaN (spurious inexact exception from y)
100 1.1 christos */
101 1.1 christos exp_x = expl(x);
102 1.1 christos sincosl(y, &s, &c);
103 1.1 christos RETURNI(CMPLXL(exp_x * c, exp_x * s));
104 1.1 christos }
105 1.1 christos }
106