PndCAMath.h

Go to the documentation of this file.

//-*- Mode: C++ -*-
// ************************************************************************
// This file is property of and copyright by the ALICE HLT Project        *
// ALICE Experiment at CERN, All rights reserved.                         *
// See cxx source for full Copyright notice                               *
//                                                                        *
//*************************************************************************

#ifndef PNDCAMATH_H
#define PNDCAMATH_H

#include "PndCADef.h"

#include <cstdlib>
#include <algorithm>
using std::abs;

#if defined(HLTCA_STANDALONE)
#include <limits>
#else
#include "TMath.h"
#endif

#ifdef __SSE__
#include <xmmintrin.h>
#endif

namespace CAMath {
template <typename T>
static inline T Min(const T &x, const T &y)
{
  return std::min(x, y);
}
template <typename T>
static inline T Max(const T &x, const T &y)
{
  return std::max(x, y);
}
template <typename T>
static inline T Sqrt(const T &x)
{
  return std::sqrt(x);
}
template <typename T>
static inline T RSqrt(const T &x)
{
  const T one = T(1.f);
  return one / std::sqrt(x);
}
template <typename T>
static inline T Abs(const T &x)
{
  return std::abs(x);
}
template <typename T>
static inline T Log(const T &x)
{
  return std::log(x);
}
template <typename T>
static inline T Log10(const T &x)
{
  return std::log10(x);
}
template <typename T>
static inline T Sin(const T &x)
{
  return std::sin(x);
}
template <typename T>
static inline T Cos(const T &x)
{
  return std::cos(x);
}
template <typename T>
static T Reciprocal(const T &x);
template <typename T>
static T ApproxSqrt(const T &x);
#ifdef USE_TBB
template <typename T>
static T AtomicMax(T volatile *addr, T val);
#endif // USE_TBB

template <typename T>
struct FiniteReturnTypeHelper {
  typedef bool R;
};
template <typename T>
static typename FiniteReturnTypeHelper<T>::R Finite(const T &x);

template <typename T>
static T Round(const T &x);

template <typename T>
static inline T Recip(const T &x)
{
  return T(1) / x;
}
template <typename T>
static T ATan2(const T &y, const T &x);
template <typename T>
static T ASin(const T &x);

float Tan(float x);
float Copysign(float x, float y);
static inline float TwoPi()
{
  return 6.283185307179586f;
}
static inline float Pi()
{
  return 3.1415926535897f;
}
int Nint(float x);

#ifdef USE_TBB
int AtomicExch(int volatile *addr, int val);
int AtomicAdd(int volatile *addr, int val);
int AtomicMin(int volatile *addr, int val);
#endif // USE_TBB
} // namespace CAMath

#if defined(HLTCA_STANDALONE)
#define choice(c1, c2, c3) c2
#else
#define choice(c1, c2, c3) c3
#endif

namespace CAMath {

template <>
inline float Reciprocal<float>(const float &x)
{
  return 1.f / x;
}

template <>
inline double Reciprocal<double>(const double &x)
{
  return 1. / x;
}

#ifdef __SSE__
template <>
inline float RSqrt<float>(const float &x)
{
  float r = x;
  __m128 tmp;
  asm("rsqrtss %0,%1\n\t"
      "movss %1,%0\n\t"
      : "+m"(r), "=x"(tmp));
  return r;
}
#endif // __SSE__

template <>
inline float ApproxSqrt<float>(const float &x)
{
  float r = x;
  asm("shr %0\n\t"
      "add $0x1fc00000,%0\n\t"
      : "+r"(r));
  return r;
}
} // namespace CAMath

inline int CAMath::Nint(float x)
{
#if defined(HLTCA_STANDALONE)
  int i;
  if (x >= 0) {
    i = int(x + 0.5f);
    if (x + 0.5f == float(i) && i & 1)
      i--;
  } else {
    i = int(x - 0.5f);
    if (x - 0.5f == float(i) && i & 1)
      i++;
  }
  return i;
#else
  return TMath::Nint(x);
#endif
}

namespace CAMath {
template <>
inline bool Finite<float>(const float &x)
{
  return choice(1, x < std::numeric_limits<float>::infinity() && -x < std::numeric_limits<float>::infinity(), TMath::Finite(x));
}

template <>
inline float Round<float>(const float &x)
{
  return static_cast<float>(Nint(x));
}

template <>
inline float ATan2<float>(const float &y, const float &x)
{
  return choice(atan2f(y, x), atan2(y, x), TMath::ATan2(y, x));
}

template <>
inline float ASin(const float &x)
{
  return choice(asinf(x), asin(x), TMath::ASin(x));
}

} // namespace CAMath

inline float CAMath::Copysign(float x, float y)
{
  x = CAMath::Abs(x);
  return (y >= 0) ? x : -x;
}

inline float CAMath::Tan(float x)
{
  return choice(tanf(x), tan(x), TMath::Tan(x));
}

#ifdef USE_TBB

#include <tbb/atomic.h>

inline int CAMath::AtomicExch(int volatile *addr, int val)
{
  tbb::atomic<int> &a = *reinterpret_cast<tbb::atomic<int> *>(const_cast<int *>(addr));
  return a.fetch_and_store(val);
}

inline int CAMath::AtomicAdd(int volatile *addr, int val)
{
  tbb::atomic<int> &a = *reinterpret_cast<tbb::atomic<int> *>(const_cast<int *>(addr));
  return a.fetch_and_add(val);
}

namespace CAMath {
template <>
inline int AtomicMax<int>(int volatile *addr, int val)
{
  tbb::atomic<int> &a = *reinterpret_cast<tbb::atomic<int> *>(const_cast<int *>(addr));
  int old = a;
  if (old < val) {
    while (old != a.compare_and_swap(val, old)) {
      old = a;
      if (old >= val) {
        break;
      }
    }
  }
  return old;
}

template <>
inline unsigned int AtomicMax<unsigned int>(unsigned int volatile *addr, unsigned int val)
{
  tbb::atomic<unsigned int> &a = *reinterpret_cast<tbb::atomic<unsigned int> *>(const_cast<unsigned int *>(addr));
  unsigned int old = a;
  if (old < val) {
    while (old != a.compare_and_swap(val, old)) {
      old = a;
      if (old >= val) {
        break;
      }
    }
  }
  return old;
}
} // namespace CAMath

inline int CAMath::AtomicMin(int volatile *addr, int val)
{
  tbb::atomic<int> &a = *reinterpret_cast<tbb::atomic<int> *>(const_cast<int *>(addr));
  int old = a;
  if (old > val) {
    while (old != a.compare_and_swap(val, old)) {
      old = a;
      if (old <= val) {
        break;
      }
    }
  }
  return old;
}
#endif // USE_TBB

#undef choice

#endif