FFT Parabolic Interpolation

Someone at the Numerix-DSP blog made the reasonable suggestion of interpolating the results of an FFT using a parabola. I did some experimenting with this method and found that it was not any worse than my other methods and it seemed less likely to turn out spurious results.

I've moved out the Maxima domain with my calculations and into Common Lisp, but to create the calculation for the parabolic interpolation, Maxima was of some help. First, we define a function for taking three points, and producing the parabola that is defined by that. 

PolySolve(A, B, C) := block([polyForm, eqn1, eqn2, eqn3], polyForm: a*x^2 + b*x + c = y, eqn1: subst([x=A[1],y=A[2]], polyForm), eqn2: subst([x=B[1],y=B[2]], polyForm), eqn3: subst([x=C[1],y=C[2]], polyForm), subst(solve([eqn1, eqn2, eqn3],[a,b,c]), polyForm) );

Then we will solve for it. We make the observation that the shape of the parabola is the same regardless of shifting it left or right along the x-axis (or index axis, which corresponds to the frequency axis when we multiply by sampling rate). This being the case, we have no need to supply the final indices but can consider our middle point as index 0, the point before it as index -1 and the point after it as index 1. This simplifies the calculation so that we can find the location of the peak relative to the index we are examining and apply the results of this calculation to it. Time to pull in some grade 12 math and differentiate that parabola and find the vertex (see Fig. 1).

Fig. 1. The green dot is approximately where the vertex of this parabola is. The result of this calculation would have x as slightly positive.

poly: PolySolve([-1,y0],[0,y1],[1,y2]); ratsimp(poly); df: diff(poly, x); freq: solve(df, [x]); \[\frac{{{x}^{2}}\, \left( \mathit{y2}-2 \mathit{y1}+\mathit{y0}\right) }{2}-\frac{x\, \left( \mathit{y0}-\mathit{y2}\right) }{2}+\mathit{y1}=y\] \[\frac{\left( {{x}^{2}}+x\right) \, \mathit{y2}+\left( 2-2 {{x}^{2}}\right) \, \mathit{y1}+\left( {{x}^{2}}-x\right) \, \mathit{y0}}{2}=y\] \[x\, \left( \mathit{y2}-2 \mathit{y1}+\mathit{y0}\right) -\frac{\mathit{y0}-\mathit{y2}}{2}=0\] \[[x=-\frac{\mathit{y2}-\mathit{y0}}{2 \mathit{y2}-4 \mathit{y1}+2 \mathit{y0}}]\]

On the basis of this we can produce a common lisp function to perform this calculation.

(defun parabola-peak-frequency (y0 y1 y2) (let* ((y2-y0 (- y0 y2)) (denom (+ (* 2 y2) (* -4 y1) (* 2 y0))) (x (/ y2-y0 denom)) (x^2 (* x x)) (amp (/ (+ (* (+ x^2 x) y2) (* 2 (- 1 x^2) y1) (* (- x^2 x) y0)) 2))) (list x amp)))

Note that if the peak is to the left of our high point, we will get a negative x value and if it is to the right, a positive value. It will be on the interval \((-1, 1)\). Some code that goes through the results of an FFT and takes the absolute value as in our previous work and finds local peaks is all we need now. Where there is a local peak of consecutive absolute values \(a\), \(b\), and \(c\) such that \(a, c < b\), do some interpolation using parabola-peak-frequency to get an interpolated amplitude and an index adjustment for the frequency index. Then the calculation for the the frequency becomes

\[(o + x) s/N\]

where o is the index of \(b\), x is the result from parabola-peak-frequency, s is the sample rate, and \(N \)is the length of elements in the FFT.

Now, maybe we can start applying some math to some real world problems, like, um, Mozart's Piano Concerto No. 13, movement II? Maybe just a snatch within the first second of it for now. 

Fig. 2. Mozart, anyone?

A general outline on how to proceed with this is to identify all of the local peaks, take the top K highest amplitude data points, interpolate those using the parabolic interpolation above, apply pitch correction (if desired) and midi note calculation. This doesn't address determining pitch durations. Another interesting problem is to try to guess which lower peaks are harmonics of lower pitches. These, and probably other confounding problems, yet await.

FFT - Local Symmetry Inferred Peak

As I was thinking about the previous post, I thought there might be a way to estimate the location of the peak. That is, to find the location in between the discrete data points where the peak probably occurs. 

I tried using realpart and imagpart and abs to see the differences and it seemed like realpart gives me the best view of the relative amplitudes when there are multiple frequencies involved. I also decided to apply an index shift since that seems to match the frequency better although I'm not satisfied with it technically yet.

Let's zoom in on the location of one of the peaks and see what it looks like:

wxplot2d([discrete, makelist((i-1)*sampleRate/N, i, 160, 170), map(lambda([u], abs(realpart(u))), firstn(rest(transformed, 159),11))]);

Fig. 1. Looks like we should be able to make a better guess than just picking one of the points.

Here are the actual points as frequency, absolute, real part value pairs:

[[427.972412109375,1.373307119161833],[430.6640625,1.783783923237471],[433.355712890625,2.525090138349324],[436.04736328125,4.273342049313825],[438.739013671875,13.47741003943344],[441.4306640625,11.94536977902466],[444.122314453125,4.166745114303215],[446.81396484375,2.532438127248608],[449.505615234375,1.822956914133206],[452.197265625,1.426082632315385],[454.888916015625,1.172306075942115]]

For a first crack at it, we might try linear interpolation on the two pairs of points on either side of the gap that contains the peak. Based on the Wikipedia article of the continuous version, we are certainly deviating from the shape of the real thing. Dauntless we press on to see what we get, because maybe we can live with a slight improvement that we know isn't perfect.

So, here we go do linear interpolation and find the line intersection using Maxima.

Line(A, B) := [x, y] . [A[2]-B[2], B[1]-A[1]] = [A[2]-B[2], B[1]-A[1]] . A; float(solve([ Line([436.04736328125,4.273342049313825],[438.739013671875,13.47741003943344]), Line([441.4306640625,11.94536977902466],[444.122314453125,4.166745114303215])], [x, y]));

[[x=439.729058165623,y=16.86285595968834]]

Let's put that point in the middle and see what it looks like:

revised: [[427.972412109375,1.373307119161833],[430.6640625,1.783783923237471],[433.355712890625,2.525090138349324], [436.04736328125,4.273342049313825],[438.739013671875,13.47741003943344], [439.729058165623,16.86285595968834], [441.4306640625,11.94536977902466],[444.122314453125,4.166745114303215], [446.81396484375,2.532438127248608],[449.505615234375,1.822956914133206],[452.197265625,1.426082632315385],[454.888916015625,1.172306075942115]]$ wxplot2d([discrete, revised]); 

Fig. 2. Hmm, closer, but looks kinda fake to me.

Somewhat predictably, this looks as fake as it really is. I think this graph makes it just how clear that linear interpolation is slightly bogus here. Not completely bogus though, it got us closer, right? So, we've ended up on the left side of the true peak which we know should happen exactly at 440 Hz. Why? The slope on the right side of the true peak is further away from the peak and therefore has a lower (absolute) slope value—it is less steep than it would be if it was closer to the peak. This is the weakness of linear interpolating this. We will end up closer to the higher of the two near-peak points than we should.

To get closer, we might want to try a different type of interpolation. Maybe a cubic spline?

FFT in Maxima

Maxima has a basic fast Fourier transform package called fft. Here's are really basic introduction to using it on music data. We define a basic sine wave form with frequency 440 Hz and amplitude 40 as our input and then run fft on it. The interpretation of the results as frequency depends on coordinating together the sample rate, the total number of samples (N), and the sample index (i).

wave(amp, freq, t) := amp*sin(2*float(%pi)*freq*t); sampleRate: 44100.0$ timeStep: 1/sampleRate$ N: 2^14; inputFFT: makelist(wave(40, 440, t*timeStep), t, 1, N)$ transformed: fft(inputFFT)$ wxplot2d([discrete, makelist(i*sampleRate/N, i, 1, 200), firstn(map(lambda([u], abs(u)), transformed), 200)]);

Fig 1. Looks like about 440 Hz

The thing that most threw me for a loop in trying to make sense of the various sources I read was that the x axis needs to be inferred in order to be comprehensible. I wasn't really clear on what that x-axis was. If you dig into the details of the output numbers, you will find the 164th value, has the maximum value. This implies 164 * 44100 / 2^14 = 441 Hz. So, a bit of error has crept into our interpretation of the result because of the discrete nature of the intermediate result. What this means in theory is that there is a theoretical peak that "should" exist somewhat before the 164th index, but because we aren't working with continuous data, we can't see that peak (at the ~163.47th index 😜--or maybe we should be considering it at the 164.47 "index" and subtract 1, interpreting it as zero based, since the peak is probably between the 164th and 165th index).

If we change the function to have two wave forms at different frequencies (340 Hz and 440 Hz), we get both of these in the result:

Fig. 2. wave(20, 340, t*timeStep) + wave(40, 440, t*timeStep)

Next step, in search of the true peak.

Packing Spheres into a Cube

In this post I compare three different alignment strategies for putting spheres into a cube. Our main variable will be \(n\), the ratio of cube side to sphere diameter. The question is, how many spheres can we get into the cube? We denote the number of spheres as \(S\).

Matrix Alignment

The simplest alignment of spheres to think about is a matrix alignment, or row, column, layer. If \(n=10\), then we can fit 1000 spheres into the cube this way. Or, in general, \(S=n^3\).

Fig. 1. Matrix alignment.

One way to summarize this alignment is that each sphere belongs to its own circumscribing cube and none of these circumscribing cubes overlap each other. For small values of \(n\), this arrangement is optimal.

Square Base Pyramid Alignment

An improvement to \(S\), for larger values of \(n\) is had by making the alignment of subsequent layers offset. Start with a layer in a square arrangement (\(n\) by \(n\)). If you look down on this first layer, it's awfully tempting to nestle the next layer of spheres in the depressions between the spheres. We will make better use of the interior space, although we won't make as good a use of the space near the edges. But if we can get enough extra layers in to make up for the losses near the outside, that would carry the day.

Fig. 2. A few red dots show where we're thinking of putting the next layer of spheres. We won't fit as many on this layer, but the second layer isn't as far removed from the first.

In the matrix alignment, consecutive layers were all a distance of 1 unit between layers (center to center). With this current alignment, we need to know the distance between layers. We need to consider a group of four spheres as a base and another sphere nestled on top and find the vertical distance from the center of the four base spheres to the center of the top sphere. (Implicitly, the flat surface the four base spheres are resting on constitute horizontal. Vertical is perpendicular to that.)

Fig. 3. Top view of a square base pyramid packing arrangement.

We can see from Fig. 3., that line segment ab is  of length \(\sqrt{2}\). We take a section parallel to this line segment, shown in Fig. 4.

Fig. 4. We recognize a 45°-45°-90° triangle, abc.

We can calculate the vertical distance between layers, line segment cd, as \(\sqrt{2}/2 \approx 0.7071\dots\). Note that the spheres in this second layer are in basically the same configuration as the bottom layer. We can see this will be the case because the top sphere in Fig. 3., has its outer circumference pass across the point where the lower spheres meet. So adjacent spheres in the second layer will be in contact. And, further, second layer spheres will follow the rectangular pattern of depressions in the bottom layer which are spaced apart the same as the centers of the spheres in the bottom layer. So we will end up with the same organization of spheres in the second layer as the bottom layer, just one fewer row and one fewer column.
We can now calculate the number of spheres that will fit in a cube using this alignment as

$$S = n^2 \lambda_1 + (n-1)^2 \lambda_2,$$

where \(\lambda_1\) is the number of full layers and \(\lambda_2\) is the number of "nestled" layers. (The third layer, which is indeed "nestled" into the second layer is also returned to the configuration of the first layer and so I refer to it as a full layer instead. The cycle repeats.) To find \(\lambda_1\) and \(\lambda_2\), we first find the total layers, \(\lambda\) as

$$\lambda = \left \lfloor{\frac{n-1}{\sqrt{2}/2}}\right \rfloor+1$$

The rationale for this formula is: find how many spaces between layers there are and compensate for the difference between layers and spaces. Imagine a layer at the very bottom and a layer at the very top (regardless of configuration). The layer at the bottom, we get for free and is the 1 term in the formula. (The top layer may or may not be tight with the top when we're done but hold it up there for now.) The distance from the center of the bottom layer to the center of the top layer is \(n-1\). The number of times we can fit our spacing interval into this number (completely), is the number of additional layers we can fit in. (Additional to the bottom layer, the top layer is the last of the layers we just found room for in the \(n-1\), or, technically, \(n-1/2\), er, whatever, do some examples.)

The break down between \(\lambda_1\) and \(\lambda_2\) will not be exactly even, with \(\lambda_1\) having one more if \(\lambda\) is odd. We can show this neatly using the ceiling operator 

$$\lambda_1 = \left \lceil{\lambda/2}\right \rceil$$

$$\lambda_2 = \lambda - \lambda_1$$

As early as \(n=4\) we get \(S=66\), which is better than the matrix arrangement by 2 spheres. I modeled this to make sure I didn't let an off by one error sneak into my formula (see Fig. 5.).

Fig. 5. Somewhat surprisingly, rearrangement of spheres (relative to matrix alignment) produces gains in volume coverage even at \(n=4\).

Here it is in Maxima:

SquareBasePacking(n) := block([lambda_1, lambda_2, lambda], lambda: floor((n-1)/(sqrt(2)/2))+1, lambda_1: ceiling(lambda/2), lambda_2: lambda - lambda_1, n^2 * lambda_1 + (n-1)^2 * lambda_2);

Tetrahedral Alignment

Another rearrangement is similar to a honeycomb. A honeycomb is an arrangement of hexagons that meet perfectly at their edges. Pictures of this are readily available online. This Wikipedia article suffices. Imagine putting a circle or a sphere in each of these honeycomb cells. The will produce a single layer which is more compact, although it will push the distance between layers apart. 

So, what do we need to know to make a formula for this one?

1. Find the distance between successive layers.
2. Formula for the number of layers.
3. Formula for the number of spheres in each layer.
1. This depends on knowing the range of layer configurations (we will only have two different layer configurations).
4. Some assembly required.

This is the same procedure we followed already for the square based pyramid alignment, but we are being more explicit about the steps.

Fig. 6. We need to find the distance of line segment ad.

We recognize we have an equilateral triangle and so we can identify a 30°-60°-90° triangle on each half of altitude ab. We have drawn all of the angle bisectors for this triangle which meet at a common point (commonly known property for triangles). We can recognize/calculate the length of bd as \(\frac{1}{2\sqrt{3}}\). We can then calculate the length of ad as \(\frac{\sqrt{3}}{2} - \frac{1}{2\sqrt{3}} = 1/\sqrt{3}\). From the Pythagorean theorem, we have a length for cd of \(\sqrt{2/3}\).
The number of spheres in the first layer will require us to use the distance for ab. We now need to find the number of rows in the bottom layer. We have
$$L_1=n \rho_1 + (n-1) \rho_2,$$
where \(\rho_1\) and \(\rho_2\) are the number of each kind of row in the layer and \(L_1\) is the number of spheres in the first layer. Since rows are spaced by \(\sqrt{3}/2\), we have a number of rows, \(\rho\), of
$$\rho = \left \lfloor{\frac{n-1}{\sqrt{3}/2}}\right \rfloor+1$$
and \(\rho_1\) and \(\rho_2\) can be broken out similar to our \(\lambda\)'s earlier.
$$\rho_1 = \left \lceil{\frac{\rho}{2}}\right \rceil$$
$$\rho_2 = \rho - \rho_1.$$

Fig. 7. We're looking for the distance of line segment cd which can be found using ac and ad.

Now, there is a little matter of how nicely (or not) the second layer will behave for us. It is noteworthy that we are not able to put a sphere into every depression. We skip entire rows of depressions. (The Wikipedia article on sphere packing referred to these locations as the "sixth sphere".) These skipped locations may play mind games on you as you imagine going to the third layer. Nevertheless, the third layer can be made identical to the first, which is what I do here.

Fig. 8. The first layer is blue. The second layer is red. Are we going to have trouble with you Mr. Hedral?

There's a couple of things that happen on the second layer. First, instead of starting with a row of size \(n\), we start with a row of size \(n-1\). The second issue is that we may not get as many total rows in. We could do the same formula for rows again but now accounting for the fact that we have lost an additional \(\frac{1}{2\sqrt{3}}\) before we start. However, we can reduce this to a question of "do we lose a row or not?" The total distance covered by the rows, \(R\), for the bottom layer is
$$R = (\rho-1) \frac{\sqrt{3}}{2} + 1$$
If \(n-R\ge\frac{1}{2\sqrt{3}}\), then we have the same number of total rows. Otherwise, we have one less row. We let \(r_L\) represent the number of rows lost in the second layer, which will either be 0 or 1. Noting that the order of the rows is now \(n-1\) followed by \(n\), we can give an expression for the number of spheres, \(L_2\), in our second layer now.
$$L_2 = (n-1)(\rho_1 - r_L) + n\rho_2$$
We have a very similar formula for the total number of spheres in the box.
$$S = L_1\lambda_1 + L_2\lambda_2,$$
where
$$\lambda_1 = \left \lceil{\frac{\lambda}{2}}\right \rceil$$
$$\lambda_2 = \lambda - \lambda_1,$$
where
$$\lambda = \left \lfloor{\frac{n-1}{\sqrt{2/3}}}\right \rfloor+1.$$
My Maxima function for this is
TetrahedralPacking(n) := block( [lambda_1, lambda_2, lambda, rho_1, rho_2, rho, r_L, L_1, L_2, R], rho: floor((n-1)/(sqrt(3)/2))+1, rho_1: ceiling(rho/2), rho_2: rho - rho_1, L_1: n*rho_1 + (n-1)*rho_2, R: (rho-1) * sqrt(3)/2 + 1, r_L: if n - R >= 1/(2*sqrt(3)) then 0 else 1, L_2: (n-1)*(rho_1 - r_L) + n*rho_2, lambda: floor((n-1)/sqrt(2/3))+1, lambda_1: ceiling(lambda/2), lambda_2: lambda - lambda_1, L_1 * lambda_1 + L_2 * lambda_2);

Comparison Between Methods

It isn't too hard to tell that one of the two second approaches is going to produce better results than the first after some value of \(n\) (Fig. 9.). What is more surprising is that our latter two alignments stay neck and neck up into very high values of \(n\). If we divide them by the volume of the containing cube, both appear to be a round about way to arrive at the square root of 2! (Fig. 10.)

Fig. 9. Both of the latter approaches give much better packing densities than the matrix alignment.

Fig. 10. Square base packing in blue, tetrahedral packing in red, both divided by total volume (# of spheres/unit volume). Unclear if there is a distinct winner and we seem to be getting closer to \(\sqrt{2}\).

Let's see which one wins more often for positive integer values of \(n\) in Maxima.

So, there's no run away winner, but if you're a betting man, bet on square base pyramid packing out of the available choices in this post. Regardless, it appears that both of these packing arrangements approach optimal packing (see Sphere packing). My density calculation (allowing the \(\sqrt{2}\) conjecture) for volume coverage comes out to

Gauss, old boy, seems to have approved of this number.

Tax Brackets Are Not So Bad

There's a common misconception that many seem to have about the way taxes work. People overthink themselves about "moving into a higher tax bracket." They think that if they just barely go over their current tax bracket into the next, the tax man will get extra hungry and they will wind up making less money than if they, for example, didn't go in for that little bit of overtime. These concerns are, in a word, unwarranted. That's not to say that anecdotal evidence in support of this idea does not abound. I'm sure you can find anecdotal evidence for anything that's popular to complain about.

In 2015, the following tax brackets applied to income taxes:

• <=$44,701, 15%
• $44,701—$89,401, 22%
• $89,401—$138,586, 26%
• >$138,586, 29%.

So, if I make $44,701 in taxable income, I pay $6705.15 in taxes on it. Suppose I make $44,702. Here's where the confusion becomes apparent. Some would try to tell you that you now pay 22% on the entire amount, or $9834.44, which would be a complaint worthy problem for someone teetering on the border of a tax bracket. (If a government actually tried to do it this way, the public outrage would produce headline news, not anecdotal evidence.) I call this false idea "folk" tax brackets.

What do I actually pay on $44,702? I pay \(15\% \times $44,701 + 22\% \times $1 = $6705.37\) in income tax.

The folk thinking taxation produces a jagged net income graph. This is the sort of thing people are implying when they speak of "making less" if they work overtime. You'd think they thought the take-home pay graph looked like this:

Fig. 1. Marginal taxes don't cause the take-home pay graph to look anything like this, but some talk as if it did.

The real graph is not jagged. It isn't smooth either (if you look close enough). It is continuous, unlike the above graph, but it does produce corners at the locations of the tax bracket cross-overs. These are not really that scary looking:

Fig. 2. Marginal taxes only increase the rate of tax "marginally". The increase in tax rate applies to the marginal amounts ("increments"), not the entire amounts. Basically you have a line with some very slight elbows in it where there is a small change in direction.

So, like you, I wouldn't mind if the whole graph was a little steeper (increased take-home pay percentage, less taxes), but it is what it is. And what it is is relatively fair, in concept at least.

Modeling Taxes in Maxima

Now, where did I get these wonderful graphs? I produced them in Maxima, of course. Here's the Maxima code to define a marginal tax formula, if given the brackets that apply:

MarginalTaxFormula(income_variable, percentList, thresholds) := block([rangeNumbers, ranges, i, u], rangeNumbers: append([0],thresholds,[inf]), ranges: makelist([percentList[i], rangeNumbers[i], rangeNumbers[i+1]],i,1,length(percentList)), lreduce("+", map(lambda([u], [p, mn, mx]: u, charfun(mn < income_variable and income_variable < mx) * (income - mn) * p + charfun(income_variable >= mx) * (mx - mn) * p), ranges))); FederalIncomeTax2015: [[0.15, 0.22, 0.26, 0.29], [44701, 89401, 138586]]; ManitobaIncomeTax2015: [[0.1080, 0.1275, 0.1740], [31000, 67000]]; marginalFederalTax(income) := ''(MarginalTaxFormula(income, FederalIncomeTax2015[1], FederalIncomeTax2015[2])); marginalManitobaTax(income) := ''(MarginalTaxFormula(income, ManitobaIncomeTax2015[1], ManitobaIncomeTax2015[2])); marginalIncomeTax(income) := marginalFederalTax(income) + marginalManitobaTax(income); afterTax(income) := income - marginalTax(income);

You can see the parts that you would need to change to produce your own formula for your marginal tax situation. FederalIncomeTax2015 is a list containing two lists. The first of these is the list of tax percentages. The next is the list of corresponding threshold values (0 is implicit as the beginning and the last percentage is assumed to apply to anything after the largest threshold amount).

Noteworthy items in the Maxima code above:

1. The characteristic function (charfun) produces the same effect as an if statement in the code. The characteristic function yields 1 if the condition is true, else 0. So, if an amount is not in the relevant range, it will not contribute. I have two different conditions to check for, but rather than branching, I just keep them in parallel. The advantage is simplicity. The disadvantage is that, in theory, there is more computation happening than really needs to happen. The simplicity carries a lot of weight in my book.
2. I append 0 and infinity to the list of boundaries for the tax brackets. Since Maxima recognizes inf and knows that any-number is less than inf, I have no case-taking to deal with the initial or final tax brackets.
3. map. When map will do the job, nothing does it more succinctly. Thanks map!
4. '' (quote-quote or double single-quote). Oh, boy. The basic meaning is, "evaluate what comes next and use the result of the evaluation instead of what's here." What this amounts to in this code is that the code that produces the function (the MarginalTaxFormula) is not re-run every time the function marginalFederalTax(income) is called with a different value of income. The MarginalTaxFormula(...) code is invoked by the '' operator and this returns the function which MarginalTaxFormual(...) produces. This is necessary because I am defining the function using the := operator which suspends evaluation until the function is invoked. Thus, you end up with a definition like this:

marginalFederalTax(income):=12788.1*charfun(income>=138586)+9834.0*charfun(income>=89401)+6705.15*charfun(income>=44701)+0.15*charfun(0

0.22*charfun(44701

Maxima in Common Lisp via ASDF

I have previously used Lisp code that ran in Maxima. This is fairly straightforward if you know Common Lisp. Basic instructions can be found in Chapter 37.1 of the Maxima help files (look for Contents link at the top of any of the help pages). But, how about the other way around—Maxima from Lisp side?

I recently looked into running Maxima inside Common Lisp using ASDF and found this. That gave me hope that it was possible and a little insight into what it was trying to do. I have sometimes wanted to create a user interface in some other environment (such as LispWorks) and call upon Maxima to do some computations for me.

To get yourself set up to access Maxima from Common Lisp you need:

1. The source, including the maxima.asd file.
2. Some idea where to put stuff and, perhaps, a few tweaks to the maxima.asd file.
3. Some file search variable set up to allow Maxima load statements in your common lisp code so you can access functionality found in the share packages or your own packages.

Tweaks

I put the source under my common-lisp directory so that a simple call to (asdf:load-system :maxima) would work.

I found that I needed to edit the maxima.asd file. I don't want to tell ASDF where to put the result or interfere with its "normal" operation. "ASDF, you're in charge!" is basically how I look at it. Accordingly, I commented out the output-files :around method. I suppose you could change the *binary-output-dir* here as well, if you would like to—perhaps to suit whichever Lisp you are compiling to.

Fig. 1. I'm pretty sure I don't need this code for my purposes, so I have it commented out.

Accessing Shared Packages

You might not realize how much stuff is in the share packages if you have been using the wxMaxima interface only. It appears that a number of widely used share packages are automatically loaded in wxMaxima, but maxima.asd only loads the core components. This is not a serious impediment—you can determine which packages to load from the Maxima help files.

But before you can use the Maxima load command, you need to set the file_search_maxima and file_search_lisp Maxima variables. I combined a call to asdf:load-system with the setting of these variables in a file under my common-lisp directory. I called this file load-max.lisp and it looks a bit like this (note how hairy it gets if you scroll right):

(asdf:load-system :maxima) ;; edit the paths in these lists to suit your own file system layout (setf maxima::$file_search_maxima '((MLIST SIMP) "C:/Users/Darren/maxima/$$$.{mac,mc}" "c:/path/common-lisp/maxima-5.37.3/share/$$$.{mac,mc}" "c:/path/common-lisp/maxima-5.37.3/share/{affine,algebra,algebra/charsets,algebra/solver,amatrix,bernstein,calculus,cobyla,cobyla/ex,cobyla/lisp,colnew,colnew/lisp,combinatorics,contrib,contrib/Eulix,contrib/Grobner,contrib/Zeilberger,contrib/alt-display,contrib/altsimp,contrib/binsplit,contrib/bitwise,contrib/boolsimp,contrib/coma,contrib/diffequations,contrib/diffequations/tests,contrib/elliptic_curves,contrib/elliptic_curves/figures,contrib/format,contrib/fresnel,contrib/gentran,contrib/gentran/man,contrib/gentran/test,contrib/gf,contrib/integration,contrib/levin,contrib/lurkmathml,contrib/maxima-odesolve,contrib/maximaMathML,contrib/mcclim,contrib/namespaces,contrib/noninteractive,contrib/odes,contrib/prim,contrib/rand,contrib/rkf45,contrib/sarag,contrib/smath,contrib/state,contrib/trigtools,contrib/unit,contrib/vector3d,descriptive,diff_form,diff_form/tests,diffequations,distrib,draw,dynamics,ezunits,finance,fourier_elim,fractals,graphs,hypergeometric,integequations,integer_sequence,integration,lapack,lapack/blas,lapack/lapack,lbfgs,linearalgebra,logic,lsquares,macro,matrix,minpack,minpack/lisp,misc,mnewton,multiadditive,numeric,numericalio,orthopoly,pdiff,physics,simplex,simplex/Tests,simplification,solve_rat_ineq,solve_rec,sound,stats,stringproc,sym,tensor,to_poly_solve,trigonometry,utils,vector,z_transform}/$$$.{mac,mc}" "c:/path-to-my-own/Useful Maxima Functions/$$$.{mac,mc}") maxima::$file_search_lisp `((MLIST SIMP) "C:/Users/Darren/maxima/$$$.{fasl,lisp,lsp}" "c:/path/common-lisp/maxima-5.37.3/share/$$$.{fasl,lisp,lsp}" "c:/path/common-lisp/maxima-5.37.3/share/{affine,algebra,algebra/charsets,algebra/solver,amatrix,bernstein,calculus,cobyla,cobyla/ex,cobyla/lisp,colnew,colnew/lisp,combinatorics,contrib,contrib/Eulix,contrib/Grobner,contrib/Zeilberger,contrib/alt-display,contrib/altsimp,contrib/binsplit,contrib/bitwise,contrib/boolsimp,contrib/coma,contrib/diffequations,contrib/diffequations/tests,contrib/elliptic_curves,contrib/elliptic_curves/figures,contrib/format,contrib/fresnel,contrib/gentran,contrib/gentran/man,contrib/gentran/test,contrib/gf,contrib/integration,contrib/levin,contrib/lurkmathml,contrib/maxima-odesolve,contrib/maximaMathML,contrib/mcclim,contrib/namespaces,contrib/noninteractive,contrib/odes,contrib/prim,contrib/rand,contrib/rkf45,contrib/sarag,contrib/smath,contrib/state,contrib/trigtools,contrib/unit,contrib/vector3d,descriptive,diff_form,diff_form/tests,diffequations,distrib,draw,dynamics,ezunits,finance,fourier_elim,fractals,graphs,hypergeometric,integequations,integer_sequence,integration,lapack,lapack/blas,lapack/lapack,lbfgs,linearalgebra,logic,lsquares,macro,matrix,minpack,minpack/lisp,misc,mnewton,multiadditive,numeric,numericalio,orthopoly,pdiff,physics,simplex,simplex/Tests,simplification,solve_rat_ineq,solve_rec,sound,stats,stringproc,sym,tensor,to_poly_solve,trigonometry,utils,vector,z_transform}/$$$.{fasl,lisp,lsp}" "c:/path/common-lisp/maxima-5.37.3/src/$$$.{fasl,lisp,lsp}" "c:/path-to-my-own/Useful Maxima Functions/$$$.{o,lisp,lsp}"))

So, if I want to use Maxima in my Lisp session, I just need to call (load "~/common-lisp/load-max.lisp"). I can move to the Maxima package by calling (in-package :maxima).

Try out the diff function:

MAXIMA> ($diff #$x^2$ #$x$)
((MTIMES SIMP) 2 $X)

Note the syntax for referring to the diff function ($diff) and the syntax for an entire maxima expression (#$expression$). There's also some lower to uppercase switching to watch out for, but see chapter 37 of the Maxima help documentation for more about this.

Suppose I want to use something in a share package or one of my own custom Maxima packages. I can use the Maxima version of load, namely, $load.

MAXIMA> ($load "stringproc")

MAXIMA> ($parse_string "x^2")
((MEXPT) $X 2)

One more example, to show how to deal with Maxima-specific features that are beyond the scope of Common Lisp. In general, use the meval function as in:

(setf n (* 3 4 5 6 7))
(maxima::meval `((maxima::$divisors) ,n))

((MAXIMA::$SET MAXIMA::SIMP) 1 2 3 4 5 6 7 8 9 10 12 14 15 18 20 21 24 28 30 35 36 40 42 45 56 60 63 70 72 84 90 105 120 126 140 168 180 210 252 280 315 360 420 504 630 840 1260 2520)

(If you're interested, see Maxima-discuss archives. HT Robert Dodier.)

Plotting a LWPOLYLINE in Maxima

In this post, I want to show how to produce a representation of LWPOLYLINEs in Maxima that is suitable for plotting to the plot2d function. The direct usefulness of such an operation is dubious. Think of this as a worked example to help build understanding of a few mathematical and computing tools.

Recap
The LWPOLYLINE is a 2d representation of polylines that represents each vertex as having x-value, y-value, and bulge value. The bulge value indicates the path of travel from the current point to the next one. Either it will go straight (bulge = 0), veer to the left and back as a CW arc (bulge < 0), or it will veer to the right and back as a CCW arc (bulge > 0).

So, it will

• bulge left (< 0),
• go straight (= 0), or
• bulge right (> 0).

More details on the mathematics of the bulge factor can be found here and here.

The Pieces

First of all, we define our problem and break its solution up into pieces. A polyline is of such a nature that is does not make a good "\(f(x)\)" type function (plotting \(y\) versus \(x\)). However, it lends itself to parametric expression well enough. So, we will have an \(x(s)\) function and a \(y(s)\) function. Ideally, I would like to specify a length along the polyline (measured from the start) and get a point \((x(s), y(s))\). (If you have a surveying background, think \(s = \) station number.)

So, I need to have parametric equations for lines and for arcs that are parameterized for arc length and initial value of the parameter. I then need to stitch these together to make both an \(x\) and a \(y\) piecewise function.

Linear

Parametric equations for lines involve a direction and a starting point. Simply normalize the direction vector to parameterize for length. This looks like

\[(x(s), y(s)) = A + \frac{(B-A)}{||B-A||}s,\]

where \(A\) is the starting point and \(B\) is the final point. We will need to adjust the initial value to account for all of the length of the polyline prior to the initial point on this segment, say, \(s_i\) which changes our equations to

\[(x(s), y(s)) = A + \frac{(B-A)}{||B-A||}(s-s_i).\]

Arc

Parametric equations for arcs use a center, radius, initial angle, direction of rotation, and total angle subtended (some variation of feasible input parameters admitted). Parameterizing for arc length isn't as obvious as it was for lines, so let's start with a basic circle:

\[x(s) = c_x + r \cos{s}\]

\[y(s) = c_y + r \sin{s},\]

where \((c_x, c_y)\) is the center and \(r\) is the radius of the circle. For our reparameterization to work we require that radian measure be used (or else require degrees and do the conversion in the function, but why make work for ourselves and introduce inefficiency?). Recall that the definition of radians is the ratio of arc length to radius (\(\theta = s/r\)). So, our parameter becomes \(s/r\) and so we have

\[x(s) = c_x + r \cos{\left(\frac{s}{r}\right)}\]

\[y(s) = c_y + r \sin{\left(\frac{s}{r}\right)}.\]

To get the right starting point for our function we need the parameter inside the trig functions to be equal to the correct initial angle when the input value of \(s_i\) (length of polyline so up to the start of the arc) is given. We need, then, the initial angle, say \(\theta_i\), plus the angle past initial that the given \(s\) corresponds to. Thus,

\[x(s) = c_x + r \cos{\left(\theta_i + \sigma \frac{s-s_i}{r}\right)}\]

\[y(s) = c_y + r \sin{\left(\theta_i + \sigma \frac{s-s_i}{r}\right)},\]

where \(\sigma\) is 1 or -1 as per the sign of the bulge factor.

Piecewise

Making a piecewise function parameterized for length requires that we keep track of the length of the polyline up to the segment we are about to address. Implementing it in Maxima also requires a bit of acquaintance with the methods available to us for doing so. It is possible to use if-then-else logic to achieve this and obtain something that works fine for plotting purposes. As an example of this, consider a really simple example:

f(x) := if x < 1 then x^2 else x^3;

plot2d(f(x),[x,-1,2]);

Fig 1. If-then-else logic in functions plots fine.

There's another, somewhat handier approach to this problem, namely, the characteristic function. The characteristic function function is a way of expressing conditions as multiplication by 1 or 0 in order to include and exclude the appropriate expression. charfun2(x,a,b) returns 1 if \(x\in [a,b)\), else 0. So, we can modify the above function to be as follows and obtain the same result. 

load(interpol);

f(x) := charfun2(x,-inf,1)*x^2 + charfun2(x,1,inf)*x^3;
plot2d(f(x),[x,-1,2]);

One nice thing about these expressions is that it is possible to create them separately and add them together later and get a relatively normal Maxima expression. (There are some hiccups if you try to differentiate or integrate them.) I can't comment on the performance of the characteristic function approach versus the if-then-else approach, but the characteristic function approach may be more manageable to produce.

The Program
There are a number of details that don't specifically relate to this problem that I will not go into detail about and other details I have given in the two previously mentioned posts on bulge factors. Here is the complete code for producing the parametric equations and plotting them followed by the result for a specific example polyline.

:lisp (defun |$MINLENGTHP| (maxima-list min-length)   (loop for p in maxima-list      for i = 1 then (1+ i)        if (> i min-length) return t))

load(interpol)$ maperror: false$ mapprint: false$ ratprint: false$ justpoint(u) := [u[1],u[2]]$ turn_right(direction) := [direction[2], -direction[1]]$ rotation_sign(dirA, dirB) := -signum(justpoint(dirA) . turn_right(dirB))$ deflection_sign(a,b,c) := rotation_sign(b-a,c-b)$ magnitude(v) := sqrt(v . v)$ polyseglength(a,b) := block([run,radius,bulge,centralAngle], run: magnitude(justpoint(b-a)), bulge: a[3], if bulge = 0 then run else ( radius: run/4*(1/bulge + bulge), centralAngle: atan(bulge)*4, centralAngle*radius))$ /* negative bulge does not cause an issue for this calc; comes out in the wash */ /* returns [center,radius,rise,run,altitude,bulge_normal] */ arc_everything(a,b) := block([dir,u,norm,bulge,s,rise,radius,alt], dir: justpoint(b-a), u: magnitude(dir), norm: turn_right(dir)/u, bulge: a[3], s: signum(bulge), rise: abs(bulge)*u/2, radius: u/4*(1/abs(bulge) + abs(bulge)), alt: radius - rise, [(a+b)/2 - s*alt*norm,radius,rise,u,alt,norm])$ parametric_line(A,B,start,parameter) := justpoint(A) + justpoint(B-A)/magnitude(justpoint(B-A)) * (parameter - start)$ parametric_arc(A,B,start,parameter) := block([center,radius,rise,run,altitude,bulge_normal,sig,theta_i,theta], [center,radius,rise,run,altitude,bulge_normal]: arc_everything(A,B), sig: signum(A[3]), [x,y]: (justpoint(A)-center)/magnitude(justpoint(A)-center), theta_i: atan2(y,x), theta: theta_i + sig * (parameter - start)/radius, center + radius*[cos(theta),sin(theta)])$ parametric_polyline(polyline,parameter) := block([cumlen,cumlen2,current,poly], cumlen: 0.0, current: polyline, poly: 0, while minlengthp(current,2) do ( term: if current[1][3] = 0 then parametric_line(current[1],current[2],cumlen,parameter) else parametric_arc(current[1],current[2],cumlen,parameter), cumlen2: cumlen + polyseglength(current[1],current[2]), poly: poly + charfun2(parameter,cumlen,cumlen2)*term, cumlen: cumlen2, current: rest(current)), [parametric,poly[1],poly[2],[parameter,0.0,cumlen-10^-8]])$ pl: [[0,0,0],[0.2805351111430803,0.2805351111430803,-0.1989123673796571], [0.4568168316928867,0.3535533905932737,0],[2.854253390593276,0.3535533905932737,-0.4142135623731028], [3.103553390593274,0.1042533905932738,0],[3.103553390593274,-3.0,-1.0],[0.0,-3.0,0.5], [1.5,-2.0,0]]$ p: parametric_polyline(pl,t)$ plot2d([p],[same_xy,true])$

Circle Fitting in Maxima

In previous posts, I've considered the problem of fitting a circle to 2D data points. Fitting problems such as this require:

1. An equation form to fit to. 
2. A function to quantify the error, which is to be minimized.

For a circle, the form of equation we want to fit to is
\[r^2=(x-a)^2 + (y-b)^2.\]
Maxima's lsquares package removes the burden of item number 2 for us. But for interest's sake, the error function we want to minimize can be given as

\[\mathrm{r}\left( a,b\right) =\frac{\sum_{i=1}^{n}\sqrt{{\left( {x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}}}{n}\]

\[\mathrm{SSE}\left( a,b\right) =\sum_{i=1}^{n}{\left( \mathrm{r}\left( a,b\right) -\sqrt{{\left( {x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}}\right) }^{2}.\]

Recalling that \(\sum_{i=1}^{n}(\bar{x}-x_i)^2 = \sum_{i=1}^{n}{x_i}^2 - (\sum_{i=1}^{n}{x_i})/n\), we see that
\[\mathrm{SSE}\left( a,b\right) =\sum_{i=1}^{n}\left({\left( {x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}\right)-\frac{{\left(\sum_{i=1}^{n}\sqrt{{\left({x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}}\right)}^{2}}{n}.\]

lsquares_estimates()

Given a set of points in matrix form (m), we need only to use the equation of a circle and indicate which are the data point variables and which are the "solve for" variables. In our problem,

lsquares_estimates(m, [x,y], (x-a)^2 + (y-b)^2 = r^2, [a,b,r]);

will suffice. The following function can be used to check the SSE.

sse(a, b, r, pts) := sum((r-sqrt((pts[i][1]-a)^2+(pts[i][2]-b)^2))^2,i,1,length(pts));

Testing for Reasonableness

It is reasonable to ask whether least squares (or other data fitting methods) will reliably produce good results with a particular quality of data for a particular application. This is true for any problem. What if the math works fine, but the likelihood of being able to "math filter" the  likely errors in the data to reach an accurate answer is not so high? In other words, how much inaccuracy in the collected data can you tolerate and still have a reasonable expectation of getting to the right answer? I'm not going to attempt to deal with probabilities (per se) or confidence intervals, but look at a simulation approach that may be able to give you an idea of the reasonableness of the circle fit based on an (intuitively) estimated amount of random error in your data.

Simulation of Circular Arc Data

We will consider the problem of data collected, using a total station, of a circular arc, with an "as-built" scenario particularly in view. As such, auto-generated sample points should simulate the expected behaviors of a rod/prism person stepping out regular sample locations from start to end of an ostensibly circular arc. This can be implemented as random normal perturbations of a circular arc.

One of the first pieces to a realistic perturbed circle is to be able to apply randomization that follows a normal distribution given an expected value and a standard deviation. Maxima has a function that behaves this way. random_normal() can be called with an expected value and standard deviation and optionally with a number of values to generate. Here is a sample use and output of the function, pictured using the histogram function:

load(distrib)$
load(draw)$

histogram(

    random_normal(50.0, 1.0, 1000),

    nclasses = 10,

    xrange = [40,60]);

Fig. 1. Normally distributed data

The standard deviation should be chosen such that roughly two-thirds of the time (≈68%), the sample value will be within one standard deviation of the mean; that is, on the interval \((\bar{x}-\sigma, \bar{x}+\sigma).\) 

Where the radius is concerned, we have the relatively predictable variations of the prism person's actions and the relatively unpredictable variations of the ostensible circular arc, which could turn out to not be an arc at all. In the absence of any sort of meaningful solution to this latter variation, we could possibly just bump up the standard deviation. But I'm not clear on whether this is a meaningful thing to do in assessing the acceptability of the data collection precision. The error should be "small" if the thing we measured conforms well to a circle. In order to check for whether this is "small" you need to produce a standard deviation from the SSE: \[\sigma_{SSE} = \sqrt{\frac{SSE}{n-1}}\] In the scenario we're dealing with here, someone would probably draw the computed arc in a CAD program and compare it to the measured data points and in an "arm-waving" sort of way decide whether it looks good or not.

Where the point sampling distances are concerned, we will assume that the rod person can step out distances (or angle sweeps) within some percentage of that intended. We will also make the pragmatic assumption that the prism person has decided a definite number of samples per given arc and will try to make them evenly spaced. The prism person will try to subdivide the remaining space at each sample point, and two thirds of the time, they will be within some percentage of the distance (or angle sweep) they intended.

Here is the perturbedArc() function with a sample usage:

perturbedArc(angle, center, radius, probable_radius_error, sample_spacing_variability, number_points) := block( [initAngle, termAngle, radii, expectedSweep, actualSweep, angles, remainingAngle, i], initAngle: random(float(2*%pi)), termAngle: float(initAngle + angle), remainingAngle: float(angle), angles: [initAngle, termAngle], radii: random_normal(radius, probable_radius_error, number_points), for i: (number_points-1) thru 2 step -1 do ( expectedSweep: remainingAngle/i, actualSweep: random_normal(expectedSweep, sample_spacing_variability*expectedSweep), remainingAngle: remainingAngle - actualSweep, angles: cons(termAngle - remainingAngle, angles) ), angles: sort(angles), map(lambda([r,theta],bfloat(center+r*[cos(theta),sin(theta)])), radii, angles) )$ pts: perturbedArc(1.5*%pi, [1000,1000], 3000, 12, 0.05, 20)$ wxplot2d([discrete,pts], [style,points], [same_xy,true]);

Fig. 2 - This is a perturbed circular arc. Doesn't look very perturbed does it? If we bumped up the standard deviation of the radius enough, we'd get something more clearly perturbed.

To produce the matrix form of the data I use a simple function I created (see here) called MakeMatrix(). Beyond that, it is a simple as

m: MakeMatrix(pts)$
lsquares_estimates(m, [x,y], (x-a)^2 + (y-b)^2 = r^2, [a,b,r]);
subst(%[1],sse(a,b,abs(r),pts));

What is the Sine of 18°?

To find the exact values of sine (and cosine) of a given angle, you need to use identities. This particular value is interesting in that it takes a fair bit of work to arrive at this value analytically, and yet the value is expressible as a fairly simple radical. The first thing to notice is that 90 = 5(18). You probably know the half angle identity, but not a fifth angle identity.  And for good reason—the fifth angle identity has multiple solutions; it's a fifth degree polynomial.

I did this calculation manually the first time through—some years ago. It goes a lot easier with a computer algebra system. Today, I'm going to use Maxima to solve the problem.

An identity which comes up in the development of the complex numbers, gives a convenient way of expressing \(\sin {n\theta}\) in terms of combinations of \(\sin \theta\) and \(\cos \theta\).

De Moivre's formula:
\[r^n (\cos \theta + i \sin \theta)^n = r^n (\cos {n\theta} + i \sin {n\theta}).\]
So, if we take \(r=1\) and \(n=5\), we can do a binomial expansion and easily isolate the \(sin \theta\) and \(\cos \theta\). In a Maxima cell, I type

eq: (cos(%theta) + %i*sin(%theta))^5 = cos(5*%theta) + %i*sin(5*%theta);
realpart(eq);
imagpart(eq);

After executing, I obtain the fifth angle identities (if you want to call them that).
\[\mathrm{cos}\left( 5\,\theta\right) =5\,\mathrm{cos}\left( \theta\right) \,{\mathrm{sin}\left( \theta\right) }^{4}-10\,{\mathrm{cos}\left( \theta\right) }^{3}\,{\mathrm{sin}\left( \theta\right) }^{2}+{\mathrm{cos}\left( \theta\right) }^{5}\]
\[\mathrm{sin}\left( 5\,\theta\right) ={\mathrm{sin}\left( \theta\right) }^{5}-10\,{\mathrm{cos}\left( \theta\right) }^{2}\,{\mathrm{sin}\left( \theta\right) }^{3}+5\,{\mathrm{cos}\left( \theta\right) }^{4}\,\mathrm{sin}\left( \theta\right) \]
To eliminate \(\cos^2 \theta\) from the sine identity we use

subst([cos(%theta)=sqrt(1-sin(%theta)^2)],%);
expand(%);

and obtain
\[\mathrm{sin}\left( 5\,\theta\right) =16\,{\mathrm{sin}\left( \theta\right) }^{5}-20\,{\mathrm{sin}\left( \theta\right) }^{3}+5\,\mathrm{sin}\left( \theta\right). \]
Taking \(\theta=18°\), using \(x\) to stand in for \(\sin 18°\), and rearranging, we have
\[0 = 16\,{x}^{5}-20\,{x}^{3}+5\,x - 1.\]
Factoring this one manually was made easier by recognizing \(x = 1\) as a root. After that I tried forming a square free polynomial and found that there was in fact a squared polynomial involved. In general, this helps in a couple ways (potentially):

1. Reduces the degree of the polynomial to be factored.
2. Sometimes produces fully reduced factors as a by-product of forming the square free polynomial.
1. This happens when each repeated factor is of different degree; that is, occurs a different number of times.

(Square free means no factors left that are raised to some power—each factor occurs exactly once.) Doing this manually involves taking the greatest common divisor between a polynomial and it's derivative. But today, let's let Maxima do it.

factor(0=16*x^5-20*x^3+5*x - 1);

\[0=\left( x-1\right) \,{\left( 4\,{x}^{2}+2\,x-1\right) }^{2}.\]

At this point, you really do have to be exceedingly lazy to not continue with the rest of the work manually. It's just a quadratic formula to be applied now (we discard \(x=1\) as an artifact). To use a computer algebra system for the rest is...well...okay, it is Saturday [at original publication], so...

solve(0=4*x^2+2*x-1,[x]);

\[[x=-\frac{\sqrt{5}+1}{4},x=\frac{\sqrt{5}-1}{4}]\]

We reject the negative solution (which is \(\sin {-54°}\) and, interestingly, is the negative of half the Fibonacci number) on the grounds that the range of the sine function on our interval is positive. Hence,

\[\sin 18° = \frac {\sqrt{5} - 1}{4}.\]

iMaxima and Slime in EMACS

First of all, I am a Windows user. Lots of (probably) very helpful documentation for how to do this is not directly applicable to my setup, because I am not a Linux user. Bummer. That being said, it worked out in the end.

If you just want to hurry up and try Common Lisp in a REPL, then try Lispbox.

I recently set out to do two things:

1. Get EMACS to run Slime (with any Common Lisp implementation)
2. Get EMACS to run iMaxima

Slime

Watch Installing Common Lisp, Emacs, Slime & Quicklisp: It's a demo of installing EMACS with SBCL, Quicklisp, and Slime.

• The primary hangup I was having until I saw this video is that I was trying to find a way to get EMACS to do something that I was supposed to do from within a running instance of a selected Common Lisp implementation. In other words, if you want to use CCL in EMACS using Slime, you should run CCL in a cmd window and use that REPL. Like many solutions, it's totally obvious when it is shown to you, which is why you may have trouble finding explicit statements about this.
• If you want to see a few of the basic environment variable modifications in slow motion (as it were), see here.

To start slime, I launch EMACS and then type ALT-X slime ENTER.

iMaxima

Visit the main site for iMaxima to see what it's about. It also has a section for installation instructions. Most of what you need to know is on that website. Below are details from my own setup and additional research that helped me get the system up and running on my computer.

After you have EMACS setup and have decided your HOME folder (see here), run EMACS to let it create a ".emacs.d" folder. Then create or edit the file init.el (el means EMACS Lisp, I believe). Here is my entire init.el (largely excepting comments), which includes commands for a REPL with SBCL and iMaxima (separately):

(load (expand-file-name "~/quicklisp/slime-helper.el")) ;; Replace "sbcl" with the path to your implementation (setq inferior-lisp-program "sbcl") (add-to-list 'load-path "C:/Program Files (x86)/Maxima-5.31.2/share/maxima/5.31.2/emacs") (load "setup-imaxima-imath.el") (autoload 'maxima-mode "maxima" "Maxima mode" t) (autoload 'imaxima "imaxima" "Frontend for maxima with Image support" t) (autoload 'maxima "maxima" "Maxima interaction" t) (autoload 'imath-mode "imath" "Imath mode for math formula input" t) (setq imaxima-use-maxima-mode-flag t) (add-to-list 'auto-mode-alist '("\\.ma[cx]" . maxima-mode))

Size of TeX output can be modified using (setq imaxima-fnt-size "xxxx"), where xxxx is one of small, normalsize, large, Large, LARGE, huge, or Huge. Also, by viewing the imaxima.el file, you can see a number of other options that are available to you.

To understand what's happening a bit better, it helps to know a few variables.

1. Go to your scratch buffer in EMACS.
2. Type inferior-lisp-program CTRL-j.
3. You will get the output "sbcl" if you used my init.el file.

I suggest the scratch buffer because it is important in debugging when you don't have a working setup of anything else. The scratch buffer is available. It runs ELISP, EMACS' own Lisp dialect.

Here are a few other variables you might benefit from looking at:

• load-path
• exec-path
• image-library-alist

I've bolded image-library-alist because it is really helpful to know what file you need to display png files. I snagged the files jpeg62.dll, libpng14-14.dll, and zlib1.dll from a couple of sources and placed them in "C:\Program Files (x86)\emacs-24.3\bin" which is where my "runemacs.exe" is. 

• Notably, I got libpng14-14.dll and zlib1.dll from "C:\Program Files (x86)\Inkscape". (You really should try Inkscape, the free SVG editor, anyways.)
• If your image-library-alist has other files in it, you either need to augment your image-library-alist (probably in your init.el file) or find one of the files already in image-library-alist and place it in the load path. You only need to do this for png files (as far as I know) in order to use imaxima.

A command you should try running (still in the scratch buffer) is

    (image-type-available-p 'png)

which will return T or NIL depending on whether the file type is available or not, respectively. If you get NIL, it means you don't have one of the files in image-library-alist in your load-path (or perhaps a dependee file is missing?).

Breqn

I had a lot of trouble with the display of the LaTeX even after I had a running imaxima session in EMACS. (Results of evaluations are usually displayed with LaTeX.) The problem stemmed from a missing package: breqn. I installed the mh package from within the MiKTeX Package Manger and it appeared to install fine, but my LaTeX was still not working properly within an imaxima session. The issue appears to have been related to dynamic loading of packages.

Here's what finally shook the problem loose: within TeXnicCenter I started a new project and included the line

    \usepackage{breqn}

along with the other \usepackage declarations that were part of that project. When I built the project I was prompted for the installation of the package and the installation appeared to go smoothly. Then I tried imaxima again and the LaTeX output worked.

To start an imaxima session, launch EMACS and then type ALT-x imaxima ENTER.

Piece of cake! :)

Additional Info (probably not necessary)

I'm still getting a warning message when I load imaxima, although it has yet to prove to be a problem. The warning is:
   (lambda (x) ...) quoted with ' rather than with #'

I also made another change, and I don't know if it was necessary or not (hopefully not), but while wrestling with something quite different I changed line 14 of maxima-build.lisp (that's "C:\Program Files (x86)\Maxima-5.31.2\share\maxima\5.31.2\src\maxima-build.lisp" on my system). I changed that line to:

(load "C:/Program Files (x86)/Maxima-5.31.2/share/maxima/5.31.2/share/lisp-utils/defsystem.lisp")

This is just an explicit file reference rather than a relative path. (In case you're wondering, I was trying, unsuccessfully, to compile maxima myself. Perhaps another time.)

Skew Cut Profiles

I recently was surprised (anew) at what happens when you cut a prism which has a uniform cross-section with a plane which is not orthogonal to the length of the prism. Especially if the plane must be described in terms of two angles. A cross-section through a square prism is (by definition) a square. If you make your cut at an angle, you will get a rectangle if you have a simple angle cut. But if you make your cut not only a miter, but also a bevel you end up with a parallelogram.

The simplest way to visualize this is to imagine a plane that intercepts the prism at a single edge first and has its steepest slope in the direction of the opposite edge:

The resulting section when viewed in the intersecting plane would be something like this:

where the blue square is the cross-section, the red parallelogram is the intersection of the above described plane through the square prism when viewed in that plane. (A few thought experiments of your own will probably serve you better than any further illustration of this phenomenon I could drum up at the moment.)

It is interesting to see the result in an I-beam shape if you use a 45° miter and 45° bevel:



As a matter of where this observation fits into mathematics generally, this is a projective or affine transformation. Such transformations preserve straight lines but they do not, in general, preserve right angles.

Aside: Convex Hulls

The fact that they preserve straight lines is a boon to the would-be-writer of a convex hull algorithm in an arbitrary plane since you can project the 3D coordinates onto a simple 2D plane and use the 2D algorithm (being careful to preserve or restore the other coordinate at the end). If you do so, you should probably choose which set of 2 coordinates to use for the algorithm, based on the direction of the plane. For example, if the normal of the plane has a larger z-coordinate than x- or y-coordinates, then the (x,y) pairs are best to use for the information to run the 2D convex hull algorithm on. There are two reasons to be picky about this:

1. If the points fed to the routine are in thy yz-plane, you will get incorrect results (if you assumed the (x,y) coordinates should always be used).
2. It may reduce the effect of rounding errors in the calculations (in terms of deciding whether points "near the line" of the hull are inside or part of the hull).

(I have lost the original code I used for this post, so I omit references to it now.)

See the Maxima project for download information for Maxima.

Round 'em up!

Over the past several months I have been dabbling in Common Lisp (LISt Processing language).  Lisp, I think, is the ugly duckling of the programming languages. The facetious "Lost In Stupid Parentheses" seems more apt as an acronym.

My first encounter with this ungainly language (Lisp) was in AutoLisp (and Visual Lisp). I was on a hiatus from programming regularly and discovered this immediately discouraging language that was the "easiest" way to interact with AutoCAD programmatically. (I don't consider AutoCAD scripts to be programming, though I would still include them under the more general category automation.) I did not make much headway at that time and it frankly amazes me that any non-programmer will step up to the challenge and learn it.

But then there was Maxima. I was happy with Maxima as a free computer algebra system. But in the documentation I found the spectre of Lisp. Specifically, Common Lisp. The computations in Maxima are all in Lisp (though there is probably some non-Lisp they use to glue together a Windows interface for Maxima). Maxima itself supports a fairly robust set of programming constructs as far as doing calculations are concerned.

• Input/output from files
• Controlled loops
• if-then-else
• functions
• plotting (via gnuplot—included in standard installs)
• formatted output
• nice math notation display
• export to latex or image

Maxima is built on top of Lisp and you really have all the abilities of Lisp in Maxima (if you know how to get at them). But, being as I am, always in search of a better way of doing something, I am also concerned with whether a technology built on top of another is "hiding" some of the features of the underlying technology. Perhaps hiding isn't exactly the right way to put it, but learning Lisp is a paradigm shift from most other (Algol-like) programming languages I was familiar with, while on the other hand, learning Maxima didn't take me through that paradigm shift.

The two most interesting features of Lisp to me are the emphasis on function references and macros. But they both help you to do the same thing: reduce the repetition in your code. Mind you, they serve this end in very different ways.

When I wanted to implement Jarvis's march [1, p. 955] to encircle a list of points in a "convex hull" it was particularly the function referencing that made the implementation of the code fairly simple.

Fig 1. The red lines join together the red + signs which indicate the convex hull for these points.

The general idea of Jarvis's march is you find an outside point to start from and then keep picking the next point so that you never have to switch directions as you walk around. I started at the left most point. In the event of a tie, I choose the lower of these. (That would place us at about (0.1, 15.0).) I decided to walk around the points in a counter-clockwise direction. To decide which point comes next, I need to find the point which requires the right-most turn. So I begin by choosing the first of the remaining points (in whatever order I received them). I go through the remaining points and test whether I would need to make a right turn or left turn to point toward them (from my vantage point on the starting point of the hull chosen at my first step). If it is a right turn, then I have a new winner which replaces the first and I continue moving forward to the points not yet checked. The algorithm is really a repeated application of a maximum function where maximum is defined according to a different function than normal. What I need are predicate functions which define the rules for choosing a winner (a "maximum") and a routine to separate the winner from the remaining points that can accept a predicate function as a parameter.  Here's the code for, remove-winner:

(defun remove-winner (test li) (let ((mx (car li)) (temp nil) (therest nil)) (setf therest (loop for i in (cdr li) collect (if (funcall test i mx) (progn (setf temp mx mx i) temp) i))) (list mx therest)))
The idea of remove-winner is to use the test function (supplied by the caller) to determine whether the next item is "better" than the current winner. If you passed in the function < (using the syntax #'<) you would get the smallest value in the list (provided it was a list of numbers). Notice the thought process in this function is not just to find the max and return it, but rather to produce a list which includes the "winner" and the rest of the values (that didn't win) separately. We don't just need to know the winner, we need to separate the winner from the rest.

When we get to the convex-hull routine, we can use the same routine (remove-winner) using two different tests: one to find the initial element and another to find the successive elements.

(defun convex-hull (points) (cond ((< (length points) 3) (list nil points)) ((= (length points) 3) (list points nil)) (t (do* ((init (remove-winner (lambda (a b) (if (= (aref a 0) (aref b 0)) (< (aref a 1) (aref b 1)) ;; pick the lower (< (aref a 0) (aref b 0)))) points)) ;; pick the lefter (start (car init)) (rest (cadr init)) (next nil) (hull (list start)) (result nil)) ((or (eq start next) (not rest)) (list hull rest)) (setf result (remove-winner (lambda (a b) (right-turn-from (vect- a (car hull)) (vect- b (car hull)))) rest)) (setf next (car result)) (if (right-turn-from (vect- start (car hull)) (vect- next (car hull))) (setf next start) ;; set terminating condition (progn ;; else (setf hull (cons next hull) rest (cadr result))))))))

The first occasion we use remove-winner, the function is of the same sort we would use to perform a two key sort. In the main body of the routine, we have a more interesting function which accepts two points, calculates their direction vectors (using the current head node of the hull as the origin) and decides whether the second is a right turn away from the first, in which case it declares the second to be the (current) winner.

Can you do this [abstract the predicate] in other languages? Yes. In many other languages you could implement the routine substantially like this. But the key is this: would you ever think to do it this way if you were using VB.NET or C# or C++? Maybe. But not because of your experience with those languages, but because of your experiences with Lisp or F# or other languages which feature a functional style of programming.

Fig. 1 was produced from running a small bit of Maxima code which loads some routines from a separate package written in Lisp (including the above routines).  You can download both files to try them out:

• points.lisp
• circlethempoints.wxm

You will need to have Maxima installed to try this example out.

Sources:

1. Cormen, T.H., Leiserson, C.E., Rivest, R.L., & Stein, C., Introduction to Algorithms, Second Ed., McGraw-Hill Book Company / The MIT Press, 2001.