<commit_before>﻿
#pragma once

#include "quantities/elementary_functions.hpp"

#include <pmmintrin.h>

#include <cmath>
#include <type_traits>

#include "quantities/si.hpp"
#include "numerics/cbrt.hpp"

namespace principia {
namespace quantities {
namespace internal_elementary_functions {

using si::Radian;

template<typename Q1, typename Q2>
Product<Q1, Q2> FusedMultiplyAdd(Q1 const& x,
                                 Q2 const& y,
                                 Product<Q1, Q2> const& z) {
  return si::Unit<Product<Q1, Q2>> * std::fma(x / si::Unit<Q1>,
                                              y / si::Unit<Q2>,
                                              z / si::Unit<Product<Q1, Q2>>);
}

template<typename Q>
FORCE_INLINE(inline) Q Abs(Q const& quantity) {
  return si::Unit<Q> * std::abs(quantity / si::Unit<Q>);
}

template<typename Q>
Q Mod(Q const& argument, Q const& modulus) {
  double const result =
      std::fmod(argument / si::Unit<Q>, modulus / si::Unit<Q>);
  if (result > 0.0) {
    return result * si::Unit<Q>;
  } else {
    return result * si::Unit<Q> + modulus;
  }
}

template<typename Q>
SquareRoot<Q> Sqrt(Q const& x) {
#if PRINCIPIA_USE_SSE3_INTRINSICS
  auto const x_128d = _mm_set_sd(x / si::Unit<Q>);
  return si::Unit<SquareRoot<Q>> * _mm_cvtsd_f64(_mm_sqrt_sd(x_128d, x_128d));
#else
  return si::Unit<SquareRoot<Q>> * std::sqrt(x / si::Unit<Q>);
#endif
}

template<typename Q>
CubeRoot<Q> Cbrt(Q const& x) {
  return si::Unit<CubeRoot<Q>> * numerics::Cbrt(x / si::Unit<Q>);
}

template<int exponent>
constexpr double Pow(double x) {
  return std::pow(x, exponent);
}

// Static specializations for frequently-used exponents, so that this gets
// turned into multiplications at compile time.

template<>
inline constexpr double Pow<-3>(double x) {
  return 1 / (x * x * x);
}

template<>
inline constexpr double Pow<-2>(double x) {<commit_msg>A faster pow using the Russian peasant method.<commit_after>﻿
#pragma once

#include "quantities/elementary_functions.hpp"

#include <pmmintrin.h>

#include <cmath>
#include <type_traits>

#include "quantities/si.hpp"
#include "numerics/cbrt.hpp"

namespace principia {
namespace quantities {
namespace internal_elementary_functions {

using si::Radian;

template<typename Q1, typename Q2>
Product<Q1, Q2> FusedMultiplyAdd(Q1 const& x,
                                 Q2 const& y,
                                 Product<Q1, Q2> const& z) {
  return si::Unit<Product<Q1, Q2>> * std::fma(x / si::Unit<Q1>,
                                              y / si::Unit<Q2>,
                                              z / si::Unit<Product<Q1, Q2>>);
}

template<typename Q>
FORCE_INLINE(inline) Q Abs(Q const& quantity) {
  return si::Unit<Q> * std::abs(quantity / si::Unit<Q>);
}

template<typename Q>
Q Mod(Q const& argument, Q const& modulus) {
  double const result =
      std::fmod(argument / si::Unit<Q>, modulus / si::Unit<Q>);
  if (result > 0.0) {
    return result * si::Unit<Q>;
  } else {
    return result * si::Unit<Q> + modulus;
  }
}

template<typename Q>
SquareRoot<Q> Sqrt(Q const& x) {
#if PRINCIPIA_USE_SSE3_INTRINSICS
  auto const x_128d = _mm_set_sd(x / si::Unit<Q>);
  return si::Unit<SquareRoot<Q>> * _mm_cvtsd_f64(_mm_sqrt_sd(x_128d, x_128d));
#else
  return si::Unit<SquareRoot<Q>> * std::sqrt(x / si::Unit<Q>);
#endif
}

template<typename Q>
CubeRoot<Q> Cbrt(Q const& x) {
  return si::Unit<CubeRoot<Q>> * numerics::Cbrt(x / si::Unit<Q>);
}

template<int exponent>
constexpr double Pow(double x) {
  // Use the Russian peasant algorithm for small exponents.
  if constexpr (exponent > 0 && exponent < 32) {
    // The end of the recursion is handled by the specializations below.
    auto const y = Pow<exponent / 2>(x);
    auto const y² = y * y;
    if constexpr (exponent % 2 == 1) {
      return y² * x;
    } else {
      return y²;
    }
  } else if constexpr (exponent < 0 && exponent > -32) {
    return 1 / Pow<-exponent>(x);
  } else {
    return std::pow(x, exponent);
  }
}

// Static specializations for frequently-used exponents, so that this gets
// turned into multiplications at compile time.

template<>
inline constexpr double Pow<-3>(double x) {
  return 1 / (x * x * x);
}

template<>
inline constexpr double Pow<-2>(double x) {
{'fraction_non_alphanumeric': 0.14279015240328252, 'fraction_numerical': 0.016881594372801876, 'mean_word_length': 3.0940499040307103, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 85, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4427733', 'n_tokens_mistral': 1557, 'n_tokens_neox': 1483, 'n_words': 386}