source: branches/2989_MovingPeaksBenchmark/HeuristicLab.ExtLibs/HeuristicLab.NativeInterpreter/0.1/NativeInterpreter-0.1/lib/vdt/tanh.h @ 16605

/*
 * tanh.h
 * The basic idea is to exploit Pade polynomials.
 * Implemented by Manuel Schiller for LHCb.
 * 
 *  Created on: Sep 23, 2017
 *      Author: Paul Seyfert, Manuel Schiller
 */

/* 
 * VDT is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef TANH_H_
#define TANH_H_

#include "vdtcore_common.h"

namespace vdt{




/// Fast tanh implementation double precision
inline double fast_tanh(double x){
  // for very large |x| > 20, tanh(x) is x/|x| anyway (at least to double
  // precision)
  //
  // NB: branch-free code takes longer to execute
  if (std::abs(x) > 20.) return std::copysign(1., x);
  // strategy for large arguments: tanh(2x) = 2 tanh(x)/(1 + tanh^2(x))
  // idea is to use this "argument halving" a couple of times, and use a
  // very short Padé approximation for the rest of the way
  const auto xx = x * 0.125;
  const auto xx2 = xx * xx;
  const auto numer = 135135 + xx2 * (17325 + xx2 * ( 378 + xx2 *  1));
  const auto denom = 135135 + xx2 * (62370 + xx2 * (3150 + xx2 * 28));

  auto tanh = xx * numer / denom;
  tanh = 2 * tanh / (tanh * tanh + 1);
  tanh = 2 * tanh / (tanh * tanh + 1);
  return 2 * tanh / (tanh * tanh + 1);
}

//------------------------------------------------------------------------------
/// Fast tanh implementation single precision
inline float fast_tanhf( float x ) {
  // same strategy as double version above, but even shorter Padé
  // approximation is sufficient for float
  //
  // NB: branch-free code takes longer to execute
  if (std::abs(x) > 9.1f) return std::copysign(1.f, x);
  const auto xx = x * 0.125f;
  const auto xx2 = xx * xx;
  auto tanh = xx * (xx2 + 15) / (6 * xx2 + 15);
  tanh = 2 * tanh / (tanh * tanh + 1);
  tanh = 2 * tanh / (tanh * tanh + 1);
  return 2 * tanh / (tanh * tanh + 1);
}

//------------------------------------------------------------------------------
// Vector signatures

}// end of vdt

#endif // end of tanh