Christian Reinsch and Roland Bulirsch both passed away recently, Reinsch on October 8 and Bulirsch on September 21. Reinsch was 88 years old and Bulirsch was 89. Both of them were retired professors of numerical analysis at the Technical University of Munich. Both of them were friends of mine. But in almost all other ways, they were very different people.... read more >>

]]>Christian Reinsch and Roland Bulirsch both passed away recently, Reinsch on October 8 and Bulirsch on September 21. Reinsch was 88 years old and Bulirsch was 89. Both of them were retired professors of numerical analysis at the Technical University of Munich. Both of them were friends of mine. But in almost all other ways, they were very different people.

The *Handbook for Automatic Computation, Volume II, Linear Algebra*, is a research monograph published in 1971 by Springer-Verlag. This *Handbook* was edited by J. H. Wilkinson and Christian Reinsch and includes dozens of Algol procedures by 19 different authors for solving systems of simultaneous equations and computing matrix eigenvalues and singular values.

Wilkinson and colleagues at the National Physical Laboratory in Teddington, England wrote around half of the Algol procedures in the *Handbook*. Reinsch authored several procedures himself and reviewed and tested most, if not all, of the entire *Handbook*.

Translations into Fortran of many of the Algol codes produced the EISPACK subroutine library that led to the first MATLAB. I think it is fair to say that without the work of Jim Wilkinson and Christian Reinsch on the *Handbook* there might never have been a MATLAB.

I would have loved to listen in to discussions between Jim and Christian during the development of the *Handbook*. Here were two very talented, proud individuals working together at the detailed level required of scientific computer programing. It must have been exciting.

One small example: what subscripts do you use for the `n-1` off-diagonal elements of the symmetric, tridiagonal matrix of order `n` produced in one procedure and passed to another? Is it `e(1:n-1)` or `e(2:n)`? (A "modern" `e(0:n-2)` is not a possibility.) The first statement of the Algol procedures for the QR algorithms is

for i := 2 step 1 until n do e[i-1] := e[i];

Was there some kind of disconnect between Teddington and Munich?

Christian Reinsch was born in Chemnitz, Germany, in 1934. He spent most of his career at the Technical University Munich, where he was a close associate of Fritz Bauer.

Christian was an intensely private person. He never married and lived alone. He visited Argonne Laboratory in the 1970s, when we were working on EISPACK, and he attended a few Gatlinburg/Householder meetings. Other than that, he rarely travelled far from Munich.

This photo, taken by his colleague Christoph Zenger, is the only photo of Reinsch that I have ever seen. And, as far as I know, it will be the first photo of him available on the Internet.

Photo credit: Christoph Zenger

Gene Golub is known as Professor SVD because he did more than anyone else to develop the algorithms for computing the decomposition and for popularizing its applications. A 1965 paper by Golub and Velvel Kahan provides the first practical method for computing the SVD. The Golub/Kahan approach is based on the eigenvalues of the *2n-by-2n* block matrix `[0 A; A' 0]` .

Peter Businger was a grad student at Stanford in the 1960's who worked with Golub on several projects. A 1967 Stanford technical report by Golub and Businger includes an Algol procedure, written by Businger, for computing the SVD that is based on bidiagonalizing `A` itself. Both Golub/Kahan and Golub/Businger suggest Sturm sequences and Givens rotations for ultimately computing the singular values.

At the same time, Christian Reinsch independently developed his own method for computing the SVD. His review work on the *Handbook* gave Reinsch access to the implicit tridiagonal QR techniques for matrix eigenvalues that Francis and Wilkinson were investigating. So, Reinsch adapted implicit QR with Wilkinson shifts to the singular value situation.

Both Gene Golub and Christian Reinsch offered SVD contributions for the *Handbook* . It turns out that, with exact arithmetic and the same shifts, Golub/Businger and Reinsch would produce the same results. But Golub and Businger never used Wilkinson shifts and Golub did not join the QR club until later. Fortunately, Fritz Bauer, editor in chief of the *Handbook*, brokered a joint authorship arrangement and the Golub-Reinsch algorithm for computing the SVD resulted.

Walter Gander has investigated the history of SVD algorithms. The slides for the talk he presented at a workshop in Lanzhou University are available at this link.

*Einfuhrung in die Numerische Mathematik*, *Introduction to Numerical Analysis*, is a classic textbook by Josef Stoer and Roland Bulirsch. The original German editions were published by Springer-Verlag in 1972 and 1976. The English translations were published in 1980 and 1993. A complete PDF of the second English edition is available at this link.

I think of Stoer & Bulirsch as one of the best theoretical textbooks in numerical analysis. It is comparable to Isaacson & Keller. There are theorems and proofs. There are algorithms, but no software. There are a few, but not many, numerical examples. There are many excellent exercises.

Roland Bulirsch was born in Liberec, in the former Czechoslovkia, in 1932. After visiting U. C. San Diego in the 1960's, he spent most of his career at the Technical University Munich. In contrast to Reinsch, Bulirsch was a very gregarious, public person. He adored his grand children and they adored him. The three-part photo collection assembled for his 80th birthday has hundreds of snapshots and portraits.

Roland was an avid body-builder. He had the broadest shoulders that I have ever seen -- and they were twice as broad as his waist. He had huge hands and huge biceps. One of my favorite stories about Roland explains the photo on his desk that was signed,

"Thanks for everything -- Arnie"

It turns out that Roland and a body-builder from Austria named Arnold Swartzenegger trained together in Munich in the 1960's. Roland and some friends from the gym took up a collection to help the ambitious young Arnie emigrate to America.

Here are a few photos from Roland's albums.

My last visit to T. U. Munich was in 2015. I gave my stump speech, The Evolution of MATLAB, in a big, new lecture hall named after Friedrich Bauer. Christian Reinsch was in the audience (skip to time stamp 13:59 in the video). He was over 70 years old at the time and had just finished his daily 50-kilometer bicycle ride.

After the talk, a few of us found our way to one of the famous Munich *Biergarten*. Christian came along but did not stay for dinner. He was not a beer drinker.

Roland Bulirsch was able to join us for dinner, however, and we managed to polish off a few of those one-liter German beer steins before closing down the place.

Get
the MATLAB code

Published with MATLAB® R2022a

A 231-character static version of colorcubes in the current MiniHack generated quite a bit of interest. So I am posting the full version, where you can change number of cubes and the view point, at this link.... read more >>

]]>A 231-character static version of colorcubes in the current MiniHack generated quite a bit of interest. So I am posting the full version, where you can change number of cubes and the view point, at this link.

What is the color of the cubelet, or cubelets, in the center?

Get
the MATLAB code

Published with MATLAB® R2022a

This nifty graphics gem started with a contribution by Paul Villain to the MATLAB 2022 Mini Hack, currently taking place on MATLAB Central.... read more >>

]]>This nifty graphics gem started with a contribution by Paul Villain to the MATLAB 2022 Mini Hack, currently taking place on MATLAB Central.

Villain's contribution is 102 mod 500 . My rewrite is `modfun`. Villain's `102` and `500` become the parameters `m` and `n`.

modfun(m,n) connects n points, z(j), equally spaced around the complex unit circle, by n+1 straight lines. The j-th line connects z(j+1) to z(mod(j*m,n)+1).

The basic code uses complex arithmetic and is only eight lines long. When the graphics is done with `line` instead of `plot`, it is not necessary to use `hold on`.

function modfun(m,n) init_fig z = exp(2i*pi*(0:n)/n); for j = 0:n zj = [z(j+1),z(mod(j*m,n)+1)]; line(real(zj),imag(zj)) end end

The initialization makes `line` possible.

function init_fig axis([-1 1 -1 1]) axis square axis off end

This animation of `modfun(105,200)` has one frame for every five lines.

A sample.

Match these calls to `modfun` to the plots in the Gallery.

modfun(88,179) modfun(89,220) modfun(99,200) modfun(101,200) modfun(111,200) modfun(113,188) modfun(126,188) modfun(126,200)

An interactive `modfun`.

Get
the MATLAB code

Published with MATLAB® R2022a

"Clever Toys" is a puzzle company in the Czech Republic. Their Web site describes five different hand-made, wooden puzzles that are related mathematically to the Rubik's Cube.... read more >>

]]>"Clever Toys" is a puzzle company in the Czech Republic. Their Web site describes five different hand-made, wooden puzzles that are related mathematically to the Rubik's Cube.

Here is the photo of the puzzle "Trio".

Here is an English translation of their description.

The task of this 2D puzzle is to assemble all the ovals according to their color into a basic assembly. By moving the balls in the grooves and turning the center wheel, you can gradually get each ball where you need it.

I am not sure that Clever Toys is still in business. They have not responded to my emails and I have not been able to purchase an actual Trio puzzle. Of course, that is all the motivation I need to make a MATLAB model.

Here is the initial configuration. There are three fixed outer, partial, discs and one inner, full, disc. Each of the outer discs has a channel containing 10 marbles. When the inner disc is positioned properly, the marbles in a channel can be rotated. Rotating the inner disc itself moves some, but not all, of the marbles, thereby scrambling the colors.

Like Rubik's Cube, the objective of the puzzle is to return to this initial "solved" state.

Mathematically, both Rubik's Cube and Trio are ultimately based upon rotation matrices. Rubik's Cube is a 3-D puzzle whose state is specified by the position and orientation of 27 cubelets; this gives Rubik's Cube about `4.3*10^19` possible positions. At each step, there are six faces and 12 possible rotations.

Trio is a 2-D puzzle whose state is specified by the colors of the marbles. There are 10 marbles of each of three colors; this implies that Trio has `30!/(10!)^3 = 5.6*10^12` possible positions. At each step, there are eight possible rotations, four discs, clockwise or counter-clockwise.

A "scramble" is an integer vector with elements between -4 and +4 that specify moves or rotations. Move 0 initializes the model. Move `d` with `d` equal to 1, 2, or 3 rotates all the marbles in the `d`-th channel counter-clockwise for one-tenth of a full rotation. Move `d` with `d` equal to -1, -2, or -3 is the corresponding clockwise rotation. Moves -4 and 4 turn the central disc clockwise or counter-clockwise for one-third of a full rotation. This rotates some, but not all, of the marbles and mixes the colors.

Here is a scramble of length 29 that provides our example.

D = [ 4 1 3 1 1 -2 3 1 4 4 2 4 4 -3 -3 -2 -3 -3 2 1 4 -1 3 -4 2]

And here is the scrambled result.

This animated gif shows the scrambling process one move at a time. The animation does not repeat automatically, so to start it over again, refresh your browser. If it still doesn't move, find another browser.

Reverse the scramble by running it backwards, changing the sign of each move. This will return the scrambled position to the initial position. I call this "unscramble"; it solves the scrambled position by a "follow the breadcrumbs" algorithm.

Both animations take a long time to run -- about 80 seconds with my browser.

I don't have any idea about how to actually solve a given position without using knowledge of how it was generated, and I don't have any idea about how to quantify the difficulty of finding a solution. This is in sharp contrast to Rubik's Cube where there are measures of difficulty and algorithms for finding optimum solutions.

The program available at this link is interactive. Click or alt-click in any one of the four discs to make a move in that disc.

Thanks to Steve Eddins and Tom Lane for help with this post.

Get
the MATLAB code

Published with MATLAB® R2022b

- How hard is it to solve a Rubik's Cube?
- When can you say that your Rubik's Cube is completely scrambled?
- Why might the answer depend on where you went to school?
- What interesting mathematical questions are involved?

- How hard is it to solve a Rubik's Cube?
- When can you say that your Rubik's Cube is completely scrambled?
- Why might the answer depend on where you went to school?
- What interesting mathematical questions are involved?

The difficulty of solving any configuration of a Rubik's Cube is the smallest number of moves required to return to the intial configuration where each face is showing a single color.

But what, exactly, is a move? The so-called *quarter-turn metric* says a move is turning any face by 90 degrees. The *half-turn metric* says turning any face by either 90 or 180 degrees is a single move. For example, using Singmaster notation and the quarter-turn metric, the sequence "L L", which turns the left face twice in the same direction, is two moves. But in the half-move metric, the sequence becomes "L2" and counts as a single move.

Tomas Rokicki, who describes himself as a programmer from Palo Alto, provides some history at cube20.org.

In the early days of cube mathematics, two camps emerged on how to measure the difficulty of a position. West coast and Stanford mathematicians, free thinkers all, tended to prefer the half-turn metric, where any twist of any face, whether 90 degrees, 180 degrees, or 270 degrees counted as a single move. The east coast crowd, including MIT, tended to prefer the rigor of the quarter-turn metric, where a half-turn counted as two moves, since of course it could be accomplished by two consecutive quarter turns.

When I began development of a Rubik's Cube simulator, `Qube`, I was unaware of this history and, even though I am a devout West-coaster, I just counted quarter-turns. Now a toggle switch in `Qube` allows use of either metric.

Let `Q` denote a cube position,

| `Q` | = minimum number of moves to solve `Q`,

` Q` = the set of all possible

`G( Q)` = maximum over

`G( Q)` is known as "God's number".

The *superflip* of Rubik's cube is a configuration where the 8 corners, the 6 face centers, and the cube center show the initial colors, but the 12 edge cubelets have the colors reversed. In 1995, Michael Reid proved that solution of the superflip requires 20 half-turn metric moves. In 2010, Tomas Rokicki and colleagues, using hundreds of computers at Google, carried out a massive computation to prove that no other configuration took more than 20 moves, cube20.org. This established that God's number for the half-turn metric is

`G( Q)` = 20

I use `Q20` to denote the superflip. Our first animation generates the superflip with 20 moves. A few rotations at the beginning and at the end are shown in more detail so that we can see the rotation matrices involved. A higher resolution video clip is available at this link: Q20.mp4

The second animation shows a solution of `Q20` in 20 moves obtained by reversing and complementing the generating moves. Reid's proof shows that any other solution requires at least 20 moves. The high resolution clip is: Q20solve.mp4

There are several other configurations that require 20 moves. Any configuration with `G(Q) = 20` can be regarded as completely shuffled in the half-turn metric.

For the quarter-turn metric, if you combine superflip and a configuration known as *fourspot* you have `Q26` . Only 8 corners and two face enters are correct. The edges, four face centers, and the cube center are all reversed. When 180 degree turns are counted as two 90 degree turns, this configuration is generated by 26 moves and solved by reversing and complementing the 26 moves. The high resolution clip is: Q26.mp4

In 2014, a massive computation at the Ohio Supercomputer Center by Rokicki and Morley Davidson proved that only `Q26` (and its two rotations) required 26 quarter-turn moves All other configurations need fewer. cube20.org. So, this established that God's number for the half-turn metric is

`G( Q)` = 26

The high resolution clip is: Q26solve.mp4

Let's compare `Q20` and `Q26` by alternating between the two.

The type 3 corner pieces are the same in `Q20` and `Q26`, and are in the correct initial position.

The type 2 edge pieces are also the same in `Q20` and `Q26`, but are reversed from their initial position.

All the action is with cubelets of type 0 and type 1. In a real, physical Rubik's Cube, this is one solid piece that holds the Cube together.

Get
the MATLAB code

Published with MATLAB® R2022b

Until recently, I knew nothing about the `polyshape` object in MATLAB. Now I can use polyshape to simulate an extraordinary puzzle.... read more >>

Until recently, I knew nothing about the `polyshape` object in MATLAB. Now I can use polyshape to simulate an extraordinary puzzle.

The Web site Art of Play, from some folks in San Diego, has well over two-hundred puzzles. Recently, one of their most popular puzzles, Mighty Cheese, caught my attention. The idea is to move the slices of plastic cheese around within the frame in order to create a hole large enough to hold the little plastic mouse. However, when I last checked, Mighty Cheese was sold out. Not wanting to be dissuaded by the unavailability of the puzzle itself, and not knowing the solution, I set out to build a simulator.

I expected Mighty Cheese to be something like the T-puzzle that I enjoyed several years ago. But, unlike the T, the geometry of cheese cannot be modelled by simple rectangular patches. I needed help. Steve Eddins answered my plea and told me about MATLAB's "polyshapes". That led to this blog post.

A "selfie" of a cheese slice provides an example of a `polyshape`. It has curved boundaries and a hole. The MATLAB documentation for `polyshape` says this a *polygon*. In my mathematical world, *polygons* can't have curved boundaries or holes. But, I guess that I am being pedantic.

I began with this photo of the Mighty Cheese puzzle.

```
puzzle = imread('Cheese_puzzle.png');
```

Two lines of code and the `L*a*b` color model find the regions in `photo` that look like cheese.

[L,a,b] = imsplit(rgb2lab(puzzle)); mask = a > 30; spy(mask)

Steve provided a function that turns `mask` into a `polyshape` object.

slice = my_polyshapes(mask);

An overloaded `plot` function then generates the selfie from `slice`.

sliceplot = plot(slice, ... 'facecolor',cheese_yellow, ... 'facealpha',1);

I didn't use AI or simulated annealing or any other modern technique to search for a solution. I just poked around while I was learning about polyshapes and programming my simulator. It was just me and a *real* mouse. No, not a *real real* mouse. You know the kind I mean.

I almost made a blog post claiming this is the solution the designers of Mighty Cheese must have had in mind. But I wasn't sure. My solution made use of the white space around the outside edges of the slices. How wide are those gaps? I didn't yet have a real physical puzzle to measure tolerances.

I am sure glad now that I didn't post that "solution". It wasn't correct, and it wasn't pretty.

This puzzle is not as easy as it looks. In fact, the solution is very hard to find and is very elegant. If you want to try to solve it yourself, stop reading this post, bookmark this spot, and come back later.

As I should have known, this puzzle is no secret to the Internet It is available for about $10 at many places and YouTube offers several solutions. But I'm glad that I didn't see any of them until a few weeks ago.

When I made more careful measurements on "Cheese_puzzle.png" and had a more accurate model to simulate, it became clear that my precious first solution just wouldn't fit.

Nobody said that the mouse had to fit neatly into a circular hole with just the right diameter. Once you realize that you might be looking for a noncircular hole, you are on the way to finding the solution.

The key to the puzzle's geometry is the purple axes in this picture. They cross at right angles at a point in the center hole at the tip of the protruding slice. The angle between the purple axes and the orange vertical/horizontal axes can be measured accurately on the photo and everything else follows. The angle is

theta = atand((y(b)-y(a))/(b-a)) = 11.62 degrees

The mask that I was using for my first solution got me started with polyshapes, but I eventually abandoned it. It is easy to make a polyshape directly from a list of points around the boundary, like this.

The three holes in the original cheese slices do not appear to serve any purpose. So, I have moved them and given them an important role in this simulation. Combine the two triangular slices in the upper left into a single slice. Rotate each of the resulting four slices by an angle of 180-theta = 168.38 degrees about its center. The four right-angled corners from the central hole go to the corners of the puzzle while the outer corners form a central square just large enough to hold the mouse.

That is really elegant.

Where did this puzzle originate? What mathematics was involved in its design?

My code is here: Mighty_Polyshape.m.

Get
the MATLAB code

Published with MATLAB® R2022a

I do not recall having seen this matrix before, but I will not be surprised to learn that somebody else already knows all about it, especially if that person's name is Nick.... read more >>

]]>I do not recall having seen this matrix before, but I will not be surprised to learn that somebody else already knows all about it, especially if that person's name is Nick.

I've been investigating the matrices generated by this elegant one-liner.

Q = @(n) (-n:n).^2 + (-n:n)'.^2;

The `Q` is for "quadratic".

The middle column contains the squares of the integers from `-n` to `n`. So does the middle row. The apostrophe summons singleton expansion. The resulting matrix has order `2*n+1`. Here is `Q(5)`.

Q5 = Q(5)

Q5 = 50 41 34 29 26 25 26 29 34 41 50 41 32 25 20 17 16 17 20 25 32 41 34 25 18 13 10 9 10 13 18 25 34 29 20 13 8 5 4 5 8 13 20 29 26 17 10 5 2 1 2 5 10 17 26 25 16 9 4 1 0 1 4 9 16 25 26 17 10 5 2 1 2 5 10 17 26 29 20 13 8 5 4 5 8 13 20 29 34 25 18 13 10 9 10 13 18 25 34 41 32 25 20 17 16 17 20 25 32 41 50 41 34 29 26 25 26 29 34 41 50

I like the contour plot.

contourf(Q(100)) axis square colorbar title('Q(100)')

For another blog post under development, I need a logical mask that carves a circular region out of graphic. This disc does the job.

D = @(n) Q(n) <= n^2;

Here is my carver.

```
spy(D(100))
title('D(100)')
```

Did you notice the digits in the count of nonzeros in `D(100)`? It happens whenever `n` is a power of 10.

fprintf('%15s %12s\n','n','nnz(D(n))') for n = 10.^(0:4) fprintf('%15d %12d\n',n, nnz(D(n))) end

n nnz(D(n)) 1 5 10 317 100 31417 1000 3141549 10000 314159053

A classic problem, described in the Online Encyclopedia of Integer Sequences, asks how many points with integer coordinates lie within a disc of increasing radius. Our nonzero count provides the answer.

fprintf('%15s %8s\n','n','a(n)') for n = [1:15 99:101 499:501 999:1001] if mod(n,100) == 99 fprintf('%15s %8s\n','-','-') end a(n) = nnz(D(n)); fprintf('%15d %8d\n',n,a(n)) end

n a(n) 1 5 2 13 3 29 4 49 5 81 6 113 7 149 8 197 9 253 10 317 11 377 12 441 13 529 14 613 15 709 - - 99 30757 100 31417 101 32017 - - 499 782197 500 785349 501 788509 - - 999 3135157 1000 3141549 1001 3147833

Taking the reciprocals of the matrix entries, and reducing the range of the anonymous index, produces a matrix that behaves a bit like the Hilbert matrix, `hilb(n)`.

R = @(n) 1./((1:n).^2 + (1:n)'.^2);

Here are the 5-by-5's.

```
format rat
R5 = R(5)
H5 = hilb(5)
```

R5 = 1/2 1/5 1/10 1/17 1/26 1/5 1/8 1/13 1/20 1/29 1/10 1/13 1/18 1/25 1/34 1/17 1/20 1/25 1/32 1/41 1/26 1/29 1/34 1/41 1/50 H5 = 1 1/2 1/3 1/4 1/5 1/2 1/3 1/4 1/5 1/6 1/3 1/4 1/5 1/6 1/7 1/4 1/5 1/6 1/7 1/8 1/5 1/6 1/7 1/8 1/9

Going away from the diagonal, the elements of `R(n)` decay more rapidly than those of `hilb(n)`, so `R(n)` is better conditioned than `hilb(n)`.

format short e fprintf('%15s %12s %12s\n','n','cond R','cond H') for n = 1:12 fprintf('%15d %12.2e %12.2e\n',n,cond(R(n)),cond(hilb(n))) end

n cond R cond H 1 1.00e+00 1.00e+00 2 1.53e+01 1.93e+01 3 2.04e+02 5.24e+02 4 2.59e+03 1.55e+04 5 3.21e+04 4.77e+05 6 3.89e+05 1.50e+07 7 4.67e+06 4.75e+08 8 5.54e+07 1.53e+10 9 6.53e+08 4.93e+11 10 7.65e+09 1.60e+13 11 8.92e+10 5.22e+14 12 1.04e+12 1.62e+16

What is the rank of `Q(n)`? Why? See a paper in SIAM Review by Strang and Moler.

Why is the table of values for `nnz(D(10^k))` so short? How might you extend this table?

Investigate `R(n)`. Is it positive definite? What are its eigenvalues? What is its inverse? What is the sign pattern of the elements of its inverse? For what values of `n` can you compute the inverse reliably using floating point arithmetic? How does all this compare with `hilb(n)` and `invhilb(n)` ?

Comments in the Comments, or in email to me, are welcome.

Get
the MATLAB code

Published with MATLAB® R2022a

The Soma Cube brings back memories.... read more >>

]]>The Soma Cube brings back memories.

Piet Hein (1905-1996) was an extraordinary Danish inventor, mathematician, poet and philosopher. He invented the Soma Cube puzzle in 1933. I wrote a blog post about Hein and some of his creations several years ago, Soma Cube 2016.

The Soma Cube puzzle has seven pieces. One of them is a V-shaped piece made from three cubelets. The other six pieces are L, T, Z, R, S, and Y with four cubelets each. That's a total of 27 cubelets, just enough to make a 3-by-3-by-3 cube. Sound familiar?

Bill McKeeman and I were buddies in grad school. He was a professor at U. C. Santa Cruz for a while, and then at the ill-fated Wang Institute of Graduate Studies in Tyngsborough, Mass. He worked for DEC in New Hampshire for a long time, taught compilers at Dartmouth, and even consulted for the MathWorks. As an exercise to learn MATLAB, he wrote the modern version of our `why` command.

Bill and I became obsessed with the Soma cube after Martin Gardiner described the puzzle in his *Scientific American* column. You may not have noticed it before, but one of Bill's programs, `soma`, is in the MATLAB `demos` directory. Bill generated all of the 240 distinctly different puzzle solutions and stored them in a 240-by-27 matrix, `demos/somasols`. His program lets you page through the solutions.

My new `Soma` code uses technology from `Qube`, the digital Rubik's Cube simulator, to plot the 240 solutions. Here are the seven Soma pieces, surrounding an animation stepping through every tenth solution.

Do you recognize the colors?

`Soma` is available from this link. You already have `somasols`, but another copy is available from this link.

I have combined my new display code and McKeeman's old program that finds all the solutions. The self-extracting archive is available at <https://blogs.mathworks.com/cleve/files/Soma_osf.m>

Get
the MATLAB code

Published with MATLAB® R2022b

The matrices in the following animations are at the heart of computer graphics. They describe objects moving in three-dimensional space and are essential to MATLAB's Handle Graphics, to Computer Added Design packages, to Computer Graphics Imagery in films, and to most popular video games. Many modern computers contain GPUs, *Graphic Processing Units*, which are optimized to compute the product of these matrices.... read more >>

The matrices in the following animations are at the heart of computer graphics. They describe objects moving in three-dimensional space and are essential to MATLAB's Handle Graphics, to Computer Added Design packages, to Computer Graphics Imagery in films, and to most popular video games. Many modern computers contain GPUs, *Graphic Processing Units*, which are optimized to compute the product of these matrices.

The *homogeneous coordinates* system used in today's computer graphics software and hardware makes it possible to describe rotations, translations and many other operations with 3-by-3 and 4-by-4 matrices. These matrices operate on vectors with the position of an object, *x*, *y* and *z* , in the first three components.

Rotations about the coordinate axes are described by three matrices. Rotations about the *x* -axis are produced by $R_x$, which rotates *y* and *z*, while leaving *x* unchanged.

$$ R_x(\theta) = \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & \cos{\theta} & -\sin{\theta} \\ 0 & \sin{\theta} & \cos{\theta} \end{array} \right) $$

Rotations about the *y* -axis are generated by

$$ R_y(\theta) = \left( \begin{array}{rrr} \cos{\theta} & 0 & -\sin{\theta} \\ 0 & 1 & 0 \\ \sin{\theta} & 0 & \cos{\theta} \end{array} \right) $$

And, rotations about *z* are provided by

$$ R_z(\theta) = \left( \begin{array}{rrr} \cos{\theta} & -\sin{\theta} & 0 \\ \sin{\theta} & \cos{\theta} & 0 \\ 0 & 0 & 1 \\ \end{array} \right) $$

Rotation angles are specified in degrees. Our MATLAB programs use the degree-based trig functions `cosd` and `sind`. Here are graphs of $\cos\theta$ and $-\sin\theta$ , evaluated with the angle $\theta$ going from `0` to `360` degrees in `10`-degree steps.

Here is another look at the same data, `cosd(theta)` and `-sind(theta)` for `theta = 0:10:360`. The columns of the rotation matrix are the coordinates of the rotating dots. The blue dot starts at (0,1) and the orange dot starts at (1,0).

If you drop the zeros from the values of `theta`, you are left with the integers from 1 to 36. This is the numbering in the international standard describing the compass direction of runways at airports. The initial position of our blue dot corresponds to runway 36 and the orange dot starts as runway 9.

**Note: Refresh your browser to synchronize these animations.**

For aircraft and space vehicles, rotation around the `x`-axis from nose to tail is known as *roll*.

Rotation about the `y`-axis parallel to the wings is *pitch*.

And, rotation about the vertical `z`-axis is *yaw*.

Our model of the Piper PA-24 Comanche has 97 *patches*. One of them is the propeller. This animation of a rotating propeller is very similar to our earlier animation of the compass.

`Qube`, our digital Rubik's Cube simulator, uses 27 copies of a single *cubelet*. This animation of a rotating cubelet shows a quarter turn clockwise about `x`, followed by a quarter turn clockwise about `y` and then a quarter turn counterclockwise about `z`. If these three rotations are repeated four times, the cubelet returns to its initial orientation. In the process, we see the traditional Rubik's colors of all six faces -- white, green and orange opposite yellow, blue and red.

**Note: Refresh your browser to synchronize these animations.**

Rubick's Cube is all about rotations. Rotating the cubelets in any face of the puzzle, while leaving the rest of the puzzle fixed, is called a "move". Like any cube, Rubik's cube has six faces. Each move rotates one of the six faces in either a clockwise or counterclockwise direction. So, after `n` moves, the puzzle is in one of `12^n` possible states. The challenge is to return the cube to its original orientation.

Here are six random rotations produced by `scramble(6)`. Because `12^6` is `2,985,984`, this is just one of almost three million six-move scrambles.

One possible way to restore any starting condition is to retrace the moves that produced it. This is the "follow the breadcrumbs" algorithm. So, I call this `unscramble`, rather than `solve`.

- 1: Which rotation matrices and what values of
`theta`are used in the animations?

- 2 (not easy): When is
`unscamble`a solution with the minimum number of moves?

The source code for `Qube` is available from this link: Qube_May18_osf.m. The `osf`, *one single file*, format is a self-extracting archive that expands into a directory of individual functions.

Get
the MATLAB code

Published with MATLAB® R2022a

Here is a link to "Qube, The Movie", a video made with `Qube`, my digital Rubik's cube simulator. Mathematically, all of the action is driven by the 3-by-3 matrices on display.

`Qube` uses a `3x3x3` array of identical *cubelets* . The coordinates `[x,y,z]` of the centers of the cubelets are all of the possible combinations of `-2`, `0` and `+2`. The vertices of an individual cubelet form a 8-by-3 matrix. The width of a cubelet is controlled by multiplying its vertices by a 3-by-3 diagonal scaling matrix.

Our video begins with the width near zero, so the cubelets appear to be points at the centers. The width is increased to a value where the cubelets almost touch each other. The small gap is maintained for its visual effect.

The primary operations with a Rubik's cube are the simultaneous rotation of all the cubelets contained in one of the six faces. For example, this video shows a slow-motion rotation of the "Left" face about the `x`-axis.

It is also possible to rotate the three interior "slices", as well as the entire cube about one of the coordinate axes.

The *type* of a cubelet is the number of nonzeros in the coordinates of its center.

type = nnz([x,y,z])

The cubelet in the center of the puzzle has `type = 0` and the six cubelets with `type = 1` are located at the center of each face. We use `0:1` to denote the set of these seven cubelets. In a real, physical Rubick's Cube, `0:1` is a single solid central core that holds the entire puzzle together.

There are twelve `type = 2` cubelets located on the edges of each face and eight `type = 3` cubelets at the corners of the puzzle. So, `0:3` denotes the entire cube and `2:3` is the corners and edges without the central core. Qube, The Movie shows the types in this order:

0:3, 0:2, 0:1, 0, 1, 2, 3, 2:3

Careful readers of this blog should recognize `2:3` as level one of the Menger sponge fractal.

Get
the MATLAB code

Published with MATLAB® R2022a