double B(double x){             /* Returns B(x) */
double d=8.333333333333333e-2;
double tmp;
double f;
double y;

if (fabs(x)<1.0e-6)
 y=1.0+x*((x*d)-0.5);
 else {
  tmp=exp(x);
  f=1.0/(tmp-1.0);
  y=x*f;
 }
return y;
}

double dB(double x) {		/* Returns dB/dx at its parameter */	
double d=8.333333333333333e-2;
double xmax=log(sqrt(DBL_MAX));	  /* DBL_MAX is the */
      /* largest possible value for a double type variable in */
      /* your system.  If the compiler does not understand this, */
      /* manually enter that value or just make it _really_ large */
double y;

    if (fabs(x)<1.0e-6)
      y=2*x*d-0.5;
    else if (x > xmax)
      y=0;
    else
      y=(exp(x)-1-x*exp(x))/((exp(x)-1)*(exp(x)-1));
    return y;
}