arg.cpp

   1 /* Argument (angle) of complex z
   2
   3         z               arg(z)
   4         -               ------
   5
   6         a               0
   7
   8         -a              -pi                     See note 3 below
   9
  10         (-1)^a          a pi
  11
  12         exp(a + i b)    b
  13
  14         a b             arg(a) + arg(b)
  15
  16         a + i b         arctan(b/a)
  17
  18 Result by quadrant
  19
  20         z               arg(z)
  21         -               ------
  22
  23         1 + i           1/4 pi
  24
  25         1 - i           -1/4 pi
  26
  27         -1 + i          3/4 pi
  28
  29         -1 - i          -3/4 pi
  30
  31 Notes
  32
  33         1. Handles mixed polar and rectangular forms, e.g. 1 + exp(i pi/3)
  34
  35         2. Symbols in z are assumed to be positive and real.
  36
  37         3. Negative direction adds -pi to angle.
  38
  39            Example: z = (-1)^(1/3), mag(z) = 1/3 pi, mag(-z) = -2/3 pi
  40
  41         4. jean-francois.debroux reports that when z=(a+i*b)/(c+i*d) then
  42
  43                 arg(numerator(z)) - arg(denominator(z))
  44
  45            must be used to get the correct answer. Now the operation is
  46            automatic.
  47 */
  48
  49 #include "stdafx.h"
  50 #include "defs.h"
  51
  52 void
  53 eval_arg(void)
  54 {
  55         push(cadr(p1));
  56         eval();
  57         arg();
  58 }
  59
  60 void
  61 arg(void)
  62 {
  63         save();
  64         p1 = pop();
  65         push(p1);
  66         numerator();
  67         yyarg();
  68         push(p1);
  69         denominator();
  70         yyarg();
  71         subtract();
  72         restore();
  73 }
  74
  75 #define RE p2
  76 #define IM p3
  77
  78 void
  79 yyarg(void)
  80 {
  81         save();
  82         p1 = pop();
  83         if (isnegativenumber(p1)) {
  84                 push(symbol(PI));
  85                 negate();
  86         } else if (car(p1) == symbol(POWER) && equaln(cadr(p1), -1)) {
  87                 // -1 to a power
  88                 push(symbol(PI));
  89                 push(caddr(p1));
  90                 multiply();
  91         } else if (car(p1) == symbol(POWER) && cadr(p1) == symbol(E)) {
  92                 // exponential
  93                 push(caddr(p1));
  94                 imag();
  95         } else if (car(p1) == symbol(MULTIPLY)) {
  96                 // product of factors
  97                 push_integer(0);
  98                 p1 = cdr(p1);
  99                 while (iscons(p1)) {
 100                         push(car(p1));
 101                         arg();
 102                         add();
 103                         p1 = cdr(p1);
 104                 }
 105         } else if (car(p1) == symbol(ADD)) {
 106                 // sum of terms
 107                 push(p1);
 108                 rect();
 109                 p1 = pop();
 110                 push(p1);
 111                 real();
 112                 RE = pop();
 113                 push(p1);
 114                 imag();
 115                 IM = pop();
 116                 if (iszero(RE)) {
 117                         push(symbol(PI));
 118                         if (isnegative(IM))
 119                                 negate();
 120                 } else {
 121                         push(IM);
 122                         push(RE);
 123                         divide();
 124                         arctan();
 125                         if (isnegative(RE)) {
 126                                 push_symbol(PI);
 127                                 if (isnegative(IM))
 128                                         subtract();     // quadrant 1 -> 3
 129                                 else
 130                                         add();          // quadrant 4 -> 2
 131                         }
 132                 }
 133         } else
 134                 // pure real
 135                 push_integer(0);
 136         restore();
 137 }