From dedde725d2b57b0d196b46fe6fec89ba84d05d23 Mon Sep 17 00:00:00 2001
From: Raph Levien <raph@google.com>
Date: Tue, 11 Jun 2013 23:41:02 -0700
Subject: [PATCH] Add tanh function

Add a fast implementation of the tanh function. This will be useful
for distortion.
---
 cpp/src/exp2.cc     | 29 +++++++++++++++++++++++++++++
 cpp/src/exp2.h      | 34 ++++++++++++++++++++++++++++++++++
 cpp/src/main.cc     | 17 +++++++++++++++++
 cpp/src/pitchenv.cc |  2 --
 4 files changed, 80 insertions(+), 2 deletions(-)

diff --git a/cpp/src/exp2.cc b/cpp/src/exp2.cc
index 7f3a93d..ac0ae43 100644
--- a/cpp/src/exp2.cc
+++ b/cpp/src/exp2.cc
@@ -19,6 +19,8 @@
 #include "synth.h"
 #include "exp2.h"
 
+#include <stdio.h>
+
 int32_t exp2tab[EXP2_N_SAMPLES << 1];
 
 void Exp2::init() {
@@ -34,3 +36,30 @@ void Exp2::init() {
   exp2tab[(EXP2_N_SAMPLES << 1) - 2] = (1U << 31) - exp2tab[(EXP2_N_SAMPLES << 1) - 1];
 }
 
+int32_t tanhtab[TANH_N_SAMPLES << 1];
+
+static double dtanh(double y) {
+  return 1 - y * y;
+}
+
+void Tanh::init() {
+  double step = 4.0 / TANH_N_SAMPLES;
+  double y = 0;
+  for (int i = 0; i < TANH_N_SAMPLES; i++) {
+    tanhtab[(i << 1) + 1] = (1 << 24) * y + 0.5;
+    //printf("%d\n", tanhtab[(i << 1) + 1]);
+    // Use a basic 4th order Runge-Kutte to compute tanh from its
+    // differential equation.
+    double k1 = dtanh(y);
+    double k2 = dtanh(y + 0.5 * step * k1);
+    double k3 = dtanh(y + 0.5 * step * k2);
+    double k4 = dtanh(y + step * k3);
+    double dy = (step / 6) * (k1 + k4 + 2 * (k2 + k3));
+    y += dy;
+  }
+  for (int i = 0; i < TANH_N_SAMPLES - 1; i++) {
+    tanhtab[i << 1] = tanhtab[(i << 1) + 3] - tanhtab[(i << 1) + 1];
+  }
+  int32_t lasty = (1 << 24) * y + 0.5;
+  tanhtab[(TANH_N_SAMPLES << 1) - 2] = lasty - tanhtab[(TANH_N_SAMPLES << 1) - 1];
+}
diff --git a/cpp/src/exp2.h b/cpp/src/exp2.h
index 2c08fb1..016ca33 100644
--- a/cpp/src/exp2.h
+++ b/cpp/src/exp2.h
@@ -44,3 +44,37 @@ int32_t Exp2::lookup(int32_t x) {
   return y >> (6 - (x >> 24));
 }
 #endif
+
+class Tanh {
+ public:
+  static void init();
+
+  // Q24 in, Q24 out
+  static int32_t lookup(int32_t x);
+};
+
+#define TANH_LG_N_SAMPLES 10
+#define TANH_N_SAMPLES (1 << TANH_LG_N_SAMPLES)
+
+extern int32_t tanhtab[TANH_N_SAMPLES << 1];
+
+inline
+int32_t Tanh::lookup(int32_t x) {
+  int32_t signum = x >> 31;
+  x ^= signum;
+  if (x >= (4 << 24)) {
+    if (x >= (17 << 23)) {
+      return signum ^ (1 << 24);
+    }
+    int32_t sx = ((int64_t)-48408812 * (int64_t)x) >> 24;
+    return signum ^ ((1 << 24) - 2 * Exp2::lookup(sx));
+  } else {
+    const int SHIFT = 26 - TANH_LG_N_SAMPLES;
+    int lowbits = x & ((1 << SHIFT) - 1);
+    int x_int = (x >> (SHIFT - 1)) & ((TANH_N_SAMPLES - 1) << 1);
+    int dy = tanhtab[x_int];
+    int y0 = tanhtab[x_int + 1];
+    int y = y0 + (((int64_t)dy * (int64_t)lowbits) >> SHIFT);
+    return y ^ signum;
+  }
+}
diff --git a/cpp/src/main.cc b/cpp/src/main.cc
index e20ca8a..ee45d59 100644
--- a/cpp/src/main.cc
+++ b/cpp/src/main.cc
@@ -81,6 +81,22 @@ void test_log_accuracy() {
   cout << "Max error: " << maxerr << endl;
 }
 
+void test_tanh_accuracy() {
+  Tanh::init();
+  double maxerr = 0;
+  for (int i = 0; i < 1000000; i++) {
+    int32_t x = (rand() & ((1 << 29) - 1)) - (1 << 28);
+    int32_t y = Tanh::lookup(x);
+    double yd = (1 << 24) * tanh(x * (1.0 / (1 << 24)));
+    double err = fabs(y - yd);
+    if (err > maxerr) {
+      cout << "x = " << x << ", y = " << y << ", yd = " << (int)yd << endl;
+      maxerr = err;
+    }
+  }
+  cout << "Max error: " << maxerr << endl;
+}
+
 void test_pure_accuracy() {
   int32_t worstfreq;
   int32_t worstphase;
@@ -264,6 +280,7 @@ int main(int argc, char **argv) {
   //FmCore::dump();
   //test_sin_accuracy();
   test_log_accuracy();
+  test_tanh_accuracy();
   //benchmark_fm_op();
   //test_pure_accuracy();
   //benchmark_sin();
diff --git a/cpp/src/pitchenv.cc b/cpp/src/pitchenv.cc
index ee91b03..7cb564c 100644
--- a/cpp/src/pitchenv.cc
+++ b/cpp/src/pitchenv.cc
@@ -17,8 +17,6 @@
 #include "synth.h"
 #include "pitchenv.h"
 
-using namespace std;
-
 int PitchEnv::unit_;
 
 void PitchEnv::init(double sample_rate) {