/* exponentiate        GCW  12/01/95

   power(x,y) defined if
       #  y is a positive integer
       #  y is a negative integer and x is a nonzero real
       #  y is a real and x is a positive real
   result has same type as x.
*/

#ifndef _BOBCONSTANTS
#include bob:Constants
#endif

#ifndef _math_power
#define _math_power 1
#endif

power(x,y)
{
 switch (typeof(y))
 {
  case REAL:
    if (x < 0.0) quit("Positive real argument needed");
    else return (x > 0.0)?(exp(y*log(x))):0.0;
    break;
  case INTEGER:
    switch (typeof(x))
    {
     case REAL:
        return (y < 0)?fastpower(1.0/x,-y):fastpower(x,y);
        break;
     case INTEGER:
        if (y < 0) quit("Positive exponent needed");
        else return fastpower(x,y);
        break;
     }
     break;
 }
}

fastpower(x,y)
{
 if (y == 0) return (typeof(x) == REAL)?1.0:1;
 if (y == 1) return x;
 if (y%2) return x*fastpower(x*x,y/2);
 else return fastpower(x*x, y/2);
}
   