ROOT  6.06/09
Reference Guide
math.cpp
Go to the documentation of this file.
1 /* This file is part of the Vc library. {{{
2 
3  Copyright (C) 2009-2012 Matthias Kretz <kretz@kde.org>
4 
5  Vc is free software: you can redistribute it and/or modify
6  it under the terms of the GNU Lesser General Public License as
7  published by the Free Software Foundation, either version 3 of
8  the License, or (at your option) any later version.
9 
10  Vc is distributed in the hope that it will be useful, but
11  WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU Lesser General Public License for more details.
14 
15  You should have received a copy of the GNU Lesser General Public
16  License along with Vc. If not, see <http://www.gnu.org/licenses/>.
17 
18 }}}*/
19 /*includes {{{*/
20 #include "unittest.h"
21 #include <iostream>
22 #include "vectormemoryhelper.h"
23 #include "const.h"
24 #include <cmath>
25 #include <algorithm>
26 #include <Vc/common/macros.h>
27 /*}}}*/
28 using namespace Vc;
29 /*fix isfinite and isnan{{{*/
30 #ifdef isfinite
31 #undef isfinite
32 #endif
33 #ifdef isnan
34 #undef isnan
35 #endif
36 /*}}}*/
37 template<typename T> struct SincosReference/*{{{*/
38 {
39  T x, s, c;
40 };
41 template<typename T> struct Reference
42 {
43  T x, ref;
44 };
45 
46 template<typename T> struct Array
47 {
48  size_t size;
49  T *data;
50  Array() : size(0), data(0) {}
51 };
52 template<typename T> struct StaticDeleter
53 {
54  T *ptr;
55  StaticDeleter(T *p) : ptr(p) {}
56  ~StaticDeleter() { delete[] ptr; }
57 };
58 
59 enum Function {
61 };
62 template<typename T, Function F> static inline const char *filename();
63 template<> inline const char *filename<float , Sincos>() { return "reference-sincos-sp.dat"; }
64 template<> inline const char *filename<double, Sincos>() { return "reference-sincos-dp.dat"; }
65 template<> inline const char *filename<float , Atan >() { return "reference-atan-sp.dat"; }
66 template<> inline const char *filename<double, Atan >() { return "reference-atan-dp.dat"; }
67 template<> inline const char *filename<float , Asin >() { return "reference-asin-sp.dat"; }
68 template<> inline const char *filename<double, Asin >() { return "reference-asin-dp.dat"; }
69 // template<> inline const char *filename<float , Acos >() { return "reference-acos-sp.dat"; }
70 // template<> inline const char *filename<double, Acos >() { return "reference-acos-dp.dat"; }
71 template<> inline const char *filename<float , Log >() { return "reference-ln-sp.dat"; }
72 template<> inline const char *filename<double, Log >() { return "reference-ln-dp.dat"; }
73 template<> inline const char *filename<float , Log2 >() { return "reference-log2-sp.dat"; }
74 template<> inline const char *filename<double, Log2 >() { return "reference-log2-dp.dat"; }
75 template<> inline const char *filename<float , Log10 >() { return "reference-log10-sp.dat"; }
76 template<> inline const char *filename<double, Log10 >() { return "reference-log10-dp.dat"; }
77 
78 template<typename T>
79 static Array<SincosReference<T> > sincosReference()
80 {
81  static Array<SincosReference<T> > data;
82  if (data.data == 0) {
83  FILE *file = fopen(filename<T, Sincos>(), "rb");
84  if (file) {
85  fseek(file, 0, SEEK_END);
86  const size_t size = ftell(file) / sizeof(SincosReference<T>);
87  rewind(file);
88  data.data = new SincosReference<T>[size];
89  static StaticDeleter<SincosReference<T> > _cleanup(data.data);
90  data.size = fread(data.data, sizeof(SincosReference<T>), size, file);
91  fclose(file);
92  } else {
93  FAIL() << "the reference data " << filename<T, Sincos>() << " does not exist in the current working directory.";
94  }
95  }
96  return data;
97 }
98 
99 template<typename T, Function Fun>
100 static Array<Reference<T> > referenceData()
101 {
102  static Array<Reference<T> > data;
103  if (data.data == 0) {
104  FILE *file = fopen(filename<T, Fun>(), "rb");
105  if (file) {
106  fseek(file, 0, SEEK_END);
107  const size_t size = ftell(file) / sizeof(Reference<T>);
108  rewind(file);
109  data.data = new Reference<T>[size];
110  static StaticDeleter<Reference<T> > _cleanup(data.data);
111  data.size = fread(data.data, sizeof(Reference<T>), size, file);
112  fclose(file);
113  } else {
114  FAIL() << "the reference data " << filename<T, Fun>() << " does not exist in the current working directory.";
115  }
116  }
117  return data;
118 }/*}}}*/
119 
120 template<typename T> struct Denormals { static T *data; };/*{{{*/
121 template<> float *Denormals<float >::data = 0;
122 template<> double *Denormals<double>::data = 0;
123 enum {
125 };
126 /*}}}*/
127 template<typename V> V apply_v(VC_ALIGNED_PARAMETER(V) x, typename V::EntryType (func)(typename V::EntryType))/*{{{*/
128 {
129  V r;
130  for (size_t i = 0; i < V::Size; ++i) {
131  r[i] = func(x[i]);
132  }
133  return r;
134 }
135 /*}}}*/
136 template<typename Vec> void testAbs()/*{{{*/
137 {
138  for (int i = 0; i < 0x7fff; ++i) {
139  Vec a(i);
140  Vec b(-i);
141  COMPARE(a, Vc::abs(a));
142  COMPARE(a, Vc::abs(b));
143  }
144 }
145 /*}}}*/
146 static inline float my_trunc(float x)/*{{{*/
147 {
148 #if __cplusplus >= 201103 /*C++11*/
149  return std::trunc(x);
150 #elif defined(_ISOC99_SOURCE)
151  return truncf(x);
152 #else
153  return x > 0 ? std::floor(x) : std::ceil(x);
154 #endif
155 }
156 
157 static inline double my_trunc(double x)
158 {
159 #if __cplusplus >= 201103 /*C++11*/
160  return std::trunc(x);
161 #elif defined(_ISOC99_SOURCE)
162  return trunc(x);
163 #else
164  return x > 0 ? std::floor(x) : std::ceil(x);
165 #endif
166 }
167 /*}}}*/
168 template<typename V> void testTrunc()/*{{{*/
169 {
170  typedef typename V::EntryType T;
171  typedef typename V::IndexType I;
172  for (size_t i = 0; i < 100000 / V::Size; ++i) {
173  V x = (V::Random() - T(0.5)) * T(100);
174  V reference = apply_v(x, my_trunc);
175  COMPARE(Vc::trunc(x), reference) << ", x = " << x << ", i = " << i;
176  }
177  V x = static_cast<V>(I::IndexesFromZero());
178  V reference = apply_v(x, my_trunc);
179  COMPARE(Vc::trunc(x), reference) << ", x = " << x;
180 }
181 /*}}}*/
182 template<typename V> void testFloor()/*{{{*/
183 {
184  typedef typename V::EntryType T;
185  typedef typename V::IndexType I;
186  for (size_t i = 0; i < 100000 / V::Size; ++i) {
187  V x = (V::Random() - T(0.5)) * T(100);
188  V reference = apply_v(x, std::floor);
189  COMPARE(Vc::floor(x), reference) << ", x = " << x << ", i = " << i;
190  }
191  V x = static_cast<V>(I::IndexesFromZero());
192  V reference = apply_v(x, std::floor);
193  COMPARE(Vc::floor(x), reference) << ", x = " << x;
194 }
195 /*}}}*/
196 template<typename V> void testCeil()/*{{{*/
197 {
198  typedef typename V::EntryType T;
199  typedef typename V::IndexType I;
200  for (size_t i = 0; i < 100000 / V::Size; ++i) {
201  V x = (V::Random() - T(0.5)) * T(100);
202  V reference = apply_v(x, std::ceil);
203  COMPARE(Vc::ceil(x), reference) << ", x = " << x << ", i = " << i;
204  }
205  V x = static_cast<V>(I::IndexesFromZero());
206  V reference = apply_v(x, std::ceil);
207  COMPARE(Vc::ceil(x), reference) << ", x = " << x;
208 }
209 /*}}}*/
210 template<typename V> void testExp()/*{{{*/
211 {
214  typedef typename V::EntryType T;
215  for (size_t i = 0; i < 100000 / V::Size; ++i) {
216  V x = (V::Random() - T(0.5)) * T(20);
217  V reference = apply_v(x, std::exp);
218  FUZZY_COMPARE(Vc::exp(x), reference) << ", x = " << x << ", i = " << i;
219  }
220  COMPARE(Vc::exp(V::Zero()), V::One());
221 }
222 /*}}}*/
223 template<typename V> void testLog()/*{{{*/
224 {
226  typedef typename V::EntryType T;
227  Array<Reference<T> > reference = referenceData<T, Log>();
228  for (size_t i = 0; i + V::Size - 1 < reference.size; i += V::Size) {
229  V x, ref;
230  for (int j = 0; j < V::Size; ++j) {
231  x[j] = reference.data[i + j].x;
232  ref[j] = reference.data[i + j].ref;
233  }
234  FUZZY_COMPARE(Vc::log(x), ref) << " x = " << x << ", i = " << i;
235  }
236 
237  COMPARE(Vc::log(V::Zero()), V(std::log(T(0))));
238  for (int i = 0; i < NDenormals; i += V::Size) {
239  V x(&Denormals<T>::data[i]);
240  V ref = apply_v(x, std::log);
241  FUZZY_COMPARE(Vc::log(x), ref) << ", x = " << x << ", i = " << i;
242  }
243 }
244 /*}}}*/
245 #if (defined(_XOPEN_SOURCE) && _XOPEN_SOURCE >= 600) || defined(_ISOC99_SOURCE) || (defined(_POSIX_C_SOURCE) && _POSIX_C_SOURCE >= 200112L)
246 static inline float my_log2(float x) { return ::log2f(x); }
247 /* I need to make sure whether the log2 that I compare against is really precise to <0.5ulp. At
248  * least I get different results when I use "double log2(double)", which is somewhat unexpected.
249  * Well, conversion from double to float goes via truncation, so if the most significant truncated
250  * mantissa bit is set the resulting float is incorrect by 1 ulp
251 
252 static inline float my_log2(float x) { return ::log2(static_cast<double>(x)); }
253 static inline float my_log2(float x) {
254  double tmp = ::log2(static_cast<double>(x));
255  int e;
256  frexp(tmp, &e); // frexp(0.5) -> e = 0
257  return tmp + ldexp(tmp < 0 ? -0.5 : 0.5, e - 24);
258 }
259  */
260 static inline double my_log2(double x) { return ::log2(x); }
261 #else
262 static inline float my_log2(float x) { return ::logf(x) / Vc::Math<float>::ln2(); }
263 static inline double my_log2(double x) { return ::log(x) / Vc::Math<double>::ln2(); }
264 #endif
265 /*}}}*/
266 template<typename V> void testLog2()/*{{{*/
267 {
268 #if defined(VC_LOG_ILP) || defined(VC_LOG_ILP2)
270 #else
272 #endif
273 #if (defined(VC_MSVC) || defined(__APPLE__)) && defined(VC_IMPL_Scalar)
275 #else
277 #endif
278  typedef typename V::EntryType T;
279  Array<Reference<T> > reference = referenceData<T, Log2>();
280  for (size_t i = 0; i + V::Size - 1 < reference.size; i += V::Size) {
281  V x, ref;
282  for (int j = 0; j < V::Size; ++j) {
283  x[j] = reference.data[i + j].x;
284  ref[j] = reference.data[i + j].ref;
285  }
286  FUZZY_COMPARE(Vc::log2(x), ref) << " x = " << x << ", i = " << i;
287  }
288 
289  COMPARE(Vc::log2(V::Zero()), V(my_log2(T(0))));
290  for (int i = 0; i < NDenormals; i += V::Size) {
291  V x(&Denormals<T>::data[i]);
292  V ref = apply_v(x, my_log2);
293  FUZZY_COMPARE(Vc::log2(x), ref) << ", x = " << x << ", i = " << i;
294  }
295 }
296 /*}}}*/
297 template<typename V> void testLog10()/*{{{*/
298 {
301  typedef typename V::EntryType T;
302  Array<Reference<T> > reference = referenceData<T, Log10>();
303  for (size_t i = 0; i + V::Size - 1 < reference.size; i += V::Size) {
304  V x, ref;
305  for (int j = 0; j < V::Size; ++j) {
306  x[j] = reference.data[i + j].x;
307  ref[j] = reference.data[i + j].ref;
308  }
309  FUZZY_COMPARE(Vc::log10(x), ref) << " x = " << x << ", i = " << i;
310  }
311 
312  COMPARE(Vc::log10(V::Zero()), V(std::log10(T(0))));
313  for (int i = 0; i < NDenormals; i += V::Size) {
314  V x(&Denormals<T>::data[i]);
315  V ref = apply_v(x, std::log10);
316  FUZZY_COMPARE(Vc::log10(x), ref) << ", x = " << x << ", i = " << i;
317  }
318 }
319 /*}}}*/
320 template<typename Vec> void testMax()/*{{{*/
321 {
322  typedef typename Vec::EntryType T;
324  T *data = mem;
325  for (int i = 0; i < Vec::Size; ++i) {
326  data[i] = i;
327  data[i + Vec::Size] = Vec::Size + 1 - i;
328  data[i + 2 * Vec::Size] = std::max(data[i], data[i + Vec::Size]);
329  }
330  Vec a(&data[0]);
331  Vec b(&data[Vec::Size]);
332  Vec c(&data[2 * Vec::Size]);
333 
334  COMPARE(Vc::max(a, b), c);
335 }
336 /*}}}*/
337  /*{{{*/
338 #define FillHelperMemory(code) \
339  typename V::Memory data; \
340  typename V::Memory reference; \
341  for (int ii = 0; ii < V::Size; ++ii) { \
342  const T i = static_cast<T>(ii); \
343  data[ii] = i; \
344  reference[ii] = code; \
345  } do {} while (false)
346 /*}}}*/
347 template<typename V> void testSqrt()/*{{{*/
348 {
349  typedef typename V::EntryType T;
351  V a(data);
352  V b(reference);
353 
354  FUZZY_COMPARE(Vc::sqrt(a), b);
355 }
356 /*}}}*/
357 template<typename V> void testRSqrt()/*{{{*/
358 {
359  typedef typename V::EntryType T;
360  for (size_t i = 0; i < 1024 / V::Size; ++i) {
361  const V x = V::Random() * T(1000);
362  // RSQRTPS is documented as having a relative error <= 1.5 * 2^-12
363  VERIFY(Vc::abs(Vc::rsqrt(x) * Vc::sqrt(x) - V::One()) < static_cast<T>(std::ldexp(1.5, -12)));
364  }
365 }
366 /*}}}*/
367 template<typename V> void testSincos()/*{{{*/
368 {
369  typedef typename V::EntryType T;
372  Array<SincosReference<T> > reference = sincosReference<T>();
373  for (size_t i = 0; i + V::Size - 1 < reference.size; i += V::Size) {
374  V x, sref, cref;
375  for (int j = 0; j < V::Size; ++j) {
376  x[j] = reference.data[i + j].x;
377  sref[j] = reference.data[i + j].s;
378  cref[j] = reference.data[i + j].c;
379  }
380  V sin, cos;
381  Vc::sincos(x, &sin, &cos);
382  FUZZY_COMPARE(sin, sref) << " x = " << x << ", i = " << i;
383  FUZZY_COMPARE(cos, cref) << " x = " << x << ", i = " << i;
384  Vc::sincos(-x, &sin, &cos);
385  FUZZY_COMPARE(sin, -sref) << " x = " << -x << ", i = " << i;
386  FUZZY_COMPARE(cos, cref) << " x = " << -x << ", i = " << i;
387  }
388 }
389 /*}}}*/
390 template<typename V> void testSin()/*{{{*/
391 {
392  typedef typename V::EntryType T;
395  Array<SincosReference<T> > reference = sincosReference<T>();
396  for (size_t i = 0; i + V::Size - 1 < reference.size; i += V::Size) {
397  V x, sref;
398  for (int j = 0; j < V::Size; ++j) {
399  x[j] = reference.data[i + j].x;
400  sref[j] = reference.data[i + j].s;
401  }
402  FUZZY_COMPARE(Vc::sin(x), sref) << " x = " << x << ", i = " << i;
403  FUZZY_COMPARE(Vc::sin(-x), -sref) << " x = " << x << ", i = " << i;
404  }
405 }
406 /*}}}*/
407 template<typename V> void testCos()/*{{{*/
408 {
409  typedef typename V::EntryType T;
412  Array<SincosReference<T> > reference = sincosReference<T>();
413  for (size_t i = 0; i + V::Size - 1 < reference.size; i += V::Size) {
414  V x, cref;
415  for (int j = 0; j < V::Size; ++j) {
416  x[j] = reference.data[i + j].x;
417  cref[j] = reference.data[i + j].c;
418  }
419  FUZZY_COMPARE(Vc::cos(x), cref) << " x = " << x << ", i = " << i;
420  FUZZY_COMPARE(Vc::cos(-x), cref) << " x = " << x << ", i = " << i;
421  }
422 }
423 /*}}}*/
424 template<typename V> void testAsin()/*{{{*/
425 {
426  typedef typename V::EntryType T;
429  Array<Reference<T> > reference = referenceData<T, Asin>();
430  for (size_t i = 0; i + V::Size - 1 < reference.size; i += V::Size) {
431  V x, ref;
432  for (int j = 0; j < V::Size; ++j) {
433  x[j] = reference.data[i + j].x;
434  ref[j] = reference.data[i + j].ref;
435  }
436  FUZZY_COMPARE(Vc::asin(x), ref) << " x = " << x << ", i = " << i;
437  FUZZY_COMPARE(Vc::asin(-x), -ref) << " -x = " << -x << ", i = " << i;
438  }
439 }
440 /*}}}*/
441 const union {
442  unsigned int hex;
443  float value;
444 } INF = { 0x7f800000 };
445 
446 #if defined(__APPLE__) && defined(VC_IMPL_Scalar)
447 #define ATAN_COMPARE FUZZY_COMPARE
448 #else
449 #define ATAN_COMPARE COMPARE
450 #endif
451 
452 template<typename V> void testAtan()/*{{{*/
453 {
454  typedef typename V::EntryType T;
457 
458  {
459  const V Pi_2 = T(Vc_buildDouble(1, 0x921fb54442d18ull, 0));
460  V nan; nan.setQnan();
461  const V inf = T(INF.value);
462 
463  VERIFY(Vc::isnan(Vc::atan(nan)));
464  ATAN_COMPARE(Vc::atan(+inf), +Pi_2);
465 #ifdef VC_MSVC
466 #pragma warning(suppress: 4756) // overflow in constant arithmetic
467 #endif
468  ATAN_COMPARE(Vc::atan(-inf), -Pi_2);
469  }
470 
471  Array<Reference<T> > reference = referenceData<T, Atan>();
472  for (size_t i = 0; i + V::Size - 1 < reference.size; i += V::Size) {
473  V x, ref;
474  for (int j = 0; j < V::Size; ++j) {
475  x[j] = reference.data[i + j].x;
476  ref[j] = reference.data[i + j].ref;
477  }
478  FUZZY_COMPARE(Vc::atan(x), ref) << " x = " << x << ", i = " << i;
479  FUZZY_COMPARE(Vc::atan(-x), -ref) << " -x = " << -x << ", i = " << i;
480  }
481 }
482 /*}}}*/
483 template<typename V> void testAtan2()/*{{{*/
484 {
485  typedef typename V::EntryType T;
488 
489  {
490  const V Pi = T(Vc_buildDouble(1, 0x921fb54442d18ull, 1));
491  const V Pi_2 = T(Vc_buildDouble(1, 0x921fb54442d18ull, 0));
492  V nan; nan.setQnan();
493  const V inf = T(INF.value);
494 
495  // If y is +0 (-0) and x is less than 0, +pi (-pi) is returned.
496  ATAN_COMPARE(Vc::atan2(V(T(+0.)), V(T(-3.))), +Pi);
497  ATAN_COMPARE(Vc::atan2(V(T(-0.)), V(T(-3.))), -Pi);
498  // If y is +0 (-0) and x is greater than 0, +0 (-0) is returned.
499  COMPARE(Vc::atan2(V(T(+0.)), V(T(+3.))), V(T(+0.)));
500  VERIFY(!Vc::atan2(V(T(+0.)), V(T(+3.))).isNegative());
501  COMPARE(Vc::atan2(V(T(-0.)), V(T(+3.))), V(T(-0.)));
502  VERIFY (Vc::atan2(V(T(-0.)), V(T(+3.))).isNegative());
503  // If y is less than 0 and x is +0 or -0, -pi/2 is returned.
504  COMPARE(Vc::atan2(V(T(-3.)), V(T(+0.))), -Pi_2);
505  COMPARE(Vc::atan2(V(T(-3.)), V(T(-0.))), -Pi_2);
506  // If y is greater than 0 and x is +0 or -0, pi/2 is returned.
507  COMPARE(Vc::atan2(V(T(+3.)), V(T(+0.))), +Pi_2);
508  COMPARE(Vc::atan2(V(T(+3.)), V(T(-0.))), +Pi_2);
509  // If either x or y is NaN, a NaN is returned.
510  VERIFY(Vc::isnan(Vc::atan2(nan, V(T(3.)))));
511  VERIFY(Vc::isnan(Vc::atan2(V(T(3.)), nan)));
512  VERIFY(Vc::isnan(Vc::atan2(nan, nan)));
513  // If y is +0 (-0) and x is -0, +pi (-pi) is returned.
514  ATAN_COMPARE(Vc::atan2(V(T(+0.)), V(T(-0.))), +Pi);
515  ATAN_COMPARE(Vc::atan2(V(T(-0.)), V(T(-0.))), -Pi);
516  // If y is +0 (-0) and x is +0, +0 (-0) is returned.
517  COMPARE(Vc::atan2(V(T(+0.)), V(T(+0.))), V(T(+0.)));
518  COMPARE(Vc::atan2(V(T(-0.)), V(T(+0.))), V(T(-0.)));
519  VERIFY(!Vc::atan2(V(T(+0.)), V(T(+0.))).isNegative());
520  VERIFY( Vc::atan2(V(T(-0.)), V(T(+0.))).isNegative());
521  // If y is a finite value greater (less) than 0, and x is negative infinity, +pi (-pi) is returned.
522  ATAN_COMPARE(Vc::atan2(V(T(+1.)), -inf), +Pi);
523  ATAN_COMPARE(Vc::atan2(V(T(-1.)), -inf), -Pi);
524  // If y is a finite value greater (less) than 0, and x is positive infinity, +0 (-0) is returned.
525  COMPARE(Vc::atan2(V(T(+3.)), +inf), V(T(+0.)));
526  VERIFY(!Vc::atan2(V(T(+3.)), +inf).isNegative());
527  COMPARE(Vc::atan2(V(T(-3.)), +inf), V(T(-0.)));
528  VERIFY (Vc::atan2(V(T(-3.)), +inf).isNegative());
529  // If y is positive infinity (negative infinity), and x is finite, pi/2 (-pi/2) is returned.
530  COMPARE(Vc::atan2(+inf, V(T(+3.))), +Pi_2);
531  COMPARE(Vc::atan2(-inf, V(T(+3.))), -Pi_2);
532  COMPARE(Vc::atan2(+inf, V(T(-3.))), +Pi_2);
533  COMPARE(Vc::atan2(-inf, V(T(-3.))), -Pi_2);
534 #ifndef _WIN32 // the Microsoft implementation of atan2 fails this test
535  const V Pi_4 = T(Vc_buildDouble(1, 0x921fb54442d18ull, -1));
536  // If y is positive infinity (negative infinity) and x is negative infinity, +3*pi/4 (-3*pi/4) is returned.
537  COMPARE(Vc::atan2(+inf, -inf), T(+3.) * Pi_4);
538  COMPARE(Vc::atan2(-inf, -inf), T(-3.) * Pi_4);
539  // If y is positive infinity (negative infinity) and x is positive infinity, +pi/4 (-pi/4) is returned.
540  COMPARE(Vc::atan2(+inf, +inf), +Pi_4);
541  COMPARE(Vc::atan2(-inf, +inf), -Pi_4);
542 #endif
543  }
544 
545  for (int xoffset = -100; xoffset < 54613; xoffset += 47 * V::Size) {
546  for (int yoffset = -100; yoffset < 54613; yoffset += 47 * V::Size) {
547  FillHelperMemory(std::atan2((i + xoffset) * T(0.15), (i + yoffset) * T(0.15)));
548  const V a(data);
549  const V b(reference);
550 
551  const V x = (a + xoffset) * T(0.15);
552  const V y = (a + yoffset) * T(0.15);
553  FUZZY_COMPARE(Vc::atan2(x, y), b) << ", x = " << x << ", y = " << y;
554  }
555  }
556 }
557 /*}}}*/
558 template<typename Vec> void testReciprocal()/*{{{*/
559 {
560  typedef typename Vec::EntryType T;
561  setFuzzyness<float>(1.258295e+07);
563  const T one = 1;
564  for (int offset = -1000; offset < 1000; offset += 10) {
565  const T scale = T(0.1);
566  typename Vec::Memory data;
567  typename Vec::Memory reference;
568  for (int ii = 0; ii < Vec::Size; ++ii) {
569  const T i = static_cast<T>(ii);
570  data[ii] = i;
571  T tmp = (i + offset) * scale;
572  reference[ii] = one / tmp;
573  }
574  Vec a(data);
575  Vec b(reference);
576 
577  FUZZY_COMPARE(Vc::reciprocal((a + offset) * scale), b);
578  }
579 }
580 /*}}}*/
581 template<typename V> void isNegative()/*{{{*/
582 {
583  typedef typename V::EntryType T;
584  VERIFY(V::One().isNegative().isEmpty());
585  VERIFY(V::Zero().isNegative().isEmpty());
586  VERIFY((-V::One()).isNegative().isFull());
587  VERIFY(V(T(-0.)).isNegative().isFull());
588 }
589 /*}}}*/
590 template<typename Vec> void testInf()/*{{{*/
591 {
592  typedef typename Vec::EntryType T;
593  const T one = 1;
594  const Vec zero(Zero);
595  VERIFY(Vc::isfinite(zero));
596  VERIFY(Vc::isfinite(Vec(one)));
597  VERIFY(!Vc::isfinite(one / zero));
598 }
599 /*}}}*/
600 template<typename Vec> void testNaN()/*{{{*/
601 {
602  typedef typename Vec::EntryType T;
603  typedef typename Vec::IndexType I;
604  typedef typename Vec::Mask M;
605  const T one = 1;
606  const Vec zero(Zero);
607  VERIFY(!Vc::isnan(zero));
608  VERIFY(!Vc::isnan(Vec(one)));
609  const Vec inf = one / zero;
610  VERIFY(Vc::isnan(Vec(inf * zero)));
611  Vec nan = Vec::Zero();
612  const M mask(I::IndexesFromZero() == I::Zero());
613  nan.setQnan(mask);
614  COMPARE(Vc::isnan(nan), mask);
615  nan.setQnan();
616  VERIFY(Vc::isnan(nan));
617 }
618 /*}}}*/
619 template<typename Vec> void testRound()/*{{{*/
620 {
621  typedef typename Vec::EntryType T;
622  enum {
623  Count = (16 + Vec::Size) / Vec::Size
624  };
625  VectorMemoryHelper<Vec> mem1(Count);
626  VectorMemoryHelper<Vec> mem2(Count);
627  T *data = mem1;
628  T *reference = mem2;
629  for (int i = 0; i < Count * Vec::Size; ++i) {
630  data[i] = i * 0.25 - 2.0;
631  reference[i] = std::floor(i * 0.25 - 2.0 + 0.5);
632  if (i % 8 == 2) {
633  reference[i] -= 1.;
634  }
635  //std::cout << reference[i] << " ";
636  }
637  //std::cout << std::endl;
638  for (int i = 0; i < Count; ++i) {
639  const Vec a(&data[i * Vec::Size]);
640  const Vec ref(&reference[i * Vec::Size]);
641  //std::cout << a << ref << std::endl;
642  COMPARE(Vc::round(a), ref);
643  }
644 }
645 /*}}}*/
646 template<typename Vec> void testReduceMin()/*{{{*/
647 {
648  typedef typename Vec::EntryType T;
649  const T one = 1;
651  T *data = mem;
652  for (int i = 0; i < Vec::Size * Vec::Size; ++i) {
653  data[i] = i % (Vec::Size + 1) + one;
654  }
655  for (int i = 0; i < Vec::Size; ++i, data += Vec::Size) {
656  const Vec a(&data[0]);
657  //std::cout << a << std::endl;
658  COMPARE(a.min(), one);
659  }
660 }
661 /*}}}*/
662 template<typename Vec> void testReduceMax()/*{{{*/
663 {
664  typedef typename Vec::EntryType T;
665  const T max = Vec::Size + 1;
667  T *data = mem;
668  for (int i = 0; i < Vec::Size * Vec::Size; ++i) {
669  data[i] = (i + Vec::Size) % (Vec::Size + 1) + 1;
670  }
671  for (int i = 0; i < Vec::Size; ++i, data += Vec::Size) {
672  const Vec a(&data[0]);
673  //std::cout << a << std::endl;
674  COMPARE(a.max(), max);
675  }
676 }
677 /*}}}*/
678 template<typename Vec> void testReduceProduct()/*{{{*/
679 {
680  enum {
681  Max = Vec::Size > 8 ? Vec::Size / 2 : Vec::Size
682  };
683  typedef typename Vec::EntryType T;
684  int _product = 1;
685  for (int i = 1; i < Vec::Size; ++i) {
686  _product *= (i % Max) + 1;
687  }
688  const T product = _product;
689  VectorMemoryHelper<Vec> mem(Vec::Size);
690  T *data = mem;
691  for (int i = 0; i < Vec::Size * Vec::Size; ++i) {
692  data[i] = ((i + (i / Vec::Size)) % Max) + 1;
693  }
694  for (int i = 0; i < Vec::Size; ++i, data += Vec::Size) {
695  const Vec a(&data[0]);
696  //std::cout << a << std::endl;
697  COMPARE(a.product(), product);
698  }
699 }
700 /*}}}*/
701 template<typename Vec> void testReduceSum()/*{{{*/
702 {
703  typedef typename Vec::EntryType T;
704  int _sum = 1;
705  for (int i = 2; i <= Vec::Size; ++i) {
706  _sum += i;
707  }
708  const T sum = _sum;
710  T *data = mem;
711  for (int i = 0; i < Vec::Size * Vec::Size; ++i) {
712  data[i] = (i + i / Vec::Size) % Vec::Size + 1;
713  }
714  for (int i = 0; i < Vec::Size; ++i, data += Vec::Size) {
715  const Vec a(&data[0]);
716  //std::cout << a << std::endl;
717  COMPARE(a.sum(), sum);
718  }
719 }
720 /*}}}*/
721 template<typename V> void testExponent()/*{{{*/
722 {
723  typedef typename V::EntryType T;
724  Vc::Memory<V, 32> input;
725  Vc::Memory<V, 32> expected;
726  input[ 0] = T(0.25); expected[ 0] = T(-2);
727  input[ 1] = T( 1); expected[ 1] = T( 0);
728  input[ 2] = T( 2); expected[ 2] = T( 1);
729  input[ 3] = T( 3); expected[ 3] = T( 1);
730  input[ 4] = T( 4); expected[ 4] = T( 2);
731  input[ 5] = T( 0.5); expected[ 5] = T(-1);
732  input[ 6] = T( 6); expected[ 6] = T( 2);
733  input[ 7] = T( 7); expected[ 7] = T( 2);
734  input[ 8] = T( 8); expected[ 8] = T( 3);
735  input[ 9] = T( 9); expected[ 9] = T( 3);
736  input[10] = T( 10); expected[10] = T( 3);
737  input[11] = T( 11); expected[11] = T( 3);
738  input[12] = T( 12); expected[12] = T( 3);
739  input[13] = T( 13); expected[13] = T( 3);
740  input[14] = T( 14); expected[14] = T( 3);
741  input[15] = T( 15); expected[15] = T( 3);
742  input[16] = T( 16); expected[16] = T( 4);
743  input[17] = T( 17); expected[17] = T( 4);
744  input[18] = T( 18); expected[18] = T( 4);
745  input[19] = T( 19); expected[19] = T( 4);
746  input[20] = T( 20); expected[20] = T( 4);
747  input[21] = T( 21); expected[21] = T( 4);
748  input[22] = T( 22); expected[22] = T( 4);
749  input[23] = T( 23); expected[23] = T( 4);
750  input[24] = T( 24); expected[24] = T( 4);
751  input[25] = T( 25); expected[25] = T( 4);
752  input[26] = T( 26); expected[26] = T( 4);
753  input[27] = T( 27); expected[27] = T( 4);
754  input[28] = T( 28); expected[28] = T( 4);
755  input[29] = T( 29); expected[29] = T( 4);
756  input[30] = T( 32); expected[30] = T( 5);
757  input[31] = T( 31); expected[31] = T( 4);
758  for (size_t i = 0; i < input.vectorsCount(); ++i) {
759  COMPARE(V(input.vector(i)).exponent(), V(expected.vector(i)));
760  }
761 }
762 /*}}}*/
763 template<typename T> struct _ExponentVector { typedef int_v Type; };
764 template<> struct _ExponentVector<sfloat_v> { typedef short_v Type; };
765 
766 template<typename V> void testFrexp()/*{{{*/
767 {
768  typedef typename V::EntryType T;
769  typedef typename _ExponentVector<V>::Type ExpV;
770  Vc::Memory<V, 32> input;
771  Vc::Memory<V, 32> expectedFraction;
772  Vc::Memory<ExpV, 32> expectedExponent;
773  input[ 0] = T(0.25); expectedFraction[ 0] = T(.5 ); expectedExponent[ 0] = -1;
774  input[ 1] = T( 1); expectedFraction[ 1] = T(.5 ); expectedExponent[ 1] = 1;
775  input[ 2] = T( 0); expectedFraction[ 2] = T(0. ); expectedExponent[ 2] = 0;
776  input[ 3] = T( 3); expectedFraction[ 3] = T(.75 ); expectedExponent[ 3] = 2;
777  input[ 4] = T( 4); expectedFraction[ 4] = T(.5 ); expectedExponent[ 4] = 3;
778  input[ 5] = T( 0.5); expectedFraction[ 5] = T(.5 ); expectedExponent[ 5] = 0;
779  input[ 6] = T( 6); expectedFraction[ 6] = T( 6./8. ); expectedExponent[ 6] = 3;
780  input[ 7] = T( 7); expectedFraction[ 7] = T( 7./8. ); expectedExponent[ 7] = 3;
781  input[ 8] = T( 8); expectedFraction[ 8] = T( 8./16.); expectedExponent[ 8] = 4;
782  input[ 9] = T( 9); expectedFraction[ 9] = T( 9./16.); expectedExponent[ 9] = 4;
783  input[10] = T( 10); expectedFraction[10] = T(10./16.); expectedExponent[10] = 4;
784  input[11] = T( 11); expectedFraction[11] = T(11./16.); expectedExponent[11] = 4;
785  input[12] = T( 12); expectedFraction[12] = T(12./16.); expectedExponent[12] = 4;
786  input[13] = T( 13); expectedFraction[13] = T(13./16.); expectedExponent[13] = 4;
787  input[14] = T( 14); expectedFraction[14] = T(14./16.); expectedExponent[14] = 4;
788  input[15] = T( 15); expectedFraction[15] = T(15./16.); expectedExponent[15] = 4;
789  input[16] = T( 16); expectedFraction[16] = T(16./32.); expectedExponent[16] = 5;
790  input[17] = T( 17); expectedFraction[17] = T(17./32.); expectedExponent[17] = 5;
791  input[18] = T( 18); expectedFraction[18] = T(18./32.); expectedExponent[18] = 5;
792  input[19] = T( 19); expectedFraction[19] = T(19./32.); expectedExponent[19] = 5;
793  input[20] = T( 20); expectedFraction[20] = T(20./32.); expectedExponent[20] = 5;
794  input[21] = T( 21); expectedFraction[21] = T(21./32.); expectedExponent[21] = 5;
795  input[22] = T( 22); expectedFraction[22] = T(22./32.); expectedExponent[22] = 5;
796  input[23] = T( 23); expectedFraction[23] = T(23./32.); expectedExponent[23] = 5;
797  input[24] = T( 24); expectedFraction[24] = T(24./32.); expectedExponent[24] = 5;
798  input[25] = T( 25); expectedFraction[25] = T(25./32.); expectedExponent[25] = 5;
799  input[26] = T( 26); expectedFraction[26] = T(26./32.); expectedExponent[26] = 5;
800  input[27] = T( 27); expectedFraction[27] = T(27./32.); expectedExponent[27] = 5;
801  input[28] = T( 28); expectedFraction[28] = T(28./32.); expectedExponent[28] = 5;
802  input[29] = T( 29); expectedFraction[29] = T(29./32.); expectedExponent[29] = 5;
803  input[30] = T( 32); expectedFraction[30] = T(32./64.); expectedExponent[30] = 6;
804  input[31] = T( 31); expectedFraction[31] = T(31./32.); expectedExponent[31] = 5;
805  for (size_t i = 0; i < input.vectorsCount(); ++i) {
806  const V v = input.vector(i);
807  ExpV exp;
808  COMPARE(frexp(v, &exp), V(expectedFraction.vector(i)));
809  if (V::Size * 2 == ExpV::Size) {
810  for (size_t j = 0; j < V::Size; ++j) {
811  COMPARE(exp[j * 2], expectedExponent[i * V::Size + j]);
812  }
813  } else {
814  COMPARE(exp, ExpV(expectedExponent.vector(i)));
815  }
816  }
817 }
818 /*}}}*/
819 template<typename V> void testLdexp()/*{{{*/
820 {
821  typedef typename V::EntryType T;
822  typedef typename _ExponentVector<V>::Type ExpV;
823  for (size_t i = 0; i < 1024 / V::Size; ++i) {
824  const V v = (V::Random() - T(0.5)) * T(1000);
825  ExpV e;
826  const V m = frexp(v, &e);
827  COMPARE(ldexp(m, e), v) << ", m = " << m << ", e = " << e;
828  }
829 }
830 /*}}}*/
831 #include "ulp.h"
832 template<typename V> void testUlpDiff()/*{{{*/
833 {
834  typedef typename V::EntryType T;
835 
839  for (size_t count = 0; count < 1024 / V::Size; ++count) {
840  const V base = (V::Random() - T(0.5)) * T(1000);
842  Vc::frexp(base, &exp);
843  const V eps = ldexp(V(std::numeric_limits<T>::epsilon()), exp - 1);
844  //std::cout << base << ", " << exp << ", " << eps << std::endl;
845  for (int i = -10000; i <= 10000; ++i) {
846  const V i_v = V(T(i));
847  const V diff = base + i_v * eps;
848 
849  // if diff and base have a different exponent then ulpDiffToReference has an uncertainty
850  // of +/-1
851  const V ulpDifference = ulpDiffToReference(diff, base);
852  const V expectedDifference = Vc::abs(i_v);
853  const V maxUncertainty = Vc::abs(abs(diff).exponent() - abs(base).exponent());
854 
855  VERIFY(Vc::abs(ulpDifference - expectedDifference) <= maxUncertainty)
856  << ", base = " << base << ", epsilon = " << eps << ", diff = " << diff;
857  for (int k = 0; k < V::Size; ++k) {
858  VERIFY(std::abs(ulpDifference[k] - expectedDifference[k]) <= maxUncertainty[k]);
859  }
860  }
861  }
862 }/*}}}*/
863 
864 int main(int argc, char **argv)/*{{{*/
865 {
866  initTest(argc, argv);
867 
868  Denormals<float>::data = Vc::malloc<float, Vc::AlignOnVector>(NDenormals);/*{{{*/
869  Denormals<float>::data[0] = std::numeric_limits<float>::denorm_min();
870  for (int i = 1; i < NDenormals; ++i) {
871  Denormals<float>::data[i] = Denormals<float>::data[i - 1] * 2.173f;
872  }
873  Denormals<double>::data = Vc::malloc<double, Vc::AlignOnVector>(NDenormals);
874  Denormals<double>::data[0] = std::numeric_limits<double>::denorm_min();
875  for (int i = 1; i < NDenormals; ++i) {
876  Denormals<double>::data[i] = Denormals<double>::data[i - 1] * 2.173;
877  }/*}}}*/
878 
882 
883  runTest(testAbs<int_v>);
884  runTest(testAbs<float_v>);
885  runTest(testAbs<double_v>);
886  runTest(testAbs<short_v>);
887  runTest(testAbs<sfloat_v>);
888 
890 
898 
899  runTest(testMax<int_v>);
900  runTest(testMax<uint_v>);
901  runTest(testMax<float_v>);
902  runTest(testMax<double_v>);
903  runTest(testMax<short_v>);
904  runTest(testMax<ushort_v>);
905  runTest(testMax<sfloat_v>);
906 
918 
919  runTest(testReduceMin<float_v>);
920  runTest(testReduceMin<sfloat_v>);
921  runTest(testReduceMin<double_v>);
922  runTest(testReduceMin<int_v>);
923  runTest(testReduceMin<uint_v>);
924  runTest(testReduceMin<short_v>);
925  runTest(testReduceMin<ushort_v>);
926 
927  runTest(testReduceMax<float_v>);
928  runTest(testReduceMax<sfloat_v>);
929  runTest(testReduceMax<double_v>);
930  runTest(testReduceMax<int_v>);
931  runTest(testReduceMax<uint_v>);
932  runTest(testReduceMax<short_v>);
933  runTest(testReduceMax<ushort_v>);
934 
935  runTest(testReduceProduct<float_v>);
936  runTest(testReduceProduct<sfloat_v>);
937  runTest(testReduceProduct<double_v>);
938  runTest(testReduceProduct<int_v>);
939  runTest(testReduceProduct<uint_v>);
940  runTest(testReduceProduct<short_v>);
941  runTest(testReduceProduct<ushort_v>);
942 
943  runTest(testReduceSum<float_v>);
944  runTest(testReduceSum<sfloat_v>);
945  runTest(testReduceSum<double_v>);
946  runTest(testReduceSum<int_v>);
947  runTest(testReduceSum<uint_v>);
948  runTest(testReduceSum<short_v>);
949  runTest(testReduceSum<ushort_v>);
950 
953 
954  return 0;
955 }/*}}}*/
956 
957 // vim: foldmethod=marker
void testSin()
Definition: math.cpp:390
void testRSqrt()
Definition: math.cpp:357
VECTOR_NAMESPACE::sfloat_v sfloat_v
Definition: vector.h:82
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
Definition: vector.h:433
unsigned int hex
Definition: math.cpp:442
V apply_v(VC_ALIGNED_PARAMETER(V) x, typename V::EntryType(func)(typename V::EntryType))
Definition: math.cpp:127
#define FUZZY_COMPARE(a, b)
Definition: unittest.h:506
static Vc_ALWAYS_INLINE float_v trunc(float_v::AsArg v)
Definition: math.h:95
void initTest(int argc, char **argv)
Definition: unittest.h:169
void testFrexp()
Definition: math.cpp:766
const char * filename< float, Log2 >()
Definition: math.cpp:73
void setFuzzyness< double >(double fuzz)
Definition: unittest.h:190
const char * filename< float, Atan >()
Definition: math.cpp:65
const char * Size
Definition: TXMLSetup.cxx:56
double T(double x)
Definition: ChebyshevPol.h:34
const char * filename< double, Atan >()
Definition: math.cpp:66
void testFloor()
Definition: math.cpp:182
void testReduceSum()
Definition: math.cpp:701
const union @200 INF
static Vc_ALWAYS_INLINE Vector< T > reciprocal(const Vector< T > &x)
Definition: vector.h:451
static const char * filename()
void testCeil()
Definition: math.cpp:196
_VC_CONSTEXPR size_t vectorsCount() const
Definition: memory.h:169
void testSqrt()
Definition: math.cpp:347
TArc * a
Definition: textangle.C:12
void testAtan()
Definition: math.cpp:452
void testTrunc()
Definition: math.cpp:168
static Vc_ALWAYS_INLINE Vector< T >::Mask isfinite(const Vector< T > &x)
Definition: vector.h:454
#define FillHelperMemory(code)
Definition: math.cpp:338
double cos(double)
void testExponent()
Definition: math.cpp:721
void testInf()
Definition: math.cpp:590
static Array< Reference< T > > referenceData()
Definition: math.cpp:100
int main(int argc, char **argv)
Definition: math.cpp:864
VECTOR_NAMESPACE::short_v short_v
Definition: vector.h:90
void testAsin()
Definition: math.cpp:424
void testAbs()
Definition: math.cpp:136
Definition: math.cpp:60
const char * filename< double, Log2 >()
Definition: math.cpp:74
double sqrt(double)
const char * filename< double, Asin >()
Definition: math.cpp:68
const char * filename< double, Log >()
Definition: math.cpp:72
static float my_log2(float x)
Definition: math.cpp:262
void testLog()
Definition: math.cpp:223
Double_t x[n]
Definition: legend1.C:17
#define COMPARE(a, b)
Definition: unittest.h:509
double_v frexp(double_v::AsArg v, int_v *e)
splits v into exponent and mantissa, the sign is kept with the mantissa
Definition: math.h:38
void testNaN()
Definition: math.cpp:600
double log10(double)
void setFuzzyness< float >(float fuzz)
Definition: unittest.h:189
Definition: math.cpp:60
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
Definition: math.cpp:60
#define testRealTypes(name)
Definition: unittest.h:51
Vc_ALWAYS_INLINE Vc_PURE VectorPointerHelper< V, AlignedFlag > vector(size_t i)
Definition: memorybase.h:310
double sin(double)
void testExp()
Definition: math.cpp:210
Vc::int_v Type
Definition: ulp.h:58
const char * filename< float, Log >()
Definition: math.cpp:71
void testLdexp()
Definition: math.cpp:819
int isnan(double)
#define VC_ALIGNED_PARAMETER(_Type)
Definition: macros.h:368
VECTOR_NAMESPACE::int_v int_v
Definition: vector.h:86
#define Vc_buildDouble(sign, mantissa, exponent)
Definition: macros.h:283
ROOT::R::TRInterface & r
Definition: Object.C:4
const char * Array
SVector< double, 2 > v
Definition: Dict.h:5
double Pi()
Mathematical constants.
Definition: Math.h:68
void testReciprocal()
Definition: math.cpp:558
Definition: math.cpp:60
const char * filename< float, Log10 >()
Definition: math.cpp:75
static Vc_ALWAYS_INLINE Vector< T > round(const Vector< T > &x)
Definition: vector.h:452
const char * filename< double, Sincos >()
Definition: math.cpp:64
TMarker * m
Definition: textangle.C:8
Definition: math.cpp:60
#define ATAN_COMPARE
Definition: math.cpp:449
double floor(double)
void testSincos()
Definition: math.cpp:367
static Vc_ALWAYS_INLINE Vector< T > rsqrt(const Vector< T > &x)
Definition: vector.h:449
double asin(double)
void testRound()
Definition: math.cpp:619
void testUlpDiff()
Definition: math.cpp:832
#define VERIFY(cond)
Definition: unittest.h:515
REAL epsilon
Definition: triangle.c:617
Type
enumeration specifying the integration types.
A helper class for fixed-size two-dimensional arrays.
Definition: memory.h:120
const char * filename< double, Log10 >()
Definition: math.cpp:76
void Random()
Definition: utils.cpp:154
#define FAIL()
Definition: unittest.h:518
void testReduceProduct()
Definition: math.cpp:678
double atan2(double, double)
void testReduceMax()
Definition: math.cpp:662
static Array< SincosReference< T > > sincosReference()
Definition: math.cpp:79
Double_t y[n]
Definition: legend1.C:17
double atan(double)
double func(double *x, double *p)
Definition: stressTF1.cxx:213
static float my_trunc(float x)
Definition: math.cpp:146
Definition: math.cpp:60
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
Vc_INTRINSIC Vc_CONST m256 exponent(param256 v)
Definition: vectorhelper.h:37
void testReduceMin()
Definition: math.cpp:646
static T ulpDiffToReference(T val, T ref)
Definition: ulp.h:34
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
double ceil(double)
void isNegative()
Definition: math.cpp:581
Definition: casts.h:28
const char * filename< float, Sincos >()
Definition: math.cpp:63
void testMax()
Definition: math.cpp:320
#define I(x, y, z)
double exp(double)
const char * filename< float, Asin >()
Definition: math.cpp:67
void testLog2()
Definition: math.cpp:266
float value
Definition: math.cpp:443
void testCos()
Definition: math.cpp:407
Definition: math.cpp:60
void testLog10()
Definition: math.cpp:297
double log(double)
double ldexp(double, int)
void testAtan2()
Definition: math.cpp:483
#define runTest(name)
Definition: unittest.h:42