Vc 1.4.5
SIMD Vector Classes for C++
 
Loading...
Searching...
No Matches
math.h
1/* This file is part of the Vc library. {{{
2Copyright © 2013-2015 Matthias Kretz <kretz@kde.org>
3
4Redistribution and use in source and binary forms, with or without
5modification, are permitted provided that the following conditions are met:
6 * Redistributions of source code must retain the above copyright
7 notice, this list of conditions and the following disclaimer.
8 * Redistributions in binary form must reproduce the above copyright
9 notice, this list of conditions and the following disclaimer in the
10 documentation and/or other materials provided with the distribution.
11 * Neither the names of contributing organizations nor the
12 names of its contributors may be used to endorse or promote products
13 derived from this software without specific prior written permission.
14
15THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER BE LIABLE FOR ANY
19DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25
26}}}*/
27
28#ifndef VC_COMMON_MATH_H_
29#define VC_COMMON_MATH_H_
30
31#define Vc_COMMON_MATH_H_INTERNAL 1
32
33#include "trigonometric.h"
34
35#include "const.h"
36#include "macros.h"
37
38namespace Vc_VERSIONED_NAMESPACE
39{
40// TODO, not vectorized:
41template <class T, class Abi>
42SimdArray<int, Vector<T, Abi>::size()> fpclassify(const Vector<T, Abi> &x)
43{
44 return SimdArray<int, Vector<T, Abi>::size()>(
45 [&](std::size_t i) { return std::fpclassify(x[i]); });
46}
47template <class T, size_t N> SimdArray<int, N> fpclassify(const SimdArray<T, N> &x)
48{
49 return SimdArray<int, N>([&](std::size_t i) { return std::fpclassify(x[i]); });
50}
51
52#ifdef Vc_IMPL_SSE
53// for SSE, AVX, and AVX2
54#include "logarithm.h"
55#include "exponential.h"
56#ifdef Vc_IMPL_AVX
57inline AVX::double_v exp(AVX::double_v _x)
58{
59 AVX::Vector<double> x = _x;
60 typedef AVX::Vector<double> V;
61 typedef V::Mask M;
62 typedef AVX::Const<double> C;
63
64 const M overflow = x > Vc::Detail::doubleConstant< 1, 0x0006232bdd7abcd2ull, 9>(); // max log
65 const M underflow = x < Vc::Detail::doubleConstant<-1, 0x0006232bdd7abcd2ull, 9>(); // min log
66
67 V px = floor(C::log2_e() * x + 0.5);
68 __m128i tmp = _mm256_cvttpd_epi32(px.data());
69 const SimdArray<int, V::Size> n = SSE::int_v{tmp};
70 x -= px * C::ln2_large(); //Vc::Detail::doubleConstant<1, 0x00062e4000000000ull, -1>(); // ln2
71 x -= px * C::ln2_small(); //Vc::Detail::doubleConstant<1, 0x0007f7d1cf79abcaull, -20>(); // ln2
72
73 const double P[] = {
74 Vc::Detail::doubleConstant<1, 0x000089cdd5e44be8ull, -13>(),
75 Vc::Detail::doubleConstant<1, 0x000f06d10cca2c7eull, -6>(),
76 Vc::Detail::doubleConstant<1, 0x0000000000000000ull, 0>()
77 };
78 const double Q[] = {
79 Vc::Detail::doubleConstant<1, 0x00092eb6bc365fa0ull, -19>(),
80 Vc::Detail::doubleConstant<1, 0x0004ae39b508b6c0ull, -9>(),
81 Vc::Detail::doubleConstant<1, 0x000d17099887e074ull, -3>(),
82 Vc::Detail::doubleConstant<1, 0x0000000000000000ull, 1>()
83 };
84 const V x2 = x * x;
85 px = x * ((P[0] * x2 + P[1]) * x2 + P[2]);
86 x = px / ((((Q[0] * x2 + Q[1]) * x2 + Q[2]) * x2 + Q[3]) - px);
87 x = V::One() + 2.0 * x;
88
89 x = ldexp(x, n); // == x * 2ⁿ
90
91 x(overflow) = std::numeric_limits<double>::infinity();
92 x.setZero(underflow);
93
94 return x;
95 }
96#endif // Vc_IMPL_AVX
97
98inline SSE::double_v exp(SSE::double_v::AsArg _x) {
99 SSE::Vector<double> x = _x;
100 typedef SSE::Vector<double> V;
101 typedef V::Mask M;
102 typedef SSE::Const<double> C;
103
104 const M overflow = x > Vc::Detail::doubleConstant< 1, 0x0006232bdd7abcd2ull, 9>(); // max log
105 const M underflow = x < Vc::Detail::doubleConstant<-1, 0x0006232bdd7abcd2ull, 9>(); // min log
106
107 V px = floor(C::log2_e() * x + 0.5);
108 SimdArray<int, V::Size> n;
109 _mm_storel_epi64(reinterpret_cast<__m128i *>(&n), _mm_cvttpd_epi32(px.data()));
110 x -= px * C::ln2_large(); //Vc::Detail::doubleConstant<1, 0x00062e4000000000ull, -1>(); // ln2
111 x -= px * C::ln2_small(); //Vc::Detail::doubleConstant<1, 0x0007f7d1cf79abcaull, -20>(); // ln2
112
113 const double P[] = {
114 Vc::Detail::doubleConstant<1, 0x000089cdd5e44be8ull, -13>(),
115 Vc::Detail::doubleConstant<1, 0x000f06d10cca2c7eull, -6>(),
116 Vc::Detail::doubleConstant<1, 0x0000000000000000ull, 0>()
117 };
118 const double Q[] = {
119 Vc::Detail::doubleConstant<1, 0x00092eb6bc365fa0ull, -19>(),
120 Vc::Detail::doubleConstant<1, 0x0004ae39b508b6c0ull, -9>(),
121 Vc::Detail::doubleConstant<1, 0x000d17099887e074ull, -3>(),
122 Vc::Detail::doubleConstant<1, 0x0000000000000000ull, 1>()
123 };
124 const V x2 = x * x;
125 px = x * ((P[0] * x2 + P[1]) * x2 + P[2]);
126 x = px / ((((Q[0] * x2 + Q[1]) * x2 + Q[2]) * x2 + Q[3]) - px);
127 x = V::One() + 2.0 * x;
128
129 x = ldexp(x, n); // == x * 2ⁿ
130
131 x(overflow) = std::numeric_limits<double>::infinity();
132 x.setZero(underflow);
133
134 return x;
135 }
136
137#endif
138} // namespace Vc
139
140#undef Vc_COMMON_MATH_H_INTERNAL
141
142#endif // VC_COMMON_MATH_H_