1 /*
2 * Copyright 2008-2013 NVIDIA Corporation
3 * Copyright 2013 Filipe RNC Maia
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*-
19 * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
20 * All rights reserved.
21 *
22 * Redistribution and use in source and binary forms, with or without
23 * modification, are permitted provided that the following conditions
24 * are met:
25 * 1. Redistributions of source code must retain the above copyright
26 * notice, this list of conditions and the following disclaimer.
27 * 2. Redistributions in binary form must reproduce the above copyright
28 * notice, this list of conditions and the following disclaimer in the
29 * documentation and/or other materials provided with the distribution.
30 *
31 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
32 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
33 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
34 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
35 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
36 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
37 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
38 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
39 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
40 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
41 * SUCH DAMAGE.
42 */
43
44 /* adapted from FreeBSD:
45 * lib/msun/src/s_cexp.c
46 * lib/msun/src/k_exp.c
47 *
48 */
49
50 #pragma once
51
52 #include <thrust/complex.h>
53 #include <thrust/detail/complex/math_private.h>
54
55 namespace thrust{
56 namespace detail{
57 namespace complex{
58 /*
59 * Compute exp(x), scaled to avoid spurious overflow. An exponent is
60 * returned separately in 'expt'.
61 *
62 * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
63 * Output: 2**1023 <= y < 2**1024
64 */
65 __host__ __device__ inline
frexp_exp(double x,int * expt)66 double frexp_exp(double x, int *expt){
67 const uint32_t k = 1799; /* constant for reduction */
68 const double kln2 = 1246.97177782734161156; /* k * ln2 */
69
70 double exp_x;
71 uint32_t hx;
72
73 /*
74 * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
75 * minimize |exp(kln2) - 2**k|. We also scale the exponent of
76 * exp_x to MAX_EXP so that the result can be multiplied by
77 * a tiny number without losing accuracy due to denormalization.
78 */
79 exp_x = exp(x - kln2);
80 get_high_word(hx, exp_x);
81 *expt = (hx >> 20) - (0x3ff + 1023) + k;
82 set_high_word(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20));
83 return (exp_x);
84 }
85
86
87 __host__ __device__ inline
ldexp_cexp(complex<double> z,int expt)88 complex<double> ldexp_cexp(complex<double> z, int expt){
89 double x, y, exp_x, scale1, scale2;
90 int ex_expt, half_expt;
91
92 x = z.real();
93 y = z.imag();
94 exp_x = frexp_exp(x, &ex_expt);
95 expt += ex_expt;
96
97 /*
98 * Arrange so that scale1 * scale2 == 2**expt. We use this to
99 * compensate for scalbn being horrendously slow.
100 */
101 half_expt = expt / 2;
102 insert_words(scale1, (0x3ff + half_expt) << 20, 0);
103 half_expt = expt - half_expt;
104 insert_words(scale2, (0x3ff + half_expt) << 20, 0);
105
106 return (complex<double>(cos(y) * exp_x * scale1 * scale2,
107 sin(y) * exp_x * scale1 * scale2));
108 }
109
110
111 __host__ __device__ inline
cexp(const complex<double> & z)112 complex<double> cexp(const complex<double>& z){
113 double x, y, exp_x;
114 uint32_t hx, hy, lx, ly;
115
116 const uint32_t
117 exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */
118 cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
119
120
121 x = z.real();
122 y = z.imag();
123
124 extract_words(hy, ly, y);
125 hy &= 0x7fffffff;
126
127 /* cexp(x + I 0) = exp(x) + I 0 */
128 if ((hy | ly) == 0)
129 return (complex<double>(exp(x), y));
130 extract_words(hx, lx, x);
131 /* cexp(0 + I y) = cos(y) + I sin(y) */
132 if (((hx & 0x7fffffff) | lx) == 0)
133 return (complex<double>(cos(y), sin(y)));
134
135 if (hy >= 0x7ff00000) {
136 if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
137 /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
138 return (complex<double>(y - y, y - y));
139 } else if (hx & 0x80000000) {
140 /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
141 return (complex<double>(0.0, 0.0));
142 } else {
143 /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
144 return (complex<double>(x, y - y));
145 }
146 }
147
148 if (hx >= exp_ovfl && hx <= cexp_ovfl) {
149 /*
150 * x is between 709.7 and 1454.3, so we must scale to avoid
151 * overflow in exp(x).
152 */
153 return (ldexp_cexp(z, 0));
154 } else {
155 /*
156 * Cases covered here:
157 * - x < exp_ovfl and exp(x) won't overflow (common case)
158 * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0
159 * - x = +-Inf (generated by exp())
160 * - x = NaN (spurious inexact exception from y)
161 */
162 exp_x = std::exp(x);
163 return (complex<double>(exp_x * cos(y), exp_x * sin(y)));
164 }
165 }
166
167 } // namespace complex
168
169 } // namespace detail
170
171 template <typename ValueType>
172 __host__ __device__
exp(const complex<ValueType> & z)173 inline complex<ValueType> exp(const complex<ValueType>& z){
174 return polar(std::exp(z.real()),z.imag());
175 }
176
177 template <>
178 __host__ __device__
exp(const complex<double> & z)179 inline complex<double> exp(const complex<double>& z){
180 return detail::complex::cexp(z);
181 }
182
183 } // namespace thrust
184