eigenmath/dpow.cpp

47 lines
760 B
C++

// power function for double precision floating point
#include "stdafx.h"
#include "defs.h"
void
dpow(void)
{
double a, b, base, expo, result, theta;
expo = pop_double();
base = pop_double();
// divide by zero?
if (base == 0.0 && expo < 0.0)
stop("divide by zero");
// nonnegative base or integer power?
if (base >= 0.0 || fmod(expo, 1.0) == 0.0) {
result = pow(base, expo);
push_double(result);
return;
}
result = pow(fabs(base), expo);
theta = M_PI * expo;
// this ensures the real part is 0.0 instead of a tiny fraction
if (fmod(expo, 0.5) == 0.0) {
a = 0.0;
b = sin(theta);
} else {
a = cos(theta);
b = sin(theta);
}
push_double(a * result);
push_double(b * result);
push(imaginaryunit);
multiply();
add();
}