private static double expm1(double x,double hiPrecOut[]){
  if (Double.isNaN(x) || x == 0.0) {
    return x;
  }
  if (x <= -1.0 || x >= 1.0) {
    double hiPrec[]=new double[2];
    exp(x,0.0,hiPrec);
    if (x > 0.0) {
      return -1.0 + hiPrec[0] + hiPrec[1];
    }
 else {
      final double ra=-1.0 + hiPrec[0];
      double rb=-(ra + 1.0 - hiPrec[0]);
      rb+=hiPrec[1];
      return ra + rb;
    }
  }
  double baseA;
  double baseB;
  double epsilon;
  boolean negative=false;
  if (x < 0.0) {
    x=-x;
    negative=true;
  }
{
    int intFrac=(int)(x * 1024.0);
    double tempA=ExpFracTable.EXP_FRAC_TABLE_A[intFrac] - 1.0;
    double tempB=ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
    double temp=tempA + tempB;
    tempB=-(temp - tempA - tempB);
    tempA=temp;
    temp=tempA * HEX_40000000;
    baseA=tempA + temp - temp;
    baseB=tempB + (tempA - baseA);
    epsilon=x - intFrac / 1024.0;
  }
  double zb=0.008336750013465571;
  zb=zb * epsilon + 0.041666663879186654;
  zb=zb * epsilon + 0.16666666666745392;
  zb=zb * epsilon + 0.49999999999999994;
  zb*=epsilon;
  zb*=epsilon;
  double za=epsilon;
  double temp=za + zb;
  zb=-(temp - za - zb);
  za=temp;
  temp=za * HEX_40000000;
  temp=za + temp - temp;
  zb+=za - temp;
  za=temp;
  double ya=za * baseA;
  temp=ya + za * baseB;
  double yb=-(temp - ya - za * baseB);
  ya=temp;
  temp=ya + zb * baseA;
  yb+=-(temp - ya - zb * baseA);
  ya=temp;
  temp=ya + zb * baseB;
  yb+=-(temp - ya - zb * baseB);
  ya=temp;
  temp=ya + baseA;
  yb+=-(temp - baseA - ya);
  ya=temp;
  temp=ya + za;
  yb+=-(temp - ya - za);
  ya=temp;
  temp=ya + baseB;
  yb+=-(temp - ya - baseB);
  ya=temp;
  temp=ya + zb;
  yb+=-(temp - ya - zb);
  ya=temp;
  if (negative) {
    double denom=1.0 + ya;
    double denomr=1.0 / denom;
    double denomb=-(denom - 1.0 - ya) + yb;
    double ratio=ya * denomr;
    temp=ratio * HEX_40000000;
    final double ra=ratio + temp - temp;
    double rb=ratio - ra;
    temp=denom * HEX_40000000;
    za=denom + temp - temp;
    zb=denom - za;
    rb+=(ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;
    rb+=yb * denomr;
    rb+=-ya * denomb * denomr* denomr;
    ya=-ra;
    yb=-rb;
  }
  if (hiPrecOut != null) {
    hiPrecOut[0]=ya;
    hiPrecOut[1]=yb;
  }
  return ya + yb;
}

private static double expm1 ( double x , double hiPrecOut [ ] ) { if ( Double . isNaN ( x ) || x == 0.0 ) { return x ; } if ( x <= - 1.0 || x >= 1.0 ) { double hiPrec [ ] = new double [ 2 ] ; exp ( x , 0.0 , hiPrec ) ; if ( x > 0.0 ) { return - 1.0 + hiPrec [ 0 ] + hiPrec [ 1 ] ; } else { final double ra = - 1.0 + hiPrec [ 0 ] ; double rb = - ( ra + 1.0 - hiPrec [ 0 ] ) ; rb += hiPrec [ 1 ] ; return ra + rb ; } } double baseA ; double baseB ; double epsilon ; boolean negative = false ; if ( x < 0.0 ) { x = - x ; negative = true ; } { int intFrac = ( int ) ( x * 1024.0 ) ; double tempA = ExpFracTable . EXP_FRAC_TABLE_A [ intFrac ] - 1.0 ; double tempB = ExpFracTable . EXP_FRAC_TABLE_B [ intFrac ] ; double temp = tempA + tempB ; tempB = - ( temp - tempA - tempB ) ; tempA = temp ; temp = tempA * HEX_40000000 ; baseA = tempA + temp - temp ; baseB = tempB + ( tempA - baseA ) ; epsilon = x - intFrac / 1024.0 ; } double zb = 0.008336750013465571 ; zb = zb * epsilon + 0.041666663879186654 ; zb = zb * epsilon + 0.16666666666745392 ; zb = zb * epsilon + 0.49999999999999994 ; zb *= epsilon ; zb *= epsilon ; double za = epsilon ; double temp = za + zb ; zb = - ( temp - za - zb ) ; za = temp ; temp = za * HEX_40000000 ; temp = za + temp - temp ; zb += za - temp ; za = temp ; double ya = za * baseA ; temp = ya + za * baseB ; double yb = - ( temp - ya - za * baseB ) ; ya = temp ; temp = ya + zb * baseA ; yb += - ( temp - ya - zb * baseA ) ; ya = temp ; temp = ya + zb * baseB ; yb += - ( temp - ya - zb * baseB ) ; ya = temp ; temp = ya + baseA ; yb += - ( temp - baseA - ya ) ; ya = temp ; temp = ya + za ; yb += - ( temp - ya - za ) ; ya = temp ; temp = ya + baseB ; yb += - ( temp - ya - baseB ) ; ya = temp ; temp = ya + zb ; yb += - ( temp - ya - zb ) ; ya = temp ; if ( negative ) { double denom = 1.0 + ya ; double denomr = 1.0 / denom ; double denomb = - ( denom - 1.0 - ya ) + yb ; double ratio = ya * denomr ; temp = ratio * HEX_40000000 ; final double ra = ratio + temp - temp ; double rb = ratio - ra ; temp = denom * HEX_40000000 ; za = denom + temp - temp ; zb = denom - za ; rb += ( ya - za * ra - za * rb - zb * ra - zb * rb ) * denomr ; rb += yb * denomr ; rb += - ya * denomb * denomr * denomr ; ya = - ra ; yb = - rb ; } if ( hiPrecOut != null ) { hiPrecOut [ 0 ] = ya ; hiPrecOut [ 1 ] = yb ; } return ya + yb ; }
internal helper method for expm1
