More Emacs Calculator Functionality

The first video demonstrates algebraic formulas and live evaluation in Emacs Calculator:

The second video covers two subjects: (1) mapping functions over vectors, and (2) using emacs to display algebraic formulas as math LaTeX (to paste into a WordPress post, for example). Please forgive the improper pronunciation of “LaTeX”, which I remembered afterwards.

Emacs Calc: Angle Between Vectors

In my geometry studies, I learned that one can get the angle between two vectors with this formula:

cos \theta = \frac{V_1 \cdot V_2}{| V_1 | | V_2 |}

I.e., the cosine of the angle equals the dot product of the two vectors over the product of their magnitudes.

Here we get about 1.05 radians or about 60.26 degrees. A cool thing about this formula is it works for vectors of any (matching) dimension, i.e., 3-D coordinates, 4-D coordinates, etc.

This is definitely doable in Emacs Calc, since we have a dot product function, called inner-product (press ‘x inner product’), But doing the angle formula involves a lot of steps, with either stack rotation or storing the vectors in variables. So I wanted to get the angle formula stored as a calc formula. Unfortunately, inner-product itself is only an interactive function, so this was problematic. However, inner-product actual calls another function, inner. So, this formula is possible:

arccos(inner(mul, add, v1, v2) / (abs(v1) abs(v2)))

How do you store this formula in Emacs? I could walk you through the steps described in section 18.4 of the Emacs Calc info manual, but the end result is that this code is stored in your ~/.emacs.d/calc.el:

(put 'calc-define 'calc-vectorsangle '(progn
 (defun calc-vectorsangle nil (interactive) (calc-wrapper (calc-enter-result 2 "vect" (cons (quote calcFunc-vectorsangle) (calc-top-list-n 2)))))
 (put 'calc-vectorsangle 'calc-user-defn 't)
 (defun calcFunc-vectorsangle (v1 v2) (math-normalize (list (quote
  calcFunc-arccos) (list (quote /) (list (quote calcFunc-inner) (quote
  (var mul var-mul)) (quote (var add var-add)) v1 v2) (list (quote *)
  (list (quote calcFunc-abs) v1) (list (quote calcFunc-abs) v2))))))
 (put 'calcFunc-vectorsangle 'calc-user-defn '(calcFunc-arccos (/
  (calcFunc-inner (var mul var-mul) (var add var-add) (var v1 var-v1)
  (var v2 var-v2)) (* (calcFunc-abs (var v1 var-v1)) (calcFunc-abs (var
  v2 var-v2))))))
 (define-key calc-mode-map "zA" 'calc-vectorsangle)
))

Then when you start calc, you can put two vectors on the stack, and enter command ‘x vectorsangle’ or ‘z A’ for short:

This will work for vectors with as many more coordinates as you want, so long as there are the same number in each vector.

Fractal with Barycentric Coordinates

This fractal idea did not originate with me, but I wrote some racket code to do the midpoint calculation using barycentric coordinates. This fractal draws a circle at the midpoint of a triangle, then subdivides the triangle and repeats:

Here is the same fractal to four iterations:

To get the midpoints, I could simple pass in the coordinates of the last triangle ABC, and then use “0.5” barycentric coordinates:

        [P1 (barycentric->complex 0.5 0.5 0 A B C)]
        [P2 (barycentric->complex 0.5 0.0 0.5 A B C)]
        [P3 (barycentric->complex 0 0.5 0.5 A B C)]

Here is the full code:

ftp://lavender.qlfiles.net/Racket/bc-fractal.7z

Barycentric Coordinates

An interesting concept I came across recently in an introductory geometry book is barycentric coordinates. The fundamental idea is defining the location of a point by its relationship to the points of a triangle. (Well, technically you can use any polygon, but we will keep things simple.)

So, each of the vertices of the triangle (A, B, and C) are given a sort of weight or pull factor. These three numbers uniquely represent the position of P. The theorem goes

If A, B, C are three non-collinear points, any vector P may be expressed in terms of the vectors A, B, C thus:

P = xA + yB +zC, where x + y + z = 1

Geometry: A Comprehensive Course, Dan Pedoe

If P is pulled all the way over to vertex A, then the coordinates become (1, 0, 0). Likewise at B it is (0, 1, 0) and for C it is (0, 0, 1). Another coordinate inside the triangle would be (⅓, ⅓, ⅓). If the point moves outside the triangle, some of those coordinates become negative.

Of course, instead of boring old vectors, you can use… complex numbers!

From the formula, it is obvious how to start with some barycentric coordinates, add a triangle, and find the point. Here is some Racket code:

(define (barycentric->complex x y z A B C)
   (+ (* x A) (* y B) (* z C)))

Note, that you can use any triangle you want with the same barycentric coordinates. This could allow you to stretch the graph of some set of coordinates in various ways.

But, how do you convert in the other direction, having a point and a triangle, and wanting the barycentric coordinates? This is a little more complicated. Fortunately, it happens that the factor for each of the vertices is equal to the proportional area of the triangle opposite that vertex. So, looking at this diagram again

x = [BCP] / [ABC], y = [CAP] / [ABC], z = [ABP] / [ABC]
where [BCP] is the area of triangle BCP, etc.

Now, a tricky part here is that the areas must be signed areas, or the formula doesn’t necessarily work correctly. Signed area means that you have to be consistent in assigning a negative or positive value to the area of the triangle, depending on the orientation of the vertices. So, if you happen to plug in the vertices in a counter-clockwise direction, you get a negative area, and in a clockwise direction, you get a positive area. Thankfully, there is an elegant formula for determining this, which I will represent by this Racket function

(define (orientation Z1 Z2 Z3)
  (let* ([rise1 (- (imag-part Z2) (imag-part Z1))]
         [rise2 (- (imag-part Z3) (imag-part Z2))]
         [run1 (- (real-part Z2) (real-part Z1))]
         [run2 (- (real-part Z3) (real-part Z2))]
         [diff (- (* rise1 run2) (* rise2 run1))])
    (if (< diff 0) 'ccw
        (if (> diff 0) 'cw 'col))))

The function returns ‘ccw for counter-clockwise, ‘cw for clockwise, and ‘col for collinear (i.e., the points are on a straight line, which is not relevant in this application).

This allows us to have a signed-area function, using the common formula for determining the area of a triangle from just the length of the sides. We get the length of the sides by subtracting the complex numbers from each other and pulling the magnitudes of the resulting complex numbers.

(define (signed-area A B C)
  (let* ([a (magnitude (- C B))]
         [b (magnitude (- A C))]
         [c (magnitude (- A B))]
         [s (* (+ a b c) 0.5)]
         [area (sqrt (* s (- s a) (- s b) (- s c)))]
         [orient (orientation A B C)])
    (if (eq? orient 'ccw)
        (* -1 area)
        area)))

And finally we determine x, y, and z:

(define (complex->barycentric A B C P)
  (let ([ABC (signed-area A B C)])
    (list (/ (signed-area B C P) ABC)
          (/ (signed-area C A P) ABC)
          (/ (signed-area A B P) ABC))))

See, simple as Pi! Heh, heh.

racket@barycentric.rkt> (complex->barycentric 0 6 2+2i 2+i)
'(0.33333333333333376 0.16666666666666657 0.5000000000000009)
racket@barycentric.rkt> (barycentric->complex 0.33333333333333376 0.16666666666666657 0.5000000000000009 0 6 2+2i)
2.0000000000000013+1.0000000000000018i

So, you can see it works, though there is some very small inaccuracy introduced by the approximations that occur during the conversion operations.

I hope to explore, in future posts, some interesting applications of barycentric coordinates.

Trigonometry Problem without Trigonometry: Clarification/Correction

In the previous post I was exulting about the coolness of solving a trigonometry problem without trigonometry, using complex numbers. Well, that isn’t quite the full truth. You can certainly do the problem described easily and elegantly with a complex number capable calculator. However, there is some hidden trigonometry in the initial rotation operation:

The complex number ei has to be raised to a power. If you raise to an integer power, you could simply multiply out the rectangular values algebraically. But how do you raise a complex number to a fractional power, like 2.5? The standard approach, near as I can tell, is to convert the complex number to polar format, with a magnitude r and an angle ϕ. Assuming complex number a+bi, pulling out the magnitude requires only the sqrt function (pythagorean theorem) but pulling out the angle requires the arctangent of b/a. Then, use the formula Zn = rn (cos nϕ + i sin nϕ). I really don’t know if this is how all calculators do it, but here is the code behind the operations in the racket interpreter:

  Scheme_Complex *cb = (Scheme_Complex *)base;
  Scheme_Complex *ce = (Scheme_Complex *)exponent;
  double a, b, c, d, bm, ba, nm, na, r1, r2;
  int d_is_zero;

  if ((ce->i == zero) && !SCHEME_FLOATP(ce->r)) {
    if (SCHEME_INTP(ce->r) || SCHEME_BIGNUMP(ce->r))
      return scheme_generic_integer_power(base, ce->r);
  }

  a = scheme_get_val_as_double(cb->r);
  b = scheme_get_val_as_double(cb->i);
  c = scheme_get_val_as_double(ce->r);
  d = scheme_get_val_as_double(ce->i);
  d_is_zero = (ce->i == zero);

  bm = sqrt(a * a + b * b);
  ba = atan2(b, a);

  /* New mag & angle */
  nm = scheme_double_expt(bm, c) * exp(-(ba * d));
  if (d_is_zero) /* precision here can avoid NaNs */
    na = ba * c;
  else
    na = log(bm) * d + ba * c;

  r1 = nm * cos(na);
  r2 = nm * sin(na);

The code there is just slightly more complicated, as it has separate angle finding conditional branches for real number and complex number exponents.

Of course, you don’t need to deal with any of that if you are just manipulating complex numbers in the interpreter or calculator.

I feel this would be a great place to dive into the interesting subject of complex roots, but that I should hold off until I have a more thorough grasp of the subject.