libf2c/libF77/pow_zi.c

   1 #include "f2c.h"
   2
   3 #ifdef KR_headers
   4 VOID pow_zi(p, a, b)    /* p = a**b  */
   5  doublecomplex *p, *a; integer *b;
   6 #else
   7 extern void z_div(doublecomplex*, doublecomplex*, doublecomplex*);
   8 void pow_zi(doublecomplex *p, doublecomplex *a, integer *b)     /* p = a**b  */
   9 #endif
  10 {
  11         integer n;
  12         unsigned long u;
  13         double t;
  14         doublecomplex q, x;
  15         static doublecomplex one = {1.0, 0.0};
  16
  17         n = *b;
  18         q.r = 1;
  19         q.i = 0;
  20
  21         if(n == 0)
  22                 goto done;
  23         if(n < 0)
  24                 {
  25                 n = -n;
  26                 z_div(&x, &one, a);
  27                 }
  28         else
  29                 {
  30                 x.r = a->r;
  31                 x.i = a->i;
  32                 }
  33
  34         for(u = n; ; )
  35                 {
  36                 if(u & 01)
  37                         {
  38                         t = q.r * x.r - q.i * x.i;
  39                         q.i = q.r * x.i + q.i * x.r;
  40                         q.r = t;
  41                         }
  42                 if(u >>= 1)
  43                         {
  44                         t = x.r * x.r - x.i * x.i;
  45                         x.i = 2 * x.r * x.i;
  46                         x.r = t;
  47                         }
  48                 else
  49                         break;
  50                 }
  51  done:
  52         p->i = q.i;
  53         p->r = q.r;
  54         }