complex numbers

The evil of COBOL: everything is in global variables

December 7, 2023 COBOL No comments , , , , , , ,

COBOL does not have stack variables.  Everything is a global variable.  There is a loose equivalent of a function, called a paragraph, which can be called using a PERFORM statement, but a paragraph does not have any input or output variables, and no return code, so if you want it to behave like a function, you have to construct some sort of complicated naming convention using your global variables. 

I’ve seen real customer COBOL programs with many thousands of global variables.  A production COBOL program is usually a giant sequence of MOVEs, MOVE A TO B, MOVE B TO C, MOVE C TO D, MOVE D TO E, … with various PERFORMs or GOTOs, or other things in between.  If you find that your variable has a bad value in it, that is probably because it has been copied from something that was copied from something, that was copied from something, that’s the output of something else, that was copied from something, 9 or 10 times.

I was toying around with the idea of coding up a COBOL implementation of 2D Euclidean geometric algebra, just as a joke, as it is surely the worst language in the world.  Yes, I work on a COBOL compiler project. The project is a lot of fun, and the people I work with are awesome, but I don’t have to like the language.

If I was to implement this simplest geometric algebra in COBOL, the logical starting place for that would be to implement complex numbers in COBOL first.  That is because we can use a pair of complex numbers to implement a 2D multivector, with one complex number for the vector part, and a complex number for the scalar and pseudoscalar parts.  That technique has been detailed on this blog previously, and also in a Mathematica module Cl20.m.

Trying to implement a couple of complex number operations in COBOL got absurd really fast.  Here’s an example.  First step was to create some complex number types.  I did that with a copybook (include file), like so:

This can be included multiple times, each time with a different name, like so:

The way that I structured all my helper functions, was with one set of global variables for input (at least one), and if appropriate, one output global variable.  Here’s an example:

So, if I want to compute and display a value, I have a whole pile of stupid MOVEs to do in and out of the appropriate global variables for each of the helper routines in question:

I wrote enough of this little complex number library that I could do conjugate, real, imaginary, multiply, inverse, and divide operations.  I can run that little program with the following JCL

//COMPLEX JOB
//A EXEC PGM=COMPLEX
//SYSOUT   DD SYSOUT=*
//STEPLIB  DD DSN=PJOOT.SAMPLE.COMPLEX,
//  DISP=SHR

and get this SYSOUT:

STEP A SYSOUT:
A                    =  .10000000000000000E 01 + ( .20000000000000000E 01) I
B                    =  .30000000000000000E 01 + ( .40000000000000000E 01) I
CONJ(A)              =  .10000000000000000E 01 + (-.20000000000000000E 01) I
RE(A)                =  .10000000000000000E 01
IM(A)                =  .20000000000000000E 01
A * B                = -.50000000000000000E 01 + ( .10000000000000000E 02) I
1/A                  =  .20000000000000000E 00 + (-.40000000000000000E 00) I
A/B                  =  .44000000000000000E 00 + ( .80000000000000000E-01) I

If you would like your eyes burned further, you can access the full program on github here. It takes almost 200 lines of code to do almost nothing.

A silly geometry problem: length of side of square in circular quadrant

September 27, 2023 math and physics play , ,

[Click here for a PDF version of this post]

Problem from:

My solution (before numerical reduction), using basic trig and complex numbers, is illustrated in fig. 1.

fig. 1. With complex numbers.

We have
\begin{equation}\label{eqn:squareInCircle:20}
\begin{aligned}
s &= x \cos\theta \\
y &= x \sin\theta \\
p &= y + x e^{i\theta} \\
q &= i s + x e^{i\theta} \\
\Abs{q} &= y + 5 \\
\Abs{p – q} &= 2.
\end{aligned}
\end{equation}
This can be reduced to
\begin{equation}\label{eqn:squareInCircle:40}
\begin{aligned}
\Abs{ x e^{i\theta} – 5 } &= 2 \\
x \Abs{ i \cos\theta + e^{i\theta} } &= x \sin\theta + 5.
\end{aligned}
\end{equation}

My wife figured out how to do it with just Pythagoras, as illustrated in fig. 2.

fig. 2. With Pythagoras.

fig. 2. With Pythagoras.

\begin{equation}\label{eqn:squareInCircle:60}
\begin{aligned}
\lr{ 5 – s }^2 + y^2 &= 4 \\
\lr{ s + y }^2 + s^2 &= \lr{ y + 5 }^2 \\
x^2 &= s^2 + y^2.
\end{aligned}
\end{equation}

Either way, the numerical solution is 4.12. The geometry looks like fig. 3.

fig. 3. Lengths to scale.

A mathematica notebook to compute the numerical part of the problem (either way) and plot the figure to scale can be found in my mathematica github repo.

A complex-pair representation of GA(2,0).

June 15, 2023 math and physics play , , , ,

[Click here for a PDF version of this post]

Motivation.

Suppose that we want to represent GA(2,0) (Euclidean) multivectors as a pair of complex numbers, with a structure like
\begin{equation}\label{eqn:bicomplexGA20:20}
M = (m_1, m_2),
\end{equation}
where
\begin{equation}\label{eqn:bicomplexGA20:40}
\begin{aligned}
\gpgrade{M}{0,2} &\sim m_1 \\
\gpgrade{M}{1} &\sim m_2.
\end{aligned}
\end{equation}
Specifically
\begin{equation}\label{eqn:bicomplexGA20:60}
\begin{aligned}
\gpgrade{M}{0} &= \textrm{Re}(m_1) \\
\gpgrade{M}{1} \cdot \Be_1 &= \textrm{Re}(m_2) \\
\gpgrade{M}{1} \cdot \Be_2 &= \textrm{Im}(m_2) \\
\gpgrade{M}{2} i^{-1} &= \textrm{Im}(m_1),
\end{aligned}
\end{equation}
where \( i \sim \Be_1 \Be_2 \).

Multivector product representation.

Let’s figure out how we can represent the various GA products, starting with the geometric product. Let
\begin{equation}\label{eqn:bicomplexGA20:80}
\begin{aligned}
M &= \gpgrade{M}{0,2} + \gpgrade{M}{1} = (m_1, m_2) = (m_{11} + m_{12} i, m_{21} + m_{22} i) \\
N &= \gpgrade{N}{0,2} + \gpgrade{N}{1} = (n_1, n_2) = (n_{11} + n_{12} i, n_{21} + n_{22} i),
\end{aligned}
\end{equation}
so
\begin{equation}\label{eqn:bicomplexGA20:200}
\begin{aligned}
M N
&= \gpgrade{M}{0,2} \gpgrade{N}{0,2} + \gpgrade{M}{1} \gpgrade{N}{1} \\
&\quad+ \gpgrade{M}{1} \gpgrade{N}{0,2} + \gpgrade{M}{0,2} \gpgrade{N}{1}
\end{aligned}
\end{equation}

The first two terms have only even grades, and the second two terms are vectors. The complete representation of the even grade components of this multivector product is
\begin{equation}\label{eqn:bicomplexGA20:240}
\gpgrade{M N}{0,2} \sim \lr{ m_1 n_1 + \textrm{Re}(m_2 n_2^\conj) – i \textrm{Im}(m_2 n_2^\conj), 0 },
\end{equation}
or
\begin{equation}\label{eqn:bicomplexGA20:260}
\begin{aligned}
\gpgrade{M N}{0} &= \textrm{Re}\lr{ m_1 n_1 + m_2 n_2^\conj } \\
\gpgrade{M N}{2} i^{-1} &= \textrm{Im}\lr{ m_1 n_1 – m_2 n_2^\conj }.
\end{aligned}
\end{equation}

For the vector components we have
\begin{equation}\label{eqn:bicomplexGA20:280}
\begin{aligned}
\gpgrade{M N}{1}
&=
\gpgrade{M}{1} \gpgrade{N}{0} + \gpgrade{M}{0} \gpgrade{N}{1}
+
\gpgrade{M}{1} \gpgrade{N}{2} + \gpgrade{M}{2} \gpgrade{N}{1} \\
&= \gpgrade{M}{1} n_{11} + m_{11} \gpgrade{N}{1} + \gpgrade{M}{1} i n_{12} + i m_{12} \gpgrade{N}{1}.
\end{aligned}
\end{equation}
For these,
\begin{equation}\label{eqn:bicomplexGA20:300}
\begin{aligned}
\gpgrade{M}{1} i
&= \lr{ m_{21} \Be_1 + m_{22} \Be_2 } \Be_{12}
&= – m_{22} \Be_1 + m_{21} \Be_2,
\end{aligned}
\end{equation}
and
\begin{equation}\label{eqn:bicomplexGA20:320}
\begin{aligned}
i \gpgrade{N}{1}
&= \Be_{12} \lr{ n_{21} \Be_1 + n_{22} \Be_2 }
&=
n_{22} \Be_1 – n_{21} \Be_2.
\end{aligned}
\end{equation}
Comparing to
\begin{equation}\label{eqn:bicomplexGA20:340}
i (a + i b)
= -b + i a,
\end{equation}
we see that
\begin{equation}\label{eqn:bicomplexGA20:360}
\gpgrade{M N}{1}
\sim
\lr{ 0, n_{11} m_2 + m_{11} n_2 + n_{12} i m_2 – m_{12} i n_2 }.
\end{equation}
If we want the vector coordinates, those are
\begin{equation}\label{eqn:bicomplexGA20:380}
\begin{aligned}
\gpgrade{M N}{1} \cdot \Be_1 &= \textrm{Re} \lr{ n_{11} m_2 + m_{11} n_2 + n_{12} i m_2 – m_{12} i n_2 } \\
\gpgrade{M N}{1} \cdot \Be_2 &= \textrm{Im} \lr{ n_{11} m_2 + m_{11} n_2 + n_{12} i m_2 – m_{12} i n_2 }.
\end{aligned}
\end{equation}

Summary.

\begin{equation}\label{eqn:bicomplexGA20:n}
M N \sim
\lr{ m_1 n_1 + \textrm{Re}(m_2 n_2^\conj) – i \textrm{Im}(m_2 n_2^\conj), n_{11} m_2 + m_{11} n_2 + n_{12} i m_2 – m_{12} i n_2 }.
\end{equation}

A sample Mathematica implementation is available, as well as an example notebook (that also doubles as a test case.)

Clarification.

I skipped a step above, showing the correspondances to the dot and wedge product.

Let \(z = a + bi\), and \(w = c + di\). Then:
\begin{equation}\label{eqn:bicomplexGA20:420}
\begin{aligned}
z w^\conj
&=
\lr{ a + b i } \lr{ c – d i } \\
&= a c + b d -i \lr{ a d – b c }.
\end{aligned}
\end{equation}
Compare that to the geometric product of two vectors \( \Bx = a \Be_1 + b \Be_2 \), and \( \By = c \Be_1 + d \Be_2 \), where we have
\begin{equation}\label{eqn:bicomplexGA20:440}
\begin{aligned}
\Bx \By &= \Bx \cdot \By + \Bx \wedge \By \\
&= \lr{ a \Be_1 + b \Be_2 } \lr{ c \Be_1 + d \Be_2 } \\
&= a c + b d + \Be_1 \Be_2 \lr{ a d – b c }.
\end{aligned}
\end{equation}
So we have
\begin{equation}\label{eqn:bicomplexGA20:460}
\begin{aligned}
a b + cd &= \Bx \cdot \By = \textrm{Re} \lr{ z w^\conj } \\
a d – b c &= \lr{ \Bx \wedge \By } \Be_{12}^{-1} = – \textrm{Im} \lr{ z w^\conj }.
\end{aligned}
\end{equation}
We see that \( \lr{z w^\conj}^\conj = z^\conj w \) can be used as a representation of the geometric product of two vectors (setting \( i = \Be_1 \Be_2 \) as usual.)

Another simplification.

We have sums of the form
\begin{equation}\label{eqn:bicomplexGA20:480}
\textrm{Re}(z) w \pm \textrm{Im}(z) i w
\end{equation}
above. Let’s see if those can be simplified. For the positive case we have
\begin{equation}\label{eqn:bicomplexGA20:500}
\begin{aligned}
\textrm{Re}(z) w + \textrm{Im}(z) i w
&=
\inv{2} \lr{ z + z^\conj } w + \inv{2} \lr{ z – z^\conj } w \\
&=
z w,
\end{aligned}
\end{equation}
and for the negative case, we have
\begin{equation}\label{eqn:bicomplexGA20:520}
\begin{aligned}
\textrm{Re}(z) w – \textrm{Im}(z) i w
&=
\inv{2} \lr{ z + z^\conj } w – \inv{2} \lr{ z – z^\conj } w \\
&=
z^\conj w.
\end{aligned}
\end{equation}
This, with the vector-vector product simplification above, means that we can represent the full multivector product in this representation as just
\begin{equation}\label{eqn:bicomplexGA20:540}
M N \sim
\lr{ m_1 n_1 + m_2^\conj n_2, m_2 n_1 + m_1^\conj n_2 }.
\end{equation}

Inscribed Triangle in circle problem

December 25, 2022 math and physics play , , , , , ,

[Click here for a PDF version of this post]

In the LinkedIn Pre-University Geometric Algebra group, James presents a problem from the MindYourDecisions youtube channel Impossible Viral Problem, as a candidate for solution using geometric algebra.

I tried this out and found a couple ways to solve it. One of those I’ll detail here. I have to admit that part of the reason that I wanted to solve this is that the figure in the beginning of the video really bugged me. The triangle that was inscribed in the circle didn’t have any of the length properties from the problem. I could do much better with a sloppy freehand sketch, but to do a good figure, you have to actually solve for the vertexes of the triangle (once you do that, the area is easy to figure out.)

Formulating the problem.

Having solved the problem, the geometry of the problem is illustrated in fig. 1.

fig. 1. Inscribed triangle in circle.

fig. 1. Inscribed triangle in circle.

I set up the problem so that the \( A,C \) triangle vertices were symmetric with respect to the x-axis, and the \(B \) vertex located elsewhere. I can describe those algebraically as
\begin{equation}\label{eqn:inscribedTriangleProblem:20}
\begin{aligned}
\BA &= r \Be_1 e^{i\theta} \\
\BC &= r \Be_1 e^{-i\theta} \\
\BB &= r \Be_1 e^{i\phi},
\end{aligned}
\end{equation}
where the radius \( r \) and two angles \( \theta, \phi \) are to be determined, and \( i = \Be_1 \Be_1\) the pseudoscalar for the \(x-y\) plane.
The vector pointing to the midpoint of the upper triangular face is given by the average of the \( \BA, \BB \) vectors, which can be seen from
\begin{equation}\label{eqn:inscribedTriangleProblem:40}
\BA + \frac{\BB – \BA}{2} = \frac{\BA + \BB}{2},
\end{equation}
and similarly, the midpoint of the lower face is found at
\begin{equation}\label{eqn:inscribedTriangleProblem:60}
\BC + \frac{\BB – \BC}{2} = \frac{\BB + \BC}{2},
\end{equation}
The problem tells us that the respective lengths of those vectors from the origin are \( r-2, r – 3\) respectively, so
\begin{equation}\label{eqn:inscribedTriangleProblem:80}
\begin{aligned}
r – 2 &= \inv{2} \Abs{ \BA + \BB } \\
r – 3 &= \inv{2} \Abs{ \BB + \BC },
\end{aligned}
\end{equation}
or
\begin{equation}\label{eqn:inscribedTriangleProblem:100}
\begin{aligned}
(r – 2)^2 &= \frac{r^2}{4} \lr{ \Be_1 e^{i\theta} + \Be_1 e^{i\phi} }^2 \\
(r – 3)^2 &= \frac{r^2}{4} \lr{ \Be_1 e^{i\phi} + \Be_1 e^{-i\theta} }^2 \\
\end{aligned}
\end{equation}
Finally, since the midpoint of the right edge is found at \( (r-1)\Be_1 \), it is clear that
\begin{equation}\label{eqn:inscribedTriangleProblem:120}
\frac{r-1}{r} = \cos\theta,
\end{equation}
or
\begin{equation}\label{eqn:inscribedTriangleProblem:140}
r = \inv{1 – \cos\theta}.
\end{equation}
This leaves us with three equations and three unknowns. Unfortunately, these are rather non-linear equations. In the video, a direct method of solving equivalent equations was demonstrated, but I picked the lazy route, and used Mathematica’s NSolve routine, solving for \( r,\theta, \phi\) numerically. Since NSolve has intrinsic complex number support, I made the following substitutions:
\begin{equation}\label{eqn:inscribedTriangleProblem:160}
\begin{aligned}
z &= e^{i\theta} \\
w &= e^{i\phi},
\end{aligned}
\end{equation}
and then plugged those into our relations above, after expanding the squares, to find
\begin{equation}\label{eqn:inscribedTriangleProblem:180}
\begin{aligned}
\lr{ \Be_1 e^{i\theta} + \Be_1 e^{i\phi} }^2
&=
2 + \Be_1 e^{i\theta} \Be_1 e^{i\phi} + \Be_1 e^{i\phi} \Be_1 e^{i\theta} \\
&=
2 + e^{-i\theta} \Be_1^2 e^{i\phi} + e^{-i\phi} \Be_1^2 e^{i\theta} \\
&=
2 + e^{-i\theta} \Be_1^2 e^{i\phi} + e^{-i\phi} \Be_1^2 e^{i\theta} \\
&=
2 + \frac{w}{z} + \frac{z}{w},
\end{aligned}
\end{equation}
and
\begin{equation}\label{eqn:inscribedTriangleProblem:200}
\begin{aligned}
\lr{ \Be_1 e^{i\phi} + \Be_1 e^{-i\theta} }^2
&=
2 + \Be_1 e^{i\phi} \Be_1 e^{-i\theta} + \Be_1 e^{-i\theta} \Be_1 e^{i\phi} \\
&=
2 + e^{-i\phi} e^{-i\theta} + e^{ i\theta} e^{i\phi} \\
&=
2 + w z + \inv{w z}.
\end{aligned}
\end{equation}
This gives us
\begin{equation}\label{eqn:inscribedTriangleProblem:220}
\begin{aligned}
4 \lr{ \frac{r – 2 }{r} }^2 &= 2 + \frac{w}{z} + \frac{z}{w} \\
4 \lr{ \frac{r – 3 }{r} }^2 &= 2 + w z + \inv{w z},
\end{aligned}
\end{equation}
where
\begin{equation}\label{eqn:inscribedTriangleProblem:240}
r = \inv{1 – \inv{2}\lr{ z + \inv{z}}}.
\end{equation}

The NSolve gave me some garbage solutions (like \(\theta = 0 \)) that must have been valid numerically, but did not encode the geometry of the problem, so I added a few additional constraints to the problem, namely
\begin{equation}\label{eqn:inscribedTriangleProblem:260}
\begin{aligned}
z \bar{z} &= 1 \\
w \bar{w} &= 1 \\
\inv{2} \lr{ z + \inv{z} } &\ne 1 \\
1/(1 – (1/2) \textrm{Re}(z + 1/z)) &> 3.
\end{aligned}
\end{equation}
This provided exactly two solutions, but when plotted, they turn out to just be mirror images of each other. After back substitution, the solution illustrated above was given by
\begin{equation}\label{eqn:inscribedTriangleProblem:280}
\begin{aligned}
r &= 3.87939 \\
\theta &= 42.078 \\
\phi &= 164.125,
\end{aligned}
\end{equation}
where these angles are in degrees, not radians.

The triangular area.

There are probably lots of formulas for the area of a triangle (that I have forgotten), but we can compute it easily by doubling the triangle, forming a parallelogram, to find
\begin{equation}\label{eqn:inscribedTriangleProblem:300}
\textrm{Area} = \inv{2} \Abs{ \lr{ \BA – \BC } \wedge {\BC – \BB } },
\end{equation}
or
\begin{equation}\label{eqn:inscribedTriangleProblem:320}
\begin{aligned}
\textrm{Area}^2
&= \frac{-1}{4} \lr{ \lr{ \BA – \BC } \wedge \lr{\BC – \BB } }^2 \\
&= \frac{-1}{4} \lr{ \BA \wedge \BC – \BA \wedge \BB + \BC \wedge \BB }^2 \\
&= \frac{-r^4}{4} \lr{\gpgradetwo{ \Be_1 e^{i\theta} \Be_1 e^{-i\theta} – \Be_1 e^{i\theta} \Be_1 e^{i\phi} + \Be_1 e^{-i\theta} \Be_1 e^{i\phi} }}^2 \\
&= \frac{-r^4}{4} \lr{\gpgradetwo{ e^{-2 i \theta} – e^{i \phi -i\theta} + e^{i\theta + i \phi} }}^2,
\end{aligned}
\end{equation}
so
\begin{equation}\label{eqn:inscribedTriangleProblem:340}
\textrm{Area} = \frac{r^2}{2} \Abs{ -\sin( 2 \theta ) – \sin(\phi- \theta) + \sin(\theta + \phi)}.
\end{equation}
Plugging in \( r, \theta, \phi \), we find
\begin{equation}\label{eqn:inscribedTriangleProblem:360}
\textrm{Area} = 17.1866.
\end{equation}
After computing this value, I then finally watched the original video to compare my answer, and was initially disturbed to find that this wasn’t even one of the possible values. However, that was because the problem itself, as originally stated, didn’t include the correct answer, and my worry that I’d made a mistake was unfounded, as the value I computed matched what was computed in the video (it also looks “about right” visually.)

Some 3D renderings of the Mandelbrot set.

February 7, 2021 math and physics play , , , , ,

As mentioned previously, using geometric algebra we can convert the iterative equation for the Mandelbrot set from complex number form
\begin{equation}
z \rightarrow z^2 + c,
\end{equation}
to an equivalent vector form
\begin{equation}
\mathbf{x} \rightarrow \mathbf{x} \mathbf{e} \mathbf{x} + \mathbf{c},
\end{equation}
where \( \mathbf{e} \) represents the x-axis (say). Geometrically, each iteration takes \( \mathbf{e} \) and reflects it about the direction of \( \mathbf{x} \), then scales that by \( \mathbf{x}^2 \) and adds \( \mathbf{c} \).

To get the usual 2D Mandelbrot set, one iterates with vectors that lie only in the x-y plane, but we can treat the Mandelbrot set as a 3D solid if we remove the x-y plane restriction.

Last time I animated slices of the 3D set, but after a whole lot of messing around I managed to save the data for all the interior points of the 3D set in netCDF format, and render the solid using Paraview. Paraview has tons of filters available, and experimenting with them is pretty time consuming, but here are some initial screenshots of the 3D Mandelbrot set:

It’s interesting that much of the characteristic detail of the Mandelbrot set is not visible in the 3D volume, but if we slice that volume, you can then you can see it again.  Here’s a slice taken close to the z=0 plane (but far enough that the “CN tower” portion of the set is not visible)

You can also see some of that detail if the opacity of the rendering is turned way down:

If you look carefully at the images above, you’ll see that the axis labels are wrong.  I think that I’ve screwed up one of the stride related parameters to my putVar call, and I end up with x+z transposed in the axes labels when visualized.

%d