1 // power function for double precision floating point
9 double a
, b
, base
, expo
, result
, theta
;
16 if (base
== 0.0 && expo
< 0.0)
17 stop("divide by zero");
19 // nonnegative base or integer power?
21 if (base
>= 0.0 || fmod(expo
, 1.0) == 0.0) {
22 result
= pow(base
, expo
);
27 result
= pow(fabs(base
), expo
);
31 // this ensures the real part is 0.0 instead of a tiny fraction
33 if (fmod(expo
, 0.5) == 0.0) {
41 push_double(a
* result
);
42 push_double(b
* result
);