# Solving Partition Problems with Linear Programming

Yesterday I promised to show how to solve partition problems with integer linear programming. Let us take the problem of yesterday: We want to split the numbers from 1 to 30 in 10 triplets (x,y,z) such that x+y+z=0 modulo 31. We will use linear problems of the following type.

$$Ax=b, \quad x_i \in \{0,1\}, \quad c^T x \to \text{max.}$$

There will be one variable in x for each possible triplet. The value of this variable determines if the triplet is in the selection (1) or not (0). So the j-th column of A refers to the triplet with number j. The i-th row of A is a linear constraint which makes sure that the number i is used only in one triplet. A contaoins only 0-1-values, and

$$a_{i,j} = 1 \quad\Leftrightarrow\quad i \in T_j = \{n_{1,j},n_{2,j},n_{3,j}\}.$$

The constraints will be

$$\sum_j a_{i,j} x_j = 1.$$

We just need any feasible point and could use any target function. Alternatively, we could use „less equal“ in the previous line and maximize the number of triplets used.

How to do that in EMT? We first store all possible triplets into one matrix.

>function makeT (n:index=31) ...
$T=[];$for i=1 to n-2;
$for j=i+1 to n-1;$       k=n-mod(i+j,n);
$if k!=i and k!=j then T=T_[i,j,k]; endif;$   end;
$end;$return T
$endfunction >makeT(7) 1 2 4 1 4 2 1 6 7 2 4 1 2 5 7 3 4 7 3 5 6 3 6 5 5 6 3 >T=makeT(31);  Then we define the matrix A. >function makeA (T) ...$v=sort(unique(flatten(T)));
$A=zeros(cols(v),rows(T));$for i=1 to cols(v);
$for j=1 to rows(T);$       if any(v[i]==T[j]) then A[i,j]=1; endif;
$end;$end;
$return A;$endfunction
>makeT(5)
1             4             5
2             3             5
>fraction makeA(makeT(5))
1         0
0         1
0         1
1         0
1         1
>A=makeA(T);
>size(A)
[31,  405]


The function makeA() can take any matrix with partitions in its rows. The first command simply selects the numbers that appear in any partition. In our case, A is simply the vector [1,…,30].

Now we need to solve the problem with integer linear programming. We use the LPSOLVE library for this (given to EMT by the developers).

>function solveP (A,T) ...
$x=ilpsolve(A,ones(rows(A))',ones(cols(A)),$         vlb=zeros(cols(A)),vub=ones(cols(A)),>max);
$return T[nonzeros(x')];$  endfunction
>solveP(A,T)
1             5            25
2             7            22
3            11            17
4             6            21
8             9            14
10            23            29
12            24            26
13            19            30
15            20            27
16            18            28


The return value of the function ilpsolve() is a 0-1-vector. We want to print the triplets which are marked by 1 in this vector. The variables vlb and vub are lower and upper bounds for the variables. Interestingly, the solution works without these restrictions. Nevertheless, more restrictions usually mean shorter calculations.

# Modulo Arithmetic with Euler Math Toolbox

Recently, I came across the problem to compute expressions in the modulo arithmetic with integer numbers. Here is an example: We like to compute expressions modulo 31. Since 31 is prime the numbers 0 to 30 with addition and multiplication modulo 31 will be a field. Each number a between 1 and 31 will have a unique multiplicative inverse b, such that ab=1 modulo 31. E.g.

>mod(27*23,31)
1


The inverse of 27 is 23 modulo 31.

How can one find the inverse element? The fastest way uses Euclid’s algorithm. If you compute the greatest common divider of 27 and 31 using Euclid’s algorithm you can derive x and y such that 27x+31y=1. Then x (taken modulo 31) will be the inverse of 27 modulo 31. Here is the computation for these numbers:

$$31=27+4, \quad 27=6 \cdot 4 + 3, \quad 4 = 3 + 1,$$

and by computing backwards

$$1 = 4 – 3 = 7 \cdot 4 – 27 = 7 \cdot 31 – 8 \cdot 27.$$

Note -8=23 modulo 31. We get our result 23.

Of course, we could program that in EMT or Maxima. Maxima contains already a function for the inverse modulo a number.

>&inv_mod(27,31)
23


We also know the small Lemma of Fermat.

$$a^{p-1} = 1 \text{ mod } p.$$

Thus

$$a \cdot a^{p-2} = 1 \text{ mod } p.$$

This can also be used to compute the inverse. To get the result in EMT we need a trick which I am going to explain below. But in Maxima we can compute very large powers, albeit very slowly.

>&mod(27^29,31)
23


For single computations or small numbers, symbolic computation is okay. But assume we need more speed. The computation above does not work in EMT, because the power is too large. Python, by the way, has infinite integers too. So the computation works there too.

>>> print 27**29
323257909929174534292273980721360271853387
>>> print 27**29 % 31
23


One of the reasons of this posting is to show how this can be done in EMT and other languages without an infinite integer arithmetic, and moreover much faster. The trick is the fact that we can as well take the modulo in each step of the multiplication. Here is the first simple code.

>mod(27^29,31) // WRONG RESULT!
12
>function map powmod (a,n,p) ...
$b=1;$  loop 1 to n; b=mod(b*a,p); end;
$return b$  endfunction
>powmod(27,29,31) // CORRECT RESULT!
23


However, this is still not the fastest way. For powers, we can use this trick.

$$x^{2y} = \left(x^y\right)^2, \quad x^{2y+1} = \left(x^y\right)^2 \cdot x.$$

This leads to the following code. By the way, the function matrixpower() of EMT uses the same technique. To understand the code, note that it is not efficient to use a mapping function recursively. Thus we call it from another function.

>function powmodrek (a,n,p) ...
$if n==0 then return 1$  elseif n==1 then return a
$elseif mod(n,2)==0 then$     b=powmodrek(a,n/2,p); return mod(b*b,p);
$else$     b=powmodrek(a,(n-1)/2,p); return mod(mod(b*b,p)*a,p);
$endfunction >function map powmod (a:number,n:nonnegative integer,p:index) ...$  return powmodrek(a,n,p);
$endfunction >powmod(27,29,31) 23 >powmod(27,1:30,31) [27, 16, 29, 8, 30, 4, 15, 2, 23, 1, 27, 16, 29, 8, 30, 4, 15, 2, 23, 1, 27, 16, 29, 8, 30, 4, 15, 2, 23, 1]  The last computation is not effective, of course. If we need to compute all powers we should use a simple loop instead. But I included that example to show that powmod() can map to vector input. As you see, the powers of 27 do not cycle through all the number 1 to 30. For other base numbers, they do. >powmod(3,1:30,31) [3, 9, 27, 19, 26, 16, 17, 20, 29, 25, 13, 8, 24, 10, 30, 28, 22, 4, 12, 5, 15, 14, 11, 2, 6, 18, 23, 7, 21, 1]  The second reason for this posting is to demonstrate the solution to another problem. We want to partition the numbers 1 to 30 into triplets with a+b+c=0 modulo 31. This is surprisingly easy. We have for all q different from 1 $$1+q+q^2 = \dfrac{q^3-1}{q-1}$$ Thus by the Lemma of Fermat $$3^0 + 3^{10} + 3^{20} = \dfrac{3^{30}-1}{3^{10}-1} = 0 \text{ mod }31.$$ Our triples can be found easily from this. Just multiply the equation with powers of 3. >for k=0 to 9; ... > powmod(3,k,31)+"+"+powmod(3,10+k,31)+"+"+powmod(3,20+k,31)+"=0", ... >end; 1+25+5=0 3+13+15=0 9+8+14=0 27+24+11=0 19+10+2=0 26+30+6=0 16+28+18=0 17+22+23=0 20+4+7=0 29+12+21=0  It is also possible to find such a partition with any desired property using linear programming and LPSOLVE in EMT. This should be the topic of another posting, however. # Euler Math Toolbox and Matlab Recently, this site got more traffic than usual. So I thought I welcome all of you in the new year (which is already three weeks old now). And I provide some new content for you to read. For me, the selected topic is important. I welcome very much any discussion about it. I am currently writing a requested paper for an on-line journal about Euler Math Toolbox, so I could benefit from your input. What is the difference between Euler Math Toolbox and Matlab? Or rather, what is the intended difference? And are there alternatives to both? First of all the following points are obvious to me. • Matlab is profiling itself in the professional market. It is called the „industry standard“, and it is advertised as a must-know to students. Matlab tries very much to serve the needs of the industry. The best example is its simulation toolbox which connects Matlab to hardware. Matlab is backed up by a group of commercial programmers. • Euler Math Toolbox is for education and research. It is not intended to be used as a professional tool in industry. EMT is backed up by open software. • Both systems fall short for specialized applications. Most professional or research software is compiled in a general programming language and does things that a one-for-all software cannot do. • We need to educate our students in general programming, not in software systems, especially not in commercial systems. With a good background in programming any software can be adapted easily. • That said, a tool like Euler Math Toolbox or Matlab can be useful for quick computations or to mathematical demos. And for this, both systems are equally well equipped with an advantage for EMT. • We should have free and open, reliable numerical libraries for our main computer languages. If you are interested in the primary differences between EMT and Matlab, here is a list. • EMT is open and free, Matlab is commercial with a considerable pricing. • The user interface of EMT is notebook oriented similar to Maple or Mathematica, Matlab is command oriented. Both have scripts, of course, and are similar in this respect. • Both systems can do symbolic computations, EMT with the free Maxima, Matlab with an additional package you need to buy. • EMT can use many open external programs like Povray, Python, C, Scilab, Latex. For Matlab, you can buy powerful external packages to do industry strength computations. • Matlab has an optional compiler. In EMT, you need to write C or Python for more speed. • Matlab is used all around the world, and you will easily find someone with the same problem as you have. EMT has a community, but it is way smaller. There might be more differences that I forgot to mention. But all in all it is clear to me that EMT deserves to be used in education much more than it is right now. Matlab, on the other hand, should be removed from general education. It should be reserved for specialized courses done by engineers who have used the program for their work. Alternatives? There are many. • R is a statistical package featuring a matrix language as well. It is free and has a big community. The programming in R is not as easy to learn as they claim, but to get statistics done it is not necessary. The interface is ugly and cumbersome to use. • Maple, Mathematica or other algebra software are commercial packages in the higher price range. They offer more advanced features than Maxima and often yield better results. But they are not designed to be used as numerical software. • Geogebra tries to be the knife for all purposes. It is a beautiful software for schools and has a large community. Clearly, it falls short in the area of numerical and other programming. Actually, it is a whole different category of software. • Python is probably the best alternative to big software packages. There are very nice packages for plotting and for numerical computations. The symbolic package falls short in comparison, however. The Sage project relies completely on Python and has a web interface which clearly is the proper way in the future. • General programming langes are the way to do software education in universities for engineers and math students. Most courses have an obligatory course in Java, and rightly so. Open for discussion. # Pure Symbolic Functions in Maxima and EMT I like to demonstrate a few tricks with Maxima and EMT. Assume you want to generate the Vandermonde matrix, and check the identity $$\det \begin{pmatrix} 1 & x_0 & \ldots & x_0^n \\ \vdots & \vdots & & \vdots \\ 1 & x_n & \ldots & x_n^n \end{pmatrix} = \prod\limits_{0 \le i<j \le n} (x_j-x_i)$$ You can generate this matrix easily in EMT for specific values. >n=3; v=0:3; V=v'^(0:3) 1 0 0 0 1 1 1 1 1 2 4 8 1 3 9 27 >det(V) 12 It is a bit more difficult to check the right hand side of the equation. Of course, you can simply write a function >function vprod (v:vector) ...$p=1; n=length(v);
$for i=1 to n;$   for j=i+1 to n;
$p=p*(v[j]-v[i]);$   end;
$end;$return p;
$endfunction >vprod(v) 12 It is also possible with Matrix tricks. We generate all differences in a matrix, square the matrix, set the diagonal to 1, take the total product of all matrix elements and take the 4-th root. >(prod(prod(setdiag((v-v')^2,0,1))'))^(1/4) 12 But that is not the general case with variables. To generate this we use purely symbolic functions in EMT. These functions evaluate in Maxima at the time when they are used, not when they are defined. >function mrow (k,v,n) &&= create_list(v[k]^m,m,0,n) m create_list(v , m, 0, n) k >function Vand (v,n) &&= apply(matrix,create_list(mrow(k,v,n),k,0,n)) apply(matrix, create_list(mrow(k, v, n), k, 0, n)) >&Vand(v,2) [ 2 ] [ 1 v v ] [ 0 0 ] [ ] [ 2 ] [ 1 v v ] [ 1 1 ] [ ] [ 2 ] [ 1 v v ] [ 2 2 ] This is almost self explaining, once you know that „matrix(a,b,c)“ generates a matrix in Maxima with rows „a,b,c“. The creation takes place for the specific symbolic v when „Vand(v,2)“ is called. Maxima prints the result with indices. But it is a string of the following form. >V &= Vand(v,2); >""+V matrix([1,v[0],v[0]^2],[1,v[1],v[1]^2],[1,v[2],v[2]^2]) Since EMT vectors start with index 1 by default, we can only evaluate this expression for zero based vectors (unless we modify the functions above a little bit). >v=0:2; zerobase v; >V() 1 0 0 1 1 1 1 2 4 Maxima can factor the determinant of any Vandermonde determinant of specific size. >&factor(determinant(Vand(v,3))) (v - v ) (v - v ) (v - v ) (v - v ) (v - v ) (v - v ) 1 0 2 0 2 1 3 0 3 1 3 2 This is in accordance with the theory. Let me try the following matrix. $$\begin{pmatrix} (x_0-x_0)^n & \ldots & (x_0-x_n)^n \\ \vdots & & \vdots \\ (x_n-x_0)^n & \ldots & (x_n-x_n)^n \end{pmatrix}$$ >function mrow (k,v,n) &&= create_list((v[k]-v[m])^n,m,0,n) n create_list((v - v ) , m, 0, n) k m >function Vand (v,n) &&= apply(matrix,create_list(mrow(k,v,n),k,0,n)) apply(matrix, create_list(mrow(k, v, n), k, 0, n)) >&Vand(v,2) [ 2 2 ] [ 0 (v - v ) (v - v ) ] [ 0 1 0 2 ] [ ] [ 2 2 ] [ (v - v ) 0 (v - v ) ] [ 1 0 1 2 ] [ ] [ 2 2 ] [ (v - v ) (v - v ) 0 ] [ 2 0 2 1 ]  I am not sure if you are aware of the following formula. >&factor(determinant(Vand(v,3))) 2 2 2 2 2 9 (v - v ) (v - v ) (v - v ) (v - v ) (v - v ) 1 0 2 0 2 1 3 0 3 1 2 (v - v ) 3 2  The factor 9 is interesting. For n=4, it is 96, and for n=5, it is 2500. Do you guesswork. The factorization of the case n=4 takes some seconds. I did not try higher degrees. Of course, it is easy to compute this numerically, even for very large matrices. >function d1(v) := det((v'-v)^(cols(v)-1)) >d1(0:10) 1.43679021499e+072  # December Problem I found the following problem on the page of the wonderful Science Blogs. It is in German. Let me translate it. Bernd writes two different numbers between 1 and 1000 on two cards of paper. Anna selects one, and after seeing it has to decide if she wants to keep it or change to the other card. The player with the higher card wins. This is a case of a two person game with the following tactics for each player. • Bernd: He can write any combination (i,j) with 1 <= i < j <= N. • Anna: She can decide to change up to a number A with 1 <= A <= N. Note that we simplified the situation for Anna. It does not make sense to change with 3, but not with a smaller number. Weather or not this simplification makes the outcome worse for Anna remains to be investigated. Now, Bernd can select any of his tactics (i,j) with probability p(i,j), and Anna can select any of her tactics with probability q(A). The problem to select the optimal strategy (i.e. the distribution of probabilities) for each player can be solved with optimization. I have done similar problems with EMT on this page. I solved the problem above for (N=10) with a bit of code. However, I want to spoil the solution here. If you want to try on your own, do not read further. • Strategy for Bernd: Bernd should select the pairs (1,2),…,(N-1,N) at random with equal probability. • Strategy for Anna: Before each game, Anna should select a number A from 1 to N-1 at random, and switch if her number is less or equal A. Let us prove that this is optimal. Assume, Anna selects tactics A. So she will change cards at any number smaller or equal A. Then Bernd will lose 1/(N-1) on average with his strategy. This is so, since the only case that makes a difference for him is the pair (A,A+1), where Anna wins all the time. So, no matter what tactics A Anna selects, the worst case is an average loss of 1/(N-1). Only if Anna would select A=N, he would do better. Now assume, Bernd selects tactics (i,j). Then Anna, with her strategy, will win in any case if i<=A<j. In all other cases, the average result is 0. We need to find the worst case for Anna based on her strategy. That is indeed the case, if Bernd selects a pair (i,i+1) and it is equal to a win of 1/(N-1). It is now clear that there cannot be a better strategy for either player. Anna cannot win more than Bernd loses on average. With a bit more care in the proof, it is also clear that it would not help Anna to add more strategies. Assume, she rejects some numbers, and accepts others, instead of rejecting the numbers 1 to A altogether for some A. Then Bernds tactics will lose less or equal 1/(N-1) on average. For an example, assume Anna rejects 1 and 3. Then (1,2) will be a loss for Bernd, (2,3) will be a win, and (3,4) will be a loss, all other cases are 0 on average. This is a very nice problem. I could not find a reference in the net, however. But it should be known. Let me add the code for EMT. >v1=flatten(dup(1:10,10)); >v2=flatten(dup((1:10)',10)); >i=nonzeros(v1<v2); v1=v1[i]; v2=v2[i]; >A=(1:10)'; >M=-((v1<=A)&&(A<v2)); >function game (R) ...$## Maximizes G such that R'.q>=G, sum q_i = 1, q_i>=0.
$## Returns {q,G}.$    m=rows(R); n=cols(R);
$M=R|-1_ones(1,n); b=zeros(m,1)_1; eq=ones(m,1)_0;$    c=zeros(1,n)|1; restr=ones(1,n)|0;
$res=simplex(M,b,c,eq,restr,>max);$    return {res[1:n]',res[-1]};
\$endfunction
>{pB,res}=game(M); fraction res
-1/9
>k=nonzeros(pB); fraction (pB[k]_v1[k]_v2[k])'
1/9         1         2
1/9         2         3
1/9         3         4
1/9         4         5
1/9         5         6
1/9         6         7
1/9         7         8
1/9         8         9
1/9         9        10
>{pA,res}=game(-M'); fraction res
1/9
>fraction pA'|(1:10)'
1/9         1
1/9         2
1/9         3
1/9         4
1/9         5
1/9         6
1/9         7
1/9         8
1/9         9
0        10

This looks a bit difficult. Here are some explanations.

The vectors v1 and v2 simply contain the pairs (i,j) with 1 <= i < j <= N = 10. The trick to get all pairs is to create matrices with 1 to 10 in each row resp. column and to „flatten“ them to vectors. Then we select only the pairs with i < j using the nonzeros() function of EMT.

The matrix M is then constructed in the following way: The element in row A and column (i,j)  is -1 (a loss for Bernd) if i <= A < j. This is the matrix of average outcomes if Anna selects A and Bernd selects (i.j).

The game() function calls the simplex algorithm for the following problem.

$$Mp \ge G, \quad \sum_k p_k = 1, \quad G \to \text{Max.}$$

The probabilities p(k) are then maximizing the worst case G, for all rows of M (tactics of Anna). The result is p and G. We print the nonzero elements of p along with the corresponding i and j.

The same is finally done with the problem seen from Anna’s side.

By the way, the fact that both results must be equal (besides the sign) is an application of duality in optimization.

Really a nice problem!

# Today in the Classroom

Me: Instead of the sum of squares we can also try to minimize the maximal distance

$$f(c) = \max_\limits{1 \le k \le n} |x_k-c|$$

What would be the result, assuming we have the points in order

$$x_1 < \ldots < x_n$$

Silence …

After two painful minutes, Me: Okay, let’s make a sketch. If c is here

What is the value of f(c)?

Silence …

I can see from some faces that they know the answers, or at least have an idea about it. But still, silence.

Only after a long while with me impatiently waiting for any response, I can hear a whisper: In the middle.

So what is this? Are they afraid to show silliness? No posing? Is this problem so hard? I cannot say.

Currently, I assume that these Bachelor students are trained to passively absorb math. It starts in school, and continues, much worse, in University. How we can break through this I do not know. There are some legal problems involved. You can no longer enforce active participation. That was an illusion anyway. But we should generate a climate where students don’t feel bad about making mistakes. After all, learning math is not flight training. How to do that against the habits of all other classes is something I would love to learn.

In a recent blog entry, I moaned about the missing logical education in schools. But what is logic actually, and why is it important?

Let us take the viewpoint of the foundations of mathematics for a moment, and sketch the main ideas behind mathematical truth very briefly. We distinguish between axioms and theorems. The theorems are derived from the axioms by logical reasoning. We best understand correct logical reasoning if we look at models of our theory. A model of a theory is a system where the axioms are fulfilled. Now, logical reasoning is correct, if in any model where the axioms are true, the derived theorems are also true. If there is one model for the axioms where the theorem is not true, then the theorem cannot be derived by correct logical steps.

As an example, we take vector spaces. There are lots of them. First, the simple Euclidean plane and space. But then are also spaces of more dimensions, and even spaces of infinite dimensions. There are also vector spaces consisting of functions. In all these spaces, the axioms for vector spaces are true, and thus are all the theorems we derive for them. E.g., every finite set spanning the space contains a basis of the space, or the theorem that all bases have the same length.

Do we need that in school? You may argue that in school we have only models, and no theories. E.g., the kids learn about the numbers 1,2,3,… and how to compute with them. There is only one model of the natural numbers, or not? We have enough troubles learning how to handle them. So why bother with logical reasoning?

There are a number of arguments against that (besides the subtle fact that there are different models of arithmetic).

First, some theorems are not obvious. As an example, take the sum of the numbers 1 to N. Now we can just learn the formula for this sum and apply it, like drilled monkeys. But I cannot imagine one young kid that is not asking „Why?“.  Once they get older some are so frustrated by math that they just do not want to know. But if you see the formula applied for the first time, there must be wow effect.

Then, some things are wrong, albeit they look right. Examples are numerous in probability theory. To be able to tell false from right should be a desirable skill for everyone. In mathematics, you have means to achieve that goal. That is why the analysis of the sequence of equations in the previous blog post is a good exercise. The job is to fill in the logical steps correctly, and point to the one that does not hold.

Finally, finding a proof is a fun game. That is, unless you are trained to hate math. You may fail, and fail again. But if you succeed the first time you will feel great. It is like solving crosswords, or, more modern, succeeding all levels of a computer game. Why not offer these challenges to the kids instead of boring application of formulas?

I do not accept the argument that there is no time for that in schools. You can save so much time omitting boring things like hand computations. I also do not accept the argument that most kids are not clever enough. Reasoning can be on very low levels. Why is 56*98=98*56? Let the small ones think about such stuff and appreciate the demonstration of flipping over a rectangle. Why is the prime number factorization unique. Do teachers even ask this question? Sadly, they don’t.

# Teachers forget the Logic

I found this in my Google+ box some time ago. The comments showed clearly the disastrous confusion that teachers produce when they refuse to teach logic. Today, blogs like this do not get a very thorough reading. But the matter is important. So let me put it in bold right at the beginning. Read aloud:

Never write a sequence of equations without the logic that connects them!

Or let me phrase it in a harder, more explicit way: Real mathematicians never omit the logic of their arguments. Only teachers do that. Have I been clear enough? As a professor of math I have to correct a lot of things the school spoiled. So get at least this right!

It is perfectly okay to set x to some value. It is also okay to multiply an equation by some value as long as you know that this is an implication only. So let us check.

$$x=(\pi+3)/2 \\ \implies 2x = \pi+3 \\ \implies 2x(\pi-3) = (\pi+3)(\pi-3) \\ \implies \ldots \\ \implies (3-x)^2=(\pi-x)^2$$

Every step of this chain of implications is perfectly valid. Note that the converse implications are only valid if we assume pi not equal to 3.  For ac=bc implies a=b only if c is not equal to 0. The problematic step is the following. I write it with the correct implication.

$$(3-x)^2 = (\pi-x)^2 \,\Longleftarrow\, 3-x = \pi – x$$

The other implication is not possible. It does not hold always. You may write

$$a^2 = b^2 \,\Longleftrightarrow\, a = \pm b$$

But then you have to clarify the meaning of the plus/minus sign on the right side. I prefer

$$a^2=b^2 \,\Longleftrightarrow\, a=b \text{ or } a=-b$$

for the sake of clarity.

The missing logic is the reason for the nonsense at the end. If you do not believe that writing down and making clear of the logic is necessary either read the comments after the post in Google+, or even better, discuss this with your own class. You might be shocked.

One final word: We may introduce the convention that writing down a sequence of equations means implication. After all, that reflects our thinking. But why not write it down explicitly? Would that really be too hard to do?

# Linear Optimization with Maxima and EMT

Maxima and EMT do of course have implementations of the Simplex algorithm. The algorithm in Maxima seems to work differently for floating point numbers and for fractions. Here is an example in fractions.

We solve a problem in standard form.

$$Ax=b, \quad c^Tx \to \text{Min.}$$

Then Maxima has the function linear_program().

>&load(simplex);
>A &:= [1,1,-1,0;2,-3,0,-1;4,-5,0,0]
1             1            -1             0
2            -3             0            -1
4            -5             0             0
>b &:= [1,1,6]
[1,  1,  6]
>c &:= [1,-2,0,0]
[1,  -2,  0,  0]
>&linear_program(A,b,c)

13     19        3
[[--, 4, --, 0], - -]
2      2         2

Since we have carefully defined the data as numerical and symbolic variables, we can use the same data for the simplex() algorithm in EMT.

>x=simplex(A,b',c,eq=0,>min); fraction x'
[13/2,  4,  19/2,  0]

Note „eq=0“, which sets all equations (in each row of A) to equalities. The default is „eq=-1“ which means inequalilites of the form „<=“. The solution would then be different.

There is another form of the command in Maxima.

>&minimize_lp(x-2*y, ...
>  [x+y-z=1,2*x-3*y-v=1,4*x-5*y=6,x>=0,y>=0,z>=0,v>=0])

3              19             13
[- -, [v = 0, z = --, y = 4, x = --]]
2              2              2

With only two variables, we can even draw the feasible set. With three variables, we can use Povray to do this. In the examples of EMT, you find the following plot.

# Povray in Euler Math Toolbox

EMT can do a lot of 3D plots, but it is not an OpenGL application. This means that it has restrictions concerning hidden surfaces and intersecting objects. But if you install Povray you can do a lot more. Note that you need to put the path to the installation into your system path, or set the „povengine“ variable as described in the tutorial about Povray in EMT. In this tutorial you will find a lot of examples to get you started. Combining EMT with Povray is a fascinating application of EMT.

The first simple function is pov3d(). This function is similar to plot3d(). Here is an example.

>load povray
Functions to generate Povray files for the open Raytracer.
>pov3d("x^2+y^3",axiscolor=red,angle=20°, ...
>  look=povlook(red,0.2),level=-1:0.5:1,zoom=3.8);

As you see, you can use an expression, set angles, zoom, and levels just as in plot3d(). The new parameter is „look“, which makes the plot a transparent red. The output will cast shadows from the light source. Moreover, there is a coordinate axis. You can control all that in the plv3d() function.

The graphics is not in the graphics window, but in the text window, just like inserted graphics from EMT. Moreover, the graphics file is in the working directory of EMT. By default, this is the directory „Euler“ in your user directory („C:\Users\Name\Euler“). By default, it is a 450×450 PNG.

But you can do much more in EMT with Povray. EMT can handle Povray is by Povray objects stored as Povray commands in strings. For some large objects files are used. If you generate a Povray scene with these objects, you need to start the generator with povstart() and finish with povend(). That command will also start the Povray program. Here is an example.

>povstart(angle=-10°,center=[0.5,0.5,0.5],zoom=7);
>x=0:0.02:1; y=x'; z=x*y; vx=-y; vy=-x; vz=1;
>mesh=povtriangles(x,y,z,"",vx,vy,vz);
>cl=povdisc([0.5,0.5,0],[1,1,0],2); ...
>ll=povdisc([0,0,1/4],[0,0,1],2);
>writeln(povdifference(mesh,povunion([cl,ll]),povlook(green)));
>writeln(povintersection([mesh,cl],povlook(red))); ...
>writeln(povintersection([mesh,ll],povlook(gray)));
>writeln(povpoint([1/2,1/2,1/4],povlook(gray),size=2*defaultpointsize));
>writeAxes(0,1,0,1,0,1,d=0.015); ...
>povend(w=600,h=600);

Let me explain the commands one by one.

First, we start the new scene with povstart(). This generates a Povray file in the user directory. It does also set the camera angle, the point to look at, and the camera zoom. It can set a lot more. You need to read the documentation for Povray in EMT for more information.

We now generate objects. The povtriangles() functions generates a file containing the description of a mesh and a string containing the Povray command for this mesh. The string „mesh“ simply contains „object { mesh1 }“.

We then generate more strings containing commands for Povray objects. E.g., cl is a disk which is later used to intersect the green surface with in height 1/2.
>cl
cylinder { <0.492929,0.492929,0>, <0.507071,0.507071,0>, 2
texture { pigment { color rgb <0.470588,0.470588,0.470588> }  }
finish { ambient 0.2 }
}

The function difference() then is a command computing the difference of the mesh „mesh“ and the disk „cl“.

>povdifference(mesh,povunion([cl,ll]),povlook(green))
difference {
object { mesh1 }
union { cylinder { <0.492929,0.492929,0>, <0.507071,0.507071,0>, 2
texture { pigment { color rgb <0.470588,0.470588,0.470588> }  }
finish { ambient 0.2 }
}
cylinder { <0,0,0.24>, <0,0,0.26>, 2
texture { pigment { color rgb <0.470588,0.470588,0.470588> }  }
finish { ambient 0.2 }
}
texture { pigment { color rgb <0.470588,0.470588,0.470588> }  }
finish { ambient 0.2 }
}
texture { pigment { color rgb <0.0627451,0.564706,0.0627451> }  }
finish { ambient 0.2 }
}

We use writeln() to write this string to the Povray file which has been opened in EMT with the open() command. That is the main trick to generate Povray scenes in EMT.

Finally, povend() generates the file. This time, we specify a larger image which can be found in the user directory of EMT.

There is a lot more you can do with this combination of EMT and Povray. In fact, you can do everything, Povray is capable of. You may need to define new functions for this. Have a look at the tutorial and the visit the links at the end for more examples.