abs.cpp

   1 // Absolute value, aka vector magnitude
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 void
   7 eval_abs(void)
   8 {
   9         push(cadr(p1));
  10         eval();
  11         absval();
  12 }
  13
  14 void
  15 absval(void)
  16 {
  17         int h;
  18         save();
  19         p1 = pop();
  20
  21         if (istensor(p1)) {
  22                 absval_tensor();
  23                 restore();
  24                 return;
  25         }
  26
  27         if (isnum(p1)) {
  28                 push(p1);
  29                 if (isnegativenumber(p1))
  30                         negate();
  31                 restore();
  32                 return;
  33         }
  34
  35         if (iscomplexnumber(p1)) {
  36                 push(p1);
  37                 push(p1);
  38                 conjugate();
  39                 multiply();
  40                 push_rational(1, 2);
  41                 power();
  42                 restore();
  43                 return;
  44         }
  45
  46         // abs(1/a) evaluates to 1/abs(a)
  47
  48         if (car(p1) == symbol(POWER) && isnegativeterm(caddr(p1))) {
  49                 push(p1);
  50                 reciprocate();
  51                 absval();
  52                 reciprocate();
  53                 restore();
  54                 return;
  55         }
  56
  57         // abs(a*b) evaluates to abs(a)*abs(b)
  58
  59         if (car(p1) == symbol(MULTIPLY)) {
  60                 h = tos;
  61                 p1 = cdr(p1);
  62                 while (iscons(p1)) {
  63                         push(car(p1));
  64                         absval();
  65                         p1 = cdr(p1);
  66                 }
  67                 multiply_all(tos - h);
  68                 restore();
  69                 return;
  70         }
  71
  72         if (isnegativeterm(p1) || (car(p1) == symbol(ADD) && isnegativeterm(cadr(p1)))) {
  73                 push(p1);
  74                 negate();
  75                 p1 = pop();
  76         }
  77
  78         push_symbol(ABS);
  79         push(p1);
  80         list(2);
  81
  82         restore();
  83 }
  84
  85 void
  86 absval_tensor(void)
  87 {
  88         if (p1->u.tensor->ndim != 1)
  89                 stop("abs(tensor) with tensor rank > 1");
  90         push(p1);
  91         push(p1);
  92         conjugate();
  93         inner();
  94         push_rational(1, 2);
  95         power();
  96         simplify();
  97         eval();
  98 }