basic math function approximations

math.c

#include <stdio.h>
#include <stdlib.h>

#define sgn(v) ((0 < (v)) - ((v) < 0))
#define abs(a) (((a) < 0) ? -(a) : (a))

static int
ilog2(int n) {
#ifdef _MSC_VER
  unsigned long msbp = 0;
  _BitScanReverse(&msbp, (unsigned long)n);
  return (int)msbp;
#elif defined(__GNUC__) || defined(__clang__)
  return (int)sizeof(unsigned long) * 8 - 1 - __builtin_clzl((unsigned long)n);
#else
  static const char tbl[256] = {
    #define lt(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
    0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,lt(4), lt(5), lt(5), lt(6), lt(6), lt(6), lt(6),
    lt(7), lt(7), lt(7), lt(7), lt(7), lt(7), lt(7), lt(7)};
  int tt, t;
  if ((tt = (n >> 16))) return (t = (tt >> 8)) ? 24 + tbl[t] : 16 + tbl[tt];
  else return (t = (n >> 8)) ? 8 + tbl[t] : tbl[n];
  #undef lt
#endif
}
static int
ipow10(int n) {
  static const int pow10[] = {
    1, 10, 100, 1000, 10000, 100000, 1000000,
    10000000, 100000000, 1000000000
  };
  return pow10[n]; 
}
static float
ipow10r(int n) {
  static const float rpow10[] = {
    1.0f, 0.1f, 0.01f, 0.001f, 0.0001f, 0.00001f, 0.000001f,
    0.0000001f, 0.00000001f, 0.000000001f, 0.0000000001f,
    0.0000000001f, 0.00000000001f, 0.000000000001f,
    0.000000000001f, 0.00000000000001f
  };
  return rpow10[n]; 
}
static int
ilog10(int v) {
  int t = (ilog2(v) + 1) * 1233 >> 12;
  return t - (v < ipow10(t));
}
static float 
log2a(float x) {
  union { float f; unsigned i; } v = { x };
  union { unsigned i; float f; } m = { (v.i & 0x007FFFFF) | 0x3f000000 };
  float y = v.i;
  y *= 1.1920928955078125e-7f;
  return y - 124.22551499f - 1.498030302f * m.f - 1.72587999f / (0.3520887068f + m.f);
}
static float
loga(float x) {
  return 0.69314718f * log2a(x);
}
static float
pow2a(float p) {
  float o = (p < 0) ? 1.0f : 0.0f;
  float c = (p < -126) ? -126.0f : p;
  int w = c;
  float z = c - w + o;
  float t = (c + 121.2740575f + 27.7280233f / (4.84252568f - z) - 1.49012907f * z);
  union { unsigned i; float f; } v = {(unsigned)((1 << 23) *  t) };
  return v.f;
}
static float
expa(float p) {
  return pow2a(1.442695040f * p);
}
static float
powa(float x, float p) {
  return pow2a(p * log2a(x));
}
static float
sinapi(float x) {
  union { float f; unsigned i; } p = { 0.20363937680730309f };
  union { float f; unsigned i; } r = { 0.015124940802184233f };
  union { float f; unsigned i; } s = { -0.0032225901625579573f };
  union { float f; unsigned i; } vx = { x };
  unsigned sign = vx.i & 0x80000000;
  vx.i = vx.i & 0x7FFFFFFF;

  float qp = 1.2732395447351627f * x - 0.40528473456935109f * x * vx.f;
  float qpsq = qp * qp;

  p.i |= sign;
  r.i |= sign;
  s.i ^= sign;
  return 0.78444488374548933f * qp + qpsq * (p.f + qpsq * (r.f + qpsq * s.f));
}
static float
sina(float x) {
  int k = x * 0.15915494309189534f;
  float half = (x < 0) ? -0.5f : 0.5f;
  return sinapi((half + k) * 6.2831853071795865f - x);
}
static float
cosapi(float x) {
  return sinapi(x + (x > 1.5707963267948966f) ? 
                -4.7123889803846899f : 1.5707963267948966f);
}
static inline float
cosa(float x) {
  return sina(x + 1.5707963267948966f);
}
static float
rsqrt(float n){
  union {float f; unsigned i;} conv = {n};
  conv.i = 0x5f375A84 - (conv.i >> 1);
  conv.f = conv.f * (1.5f - ((n * 0.5f) * conv.f * conv.f));
  return conv.f;
}
static float
sqrta(float x) {
  return x * rsqrt(x);
}
int main() {
  // only works for integer representable float values
  float v = 25.65f;
  if (v <= -1.0f || v >= 1.0f) {
    int e = ilog10(abs(v));
    float m = v * ipow10r(e);
    printf("scientific notation: %fe%d", m, e);  
  } else {
    float s = sgn(v), n = abs(v);
    int e = ilog10(1.0f/n) + 1;
    float m = s * n * ipow10(e);
    printf("scientific notation: %fe%d", m, -e);  
  }
  return 0;
}