e_pow.c

/***************************************************************************/
/* An ultimate power routine. Given two IEEE double machine numbers y,x    */
/* it computes the correctly rounded (to nearest) value of X^y.            */
/***************************************************************************/
double
SECTION
__ieee754_pow(double x, double y) {
  double z,a,aa,error, t,a1,a2,y1,y2;
#if 0
  double gor=1.0;
#endif
  mynumber u,v;
  int k;
  int4 qx,qy;
  v.x=y;
  u.x=x;
  if (v.i[LOW_HALF] == 0) { /* of y */
    qx = u.i[HIGH_HALF]&0x7fffffff;
    /* Checking  if x is not too small to compute */
    if (((qx==0x7ff00000)&&(u.i[LOW_HALF]!=0))||(qx>0x7ff00000)) return NaNQ.x;
    if (y == 1.0) return x;
    if (y == 2.0) return x*x;
    if (y == -1.0) return 1.0/x;
    if (y == 0) return 1.0;
  }
  /* else */
  if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)||        /* x>0 and not x->0 */
       (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0))  &&
				      /*   2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
      (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) {              /* if y<-1 or y>1   */
    double retval;

    SET_RESTORE_ROUND (FE_TONEAREST);

    z = log1(x,&aa,&error);                                 /* x^y  =e^(y log (X)) */
    t = y*134217729.0;
    y1 = t - (t-y);
    y2 = y - y1;
    t = z*134217729.0;
    a1 = t - (t-z);
    a2 = (z - a1)+aa;
    a = y1*a1;
    aa = y2*a1 + y*a2;
    a1 = a+aa;
    a2 = (a-a1)+aa;
    error = error*ABS(y);
    t = __exp1(a1,a2,1.9e16*error);     /* return -10 or 0 if wasn't computed exactly */
    retval = (t>0)?t:power1(x,y);

    return retval;
  }

  if (x == 0) {
    if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
	|| (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
      return y;
    if (ABS(y) > 1.0e20) return (y>0)?0:1.0/0.0;
    k = checkint(y);
    if (k == -1)
      return y < 0 ? 1.0/x : x;
    else
      return y < 0 ? 1.0/0.0 : 0.0;                               /* return 0 */
  }

  qx = u.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
  qy = v.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */

  if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0)) return NaNQ.x;
  if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))
    return x == 1.0 ? 1.0 : NaNQ.x;

  /* if x<0 */
  if (u.i[HIGH_HALF] < 0) {
    k = checkint(y);
    if (k==0) {
      if (qy == 0x7ff00000) {
	if (x == -1.0) return 1.0;
	else if (x > -1.0) return v.i[HIGH_HALF] < 0 ? INF.x : 0.0;
	else return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x;
      }
      else if (qx == 0x7ff00000)
	return y < 0 ? 0.0 : INF.x;
      return NaNQ.x;                              /* y not integer and x<0 */
    }
    else if (qx == 0x7ff00000)
      {
	if (k < 0)
	  return y < 0 ? nZERO.x : nINF.x;
	else
	  return y < 0 ? 0.0 : INF.x;
      }
    return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */
  }
  /* x>0 */

  if (qx == 0x7ff00000)                              /* x= 2^-0x3ff */
    {if (y == 0) return NaNQ.x;
    return (y>0)?x:0; }

  if (qy > 0x45f00000 && qy < 0x7ff00000) {
    if (x == 1.0) return 1.0;
    if (y>0) return (x>1.0)?huge*huge:tiny*tiny;
    if (y<0) return (x<1.0)?huge*huge:tiny*tiny;
  }

  if (x == 1.0) return 1.0;
  if (y>0) return (x>1.0)?INF.x:0;
  if (y<0) return (x<1.0)?INF.x:0;
  return 0;     /* unreachable, to make the compiler happy */
}
#ifndef __ieee754_pow
strong_alias (__ieee754_pow, __pow_finite)
#endif

/**************************************************************************/
/* Computing x^y using more accurate but more slow log routine            */
/**************************************************************************/
static double
SECTION
power1(double x, double y) {
  double z,a,aa,error, t,a1,a2,y1,y2;
  z = my_log2(x,&aa,&error);
  t = y*134217729.0;
  y1 = t - (t-y);
  y2 = y - y1;
  t = z*134217729.0;
  a1 = t - (t-z);
  a2 = z - a1;
  a = y*z;
  aa = ((y1*a1-a)+y1*a2+y2*a1)+y2*a2+aa*y;
  a1 = a+aa;
  a2 = (a-a1)+aa;
  error = error*ABS(y);
  t = __exp1(a1,a2,1.9e16*error);
  return (t >= 0)?t:__slowpow(x,y,z);
}