function [bp,bn]=bern(x)xlim=1e-2; ax=abs(x);if (ax == 0)   bp=1.;    bn=1.; returnendif (ax > 80),   if (x > 0)     bp=0.;      bn=x;       return   else           bp=-x;      bn=0.;      return     endendif (ax > xlim)   bp=x/(exp(x)-1);    bn=x+bp;    returnelse   ii=1;    fp=1.;   fn=1.;    df=1.;    s=1.;   while (abs(df) > eps),     ii=ii+1;       s=-s;      df=df*x/ii;     fp=fp+df;      fn=fn+s*df;     bp=1./fp;      bn=1./fn;   end   returnendreturn