#include <math.h>

int __fpclassifyd(double x) {
  union fp {
    double d;
    long i;
  } f;

  f.d = x;

  if((f.i & 0x7fffffffffffffff) == 0) {
    return FP_ZERO;
  }

  if((f.i & 0x7ff0000000000000) == 0 && (f.i & 0xfffffffffffff != 0)) {
    return FP_SUBNORMAL;
  }

  if((f.i & 0x7ff0000000000000) == 0x7ff0000000000000) {
    if((f.i & 0xfffffffffffff) == 0) {
      return FP_INFINITE;
    }
    else {
      return FP_NAN;
    }
  }

  return FP_NORMAL;
}

int __fpclassifyf(float x) {
  union fp {
    float f;
    int i;
  } f;

  f.f = x;

  if((f.i & 0x7fffffff) == 0) {
    return FP_ZERO;
  }

  if((f.i & 0x7f800000) == 0 && (f.i & 0x7fffff != 0)) {
    return FP_SUBNORMAL;
  }

  if((f.i & 0x7f800000) == 0x7f800000) {
    if((f.i & 0x7fffff) == 0) {
      return FP_INFINITE;
    }
    else {
      return FP_NAN;
    }
  }

  return FP_NORMAL;
}

/*
  This is based off IEEE 754 binary128 (A.K.A quad-precision) FP numbers. 
  Apparently, C tends to actually use 80-bit extended precision, which would
  mean this doesn't actually work. We'll worry about that later, though.
  -Kat
*/
int __fpclassifyl(long double x) {
  union fp {
    long double d;
    long long i;
  } f;

  f.d = x;

  if((f.i & 0x7fffffffffffffffffffffffffffffff) == 0) {
    return FP_ZERO;
  }

  if((f.i & 0x7fff0000000000000000000000000000) == 0 && (f.i & 0xffffffffffffffffffffffffffff != 0)) {
    return FP_SUBNORMAL;
  }

  if((f.i & 0x7fff0000000000000000000000000000) == 0x7fff0000000000000000000000000000) {
    if((f.i & 0xffffffffffffffffffffffffffff) == 0) {
      return FP_INFINITE;
    }
    else {
      return FP_NAN;
    }
  }

  return FP_NORMAL;
}