fast_exp.nim

# Laser
# Copyright (c) 2018 Mamy André-Ratsimbazafy
# Distributed under the Apache v2 License (license terms are at http://www.apache.org/licenses/LICENSE-2.0).
# This file may not be copied, modified, or distributed except according to those terms.

import
  math

const
  bits = 10
  PowBits = 1 shl bits
  ln2 = ln(2'f32)
  ExpMax = 88'i32
  ExpMin = -88'i32
  ExpA = PowBits.float32 / ln2
  ExpB = ln2 / PowBits.float32


func initExpLUT(): array[PowBits, int32] =
  for i, val in result.mpairs:
    let y = pow(2'f32, i.float32 / PowBits.float32)
    val = cast[int32](y) and ((1 shl 23) - 1)

let ExpLUT = initExpLUT()

proc fast_exp(x: float32): float32 =

  static: assert cpuEndian == littleEndian, "Fast exp has not been tested on big-endian architecture"

  let r = int32(x * ExpA) # ->f32->i32->f32 Rounding is important
  let t = x - r.float32 * ExpB + 1

  let v = r and ((1 shl 10) - 1)
  var u = r + (127 shl 10)
  u = u shr bits
  u = u shl 23

  let ti = ExpLut[v] or u
  result = t * cast[float32](ti)

when isMainModule:

  echo exp(0.5'f32)
  echo fast_exp(0.5'f32)

  echo exp(-0.5'f32)
  echo fast_exp(-0.5'f32)