// Absolute value, aka vector magnitude #include "stdafx.h" #include "defs.h" void eval_abs(void) { push(cadr(p1)); eval(); absval(); p1 = pop(); push(p1); if (car(p1) != symbol(ABS)) { // to prevent circular reference simplify(); eval(); // to normalize } } void absval(void) { save(); yyabsval(); restore(); } void yyabsval(void) { int h; p1 = pop(); if (istensor(p1)) { absval_tensor(); return; } if (isnum(p1)) { push(p1); if (isnegativenumber(p1)) negate(); return; } if (iscomplexnumber(p1)) { push(p1); push(p1); conjugate(); multiply(); push_rational(1, 2); power(); return; } // abs(1/a) evaluates to 1/abs(a) if (car(p1) == symbol(POWER) && isnegativeterm(caddr(p1))) { push(p1); reciprocate(); absval(); reciprocate(); return; } // abs(a*b) evaluates to abs(a)*abs(b) if (car(p1) == symbol(MULTIPLY)) { h = tos; p1 = cdr(p1); while (iscons(p1)) { push(car(p1)); absval(); p1 = cdr(p1); } multiply_all(tos - h); return; } if (isnegativeterm(p1) || (car(p1) == symbol(ADD) && isnegativeterm(cadr(p1)))) { push(p1); negate(); p1 = pop(); } push_symbol(ABS); push(p1); list(2); } void absval_tensor(void) { if (p1->u.tensor->ndim != 1) stop("abs(tensor) with tensor rank > 1"); push(p1); push(p1); conjugate(); inner(); push_rational(1, 2); power(); } #if SELFTEST static char *s[] = { "abs(2)", "2", "abs(2.0)", "2", "abs(-2)", "2", "abs(-2.0)", "2", "abs(a)", "abs(a)", "abs(-a)", "abs(a)", "abs(2*a)", "2*abs(a)", "abs(-2*a)", "2*abs(a)", "abs(2.0*a)", "2*abs(a)", "abs(-2.0*a)", "2*abs(a)", "abs(a-b)+abs(b-a)", "2*abs(a-b)", "abs(3 + 4 i)", "5", "abs((2,3,4))", "29^(1/2)", "abs(a*b)", "abs(a)*abs(b)", "abs(a/b)", "abs(a)/abs(b)", "abs(1/a^b)", "1/(abs(a^b))", // this test makes sure the result of abs is simplified "P=(u*cos(v),u*sin(v),v)", "", "abs(cross(d(P,u),d(P,v)))", "(1+u^2)^(1/2)", }; void test_abs(void) { test(__FILE__, s, sizeof s / sizeof (char *)); } #endif