Approximationa rough calculation without math.h （for my homework）

aprrox.c

#include <stdio.h>
#define LOG2   0.693147180559945309417  //log(2)
#define chk_nan(x) x!=x

double log_fast(double x)
{
  static double table[17]={
    .0                       ,  // log( 16 /16)
    .0606246218164348425806  ,  // log( 17 /16)
    .1177830356563834545387  ,  // log( 18 /16)
    .17185025692665922234    ,  // log( 19 /16)
    .2231435513142097557662  ,  // log( 20 /16)
    .2719337154836417588316  ,  // log( 21 /16)
    .3184537311185346158102  ,  // log( 22 /16)
    .3629054936893684531378  ,  // log( 23 /16)
    .405465108108164381978   ,  // log( 24 /16)
    .4462871026284195115325  ,  // log( 25 /16)
    .4855078157817008078017  ,  // log( 26 /16)
    .5232481437645478365168  ,  // log( 27 /16)
    .5596157879354226862708  ,  // log( 28 /16)
    .5947071077466927895143  ,  // log( 29 /16)
    .6286086594223741377443  ,  // log( 30 /16)
    .6613984822453650082602  ,  // log( 31 /16)
    .6931471805599453094172  ,  // log( 32 /16)
  };   
  unsigned long long w,mantissa16;
  int q;
  double y,h,z;
  w=*(unsigned long long*)&x;
  q=(((int)(w>>47)&0x1F)+1)>>1;
  mantissa16=(w & 0xFFFFFFFFFFFFFULL)^0x4030000000000000ULL;//仮数*16  16<=mantissa16<32 
  h=*(double*)&mantissa16;
  z=(double)(q+16);
  h=(h-z)/(h+z);
  z=h*h;
  y=(2.0/9)*z+2.0/7;
  y=y*z+2.0/5;
  y=y*z+2.0/3;
  y=y*z+2.0;
  y=y*h;
  return ((int)(w>>52)-1023)*LOG2+table[q]+y;
}


double exp_fast(double x)
{
  static double devisor[7]={1.0/(2*3*4*5*6),1.0/(2*3*4*5),1.0/(2*3*4),1.0/(2*3),1.0/2,1.0,1.0};
  static unsigned long long table[16]={
    0x059B0D3158574ull,  // 2^( 1 /32)-1
    0x11301D0125B51ull,  // 2^( 3 /32)-1
    0x1D4873168B9AAull,  // 2^( 5 /32)-1
    0x29E9DF51FDEE1ull,  // 2^( 7 /32)-1
    0x371A7373AA9CBull,  // 2^( 9 /32)-1
    0x44E086061892Dull,  // 2^( 11 /32)-1
    0x5342B569D4F82ull,  // 2^( 13 /32)-1
    0x6247EB03A5585ull,  // 2^( 15 /32)-1
    0x71F75E8EC5F74ull,  // 2^( 17 /32)-1
    0x82589994CCE13ull,  // 2^( 19 /32)-1
    0x93737B0CDC5E5ull,  // 2^( 21 /32)-1
    0xA5503B23E255Dull,  // 2^( 23 /32)-1
    0xB7F76F2FB5E47ull,  // 2^( 25 /32)-1
    0xCB720DCEF9069ull,  // 2^( 27 /32)-1
    0xDFC97337B9B5Full,  // 2^( 29 /32)-1
    0xF50765B6E4540ull,  // 2^( 31 /32)-1
  };
  double y=1.0/(2*3*4*5*6*7),*p=devisor,z,r;
  int q;
  unsigned long long w;
  z=x*(16.0/LOG2);
  q=(int)z-(x<0);
  r=x-((q<<1)+1)*(LOG2/32.0);
  w=(unsigned long long)(1023+(q>>4))<<52 ^ table[q & 0xF];
  z=*(double*)&w;
  do{
    y=y*r+(*p);
    p++;
  }while (p<devisor+7);
  return y*z;
}

double powf_fast(double x, double y){
  return exp_fast(log_fast(x)*y);
}


long double calc_napier(long int count){
  long double ans = 1;
  long long int m = 1;
  for (long int i = 1; i < count; i++)
  {
    long long int tmp =  (long long int)(m *= i);
    if(tmp <= 0){
      printf("Overflow! Stop at %ld times\n",i);
      return ans;
    }
    ans += (long double)1.0 / tmp;
  }
  return ans;
}

long double calc_pi_wallisformula(long int count){
  long double pi = 1;
  for(long int i = 1; i <= count; i++){
    long int n = 4 * i * i;
    pi *= (long double)n / (n-1);
  }
  return pi * 2;
}

long double calc_pi_arctan(long int count){
  long double ans = 0;
  long double tmp = 0;
  const long double ROOT3 = 1.73205080757; 
  for (long int n = 1; n < count; n++)
  {
    if (n % 2 == 1){
      tmp += powf_fast(ROOT3 / 3, (2 * n - 1)) * (1.0 / (2 * n - 1));
    }
    else{
      tmp -= powf_fast(ROOT3 / 3, (2 * n - 1)) * (1.0 / (2 * n - 1));
    }

    if(chk_nan(tmp)) 
    {
      printf("Overflow! Stop at %ld times\n",n);
      return ans * 6;
    }
    ans = tmp;
  }
  return ans * 6;
}

long double calc_pi_gauss_legendre(long int count){
  // b0 = 1/sqrt(2)
  long double a=1,b=0.70710678118654752440084436210485,t=0.25,p=1;
  for (long int i = 0; i < count; i++)
  {
    long double tmpa = a,tmpb = b,tmpt = t, tmpp = p;
    a = (a + b) * 0.5;
    b = powf_fast(tmpa*b,0.5);
    t -= p * (a - tmpa) * (a - tmpa);
    p *= 2;
    if(chk_nan(a) || chk_nan(b) || chk_nan(t) || chk_nan(p)){
      printf("Overflow! Stop at %ld times\n",i);
      return (tmpa + tmpb) * (tmpa + tmpb) / (4 * tmpt);
    }
  }
  return (a + b) * (a + b) / (4 * t);
}

int main(int argc, char *argv[]){

  int select_num = 0;
  long int count = 0;
  long double ans;

  printf("Select Number \n"
      "1: napier's number (taylor) \n"
      "2: pi (arcTan taylor) \n"
      "3: pi (wallis' formula) \n"
      "4: pi (gauss–legendre algorithm) \n"
      " > ");
  scanf("%d",&select_num);
  printf("Number of trials > ");
  scanf("%ld",&count);

  if (count <= 0)
  {
    printf("%ld is invalid number of trials.\n",count);
    return 0;
  }

  switch (select_num)
  {
    case 1:
      ans = calc_napier(count);
      break;
    case 2:
      ans = calc_pi_arctan(count);
      break;
    case 3:
      ans = calc_pi_wallisformula(count);
      break;
    case 4:
      ans = calc_pi_gauss_legendre(count);
      break;
    default:
      printf("%d's method is not found.\n",select_num);
      return 0;
      break;
  }

  printf("RESULT: %14.14Lf\n", ans);
  return 0;
}