fixed_point_arithmetic.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import math\n",
    "\n",
    "from quantization import Tensor\n",
    "from quantization_utils import *\n",
    "from data_utils import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fixed Point\n",
    "\n",
    "---\n",
    "## What's Fixed Point?\n",
    "\n",
    "You may refer [here](https://stackoverflow.com/questions/7524838/fixed-point-vs-floating-point-number) for the difference between fixed point & floating point.\n",
    "\n",
    "We usually use [Qm.n notation](https://en.wikipedia.org/wiki/Q_%28number_format%29) to represent the fixed point.\n",
    "\n",
    "For example Q3.12 means we use 3 bits to represent the integer parts and 12 bits to represent the fractional parts (and one bit to sign).\n",
    "\n",
    "So the range for Qm.n is [-2^m, 2^(m+1) - 2^-n]\n",
    "\n",
    "---\n",
    "## Fixed Point <-> Floating point?\n",
    "\n",
    "In this case, we will represent all the fixed point using integer values (int16, int32, etc.), and the floating value is simply:\n",
    "\n",
    "$$value_{float} = value_{int} * 2^{-frac_{bits}}$$\n",
    "\n",
    "where $$frac_{bits} = total_{bits} - integer_{bits} - 1$$\n",
    "\n",
    "So convert floating point value x to Qm.n format is simply:\n",
    "\n",
    "$$value_{int} = value_{float} * 2^{n}$$\n",
    "\n",
    "---\n",
    "## Why Fixed Point matters?\n",
    "\n",
    "Fixed point arithmetics can be very useful for division, tanh, exp and others.\n",
    "\n",
    "This codelab will focus on division mainly with brief introduction to others."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fixed Point Arithmetic\n",
    "\n",
    "Let's just focus on Qm.n case first. It's easier to convert Qm.n into Qm1.n1 case (where just bits shifting).\n",
    "\n",
    "---\n",
    "## Add\n",
    "\n",
    "So we have\n",
    "\n",
    "$$lhs_{float} = lhs_{int} * 2^{-n}$$\n",
    "\n",
    "and\n",
    "\n",
    "$$rhs_{float} = rhs_{int} * 2^{-n}$$\n",
    "\n",
    "Our desired output is:\n",
    "\n",
    "$$out_{int} = (lhs_{float} + rhs_{float}) * 2^{n}$$\n",
    "\n",
    "=>\n",
    "\n",
    "$$out_{int} = (lhs_{int} * 2^{-n} + rhs_{int} * 2^{-n}) * 2^{n}$$\n",
    "\n",
    "=>\n",
    "\n",
    "$$out_{int} = (lhs_{int} + rhs_{int}) * 2^{-n} * 2^{n}$$\n",
    "\n",
    "=>\n",
    "\n",
    "$$out_{int} = lhs_{int} + rhs_{int}$$\n",
    "\n",
    "So, the output is simply the add of the int value of lhs and rhs, exactly like the int math! Very simple, right?\n",
    "\n",
    "---\n",
    "## Sub\n",
    "\n",
    "Sub is actually very similar to add, so we can just skip the deduction part :)\n",
    "\n",
    "We will have\n",
    "\n",
    "$$out_{int} = lhs_{int} - rhs_{int}$$\n",
    "\n",
    "for the sub.\n",
    "\n",
    "---\n",
    "## Mul\n",
    "\n",
    "Let's take a look at the mul:\n",
    "\n",
    "Still we have:\n",
    "\n",
    "$$lhs_{float} = lhs_{int} * 2^{-n}$$\n",
    "\n",
    "and\n",
    "\n",
    "$$rhs_{float} = rhs_{int} * 2^{-n}$$\n",
    "\n",
    "Our desired output is:\n",
    "\n",
    "$$out_{int} = (lhs_{float} * rhs_{float}) * 2^{n}$$\n",
    "\n",
    "=>\n",
    "\n",
    "$$out_{int} = (lhs_{int} * 2^{-n} * rhs_{int} * 2^{-n}) * 2^{n}$$\n",
    "\n",
    "=>\n",
    "\n",
    "$$out_{int} = (lhs_{int} * rhs_{int}) * 2^{-n} * 2^{-n} * 2^{n}$$\n",
    "\n",
    "=>\n",
    "\n",
    "$$out_{int} = (lhs_{int} * rhs_{int}) * 2^{-n}$$\n",
    "\n",
    "Given the reality that n is always >= 0, so we essentially have a normal integer mul with a right shift of n bits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Div\n",
    "\n",
    "---\n",
    "\n",
    "Div is indeed very complicated as we're limited by two constraints: 1) we're dealing with integer-only arithmetic; 2) most of the time, we don't have direct and accurate 'div' instruction available, so we really need to emulate the 'div' behavior.\n",
    "\n",
    "How are we going to do that exactly?\n",
    "\n",
    "---\n",
    "\n",
    "### Newton-Ralpson Division\n",
    "\n",
    "See the detailed reference [here](https://en.wikipedia.org/wiki/Division_algorithm#Newton.E2.80.93Raphson_division)\n",
    "\n",
    "The first observation is we can transform:\n",
    "\n",
    "$$out = lhs \\div rhs$$\n",
    "\n",
    "into $$out = lhs * (\\frac{1}{rhs})$$\n",
    "\n",
    "and for x in range of [0.5, 1], we can use the newton-ralpson method to estimate the value.\n",
    "\n",
    "Initial value is:\n",
    "\n",
    "$$x_0 = (\\frac{48}{17}) - (\\frac{32}{17}) * x$$\n",
    "\n",
    "Then we can iterate using:\n",
    "\n",
    "$$x_{new} = x_{pre} + x_{pre} * (1 - x * x_{pre})$$\n",
    "\n",
    "Normally we can get good accuracy with 2 or 3 iterations, let's just verify it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg Diff for 2 iterations is 5.553048509880476e-06\n",
      "Avg Diff for 3 iterations is 4.685110077673471e-11\n",
      "Max Diff for 2 iterations is 2.394607344258226e-05\n",
      "Max Diff for 3 iterations is 2.867073245482743e-10\n"
     ]
    }
   ],
   "source": [
    "def newton_ralpson_reci(x, iteration=2):\n",
    "  a = 48.0 / 17.0\n",
    "  b = 32.0 / 17.0\n",
    "  result = a - b * x\n",
    "  for i in range(iteration):\n",
    "    result = result + result * (1 - x * result)\n",
    "  return result\n",
    "\n",
    "test_data = np.arange(0.5, 1, 0.001)\n",
    "\n",
    "expected_result = 1.0 / test_data\n",
    "\n",
    "result_iter2 = newton_ralpson_reci(test_data, 2)\n",
    "result_iter3 = newton_ralpson_reci(test_data, 3)\n",
    "\n",
    "print(\"Avg Diff for 2 iterations is %s\" % (diff(expected_result, result_iter2)))\n",
    "print(\"Avg Diff for 3 iterations is %s\" % (diff(expected_result, result_iter3)))\n",
    "print(\"Max Diff for 2 iterations is %s\" % (max_diff(expected_result, result_iter2)))\n",
    "print(\"Max Diff for 3 iterations is %s\" % (max_diff(expected_result, result_iter3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reciprocal\n",
    "---\n",
    "\n",
    "So far so good, but wait a second, isn't newton_ralpson_reci only works for [0.5, 1]? \n",
    "\n",
    "Well, we can always \"canonicalize\" the value to make it between [0.5, 1] with a proper shift by counting the leading zeros (in arm, it's just instruction 'CLS'). And yes, that only applies to positive numbers, but it should be very cheap for negative values to apply 'NEG'.\n",
    "\n",
    "---\n",
    "#### Canonicalize\n",
    "\n",
    "So for Q0.n+m case (all bits go to the fractional bits except the sign bits), we just to make sure the number of leading zeros is 1 (which is the sign bit).\n",
    "\n",
    "So we have:\n",
    "\n",
    "$$value_{FloatQ0.n+m} = value_{IntQm.n} * 2^{m} * 2^{-(m+n)}$$\n",
    "\n",
    "It's easily to get\n",
    "\n",
    "$$value_{IntQ0.n+m} = value_{IntQm.n} * 2^{m}$$\n",
    "\n",
    "In reality, simply shifting Qm.n int value left by m bits can easily result in overflow, so we don't shift it just yet, but just do that at the end.\n",
    "\n",
    "Moving forward:\n",
    "\n",
    "We know\n",
    "\n",
    "$$value_{IntCano} = value_{IntQ0.n+m} * 2^{zeros_{leading} - 1}$$\n",
    "\n",
    "so we have\n",
    "\n",
    "$$value_{FloatCano} = value_{IntQ0.n+m} * 2^{zeros_{leading} - 1} * 2^{-(zeros_{leading} - 1)} * 2^{-(m+n)}$$\n",
    "\n",
    "=>\n",
    "\n",
    "$$value_{FloatCano} = value_{IntQm.n} * 2^{m} * 2^{zeros_{leading} - 1} * 2^{-(zeros_{leading} - 1)} * 2^{-(m+n)}$$\n",
    "\n",
    "=>\n",
    "\n",
    "$$value_{FloatCano} = (value_{IntQm.n} * 2^{zeros_{leading} - 1} * 2^{-(m+n)}) * 2^{m-(zeros_{leading} - 1)}$$\n",
    "\n",
    "Let \n",
    "\n",
    "$$value_{shifted} = value_{IntQm.n} * 2^{zeros_{leading} - 1} * 2^{-(m+n)}$$\n",
    "\n",
    "So\n",
    "\n",
    "$$\\frac{1}{value_{FloatCano}} = f(value_{shifted}) * 2^{(zeros_{leading} - 1) - m}$$\n",
    "\n",
    "and $$value_{shifted}$$ can be interpreted as the canonicalize representation for Q0.m+n format.\n",
    "\n",
    "where function f is the newton-ralpson method.\n",
    "\n",
    "---\n",
    "\n",
    "Another view of this problem, when we transform Qm.n into Q0.m+n format, we need to shift left by m bits, but we don't shift that just yet, we can shift leading_zeros - 1, so we will still need to shift left by m + 1 - leading_zeros bits.\n",
    "\n",
    "---\n",
    "\n",
    "In reality, since for newton-ralpson method we need to represent 48/17, which will require Q2.m+n-2 representation, note, from Q0.m+n to Q2.m+n-2 is just to shift right by two bits or shift left by -2 bits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clz(value):\n",
    "  assert value >= 0\n",
    "  if value.dtype == np.int32:\n",
    "    total_bits = 32\n",
    "  elif value.dtype == np.int16:\n",
    "    total_bits = 16\n",
    "  else:\n",
    "    assert False\n",
    "    \n",
    "  return total_bits - (len(bin(value)) - 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FixedPoint:\n",
    "  def __init__(self, integer, fractional, datatype=np.int32, scale=None, zero_point=None):\n",
    "    # datatype is either np.int16 or np.int32.\n",
    "    self.integer = integer\n",
    "    self.fractional = fractional\n",
    "    self.datatype = datatype\n",
    "    self.value = None\n",
    "    self.scale = scale\n",
    "    self.zero_point = zero_point\n",
    "    if self.scale:\n",
    "      self.multiplier = self.from_float(self.scale)\n",
    "    \n",
    "  def set_raw_value(self, value):\n",
    "    self.value = self.datatype(value)\n",
    "    \n",
    "  def get_raw_value(self):\n",
    "    return self.value\n",
    "    \n",
    "  def to_float(self):\n",
    "    return self.value * (2 ** (-self.fractional))\n",
    "\n",
    "  def from_float(self, value):\n",
    "    temp = value * (2 ** self.fractional)\n",
    "    if self.datatype == np.int32:\n",
    "      max_value = 2 ** 31 - 1\n",
    "      min_value = - 2 ** 31\n",
    "    elif self.datatype == np.int16:\n",
    "      max_value = 2 ** 15 - 1\n",
    "      min_value = -2 ** 15\n",
    "    else:\n",
    "      # Should not reached\n",
    "      assert False\n",
    "    \n",
    "    temp = max(min(temp, max_value), min_value)\n",
    "    return self.datatype(temp)\n",
    "    \n",
    "  def from_q_tensor(self, value):\n",
    "    # This may result in overflow\n",
    "    return (value - self.zero_point) * self.multiplier\n",
    "\n",
    "  # currently we only support same type arithmetics.\n",
    "  def check(self, another):\n",
    "    assert self.integer == another.integer\n",
    "    assert self.fractional == another.fractional\n",
    "    assert self.datatype == another.datatype\n",
    "    \n",
    "  def add(self, another):\n",
    "    self.check(another)\n",
    "    result = FixedPoint(self.integer, self.fractional, self.datatype, self.scale, self.zero_point)\n",
    "    # We need to be careful about overflow.\n",
    "    result.set_raw_value(self.get_raw_value() + another.get_raw_value())\n",
    "    return result\n",
    "\n",
    "  def sub(self, another):\n",
    "    self.check(another)\n",
    "    result = FixedPoint(self.integer, self.fractional, self.datatype, self.scale, self.zero_point)\n",
    "    # We need to be careful about overflow.\n",
    "    result.set_raw_value(self.get_raw_value() - another.get_raw_value())\n",
    "    return result\n",
    "\n",
    "  def mul(self, another):\n",
    "    self.check(another)\n",
    "    result = FixedPoint(self.integer, self.fractional, self.datatype, self.scale, self.zero_point)\n",
    "    # We need to be careful about overflow.\n",
    "    temp = np.int64(self.get_raw_value()) * np.int64(another.get_raw_value())\n",
    "    temp = rounding_divide_by_POT(temp, self.fractional)                        \n",
    "    result.set_raw_value(temp)\n",
    "    return result\n",
    "\n",
    "  def neg(self):\n",
    "    result = FixedPoint(self.integer, self.fractional, self.datatype, self.scale, self.zero_point)\n",
    "    assert self.get_raw_value() != None\n",
    "    result.set_raw_value(-self.get_raw_value())\n",
    "    return result\n",
    "\n",
    "  # this only handles cannonical format of Q2.x.\n",
    "  def _reciprocal_for_cano(self, iteration=2):\n",
    "    a = FixedPoint(2, 29, np.int32)\n",
    "    a.set_raw_value(a.from_float(48.0 / 17.0))\n",
    "    b = FixedPoint(2, 29, np.int32)\n",
    "    b.set_raw_value(a.from_float(32.0 / 17.0))\n",
    "    one = FixedPoint(2, 29, np.int32)\n",
    "    one.set_raw_value(one.from_float(1.0))\n",
    "    result = a.sub(self.mul(b))\n",
    "    for i in range(iteration):\n",
    "      temp = one.sub(self.mul(result))\n",
    "      temp = temp.mul(result)\n",
    "      result = result.add(temp)\n",
    "    return result\n",
    "\n",
    "  def reciprocal(self, iteration=2):\n",
    "    # We should also check hte raw value should not be 0.\n",
    "    if self.get_raw_value() < 0:\n",
    "      return self.neg().reciprocal().neg()\n",
    "\n",
    "    leading_zeros = clz(self.get_raw_value())\n",
    "    # We need to shift left leading_zeros - 1 bits so the raw value is in the canonical format for Q0.n.\n",
    "    # We also need to shift right by 2 bits so the raw_value is in the Q2.n case.\n",
    "    shift_left_bits = leading_zeros - 1 - 2\n",
    "    raw_value = self.get_raw_value() * (2 ** shift_left_bits)\n",
    "    cano = FixedPoint(2, 29, self.datatype)\n",
    "    cano.set_raw_value(raw_value)\n",
    "    # Note the return is actually in Q2.n format.\n",
    "    # we still need to shift left by leading_zeros bits -1 - m.\n",
    "    cano_reci = cano._reciprocal_for_cano(iteration)\n",
    "    final_shift_left_bits = leading_zeros - 1 - self.integer\n",
    "    final_raw_value = cano_reci.get_raw_value()\n",
    "    # Since we're using Q2.n, it's possible we run into shift left < -2 scenario,\n",
    "    # in that sense, we need to shift right the value.\n",
    "    if final_shift_left_bits < -2:\n",
    "      value_shift_right_bits = -2 - final_shift_left_bits\n",
    "      final_raw_value = rounding_divide_by_POT(final_raw_value, value_shift_right_bits)\n",
    "      final_shift_left_bits = -2\n",
    "    result = FixedPoint(cano_reci.integer + final_shift_left_bits,\n",
    "                        cano_reci.fractional - final_shift_left_bits,\n",
    "                        cano_reci.datatype)\n",
    "    result.set_raw_value(final_raw_value)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg Diff is 5.491479510641723e-06\n",
      "Avg Diff is 5.491479510641723e-06\n",
      "Max Diff is 2.3946166038513184e-05\n",
      "Max Diff is 2.3946166038513184e-05\n"
     ]
    }
   ],
   "source": [
    "# Let's verify reciprocal for the canonical format first.\n",
    "cano_test_data = np.arange(0.5, 0.99, 0.001)\n",
    "\n",
    "cano_expect_result = []\n",
    "cano_true_result = []\n",
    "for item in cano_test_data:\n",
    "  expect_value = 1.0 / item\n",
    "  test = FixedPoint(2, 29, np.int32)\n",
    "  test.set_raw_value(test.from_float(item))\n",
    "  result = test.reciprocal().to_float()\n",
    "  cano_expect_result.append(expect_value)\n",
    "  cano_true_result.append(result)\n",
    "    \n",
    "cano_expect_result = np.array(cano_expect_result)\n",
    "cano_true_result = np.array(cano_true_result)\n",
    "\n",
    "print(\"Avg Diff is %s\" % (diff(cano_expect_result, cano_true_result)))\n",
    "print(\"Avg Diff is %s\" % (diff(cano_expect_result, cano_true_result)))\n",
    "print(\"Max Diff is %s\" % (max_diff(cano_expect_result, cano_true_result)))\n",
    "print(\"Max Diff is %s\" % (max_diff(cano_expect_result, cano_true_result)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg Diff is 3.077560695339107e-06\n",
      "Max Diff is 8.51750411499097e-05\n"
     ]
    }
   ],
   "source": [
    "# Let's also verify reciprocal for a wider range (including negative values.)\n",
    "test_data = np.arange(-10.5, 10.5, 0.2)\n",
    "\n",
    "expect_result = []\n",
    "true_result = []\n",
    "for item in test_data:\n",
    "  expect_value = 1.0 / item\n",
    "  test = FixedPoint(4, 27, np.int32)\n",
    "  test.set_raw_value(test.from_float(item))\n",
    "  test = test.reciprocal()\n",
    "  result = test.to_float()\n",
    "  expect_result.append(expect_value)\n",
    "  true_result.append(result)\n",
    "\n",
    "expect_result = np.array(expect_result)\n",
    "true_result = np.array(true_result)\n",
    "\n",
    "print(\"Avg Diff is %s\" % (diff(expect_result, true_result)))\n",
    "print(\"Max Diff is %s\" % (max_diff(expect_result, true_result)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thanks to my college tianlin, I decided to expand this codelab to cover sin. (well, if you know how to handle sin, cos should be just trivial. :) )\n",
    "\n",
    "---\n",
    "### Sin\n",
    "\n",
    "For reference, please see details [here](https://en.wikipedia.org/wiki/Small-angle_approximation).\n",
    "\n",
    "$$\\sin \\theta = \\sum^{\\infty}_{n=0} \\frac{(-1)^n}{(2n+1)!} \\theta^{2n+1}$$\n",
    "\n",
    "And I found we can get good approximation with a few expansions (6~7), See the followig code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg Diff is 3.23777951746771e-05\n",
      "Max Diff is 0.0004422558761439342\n"
     ]
    }
   ],
   "source": [
    "def sin(radian, expansion=6):\n",
    "  sign = 1\n",
    "  divi = 1\n",
    "  result = radian\n",
    "  theta = radian\n",
    "  for i in range(1, expansion):\n",
    "    divi *= (2 * i) * (2 * i + 1)\n",
    "    theta *= radian * radian\n",
    "    sign *= -1\n",
    "    result += sign * theta / divi\n",
    "  return result\n",
    "\n",
    "test_data = np.arange(0, math.pi, 0.01)\n",
    "\n",
    "expect_result = np.sin(test_data)\n",
    "\n",
    "true_result = []\n",
    "\n",
    "for item in test_data:\n",
    "  true_result.append(sin(item))\n",
    "\n",
    "true_result = np.array(true_result)\n",
    "\n",
    "print(\"Avg Diff is %s\" % (diff(expect_result, true_result)))\n",
    "print(\"Max Diff is %s\" % (max_diff(expect_result, true_result)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg Diff is 0.0007714671227579742\n",
      "Max Diff is 0.014670956879854202\n"
     ]
    }
   ],
   "source": [
    "# Let's implement sin with fixed point math!\n",
    "\n",
    "# Let's just assume radian is in Q2.29.\n",
    "# Also, Let's just assume all the value are within [-pi/2, pi/2]\n",
    "def sin_fixed_point(radian):\n",
    "  multiplier_1 = FixedPoint(2, 29, np.int32)\n",
    "  multiplier_2 = FixedPoint(2, 29, np.int32)\n",
    "  multiplier_3 = FixedPoint(2, 29, np.int32)\n",
    "  multiplier_4 = FixedPoint(2, 29, np.int32)\n",
    "  multiplier_5 = FixedPoint(2, 29, np.int32)\n",
    "\n",
    "  multiplier_1.set_raw_value(-89478485)\n",
    "  multiplier_2.set_raw_value(4473924)\n",
    "  multiplier_3.set_raw_value(-106522)\n",
    "  multiplier_4.set_raw_value(1479)\n",
    "  multiplier_5.set_raw_value(-13)\n",
    "  \n",
    "  result = radian\n",
    "  radian_square = radian.mul(radian)\n",
    "  \n",
    "  # We're trying to avoid the overflow issupackste by afford some computation cost.\n",
    "  multiplier_1 = multiplier_1.mul(radian_square)\n",
    "  result = result.add(radian.mul(multiplier_1))\n",
    "  multiplier_2 = multiplier_2.mul(radian_square).mul(radian_square)\n",
    "  result = result.add(radian.mul(multiplier_2))\n",
    "  multiplier_3 = multiplier_3.mul(radian_square).mul(radian_square).mul(radian_square)\n",
    "  result = result.add(radian.mul(multiplier_3))\n",
    "  multiplier_4 = multiplier_4.mul(radian_square).mul(radian_square).mul(radian_square).mul(radian_square)\n",
    "  result = result.add(radian.mul(multiplier_4))\n",
    "  multiplier_5 = multiplier_4.mul(radian_square).mul(radian_square).mul(radian_square).mul(radian_square).mul(radian_square)\n",
    "  result = result.add(radian.mul(multiplier_5))\n",
    "  return result\n",
    "\n",
    "fixed_point_test_data = np.arange(-math.pi / 2.0, math.pi / 2.0, 0.01)\n",
    "\n",
    "fixed_point_expect_result = np.sin(fixed_point_test_data)\n",
    "\n",
    "fixed_point_true_result = []\n",
    "\n",
    "for item in fixed_point_test_data:\n",
    "  test = FixedPoint(2, 29, np.int32)\n",
    "  test.set_raw_value(test.from_float(item))\n",
    "  result = sin_fixed_point(test)\n",
    "  fixed_point_true_result.append(result.to_float())\n",
    "\n",
    "fixed_point_true_result = np.array(fixed_point_true_result)\n",
    "\n",
    "print(\"Avg Diff is %s\" % (diff(fixed_point_expect_result, fixed_point_true_result)))\n",
    "print(\"Max Diff is %s\" % (max_diff(fixed_point_expect_result, fixed_point_true_result)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thanks to my colleague Feng Liu, I'm adding Tanh as well.\n",
    "\n",
    "The implementation is also referencing benoitjacob's GEMMLOWP. :)\n",
    "\n",
    "---\n",
    "\n",
    "## Tanh\n",
    "\n",
    "---\n",
    "\n",
    "### Intro\n",
    "\n",
    "Tanh is properly the most sophisticated quantized kernel to implement, in practice, we found Tanh is more performant and much easier with the Lookup Table approach, but I think it's still worth mentioning the implementation of interger-only Tanh.\n",
    "\n",
    "Let's take a look at what is Tanh first:\n",
    "\n",
    "$$tanh(x) = \\frac{e^x - e^{-x}}{e^x + e^{-x}}$$\n",
    "\n",
    "It's very easier for us to observe that $$tanh(-x) = -tanh(x)$$\n",
    "\n",
    "So let's just focus on tanh(x) for x > 0.\n",
    "\n",
    "---\n",
    "\n",
    "### Transform\n",
    "\n",
    "Actually from $$tanh(x) = \\frac{e^x - e^{-x}}{e^x + e^{-x}}$$ we have two transforms:\n",
    "\n",
    "Option 1):\n",
    "\n",
    "$$tanh(x) = \\frac{e^{2x} - 1}{e^{2x} + 1}$$\n",
    "\n",
    "Option 2):\n",
    "\n",
    "$$tanh(x) = \\frac{1 - e^{-2x}}{1 + e^{-2x}}$$\n",
    "\n",
    "And we perfer the second option.\n",
    "\n",
    "Why?\n",
    "\n",
    "Well, it's easy to tell that for\n",
    "\n",
    "$$x > 0,$$\n",
    "$$e^{2x} \\in (1, +\\infty)$$\n",
    "$$e^{-2x} \\in (0, 1)$$\n",
    "\n",
    "So e^-2x is bounded and easier for our fixed point representation.\n",
    "\n",
    "So Let $$F(x) = e^{-2x},$$\n",
    "\n",
    "We will have $$tanh(x) = \\frac{1 - F(x)}{1 + F(x)}$$,\n",
    "\n",
    "In previous sections, we have already covered how to get division, so there's only one missing piece for us to handle.\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "### Exp\n",
    "\n",
    "For exp, we can actually use Taylor Series expansion to estimate the result for small x (< 0.25):\n",
    "\n",
    "$$ e^x = 1 + \\sum^{\\infty}_{n=1}\\frac{x^n}{n!}$$\n",
    "\n",
    "Let's see the following test example, note we're focusing on exp(-2x) for x > 0 case, so essentialy we care for exp(t) for t < 0 case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg Diff is 1.5972655070117802e-06\n",
      "Max Diff is 7.757713531231225e-06\n"
     ]
    }
   ],
   "source": [
    "def exp_estimate(x, expansion=4):\n",
    "  temp = 1\n",
    "  divi = 1\n",
    "  x_exp = 1\n",
    "  for i in range(1, expansion + 1):\n",
    "    divi *= i\n",
    "    x_exp *= x\n",
    "    temp += x_exp / float(divi)\n",
    "  return temp\n",
    "\n",
    "test_data = np.random.random(50) * -0.25\n",
    "\n",
    "expect_result = np.exp(test_data)\n",
    "\n",
    "true_result = []\n",
    "\n",
    "for item in test_data:\n",
    "  true_result.append(exp_estimate(item))\n",
    "\n",
    "true_result = np.array(true_result)\n",
    "\n",
    "print(\"Avg Diff is %s\" % (diff(expect_result, true_result)))\n",
    "print(\"Max Diff is %s\" % (max_diff(expect_result, true_result)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exp Continuous\n",
    "\n",
    "Since we can handle -0.25 < x < 0 case, how about other cases?\n",
    "\n",
    "Let\n",
    "\n",
    "$$x = 2^n * x_{cano}$$\n",
    "\n",
    "where $$x_{cano} \\in (-0.25, 0)$$\n",
    "\n",
    "So $$exp(x) = exp(2^n * x_{cano})$$\n",
    "\n",
    "=>\n",
    "\n",
    "$$exp(x) = \\prod_{k=1}^{2^n}exp(x_{cano})$$\n",
    "\n",
    "Great, so it's essentially computing the square of exp(x_cano) n times.\n",
    "\n",
    "We will be computing the exp(x_cano) in Q0.31 format, so for Qm.31-m case, \"n\" is essentially can be computed as\n",
    "\n",
    "(3 - couting_leading_zeros(x)) + m.\n",
    "\n",
    "Let's dive deeper into how to really implement this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg Diff is 1.2523731640690094e-06\n",
      "Max Diff is 6.627994739161203e-06\n"
     ]
    }
   ],
   "source": [
    "def exp_cano(x):\n",
    "  # Support exp for -0.25 < x < 0 case.\n",
    "  # Also assume data is in Q0.31 format.\n",
    "  result = FixedPoint(0, 31, np.int32)\n",
    "  result.set_raw_value((1 << 31) - 1)  ## This is actually 1, minus one to avoid overflow.\n",
    "\n",
    "  multiplier_2 = FixedPoint(0, 31, np.int32)\n",
    "  multiplier_3 = FixedPoint(0, 31, np.int32)\n",
    "  multiplier_4 = FixedPoint(0, 31, np.int32)\n",
    "  \n",
    "  multiplier_2.set_raw_value(1073741824) ## 1 / 2!\n",
    "  multiplier_3.set_raw_value(357913941)  ## 1 / 3!\n",
    "  multiplier_4.set_raw_value(89478485)  ## 1 / 4!\n",
    "\n",
    "  x_exp = x\n",
    "  result = result.add(x_exp)\n",
    "  x_exp = x_exp.mul(x)\n",
    "  result = result.add(x_exp.mul(multiplier_2))\n",
    "  x_exp = x_exp.mul(x)\n",
    "  result = result.add(x_exp.mul(multiplier_3))\n",
    "  x_exp = x_exp.mul(x)\n",
    "  result = result.add(x_exp.mul(multiplier_4))\n",
    "  return result\n",
    "\n",
    "# Let's test it.\n",
    "fixed_point_test_data = np.random.random(50) * -0.25\n",
    "\n",
    "fixed_point_expect_result = np.exp(fixed_point_test_data)\n",
    "\n",
    "fixed_point_true_result = []\n",
    "for item in fixed_point_test_data:\n",
    "  test = FixedPoint(0, 31, np.int32)\n",
    "  test.set_raw_value(test.from_float(item))\n",
    "  result = exp_cano(test)\n",
    "  fixed_point_true_result.append(result.to_float())\n",
    "    \n",
    "fixed_point_true_result = np.array(fixed_point_true_result)\n",
    "\n",
    "print(\"Avg Diff is %s\" % (diff(fixed_point_expect_result, fixed_point_true_result)))\n",
    "print(\"Max Diff is %s\" % (max_diff(fixed_point_expect_result, fixed_point_true_result)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg Diff is 2.534881819946447e-06\n",
      "Max Diff is 1.2134654507722775e-05\n",
      "Avg Diff is 2.533447583181539e-06\n",
      "Max Diff is 1.2132326201286237e-05\n"
     ]
    }
   ],
   "source": [
    "# Great, let's support for x < 0 case.\n",
    "def exp_negative_x(x):\n",
    "  leading_zeros = clz(x.neg().get_raw_value())\n",
    "  n = 3 - leading_zeros + x.integer\n",
    "  # We will represent x_cano in Q0.31 format.\n",
    "  x_cano = FixedPoint(0, 31, np.int32)\n",
    "  if n <= 0:\n",
    "    x_cano.set_raw_value(x.get_raw_value() * (2 ** x.integer))\n",
    "    n = 0\n",
    "  else:\n",
    "    # As you can see the following it's redudant: we shift right then shift left.\n",
    "    # It's just to make things clear.\n",
    "    # We first get the x_cano in Qm.31-m format.\n",
    "    # Then we rescale to Q0.31 format.\n",
    "    x_cano_raw = rounding_divide_by_POT(x.get_raw_value(), n)\n",
    "    x_cano.set_raw_value(x_cano_raw * (2 ** x.integer))\n",
    "  result = exp_cano(x_cano)\n",
    "  for i in range(n):\n",
    "    result = result.mul(result)\n",
    "  return result\n",
    "\n",
    "# Let's test it.\n",
    "fixed_point_test_data = np.random.random(50) * -5\n",
    "\n",
    "fixed_point_expect_result = np.exp(fixed_point_test_data)\n",
    "\n",
    "fixed_point_true_result_1 = []\n",
    "fixed_point_true_result_2 = []\n",
    "for item in fixed_point_test_data:\n",
    "  test_1 = FixedPoint(3, 28, np.int32)\n",
    "  test_2 = FixedPoint(4, 27, np.int32)\n",
    "  test_1.set_raw_value(test_1.from_float(item))\n",
    "  test_2.set_raw_value(test_2.from_float(item))\n",
    "  result_1 = exp_negative_x(test_1)\n",
    "  result_2 = exp_negative_x(test_2)\n",
    "  fixed_point_true_result_1.append(result_1.to_float())\n",
    "  fixed_point_true_result_2.append(result_2.to_float())\n",
    "    \n",
    "fixed_point_true_result_1 = np.array(fixed_point_true_result_1)\n",
    "fixed_point_true_result_2 = np.array(fixed_point_true_result_2)\n",
    "\n",
    "print(\"Avg Diff is %s\" % (diff(fixed_point_expect_result, fixed_point_true_result_1)))\n",
    "print(\"Max Diff is %s\" % (max_diff(fixed_point_expect_result, fixed_point_true_result_1)))\n",
    "print(\"Avg Diff is %s\" % (diff(fixed_point_expect_result, fixed_point_true_result_2)))\n",
    "print(\"Max Diff is %s\" % (max_diff(fixed_point_expect_result, fixed_point_true_result_2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg Diff is 8.401370436255383e-06\n",
      "Max Diff is 1.2965664681119371e-05\n"
     ]
    }
   ],
   "source": [
    "# Let's put it together for tanh.\n",
    "# Assume x use np.int32.\n",
    "def tanh(x):\n",
    "  # Handle x < 0 case.\n",
    "  if x.get_raw_value() < 0:\n",
    "    return tanh(x.neg()).neg()\n",
    "\n",
    "  minus_two_x = FixedPoint(x.integer + 1, x.fractional - 1, np.int32)\n",
    "  minus_two_x.set_raw_value(-x.get_raw_value())\n",
    "  t = exp_negative_x(minus_two_x)\n",
    "\n",
    "  # Note: t is in Q0.31 format, we need it in Q1.30.\n",
    "  new_t = FixedPoint(1, 30, np.int32)\n",
    "  new_t.set_raw_value(rounding_divide_by_POT(t.get_raw_value(), 1))\n",
    "  \n",
    "  numerator = FixedPoint(1, 30, np.int32)\n",
    "  numerator.set_raw_value(numerator.from_float(1.0))\n",
    "  denominator = FixedPoint(1, 30, np.int32)\n",
    "  denominator.set_raw_value(denominator.from_float(1.0))\n",
    "\n",
    "  numerator = numerator.sub(new_t)\n",
    "  denominator = denominator.add(new_t)\n",
    "  \n",
    "  denominator = denominator.reciprocal()\n",
    "  return numerator.mul(denominator)\n",
    "\n",
    "# Let's test it.\n",
    "fixed_point_test_data = (np.random.random(50) - 0.5) * 8\n",
    "\n",
    "fixed_point_expect_result = np.tanh(fixed_point_test_data)\n",
    "\n",
    "fixed_point_true_result = []\n",
    "for item in fixed_point_test_data:\n",
    "  test = FixedPoint(3, 28, np.int32)\n",
    "  test.set_raw_value(test.from_float(item))\n",
    "  result = tanh(test)\n",
    "  fixed_point_true_result.append(result.to_float())\n",
    "    \n",
    "fixed_point_true_result = np.array(fixed_point_true_result)\n",
    "\n",
    "print(\"Avg Diff is %s\" % (diff(fixed_point_expect_result, fixed_point_true_result)))\n",
    "print(\"Max Diff is %s\" % (max_diff(fixed_point_expect_result, fixed_point_true_result)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}