dpow.cpp

   1 // power function for double precision floating point
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 void
   7 dpow(void)
   8 {
   9         double a, b, base, expo, result, theta;
  10
  11         expo = pop_double();
  12         base = pop_double();
  13
  14         // divide by zero?
  15
  16         if (base == 0.0 && expo < 0.0)
  17                 stop("divide by zero");
  18
  19         // nonnegative base or integer power?
  20
  21         if (base >= 0.0 || fmod(expo, 1.0) == 0.0) {
  22                 result = pow(base, expo);
  23                 push_double(result);
  24                 return;
  25         }
  26
  27         result = pow(fabs(base), expo);
  28
  29         theta = M_PI * expo;
  30
  31         // this ensures the real part is 0.0 instead of a tiny fraction
  32
  33         if (fmod(expo, 0.5) == 0.0) {
  34                 a = 0.0;
  35                 b = sin(theta);
  36         } else {
  37                 a = cos(theta);
  38                 b = sin(theta);
  39         }
  40
  41         push_double(a * result);
  42         push_double(b * result);
  43         push(imaginaryunit);
  44         multiply();
  45         add();
  46 }