libf2c/libF77/pow_zi.c

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