## Floating Point Arithmetic – Youtube Video

I demonstrate floating point arithmetic and big floats in this video on Youtube. I will continue this in a series with videos about EMT.

## Let’s talk about Logic

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

Of course, the simplex() command and the Maxima commands can do a lot more. EMT does also have efficient packages for integer programming. Start with the tutorial to learn more about this.

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.

## Random Dices

This is just an answer to a question that came up in a discussion. Can you simulate 2 dice throws and display the frequencies in EMT?

Actually, this is quite simple. We demonstrate it with 20 repetitions first. The trick in a matrix language is always to do all random stuff beforehand, and then to apply matrix functions to get the desired statistics. In our case sum(W)‘ does the job. It sums all rows of W and returns a vector with of sums. Then getmultiplicities() extracts the statistics. In our case, we know that the numbers 2 to 12 can appear in the sum. If we don’t know that, sort(unique()) extracts a sorted set of unique elements in the data.

>n=20; >W=intrandom(n,2,6) 6 2 4 1 5 3 1 3 5 6 3 3 1 2 6 4 5 6 3 2 1 3 5 4 6 3 3 3 6 1 3 3 6 2 6 6 4 4 5 3 >sum(W)' [8, 5, 8, 4, 11, 6, 3, 10, 11, 5, 4, 9, 9, 6, 7, 6, 8, 12, 8, 8] >f=getmultiplicities(2:12,sum(W)') [0, 1, 2, 2, 3, 1, 5, 2, 1, 2, 1]

We change n to 20000, and see the expected distribution.

>columnsplot(f,lab=2:12):

## About Multi-Variate Functions in EMT

You define a vector valued function of two variables in EMT as follows.

>function f([x,y]) := [x^2+y^2-10,x*y-1] >broyden("f",[1,2]) [0.317837, 3.14626]

Then f can be used in the form f(x,y) or f(v) with a 1×2 vector v. The Broyden algorithm needs the second form. This can also be done for a symbolic function.

>function f([x,y]) &= [x^2+y^2-10,x*y-1] 2 2 [y + x - 10, x y - 1] >function Df([x,y]) &= jacobian(f(x,y),[x,y]) [ 2 x 2 y ] [ ] [ y x ] >newton2("f","Df",[1,2]) [0.317837, 3.14626]

The function f can now work numerically just as above. As a symbolic function, it is defined in Maxima too. There, it can only be used as f(x,y). Likewise for the derivative of f defined in Df. The Newton algorithm needs this derivative matrix to do its job.

Note, that the Broyden and the Newton algorithms work for row and for column vectors alike. It depends on the start vector. E.g., if you start with a row vector, the algorithms expect the functions to return row vectors.

You can also use the following trick.

>function f(v) &= [v[1]^2+v[2]^2-10,v[1]*v[2]-1] 2 2 [v + v - 10, v v - 1] 2 1 1 2 >function Df(v) &= jacobian(f(v),[v[1],v[2]]) [ 2 v 2 v ] [ 1 2 ] [ ] [ v v ] [ 2 1 ] >newton2("f","Df",[1,2]) [0.317837, 3.14626]

The first definition is a valid definition of a numerical function. But for Maxima, the symbolic function can also be defined this way. Maxima interprets v[1] as a variable, as you see in the definition of Df. Of course, you can no longer write f(x,y).

## Minimize the Surface of a Cubuid

The problem I discuss in this blog post is from the nice page by Dan Meyer. The blog is interesting anyway because of the way Dan thinks about teaching. The concept „less helpful“ says it all. I was sometimes asked by a newbie teacher how to plan a course. My answer was: „Make them work!“. And with „work“ I am not talking about stupid fact learning. It is all about questions, discoveries, and the satisfaction to have found an answer, even it is not the best one. They will appreciate the solution only after having appreciated the problem.

Sorry about this digression. The problem is to minimize the surface

\(S = 2(ab+bc+cd)\)

given the volume

\(V = abc.\)

With no further restrictions the answer is …

That is one of the points where Dan would allow guesswork. Even better, he would allow you to discover the problem itself from real world situations. And to see, why it is a problem at all.

You can surf the answer, or find it using Lagrange’s method or any other method involving Calculus. The result is no surprise.

\(a=b=c=V^{1/3}.\)

If students do this I often have to ask the questions they do not ask. E.g., why is the solution of the Lagrange equations a minimum? Couldn’t it be a maximum? What happens if a, b, or c tend to 0? Can you find a proof without derivatives?

But the cited page does pose the constraint that a,b,c should be integer. This changes the problem and puts it to a higher level. In fact, I am quite sure that there is no easier algorithm than trying all triplets of divisors of V. There are some shortcuts, for sure, and to find them is a nice job. For a start, try simple numbers like 24.

The twist on the page is something deeper, however. Dan asks the readers for easy algorithms, and then puts up a contest to break these algorithms. This leads to endless trials and errors, discovering nice and useful math along the way. I like this very much.

Let me give you an example. The most obvious candidate for a solution is to take divisors close the cube root of V, unless the remaining divisor is a large prime. You need to go a long way to find counterexamples to this. I had to write a program to find one, 5850. The minimal solution does not contain the divisors closest to the cube root, and none is a prime. But read the comments on that page for yourself.

How this all translates into teaching at a university, I do not know. But I know that it should.

## Visualizing the Newton Algorithm in the Plane

Continuing my last blog entry, let us try the Newton Algorithm with two variables. First we define the function f(x,y). We use an example where we know the zeros exactly.

\(f(x,y) = (x+iy)^3-1.\)

The zeros are the three roots of unity. In EMT, we define the real and the imaginary parts of this function as f1(x,y) and f2(x,y).

>fe &= z^3-1 with z=x+I*y 3 (I y + x) - 1 >function f1(x,y) &= realpart(fe) 2 3 - 3 x y + x - 1 >function f2(x,y) &= imagpart(fe) 2 3 3 x y - y

We can now plot the function.

>plot2d("abs(f1(x,y)+f2(x,y)*I)",r=1.2,>contour,>hue,>spectral,n=100);

We can now define the function f([x,y]) as a hybrid function in EMT which can be used as f(x,y) or f(v) numerically, and likewise the derivative matrix of this function, courtesy by Maxima. Furthermore, we define the Newton step.

\(f(v) = v – Df(v)^{-1} \cdot f(v).\)

>function f([x,y]) &= [f1(x,y),f2(x,y)] 2 3 2 3 [- 3 x y + x - 1, 3 x y - y ] >function Df([x,y]) &= jacobian(f(x,y),[x,y]) [ 2 2 ] [ 3 x - 3 y - 6 x y ] [ ] [ 2 2 ] [ 6 x y 3 x - 3 y ] >function newtonstep ([x,y]) &= transpose([x,y]-invert(Df(x,y)).f(x,y));

The function is the following.

\begin{pmatrix}\frac{2\,x\,y^4+4\,x^3\,y^2-y^2+2\,x^5+x^2}{3\, \left(y^2+x^2\right)^2} & \frac{2\,y\,\left(y^4+2\,x^2\,y^2+x^4-x \right)}{3\,\left(y^2+x^2\right)^2} \\ \end{pmatrix}

To get the Latex code for this, use the tex() command of Maxima.

>tex(&factor(newtonstep(x,y))) \begin{pmatrix}\frac{2\,x\,y^4+4\,x^3\,y^2-y^2+2\,x^5+x^2}{3\, \left(y^2+x^2\right)^2} & \frac{2\,y\,\left(y^4+2\,x^2\,y^2+x^4-x \right)}{3\,\left(y^2+x^2\right)^2} \\ \end{pmatrix}

Let us test.

>v=[1.5,0.4]; loop 1 to 6; v=newtonstep(v), end; [1.11995, 0.197797] [0.988758, 0.043586] [0.998147, -0.000885968] [1, 3.29448e-006] [1, 1.74668e-011] [1, 0]

Okay, [1,0] is one of the zeros. Finally, we plot one path of the Newton Algorithm.

>v=[0.5,0.3]; w=v; loop 1 to 10; v=newtonstep(v); w=w_v; end; w, 0.5 0.3 0.794694 -0.665052 0.584502 -0.137812 1.21661 0.321028 0.994208 0.110139 0.987868 0.000509852 1.00015 -1.26776e-005 1 -3.78536e-009 1 0 1 0 1 0 >plot2d(w'[1],w'[2],color=red,thickness=2,>add):

We see that the algorithm can take quite erratic paths. So, it is interesting from which points the algorithm converges to which zeros. We simulate that.

>function map res(x,y) ... $v=[x,y]; $loop 1 to 100; $ if all(v~=0) then return 1; endif; $ vnew=newtonstep(v); $ if all(vnew~=v) then $ if vnew[1]>0 then return 3; $ elseif vnew[2]>0 then return 4; $ else return 5; $ endif; $ endif; $ v=vnew; $end; $return 1; $endfunction >res(0.5,0.3) 3

Since we used the „map“ keyword we can apply the code to many pairs [x,y]. It will vectorize to a matrix of X and Y values too. We use it to compute the results of 500×500 points. On my computer, this takes a few seconds.

>X=linspace(-1.2,1.2,500); Y=X'; >R=res(X,Y); >plotrgb(R):

As you see there are three main regions of convergence. The image is a nice fractal between these main regions.

## Visualizing the Newton Method in EMT

In a post by squareCircleZ I found a visualization of the Newton method with JavaScript. This looks like a modern and sleek way of providing demos to students. I have some doubts about the depth of the insights gained by this click and try demos. But, for a first clue, they seem to be okay.

However, I’d much prefer this to be done with a little bit of programming in EMT. Here is the code. I am going to explain it below.

>function f(x) &= x^4+x-1 4 x + x - 1 >function df(x) &= diff(f(x),x) 3 4 x + 1 >function newstep (x) &= x - f(x)/df(x) 4 x + x - 1 x - ---------- 3 4 x + 1 >&ratsimp(newstep(x)) 4 3 x + 1 -------- 3 4 x + 1 >plot2d("f",0,2); >function vstep (x) ... $ plot2d([x,x],[0,f(x)],style="--",color=blue,>addpoints,>add); $ x1=newstep(x); $ plot2d([x,x1],[f(x),0],style="-",color=blue,>add); $ return x1; $endfunction >x=1.9; >loop 1 to 8; x=vstep(x); end:

It is all quit straightforward. Using Maxima, we compute the derivative or our symbolic function f(x), and then the Newton step. The simplification with ratsimp() is only to show that the original form can be simplified quite substantially. We plot the function, and define a function vstep() which draws two lines (one with endpoints) and returns the new value of the Newton algorithm. Then we do this 8 times starting from x=1.9. Note the colon : after the last statement, which inserts the graphics into the notebook.

If you save all that in an EMT notebook it should provide a nice environment for own experiments by the students. With a little bit of thought, we can do a lot more now. E.g., we could try to experimentally prove the the order of the algorithm is quadratic. This involves logarithms and quite a bit of thought.

By the way, Maxima can also give us the formula for the root of this particular function. As you will see, it is not very useful.

>$solve(f(x),x)[2]