lang/fortran/lib/libF77/pow_zi.c

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