#include <stdio.h>
#include <math.h>

/* Four ways of computing (e^z - 1)/z */

double expm1_a(double z)
{
  return (exp(z)-1)/z;
  
}

double expm1_b(double z)
{
  double w;
  w = exp(z);
  if ( w == 0 ) return -1/z;
  if ( w == 1 ) return 1.0;
  return (w-1)/log(w);
  
}

double expm1_c(double z)
{
  double u, w;
  w = exp(z);  u = 1 + z;
  if ( w == 1 || u == 1 ) return 1;
  return (w-1)/(u-1);
  
}

#define N 20
double expm1_d(double z)
{
  double coeffs[N], sum;
  int    i;

  if ( fabs(z) >= 1 )
    return (exp(z)-1)/z;

  coeffs[0] = 1;
  for ( i = 2; i <= N; i++ )
    coeffs[i-1] = coeffs[i-2] / (double)i;

  /* use Horner's rule to evaluate power series */
  sum = coeffs[N-1];
  for ( i = N-2; i >= 0; i-- )
    sum = coeffs[i] + z*sum;

  return sum;
}

int main(int argc, char *argv[])
{
  double z;
  int i;

  for ( i = -2; i <= 16; i++ )
    {
      z = 3*pow(10.0,-(double)i);
      printf("%11.6g  %11.6g  %11.6g  %11.6g\n",
             z, expm1_a(z)-expm1_d(z), 
             expm1_b(z)-expm1_d(z), expm1_c(z)-expm1_d(z));
      z = pow(10.0,-(double)i);
      printf("%11.6g  %11.6g  %11.6g  %11.6g\n",
             z, expm1_a(z)-expm1_d(z), 
             expm1_b(z)-expm1_d(z), expm1_c(z)-expm1_d(z));
    }
}